[gmx-users] Re: problem with free energy calculations

[Berk Hess]

>From: Michael Brunsteiner <mbx0009 at yahoo.com>
>Reply-To: Discussion list for GROMACS users <gmx-users at gromacs.org>
>To: gmx users <gmx-users at gromacs.org>
>Subject: [gmx-users] Re: problem with free energy calculations
>Date: Tue, 6 Feb 2007 12:55:06 -0800 (PST)
>
>
>Thanks for your comments!
>
>Berk Hess wrote:
>
> > One "obvious" problem is that using two thermostats is incorrect.
> > In this way the com of the water with respect to the com of the rest
> > will not be sampled correctly (both are fixed).
>
>I am not sure I understand what you mean, is the COM reset
>for each temperature group seperately or not ??
>(I used nstcomm = 100) ... when the systems are still
>interacting, i.e., when lambe < 1 I wouldn't expext
>this to be a problem, while when the water is completely
>de-coupled (lambda=1) then I actually NEED to use two
>thermostates plus a finite nstcomm, since otherwise I
>might end up with a water thats either much too "hot"
>or too "cold", or a tumbling ice-cube ...
>

You indeed need to couple the two groups separately
when approaching the decoupled state.
But you also want to sample the motion of the one
water with respect to the rest of the system correctly.
For that you would need a third thermostat for these
three degrees of freedom. Gromacs does not have such
a feature and I guess neither do many other packages.
In a recent paper I found out that for amino-acid analog
solvation free energies this gives an error of roughly 2 kJ/mol.
So it is probably not the explanation for your problem.

> > The solution is using the sd integrator.
>
>There must be a way to do that with a normal verlet algorithm
>or the like ... the sd integrator, i.e. the Langevin dynamics, was
>meant to be used with implicit solvent, after all. There is
>this friction term that's meant to mimic the effect of the
>explicit water, and even if you let the friction constant zeta
>be very large, i.e., if you do away with the friction term, you
>just end up with a primitive kind of thermostate, and i can't
>see a reason why this should be any better than Nose Hoover !?

SD is not a primitive kind of thermostat. It is a pefectly valid thermostat
that gives you a guaranteed Boltzmann ensemble, without the nasty
problems that global thermostats have for systems with nearly decoupled
modes. When you take the friction coefficient smaller than the viscosity of 
water
the effect of friction is negligible and anyhow for free-energy calculations
you don't care about the dynamics.

>
> > Another issue is that since water is strongly hydrogen bonding
> > and the hydrogen bond interaction is rapidly soft-cored,
> > you probably need much closer spaced lambda point at lambda=0
> > (i.e. use an inhomogeneous spacing).
>
>that would mean to forgo one of the advantages of softcore
>potentials altogether, since this is exactly why we use
>them: to spare ourselves the extensive sampling required
>when approaching lambda=1 (where the interactions between
>particle and the rest approach zero, and dA/dlambda
>diverges)

No.
You need the soft-core to avoid the singularities at r=0.
The problem is that water is very strongly hydrogen bonding
and you need to brake these hydrogen bonds.
How this happen along your lambda path depends very
much on the choice of interpolation and also on the water model.

>
>I've simulated the very same system in CHARMM, with
>soft core potentials plus NH thermostate, and 21 equally
>spaced windows, and got perfectly reasonable results.
>can it be that the small LJ interaction parameters
>that CHARMM gives the TIP3P water hydrogens (for reasons
>related to a completely different problem, and in fact
>irrelevant for the resulting (free) energies) makes
>all the difference ??

CHARMM might use a very different way of interpolating
between lambda=0 and lambda=1. But I would not know
what it does.

>From your lambda curve you should be able to judge
how if interpolation is reasoable.

A much more robust way of calculating free energies
in hydrogens bonding systems is doing it in two steps:
first make the Coulomb interactions dissappear and
then the LJ interactions. This can be done with
less lambda points.

Berk.

_________________________________________________________________
Play online games with your friends with Messenger 
http://www.join.msn.com/messenger/overview