[gmx-users] problem with free energy calculations

[David Mobley]

Michael,

> 2) i get the same results when using md+NH thermostate
> with two thermostates, or langevin dynamics (sd), the
> error bars are somewhat larger with sd though.

Well, that's fine -- but it doesn't mean that everything is OK. NH has
known problems. There is also a bugzilla relating to this
implementation of NH not sampling the correct distribution of
temperatures. I guess your result means that it doesn't matter (at
least in this case) for free energy calculations; I'm just not
convinced that this is never something to be worried about.

> 3) i get the same results when doing slectrostatics
> and VdW seperately or in one go, (provided i use p=2).

That's useful to know. Have you compared efficiency?

> system/protocols: 462 TIP3P water in a cubic box, 12 A cut-off,
> PME with nk=32 or 64. free energy from trapecoidal integration of
> <dU/dlam> from 21 windows, at least 100 ps per window,
> sc_alpha = 0.5-2.5, sc_p= 1 or 2, sc_sigma=0.3.
> error-bars taken to be the RMSD of dGs obtained
> from 10ps/window subsets.

Again, we did a lot of work optimizing the choice of soft core
parameters. For LJ it is optimal to use sc_power=1 with sc_alpha=0.5.
This actually makes a *big* difference; it's not a slight improvement.
Even sc_power=1 and sc_alpha=0.7, say, can be a lot worse... (as can
sc_power=2). It makes a big difference in how much structure there is
in the dV/dlambda curve, and thus how many lambda points you need to
get your integration error sufficiently small. It has similar effects
for BAR (which is what we use, since it doesn't make integration
error).

If you're really set on finding the most efficient way of doing
things, here's the test I would do. Compare accuracy for the following
two scenarios:
1) Doing Coulomb and LJ separately, using sc_power=0 for the coulomb
part and sc_power=1 with sc-alpha=0.5 for the LJ part.
2) Doing Coulomb and LJ together in whatever manner you find works best.

My feeling is that (1) is more efficient (in that you can get lower
uncertainties with fewer lambda values or less simulation at each
lambda value) but I must admit I never really put this to the test,
because I never found a soft core set of parameters that would work
reasonably efficiently for (2). If you try this out, I'd be really
interested in the outcome...

It's probably also worth pointing out that if you do use BAR, you can
typically get by with wider lambda spacings for the same amount of
uncertainty than you can with TI (plus you don't make the integration
error). Unfortunately BAR isn't implemented yet so you have to do it
by storing trajectories and then reprocessing them, although I think
Erik is working on this.

Speaking of integration error, that reminds me of when I did
discharging of an SPC water in SPC water. This is a fairly smooth
curve, and I can do it fine with 4 or 5 lambda value s with BAR, but
TI (with the trapezoid method) makes an integration error of roughly
0.3 kcal/mol with 5 lambda values (0.3 out of about 8 kcal/mol). You
have to go up to 9 or 10 lambda values to get this down to what you
get with BAR for 5 lambda values, because the curve has a uniform
curvature that you botch with the trapezoid method.

David


> for what its worth ...
>
> mic
>
>
>
>
>
> ____________________________________________________________________________________
> Never Miss an Email
> Stay connected with Yahoo! Mail on your mobile.  Get started!
> http://mobile.yahoo.com/services?promote=mail
> _______________________________________________
> gmx-users mailing list    gmx-users at gromacs.org
> http://www.gromacs.org/mailman/listinfo/gmx-users
> 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
>