[gmx-users] Kinetic Energy in FEP calculations

[Connie Wang]

Hello,

 I'm attempting to calculate a relative free energy of hydration between
Argon and Neon. However when perturbing from Argon(aq)->Neon(aq) I was
surprised to compute a change in free energy of 1.2kJ/mol, where the
literature values suggest that the hydration energy of Argon is ~-9kJ/mol
and neon is ~-8kJ/mol, which means that I should have seen a change in free
energy of ~-1kJ/mol.

 Interestingly, I noticed that this 2kJ/mol difference could be resolved if
I subtract the dEkin/dlambda value obtained from the edr file. A clumsy
look into the source code seemed to confirm that there is a dEkin/dlambda
term included in the dH/dlambda calculation which results from the change
in mass during the mutation.  However it seems in the literature that most
people calculate the free energy between two states as the difference of
the configurational partition functions, and thus should be a function of
the positions and independent of the mass.

It seems tempting to me to resolve my problem by simply subtracting off the
dEkin/dlambda term, since I am in the NPT ensemble (and at constant
temperature the Ekin shouldn't change). Does this seem reasonable? Is there
a good reason for including this dEkin/dlambda term in FEP calculations? If
I'm interested only in the difference between the configurational partition
functions in the NPT ensemble, should the dEkin/dlambda term be included in
the first place?



In case it helps, my topology file and mdp file:

*Topology File*
; Argon parameters from Maitland, G. C.; Rigby, M.; Smith, E. B.; Wakeham,
W. A.
;Intermolecular Forces: Their Origin and Determination;

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               2               yes             1.0    1.0

[ atomtypes ]
  Ar    Ar      39.948     0.000        A     3.41000e-01   9.96000e-01
  Ne    Ne      20.180     0.000        A     2.72000e-01   3.91000e-01
  DNe   DNe     20.180     0.000        A     0.00000e-01   0.00000e-01
  DAr   DAr     39.948     0.000        A     0.00000e-01   0.00000e-01

; Include water topology
[ atomtypes ]
OW      OW     16.00   0.0000  A   3.15061e-01  6.36386e-01
HW      HW     1.008   0.0000  A   0.00000e+00  0.00000e+00
[ moleculetype ]
; molname       nrexcl
SOL             1

[ atoms ]
; id  at type     res nr  res name  at name  cg nr  charge    mass
  1   OW          1       SOL       OW       1      -0.834    16.00000
  2   HW          1       SOL       HW1      1       0.417     1.00800
  3   HW          1       SOL       HW2      1       0.417     1.00800

[ settles ]
; i     j       funct   length
1       1       0.09572 0.15139

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


[ moleculetype ]
; Name            nrexcl
 Argon           1

[ atoms ]
; residue   3 LEU rtp CLEU q -1.0
;   nr       type   resnr residue atom   cgnr   charge       mass    typeB
chargeB massB
    1         Ar      1    Ar      Ar     1     0.0000     39.948    Ne
 0.0000  20.180

[ system ]
; Name
Argon Mutate Neon

[ molecules ]
; Compound        #mols
Argon               1
SOL                215

*MDP File*
; Run control
integrator               = sd       ; Langevin dynamics
dt                       = 0.002
nsteps                   = 2500000  ; 5 ns

; Neighborsearching and short-range nonbonded interactions
nstlist                  = 10
ns_type                  = grid
pbc                      = xyz
rlist                    = 0.9
; Electrostatics
coulombtype              = PME
rcoulomb                 = 0.9

; van der Waals
vdw-type                 = switch
rvdw-switch              = 0.7
rvdw                     = 0.8

; Apply long range dispersion corrections for Energy and Pressure
DispCorr                  = EnerPres

; Spacing for the PME/PPPM FFT grid
fourierspacing           = 0.1
; EWALD/PME/PPPM parameters
pme_order                = 6
ewald_rtol               = 1e-06
epsilon_surface          = 0
optimize_fft             = no

; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tc_grps                  = system
tau_t                    = 1.0
ref_t                    = 298

; Pressure coupling is on for NPT
Pcoupl                   = Parrinello-Rahman
tau_p                    = 0.5
compressibility          = 4.5e-05
ref_p                    = 1.0

; Free energy control stuff
free_energy              = yes
init_lambda              = 0.0
delta_lambda             = 0
foreign_lambda           = 0.05
sc-alpha                 = 0.5
sc-power                 = 1.0
sc-sigma                 = 0.3
; options for bonds
constraints              = h-bonds  ; we only have C-H bonds here
; Type of constraint algorithm
constraint-algorithm     = lincs
lincs-order              = 12

I used 10 windows for the perturbation and ran each window for 3ns. The
results were then combined using the g_bar utility

Thanks!
-Connie
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://maillist.sys.kth.se/pipermail/gromacs.org_gmx-users/attachments/20111101/38031bdc/attachment.html>