[gmx-users] Free energy calculation problems (bug?)

[Олег Титов]

Dear GROMACS users!

I've encountered a problem in thermodynamic integration with GROMACS. I'm
interested in calculating of solvation free energy differences between
holbenzenes and benzene. To calculate this values I use thermodynamic
integration approach (5 point formula).

Half a year ago I did my calculation with GROMACS 4.6.1. The results looked
good and were consistent with literature data. Recently, I've discovered
that versions 4.6.0 and 4.6.1 had a bug in free energy code, so I decided
to recalculate everything with 4.6.5 version of GROMACS. New results look
not so good and, moreover, are not consistent with literature.

I use GAFF forcefield with TIP3P water models. To obtain GROMACS topology
and coordinate files I convert AMBER topologies (generated with tleap) with
amb2gmx.pl script and this topologies look fine.

 For phenylbromide molecule calculated solvation free difference should be
near -0.1 - -0.2 kcal/mol (Experimental data is 0.59 kcal/mol).
4.6.1 version gave me -0.1+-0.13 kcal/mol - Fine!
I've performed the same calculation with AMBER11  - -0.27 +-0.16 kcal/mol -
also fine.
4.6.5 version for the same system gives me 1.1 - 1.2 kcal/mol

I've played with some options (tried adding soft-core potentials and turned
on/off couple-intramol option). The results are:
no_sc no_intramol 1.27 +- 0.18
no_sc intramol 1.20 +- 0.20
sc no_intramol 1.35 +- 0.21
sc intramol 1.22 +- 0.19

So now I don't know what to do. Old buggy version produced good results
consistent with literature and other MD packages and new fixed version
produces different results. The dH/dl curve for two versions of GROMACS
looks similar for lambdas 0-0.7 (mutation is PhBr -> PhH ) but becomes
different at lambda = 0.9 So something changed there... Moreover different
elecrostatic approaches for PhBr molecule are consistent with each other
(dddG are fine), but new version of GROMACS gives me overestimated ddG as
shown above.

Can anyone give me some advice what to do with it? May be there is some bug
in new version? May be I'm doing something wrong?

Attaching sample .mdp (for solvated and vacuum system) and .top files.

Thanks for the help,
Oleg
-------------- next part --------------
; VARIOUS PREPROCESSING OPTIONS
title                    = Production run

; RUN CONTROL PARAMETERS
integrator               = sd
dt                       = 0.001
nsteps                   = 9500000

; OUTPUT CONTROL OPTIONS
nstxout                  = 10000
nstvout                  = 0
nstfout                  = 0
nstlog                   = 10000
nstenergy                = 100
nstxtcout		 = 10000
nstdhdl			 = 100
xtc_precision            = 2000
xtc-grps                 = Other
energygrps               = Other

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 20
ns-type                  = Grid
rlist                    = 0.95

; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype              = cut-off
rcoulomb                 = 0.95
vdw-type                 = Cut-off
rvdw                     = 0.95
cutoff-scheme		 = group

; Temperature coupling  
tc-grps                  = System 
tau_t                    =  0.2    
ref_t                    = 300   

; Pressure coupling     
;Pcoupl                   = Berendsen
;Pcoupltype               = Isotropic
;tau_p                    = 1.0
;compressibility          = 4.5e-5
;ref_p                    = 1.0

; OPTIONS FOR BONDS   
constraints              = none

; FREE ENERGY OPTIONS
free_energy 		 = yes
init_lambda		 = 0.5
delta_lambda		 = 0	
couple-intramol          = no
sc-alpha		 = 0
-------------- next part --------------
; VARIOUS PREPROCESSING OPTIONS
title                    = Production run

; RUN CONTROL PARAMETERS
integrator               = sd
dt                       = 0.001
nsteps                   = 9500000

; OUTPUT CONTROL OPTIONS
nstxout                  = 10000
nstvout                  = 0
nstfout                  = 0
nstlog                   = 10000
nstenergy                = 100
nstxtcout		 = 10000
nstdhdl			 = 100
xtc_precision            = 2000
xtc-grps                 = Other
energygrps               = System

; NEIGHBORSEARCHING PARAMETERS
nstlist                  = 20
ns-type                  = Grid
pbc                      = xyz
rlist                    = 0.95

; OPTIONS FOR ELECTROSTATICS AND VDW
coulombtype              = pme
rcoulomb                 = 0.95
vdw-type                 = Cut-off
rvdw                     = 0.95
cutoff-scheme		 = group

; Temperature coupling  
tc-grps                  = System 
tau_t                    = 0.2    
ref_t                    = 300   

; Pressure coupling     
Pcoupl                   = berendsen
Pcoupltype               = Isotropic
tau_p                    = 1.0
compressibility          = 4.5e-5
ref_p                    = 1.0


; OPTIONS FOR BONDS   
constraints              = none

; FREE ENERGY OPTIONS
free_energy 		 = yes
init_lambda		 = 0.5
couple-intramol          = no
delta_lambda		 = 0	
sc-alpha		 = 0