[gmx-users] TI calculation in vacuum gives number which is not zero for small molecule.

[Itamar Kass]

  Shalom all,

I wish to calculate the free energy of solvation of acetic acid using 
thermodynamics integration. So I built the topology:

; Include forcefield parameters
#include "ffG53a6.itp"

[ moleculetype ]
; Name            nrexcl
Protein             3

[ atoms ]
;   nr       type  resnr residue  atom   cgnr     charge       mass  
typeB    chargeB      massB
      1        CH3      1   ASPH     CB      2          0     15.035   
DUM        0        15.035
      2          C      1   ASPH     CG      3       0.33     12.011   
DUM        0        12.011
      3          O      1   ASPH    OD1      3      -0.45    15.9994   
DUM        0        15.9994
      4         OA      1   ASPH    OD2      3     -0.288    15.9994   
DUM        0        15.9994
      5          H      1   ASPH    HD2      3      0.408      1.008   
DUM        0        1.008
[ bonds ]
;  ai    aj funct            c0            c1            c2            c3
     1     2     2    gb_27
     2     3     2    gb_5
     2     4     2    gb_13
     4     5     2    gb_1

[ pairs ]
;  ai    aj funct            c0            c1            c2            c3
     1     5     1

[ angles ]
;  ai    aj    ak funct            c0            c1            
c2            c3
     1     2     3     2    ga_30
     1     2     4     2    ga_19
     3     2     4     2    ga_33
     2     4     5     2    ga_12

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            
c2            c3            c4            c5
     1     2     4     5     1    gd_12

[ dihedrals ]
;  ai    aj    ak    al funct            c0            c1            
c2            c3
     1     4     3     2     2    gi_1

and put the molecule in an empty box at the size of 5.32 nm^3 and 
simulated using the below parameters:

integrator                =  sd
ld_seed                  =  -1
dt                        =  0.002    ; ps !
nsteps                    =  550000    ; total 1.1ns.

; OPTIONS FOR BONDS
; Constrain control
constraints             =  all-bonds
; Do not constrain the start configuration
continuation         = no
; Type of constraint algorithm
constraint-algorithm     =  lincs

; OUTPUT CONTROL OPTIONS
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout                   =  10000
nstvout                   =  10000
nstfout                   =  10000
; Output frequency and precision for xtc file
nstxtcout                 =  500
xtc-precision             =  1000
; Energy monitoring
energygrps                =  protein
nstenergy             =  500

; NEIGHBORSEARCHING PARAMETERS
; nblist update frequency
nstlist                   =  5
; ns algorithm (simple or grid)
ns-type                   =  Grid
; nblist cut-off
rlist                     =  0.8

; OPTIONS FOR ELECTROSTATICS AND VDW
; Method for doing electrostatics
coulombtype              = Generalized-Reaction-Field
rcoulomb                 = 1.4
epsilon_rf               = 62
; Method for doing Van der Waals
vdw-type                 = Cut-off
; cut-off lengths
rvdw-switch              = 0
rvdw                     = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr                 = no

; OPTIONS FOR WEAK COUPLING ALGORITHMS
; Temperature coupling
tcoupl                   = no
; Groups to couple separately, time constant (ps) and reference 
temperature (K)
tc-grps                  = system
tau-t                    = 0.2
ref-t                    = 298
; Pressure coupling
Pcoupl                   = no
; Generate velocites is on at 290K - do not get velocity from gro file.
gen_vel             =  yes
gen_temp            =  290
gen-seed            =  -1

; Free energy control stuff
free-energy              = yes
init-lambda              = 0 ; Topology A (lambda=0) to topology B 
(lambda=1).
delta-lambda             = 0
sc-alpha                 = 0.5 ; soft core potantial
sc-power                 = 1
sc-sigma                 = 0.3

; Center of mass control
nstcomm              =  1000
; Periodic boundary conditions
pbc                  =  xyz
; Mode for center of mass motion removal
comm_mode            =  Linear
; Groups for center of mass motion removal
comm-grps            =  system

The FE value out of the in vacuum simulation should be 0, as I don't 
have any interactions larger then 1-4, but what I get is -0.0956224 
kJ/mol. Some more technical data - I have done 40 (delta_lambda = 
0.025), for each lambda value I had calculated the dg/dl for the last 
1ns (out of 1.1ns).  Any idea why this happen?

All the best,
Itamar

-- 


"In theory, there is no difference between theory and practice. But, in practice, there is." - Jan L.A. van de Snepscheut

===========================================
| Itamar Kass, Ph.D.
| Postdoctoral Research Fellow
|
| Department of Biochemistry and Molecular Biology
| Building 77 Clayton Campus
| Wellington Road
| Monash University,
| Victoria 3800
| Australia
|
| Tel: +61 3 9902 9376
| Fax: +61 3 9902 9500
| E-mail: Itamar.Kass at monash.edu
============================================