Compute TI and bonded interactions

Dear Anna,

Thanks for your email. I have a coarse grained polymer system with FENE
bonds and LJ pair-wise repulsion (self avoiding chain). What I did is to
look at the free energy profile based on thermodynamic integration, in
which case I could look at the free energy profile when lj component is
turned off and on, but turning off FENE bonds cannot really be turned off,
and therefore I am interested whether there is a way to see the free
energy contribution of the FENE bonds?

TI provides differences in free energies e.g. from a reference state, but why
would you want to compare the free energy difference between a bonded
and non-bonded system?

Do not forget that nature does not distinguish between pair and bonded
interactions, in reality bonded interactions are just pair interactions with
such a high barrier, that we implicitly coarse-grain them into permanent
bonds. This is computationally very convenient, but does not necessarily
have anything to do with nature.

Perhaps a better reference state for TI would be to take a random chain
conformation, and replace bond interactions by self-springs fixing the
monomers to their instantaneous position in space and perhaps fix the
self-spring constant to reproduce the same monomer fluctuations as in
the polymer?

Anyways, if you take a step back (ignoring pair interactions (and units)
for simplicity), then in the absence of bonds you have an ideal gas, the
partition function is Z ~ V^N where N is the number of particles, and V
the volume.

If you link identical monomers into a RW polymer the partition function
becomes Z ~ V (4pi)^(N-1). Here I assume you integrate over the position
of the first monomer, and bond lengths are fixed so the remaining (N-1)
integrals are just over the solid angles specifying the chain conformation.

In both cases F= kT ln(Z) ~ kT N * const1 + const2 for large N, with
different constants for the free and bonded cases. For more complicated
polymer models, I would expect the same generic behaviour with
suitably renormalised constants.

In the case where you replace bonded interactions with self-springs, the
partition function should still be Z ~Veff^N, where Veff ~ sqrt( k/kT )^3
is the effective volume of monomer fluctuations and k the spring. This
follows since 0.5 k <X^2> = 0.5 kT by the equipartition theorem. But if
you now gradually switch on pair AND bonded interactions, while letting
the self-spring constant go to zero, that is perhaps an interesting free
energy difference, and I would expect that this can be done with LAMMPS?
Note you can still derive the exact free energy of the reference state
analytically.

Note that the self-springed state is fixed in space, while the final state
is free to move, hence you acquire a -kT Log(V) additional free energy
due to the translational entropy of the final state. I would subtract
this off the TI result, to get a number that is reflects the free energy
contribution due to the specific molecular interactions. Secondly note
that this should be averaged over several initial states.

Dear Carsten,

Thank you so very much for your insightful email, which includes both the understanding of theoretical and simulation concepts and beautifully combines them together; I started thinking about the appropriate reference state as well, but I think your thoughts are justified and clear and I shall adapt them to the problem I am working on. Thank you!

With best wishes,
Anna