[gmx-users] problem calculating PMF with distant constraints

[rzangi at att.net]

> >>> > >>>I am trying to calculate potential of mean force between two atoms 
> >>>(neon,
> >>> > >>
> >>> > >>methane) by running multiples simulations with distances that are 
> >>>constrained
> >>> > >>from lambda=0 to lambda=1 and then integrating <dHdl>. However, this 
> >>>doesn't
> >>> > >>seem to give the right curve as compare to calculations obtained 
> >>>with other
> >>> > >>methods. There appears to be a serious discrepancy.
> >>> > >>
> >>> > >>the dHdl curve is meaningless, only the integral is physcially
> >>> > >>meaningfull. Of course you can compare to other simulations done in 
> >>>an
> >>> > >>identical fashion.
> >>> > >
> >>> > >
> >>> > > We compared difference free energies (i.e. the PMF curves 
> >>>supperimposed on
> >>> > each other) with identical settings and interaction parameters.
> >>> > >
> >>> > There are many hidden parameters, not only in gromacs but also in 
> >>>other
> >>> > programs, like combination rules.
> >>> > By the way you still don't say how you compute it.
> >>>
> >>>We used the same combination rule as the reference and also used the 
> >>>second type of combination rule.
> >>>
> >>>Below is one of the topology files. Similar results we got also with SPC, 
> >>>SPCE and tip5p, and also with different LJ particles.
> >>>
> >>>The locations of the minimums of the PMF are good but the depth at 
> >>>contact is higher than the state at large distances, where is should be 
> >>>lower.
> >>
> >>
> >>This sounds like you might have forgotten about the trivial volume or 
> >>entropy term.
> >>You can not directly integrate the constraint force.
> >>You need to correct for the fact that the volume sampled by the two 
> >>molecules
> >>increases with the constraint length.
> >>The free-energy volume contribution is:
> >>-2 k T log(r)
> >>The force thus is:
> >>2 k T / r

This correction term (2kT/r) comes to be exactly the same as the thermal average of the correction to the froces due the centripetal force [JCP 100 9129(1994)]. We already did this correction but still it didn't give the PMF of the reference system. At short distances the corrections was better than at larger distances (as one would expect from better rotational sampling at shoter constraint length) but still it was very different (for example the difference between the first minimum and the first maximum). At large distances it actually made the PMF worse. At each lambda point (51 in total) we ran for 6ns simulations and at some points for even longer. The value of <dgdl> seems to converage (as indicated by block analysis - g_analyze).


> >
> >Does this apply to any PMF calculation, even when no constraints are 
> >involved (i.e. measuring the force at discrete distances)?
> 
> Yes.
> This applies to any PMF calculation where you look at a distance,
> but the particles are free to rotate in 1 or 2 other directions.
> 
> The volume or entropy contribution to the PMF is:
> -k T log(V) = -k T log(c r^(D-1)) = -(D-1) k T (log(r) + log(c))
> where D is the dimensionality of the space the particles move in
> and c is a constant that drops out in the force.
> 
> This correction is equivalent to the 4 pi r^2 correction in an rdf.
> 
> Not that in some cases one might actually want this term in the PMF.
> But not if you want effective interactions between particles.
> The term always makes the PMF go down as -log(r) at large distances.
> You can find many artcicles in the literature with such incorrect
> effective interactions.
> 

As for David question:
The reference system was with 0.95nm cut-off, and that is what we took when we tried to it compare with. So the box is big enough for the cut-off.

I am trying now to calculate the total forces on the `frozen' atoms.

Groetjes.