[gmx-users] rcoulomb, rlist

[David Mobley]

All,

I need some help understanding the behavior of GROMACS 3.3. (
Apologies on the reference to the old version, but I published some
data that now others are having trouble reproducing, which
necessitates me using the version I originally ran the calculations
in).

I am computing the hydration free energy of acetamide in TIP4P-Ew
water using the FFAMBER ports, and I am having a problem in the
component of the calculation where the A state has acetamide in water
with no partial charges on it, and the B state has acetamide in water
with all of its interactions with water turned off (but intramolecular
interactions on) which is accomplished through some topology file
trickery (modifications to the pairs section of the topology and
explicit specification of intramolecular nonbonded interactions).

Anyway, the expected value for this component calculation is something
around -1.8 kcal/mol. My published result (in GROMACS 3.3) was -2.1
kcal/mol (+/- 0.05 or so). I tried to reproduce the result of -2.1
kcal/mol to confirm that the difference was legitimate, and found that
I instead got -4.1 kcal/mol roughly.

After some testing, I established that the difference was due to a
single setting: In my original calculations, I had used PME, with
rcoulomb = 0.9 nm and rlist = 1.0 nm. When I tried to reproduce these,
I used PME, with rcoulomb=1.0 nm and rlist=1.0 nm.

Can someone help me understand what exactly is going on here with this
(seemingly small) change in the mdp files making such a large
difference in the computed hydration free energies? I realize that in
later (3.3.1 and after) versions of gromacs, the code insists that
rcoulomb=rlist with pme, essentially (which still seems strange to
me), but I still am not clear on why this is (and what exactly is
going wrong in the calculations).

For what it's worth, I still get the correct density of water with
both sets of parameters. Yet this parameter change seems to cause a
*systematic* difference of some sort -- for example, if I look at the
plot of dV/dlambda values
(http://www.dillgroup.ucsf.edu/~dmobley/uploads/novdw_comparison.pdf)
the differences between the two cases are entirely within a specific
range of lambda values, and are in a consistent direction. Yet I find
no systematic differences in system volume or density at these lambda
values.

That's my question -- any help will be appreciated.

The astute will notice I now also have a separate problem: If the
"correct" input parameters to use are rcoulomb=rlist=1.0 nm, then the
correct result is -4.1 kcal/mol, which is even further from the
expected answer (-1.8 kcal/mol) than I started off with. This is the
expected answer, in that it is the answer I arrive at using other
simulation packages (DESMOND, for example) with the same input
parameters, and the answer my collaborators (using several other
packages) are getting with the same input parameters. Obviously this
is a cause for some concern and suggests there may be some sort of an
underlying problem here. On my list of things to do is to see whether
I still get the same results in GROMACS 4.

Thanks for your help,
David Mobley