[gmx-users] problem with free energy calculations

[Michael Brunsteiner]

Dear all,

I seem to have a problem with my free energy TI calculations 
using Gromacs, I'd be grateful for any suggestions about what
could have gone wrong there...

I calculated the free energy of solvation of water
(basically its chemical potential). The literature value for TIP3P 
water, which I used, is about -6.3 kcal/mol - 
I get something like +2.5 (!) ...

my calculations are fairly standard, I use PME with 
conservative parameters, soft core potentials,
do 21 windows, each one lasting 500 pico seconds, then 
after discarding the first 100 ps of each dgdl*xvg file
i use the averages and do a numerical integration, 
the works ...
I also use two thermostates (one for the water that disappears,
one for the rest)

I see that the fluctuations of my dU/dlambda values are fairly
large ... can it be that the soft core parameters i use (sc_alpha=0.7,
sc_power=1, sc_sigma=0.3) are inappropriate if, as i do here, 
i turn on/off the VdW and Coulomb interactions simultanioulsy ??

mdp plus topology files and the results are included below,

thanks in advance for any suggestions!

cheers,
michael


================ md.mdp ==================================

title                    = water mu
cpp                      = /usr/bin/cpp
;
integrator               = md
nsteps                   = 250000
dt                       = 0.002
;
nstxtcout                = 20
nstenergy                = 20
nst_log                  = 20
nstxout                  = 10000
nstvout                  = 10000
;
nstlist                  = 12
nstype                   = grid
pbc                      = xyz
nstcomm                  = 100
rlist                    = 1.2
;
coulombtype              = PME
rcoulomb                 = 1.2
vdwtype                  = shift
rvdw                     = 1.2
;
Tcoupl                   = nose-hoover
tc_grps                  = w1 rest
ref_t                    = 300.0 300.0
tau_t                    = 0.1 0.1
gen_vel                  = yes
;
Pcoupl                   = no
;
pme_order                = 4
ewald_rtol               = 1e-05
ewald_geometry           = 3d
epsilon_surface          = 0
optimize_fft             = no
fourier_nx               = 64
fourier_ny               = 64
fourier_nz               = 64
;
constraints              = hbonds
unconstrained-start      = no
lincs_iter               = 4
;
free_energy              = yes    
init_lambda              = 0.00 ; (0.05, 0.10, etc in the other runs)
delta_lambda             = 0.0    
sc_alpha                 = 0.7    
sc_power                 = 1    
sc_sigma                 = 0.3    


=== topol.top =====================================================

; Include forcefield parameters
#include "ffamber03.itp"

; a TIP3P water with zero charge and no VdW interactions in the B
; topology the dummy atom types dumow and dumhw are defined in the
; ff*.atp and ff*nb.itp files used here                                        
                       

[ moleculetype ]
; molname       nrexcl
wdummy          1

[ atoms ]
;   nr  type   resnr residue atom cgnr  charge  mass typeB  chargeB   massB
 1  amber99_42   1  SOL     OW  1   0.000   16.00000  dumow  0.00   16.00000
 2  amber99_27   1  SOL    HW1  1   0.000    1.00800  dumhw  0.00    1.00800
 3  amber99_27   1  SOL    HW2  1   0.000    1.00800  dumhw  0.00    1.00800

[ settles ]
; OW    funct   doh     dhh
1       1       0.09572 0.15139

[ exclusions ]
1       2       3
2       1       3
3       1       2


; the other waters 

; Include water topology
#include "ffamber_tip3p.itp"

; the system 

[ system ]
; Name
462 tip3p water

[ molecules ]
; Compound        #mols
wdummy              1
SOL               461


=== the results (averages from dgdl*xvg files) ===

# lamda <dU/dl>   rmsd
0.00  -4.85       14.04
0.05  -5.44       14.42
0.10  -5.66       14.40
0.15  -6.77       14.93
0.20  -7.98       17.11
0.25  -9.17       17.85
0.30 -11.09       19.52
0.35 -14.23       21.80
0.40 -18.69       27.28
0.45 -27.41       29.40
0.50 -27.75       22.81
0.55 -26.55       17.26
0.60 -20.99       11.85
0.65 -15.47        8.66
0.70 -10.12        6.01
0.75  -6.54        4.53
0.80  -2.75        3.30
0.85  -0.02        2.50
0.90   2.37        1.93
0.95   3.78        1.45
1.00   5.07        1.15





 
____________________________________________________________________________________
Get your own web address.  
Have a HUGE year through Yahoo! Small Business.
http://smallbusiness.yahoo.com/domains/?p=BESTDEAL