[gmx-users] Deviations in free energies with slow growth (single and 3-step process)

Fri Feb 8 21:48:27 CET 2008

Maik,

> Thanks for your answer.
> I'm aware of that problem, but the idea was, that such a small system is
> very close to equilibrium on a 10ns timescale (maybe to optimistic).

I am not sure the issue with slow growth is as simple as "what
timescale does it take for the system to equilibrate?". In particular,
with repeated realizations of the "same" slow growth simulation (i.e.
from different starting configurations) you should get a distribution
of work values; the width of the distribution will be determined by
how slowly you do the calculations. But I don't see any reason why
this distribution should be that narrow even on 10 ns timescales.

> Actually, the thought behind it was to compare the results of the
> different free energy calculation methods.

If you want to do that, instead of doing slow growth you may want to
use the Jarzynski method to analyze output of slow growth simulations.

> Now, if I assume slow growth
> has a systematic error of X because of non-equilibrium,

That's just the thing -- it won't have a systematic error. It will
give you a distribution of different answers if you do it again and
again, since it's nonequilibrium. For any particular simulation, you
have no idea how the computed "free energy" will relate to the true
free energy you would get if you were able to evaluate the free energy
using i.e. the Jarzynski relation. Hence, it's not a systematic error,
but rather a random error (from a particular distribution).

>I expect this
> error to occur more or less in the 1_step and in the 3_step calculation.
> However, I didn't mention, but I also did non-equilibrium tests (e.g.
> BAR) with this system and the deviation is exactly the same;

Sorry, I'm unclear what you mean here. BAR is usually used to analyze
results of equilibrium simulations (i.e., simulations run in
equilibrium at distinct lambda values). There is of course the Crooks
equation (generalization of Jarzynski) that looks like the BAR
equation but is for nonequilibrium switching; is that what you're
talking about?

And I'm not clear what you mean when you say, "I expect this error to
occur more or less in the 1_step and in the 3_step calculation."

> means:
> BAR gives a difference of 7 kJ/mol for the ethane to methanol in 1_step
> hardcore and 3_step hc/sc/hc for the first and third being coulomb and
> the second being Lennard Jones.

Oh, you mean using the Crooks relationship for computing free energies
using nonequilibrium switching from ethane-> methanol and methanol->
ethane? Of the nonequilibrium approaches you've described that seems
most reasonable.

I think if you want me to comment further I'll probably need a table
showing what exact calculations you've done, with a brief description
of the protocols, and the results -- I'm getting a bit jumbled up
about all of the different things you've tried. Is this supposed to be
a methods comparison paper or something?

> This actually puzzles me a bit, cause in the case of some tripeptide
> sidechain morphing (SER to ALA), the results of 1_step and 3_step
> perfectly match. I don't have to tell you, that the cycle together with
> the sim in gas phase gives a difference in free energy of solvation
> which is somewhat 15 kJ/mol away from the experiment, though. ;)

Again, if you do a particular slow growth simulation (what you're
calling 1_step, I think?) and compare results with equilibrium
simulations, or slow-growth broken into components, and manage to get
perfect agreement, you're probably just lucky (or unlucky, depending
on how you look at it). You could even try repeating the slow growth
simulations from a different starting config -- I bet you'd get a
different result, indicating agreement is coincidental.

> For me the question arises, how to publish such "shitty" results.
> Do you think, 7 kJ/mol lies within the usual error of free energy
> calculations? If, one could never resolve reliably smaller free energy
> differences and the whole method won't be applicable for say drug
> screening or so.

This is of course why error analysis is essential. When I do these
sorts of calculations  I make sure the statistical uncertainty is much
lower than 7 kJ/mol. As far as your results go, since I don't have
enough detail on your protocol I just don't know.

David

