[gmx-users] Free energy of VDW transformation from state A to State B

[deepak bapat]

Dear gmx-users

I am trying to calculate free energy of VDW transformation. The system is
very simple as I want to test the methodology. Single neutral LJ particle
is transformed from state A to state B in TIP4P water.

I have tried simulations with different state B LJ parameters without
changing state A LJ parameters.

I find that at init_lambda = 0.0 dhbydl values are dependent on state B
parameters. And these values vary considerably depending on state B
parameters.

e.g.

1) transformation 1: (TIP4P OW like LJ)

9e-02 & 6.50194e-02        TO      3.15365e-01 & 6.48520e-01 : <dh/dlambda>
= 5061 +- 95 kJ/mol

2) transformation 2: (Neon like LJ)

9e-02 & 6.50194e-02        TO      2.78000e-01 & 2.88696e-01 : <dh/dlambda>
= 1632 +- 30 kJ/mol

3) transformation 3: (All atom methane like LJ)

9e-02 & 6.50194e-02      TO     3.73000e-01 & 1.23010e+00 : <dh/dlambda> =
19062 +- 417 kJ/mol

Is this behavior normal if initial state is similar in all three
simulations at init_lambda = 0.0

Shouldn't results be equivalent irrespective state B LJ parameters ? Or am
I doing something wrong?

topology file reads something like this


;begin top file

#include "forcefield.itp"

#include "nonbonded.itp"

#include "tip4p.itp"

#include "TPI.itp"

[ system ]

TPI in water

[ molecules ]

SOL 907

TPI 1

; end of top file


I have forcefield.itp, nonbonded.itp, tip4p.itp and TPI.itp in my working
folder.

TPI.itp is as follows


;begin TPI file

[ moleculetype ]

; molname nrexcl

TPI 1

[ atoms ]

; id at type res nr residu name at name cg nr charge mass

1 opls_test 1 TPI OW 2 0.0 20.17970 opls_test_transform 0.0 20.17970

; end TPI file


and nonbonded.itp is shown below


;begin nonbonded

[ atomtypes ]

; TPI

opls_test OW 10 20.17970 0.000 A 9e-02 6.50194e-02 ; state A

opls_test_transform OW 10 20.17970 0.000 A 3.73000e-01 1.23010e+00 ;
transformed to state B

;TIP4P

opls_113 OW 8 15.99940 0.000 A 3.15365e-01 6.48520e-01

opls_114 HW 1 1.00800 0.520 A 0.00000e+00 0.00000e+00

opls_115 MW 0 0.00000 -1.040 D 0.00000e+00 0.00000e+00

;end nonbonded


some relevant parameters in md.mdp file are

;begin mdp

integrator = md

tinit = 0

dt = 0.001

nsteps = 1000000

nstcomm = 100

nstxout = 0

nstvout = 00

nstfout = 0

nstlog = 5000

nstenergy = 5000

nstxtcout = 5000

cutoff-scheme = verlet

nstlist = 10

ns_type = grid

pbc = xyz

rlist = 1.0

coulombtype = PME

rcoulomb = 1.0

vdwtype = cutoff

rvdw = 1.0

DispCorr = EnerPres

fourierspacing = 0.12

pme_order = 6

epsilon_surface = 0

tc_grps = system

tcoupl = nose-hoover

tau_t = 0.5

ref_t = 300

Pcoupl = Parrinello-Rahman

tau_p = 1.0

compressibility = 4.5e-05

ref_p = 1.0

; Free energy variables

free_energy = yes

couple-moltype =

couple-lambda0 = vdw-q

couple-lambda1 = vdw-q

couple-intramol = no

init_lambda = 0.0

delta_lambda = 0

nstdhdl = 10

dhdl-print-energy = no

sc-alpha = 0

sc-power = 1

sc-r-power = 6

sc-sigma = 0.3

sc-coul = no

separate-dhdl-file = yes

dhdl-derivatives = yes

dh_hist_size = 0

dh_hist_spacing = 0.1

; end mdp file


Regards

-- 
Deepak