[gmx-users] Trying to replicate Aqvist's results (solvation free energy).

[Francis de Lasalle]

Hello,

As we're trying to get familiar with Gromacs, in order to be able to 
perform solvation free energy calculations of ions, we're trying to 
replicate the results obtained by Aqvist (the author of the article from 
which were obtained Lennard-Jones parameters for alkaline and alkaline 
earth ions, for the OPLS-AA force field). The solvation free energies 
we're computing are all about 50 kJ/mol different from the values 
reported by Aqvist. We verified that the Lennard-Jones parameters for 
our ions and for the water are identical to that of Aqvist and we're 
also using SPC water. We're wondering why we can't reproduce these 
results and any input would be greatly appreciated. Here's the article :

http://pubs.acs.org/doi/pdf/10.1021/j100384a009

We're using 200 lambda states with 2000000 steps incremented by 0.002 ps 
(4000 ps). Here's our run file for Li+ :

integrator               = sd
nsteps                   = 2000000
dt                       = 0.002
nstenergy                = 1
nstlog                   = 1
nstxout                  = 0
nstvout                  = 0
cutoff-scheme            = group
rlist                    = 1.0
dispcorr                 = EnerPres
vdw-type                 = cut-off
rvdw                     = 1.0
coulombtype              = pme
rcoulomb                 = 1.0
fourierspacing           = 0.12
constraint-algorithm     = SHAKE
continuation             = no
shake-tol                = 0.00001
constraints              = all-angles
tcoupl                   = v-rescale
tc-grps                  = system
tau-t                    = 0.2
ref-t                    = 298
ref-p                    = 1
compressibility          = 4.5e-5
tau-p                    = 5
free-energy              = yes
couple-moltype           = lithium
sc-power                 = 1
sc-sigma                 = 0.3
sc-alpha                 = 1.0
couple-intramol          = no
couple-lambda1           = vdwq
couple-lambda0           = none
init-lambda-state        = $LAMBDA
fep-lambdas              = $LAMBDAS

The box has been equilibrated for 1 ns at 298 K. The box is a 
dodecahedron with 240 water molecules.

Thanks a lot.

Francis