[gmx-users] problem calculating PMF with distant constraints

[Berk Hess]

>From: rzangi at att.net
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>Subject: Re: [gmx-users] problem calculating PMF with distant constraints
>Date: Fri, 16 Dec 2005 15:25:52 +0000
>
> > >>> > >>>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).

You could indeed also call this term a centripetal correction.

I get perfect PMF's between two ions in water up to 1.2 nm, for distances
between 1 and 1.2 nm I used 8 ns.
Colleagues of mine used 10 ns for similarly long distances between 
hydrophobic
molecules in water.
So your 6 ns should be enough.

But the results are very sensitive to the simulation setup.
You will for instance see differences between RF with a 1.4 nm cut-off and 
PME.

Berk.