[gmx-users] error estimates of free energy calculations

[David Mobley]

Hi,

Berk may weigh in on this, but...

> To test my protocol, I have been trying to reproduce some of Berk
> Hess's hydration free energy results from his 2006 paper (J. Phys.
> Chem. B 110, 17616-17626, 2006). Although my obtained average free
> energies are quite comparable, my error estimates are much higher
> than the errors reported in that paper. I'd like to figure out what's
> causing this difference in accuracy, because I really need optimal
> accuracy for my future calculations.
>
> I'm using GROMACS 3.3.1, the G53a6 force field and 3.0 nm
> dodecahedron simulation boxes with ~631 SPC waters. After
> mimimization and 500 ps of equilibration for each of the 21 equally
> spaced lambda-values, I simulated 5 ns for the propane, phenylalanine
> and asparagine analogues. Theoretically, this should give me a
> significantly better accuracy than the paper's 500 ps production runs
> (approx. 1/sqrt(10) times better).

Looks to me like Berk is using 42 lambda values in solvent, as he
spaces them more tightly around the feature in the curve, and 29 in
vacuum. This should give lower uncertainties.

Without digging too much into details on my part (I could check these
things), here are some other things I'd check if I were you:
1) What thermostat did he use? Same friction?
2) You're constraining all bonds? Did he?
3) Is his method for computing the uncertainties the same as yours?
4) Barostat differences could be significant... What barostat are you
using versus his?
5) Are your equilibration protocols the same?

> My accuracy is also significantly worse than the ones reported in A.
> Villa, J Comput Chem 23, 548-553, 2002 and J. Maccallum, J Comput
> Chem 24, 1930-1935, 2003. And it doesn't even come near to the
> extraordinary accuracy reported by Michael Shirts in his papers (who
> btw also uses 5 ns runs).

If I remember correctly, Michael averaged over multiple simulations at
each lambda value (or rather, ran multiple simulations at each lambda
value and averaged the total free energies). He was also using Bennett
acceptance ratio (BAR) which he had in his own implementation of
GROMACS; this gives higher accuracies with the same amount of data,
generally speaking. I use BAR, too (see my JPCB paper from this year
on solvation free energies) but I do it by reprocessing (if you want
details I can explain).

Supposedly Erik Lindahl *almost* has implemented the ability to write
out the potential energies for BAR into GROMACS, but I haven't heard
anything on this in a month or two.

Best wishes,
David Mobley
UCSF