[gmx-users] TI, sampling, sc_power, and sc_alpha

[mernst at tricity.wsu.edu]

Greetings once again fellow Gromacs-users,

I have seen a substantial and disturbing change in calculated free energy differences
between systems run with Gromacs 3.3 and Gromacs 3.3.1. I am trying to determine if
these differences are due to the new software, new parameters, insufficient sampling, or
something else, and I don't have enough computer time to investigate every possible
combination. I turn to the list for insight.

With Gromacs 3.3, I was using thermodynamic integration over 21 equally spaced and
independent (delta_lambda=0) windows with the AMBER99 force field ported by the Pande
group. I was running two slightly different DNA systems and using thermodynamic cycles
to find the quantity of interest by comparing free energy differences. I was using soft
cores and incorrectly using sc_alpha=1.51 because I did not realize that the default
sc_power was 1 instead of 2 as in the original soft core method. Nonetheless, my
calculations seemed to be in reasonable agreement with experiment and with previous
calculations I had performed using NWChem. Recent calculations performed via a slightly
different method using AMBER99 within AMBER itself are also in reasonable harmony with
these results.

When Gromacs 3.3.1 was released, I was chastised by the software for not explictly
setting sc_power, and that's when I realized that the default was 1 and that I was using
an incorrect value of sc_alpha. I undertook to rerun my calculations with 3.3.1, the
same sc_power=1, and the corrected sc_alpha=0.7. I am now getting qualitatively wrong
results, though. According to experiment and previous calculations, I should be seeing a
difference between my two systems of about -10 kj/mol. In my current simulations I am
seeing a difference of about +20 kj/mol. This is after 300 ps of equilibration in each
window followed by 500 ps of energy statistics collection. The starting structure for
each window was taken from the final frame of a 2 ns ordinary MD production run with
free_energy=no and init_lambda=0.

The obvious possibility seems to be that the difference in soft cores has nudged the
trajectories so that the luck of the draw has me collecting less representative
statistics. Is such a large change between free energy cycles credibly attributable to
slight sampling differences, though? 30 kj/mol seems pretty substantial. On the other
hand, my systems do involve a vast number of degrees of freedom (350-700 DNA atoms, ~20
of which have altered B-states, ~40,000 water atoms, 11 or 22 Na+) and the absolute free
energy changes within each system are on the order of 150 kj/mol.

Is it possible that the change in sc_alpha alone is responsible for the difference seen
with my older runs and newer ones? I have looked at the log files and parameters don't
seem to be different between the old systems and the new ones. I know that there was
some change in how B-state parameters were assigned with 3.3.1, but for a while I ran a
CVS version of Gromacs that included this B-state change and I didn't see any alteration
in free energy calculations. Has there been any other change to free energy calculations
in 3.3.1?

If there are no obvious things to check I will probably be trying fewer windows, with
nonlinear spacing, and see if that gives me more believable results that converge
faster.

Matt Ernst
Washington State University