> Regards
>
> Maik Goette, Dipl. Biol.
> Max Planck Institute for Biophysical Chemistry
> Theoretical & computational biophysics department
> Am Fassberg 11
> 37077 Goettingen
> Germany
> Tel.  : ++49 551 201 2310
> Fax   : ++49 551 201 2302
> Email : mgoette[at]mpi-bpc.mpg.de
>          mgoette2[at]gwdg.de
> WWW   : http://www.mpibpc.gwdg.de/groups/grubmueller/
>
>
>
> David Mobley wrote:
> > Maik,
> >
> >> I simulated a switching process (slow growth TI) over 10ns of ethane to
> >> methanol with hardcore slow-growth in water (365 TIP4P waters,PME) and
> >> in vacuum. The thermodynamic cycle of this calculation yields a DeltaG,
> >> which is in perfect agreement with the experiment and other calculations.
> >
> > Slow growth is notoriously problematic, because it only gives the true
> > free energy in the limit that you do it infinitely slowly. Otherwise
> > you're really measuring nonequilibrium work values, which are
> > connected to free energies but only in an average way (see the
> > Jarzynski relationship for a more recent discussion of this).
> >
> > This could be a classic example of "just because a method gives
> > perfect agreement with experiment, you can't be sure it's giving the
> > correct value for the force field." In other words, since force fields
> > have their own limitations, we can't necessarily expect that when we
> > do things right, we'll get optimal agreement with experiment. In some
> > cases one may get better agreement with experiment by doing things
> > *wrong* than by doing them right. i.e., suppose the correct value for
> > the force field is actually off by 2 kJ/mol from the experimental
> > value. If you make a protocol mistake, there's roughly a 50% chance
> > (assuming the errors are randomly distributed!) that it will give you
> > a value that's closer to the experimental value than your original
> > value was.
> >
> >> Now, when splitting this simulation into 3 steps (3x 3.2ns), where I
> >> switch the charges to zero in the first step (hardcore), morph the LJ
> >> and bonded in the second (softcore or hardcore) and switch the charges
> >> back on from zero in the third step (hardcore), all in solvent, the sum
> >> of the single contributions does NOT yield the same number as the
> >> "single-step switching".
> >> In vacuum, the work values I get from the single step process and from
> >> adding up the values from the 3-step process perfectly match.
> >>
> >> The topologies for all systems are the same (except +/- water) for the
> >> 1-step and the 3-step simulations.
> >>
> >> Now I'm a bit puzzled and can't get an idea, where this difference in
> >> the solvated system may come from. I exclude a problem due to softcore,
> >> cause I also simulated the 3-step process all hardcore...
> >>
> >> Any ideas?
> >
> > I'm betting your problems are largely due to the fact you're doing
> > slow growth. Again, see the Jarzynski relationship; I especially
> > recommend his Phys. Rev. E. paper on convergence (2001?). It's also
> > worth noting that slow growth transformations in different directions
> > have different convergence properties, and with slow growth
> > (interestingly enough) particle insertion actually converges more
> > quickly than particle deletion. So you could try reversing the
> > direction of the transformations and see what happens. But again, slow
> > growth isn't recommended these days -- you should probably make sure
> > you have a very good reason for doing it and are aware of the
> > limitations.
> >
> > Best,
> > David
> >
> >
> >> Regards
> >>
> >> --
> >> Maik Goette, Dipl. Biol.
> >> Max Planck Institute for Biophysical Chemistry
> >> Theoretical & computational biophysics department
> >> Am Fassberg 11
> >> 37077 Goettingen
> >> Germany
> >> Tel.  : ++49 551 201 2310
> >> Fax   : ++49 551 201 2302
> >> Email : mgoette[at]mpi-bpc.mpg.de
> >>          mgoette2[at]gwdg.de
> >> WWW   : http://www.mpibpc.gwdg.de/groups/grubmueller/
> >> _______________________________________________
> >> gmx-users mailing list    gmx-users at gromacs.org
> >> http://www.gromacs.org/mailman/listinfo/gmx-users
> >> Please search the archive at http://www.gromacs.org/search before posting!
> >> Please don't post (un)subscribe requests to the list. Use the
> >> www interface or send it to gmx-users-request at gromacs.org.
> >> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> >>
> > _______________________________________________
> > gmx-users mailing list    gmx-users at gromacs.org
> > http://www.gromacs.org/mailman/listinfo/gmx-users
> > Please search the archive at http://www.gromacs.org/search before posting!
> > Please don't post (un)subscribe requests to the list. Use the
> > www interface or send it to gmx-users-request at gromacs.org.
> > Can't post? Read http://www.gromacs.org/mailing_lists/users.php
> >
> > .
> >
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> Please search the archive at http://www.gromacs.org/search before posting!
> Please don't post (un)subscribe requests to the list. Use the
> www interface or send it to gmx-users-request at gromacs.org.
> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>