# [gmx-users] Free Energy Calculations

Justin A. Lemkul jalemkul at vt.edu
Sat Feb 16 23:11:17 CET 2008

```Hello all,

I'll preface this message by saying that this is my first foray into free energy
calculations, and I am hoping that someone more well-versed in the matter may be
able to help me sort out some strange results.  I have done a good bit of
reading in terms of tutorials and in the archives, so I decided to give my
system a try.  I am attempting to determine DeltaG of solvation for a
polyphenolic compound.  My approach is to turn off charges in one step, and to
use soft-core to turn off Lennard-Jones interactions in a second step, using
thermodynamic integration to determine my DeltaG values.  Bonded parameters
remain untouched.  I follow closely the procedure posted by David Mobley at the
DillGroup Wiki site, changing only a few parameters for my simulation, as I am
using Gromos96, not OPLS-AA.

Relevant snapshots of my topologies and .mdp files follow this message.  The
problem I am having (well, at least I think it's a problem), is that I get
different results for the forward and reverse simulations.  I define "forward"
as turning off the charges/LJ parameters, and "reverse" as turning said
parameters back on.

For example, for an in vacuo simulation:

dG(forward, LJ component) = -113.276 kJ/mol
dG(forward, Coulombic component) = 37.9936 kJ/mol
dG(forward, total) = -75.282 kJ/mol

dG(reverse, LJ component) = 112.695 kJ/mol
dG(reverse, Coul. component) = -65.3451 kJ/mol
dG(reverse, total) = 47.350

The Lennard-Jones components seem to come out in good agreement, but the charge
component shows a substantial difference.

Thanks in advance for any ideas.

-Justin

======topol_forward-QQ.top======
[ atoms ]
;   nr      type  resnr resid  atom  cgnr   charge     mass     typeB chargeB
massB
1       CR1     1  MYR     CAI     1   -0.100  12.0110     CR1     0.000
12.0110
2        HC     1  MYR     HAI     1    0.100   1.0080      HC     0.000
1.0080
3         C     1  MYR     CAN     2    0.150  12.0110       C     0.000
12.0110
4        OA     1  MYR     OAC     2   -0.548  15.9994      OA     0.000
15.9994
5         H     1  MYR     HAC     2    0.398   1.0080       H     0.000
1.0080

======topol_reverse-QQ.top======
[ atoms ]
;   nr      type  resnr resid  atom  cgnr   charge     mass     typeB chargeB
massB
1       CR1     1  MYR     CAI     1    0.000  12.0110     CR1    -0.100
12.0110
2        HC     1  MYR     HAI     1    0.000   1.0080      HC     0.100
1.0080
3         C     1  MYR     CAN     2    0.000  12.0110       C     0.150
12.0110
4        OA     1  MYR     OAC     2    0.000  15.9994      OA    -0.548
15.9994
5         H     1  MYR     HAC     2    0.000   1.0080       H     0.398
1.0080

======topol_forward-LJ.top=======
[ atoms ]
;   nr      type  resnr resid  atom  cgnr   charge     mass     typeB chargeB
massB
1       CR1     1  MYR     CAI     1    0.000  12.0110     DUM     0.000
12.0110
2        HC     1  MYR     HAI     1    0.000   1.0080     DUM     0.000
1.0080
3         C     1  MYR     CAN     2    0.000  12.0110     DUM     0.000
12.0110
4        OA     1  MYR     OAC     2    0.000  15.9994     DUM     0.000
15.9994
5         H     1  MYR     HAC     2    0.000   1.0080     DUM     0.000
1.0080

======topol_reverse-LJ.top======
[ atoms ]
;   nr      type  resnr resid  atom  cgnr   charge     mass     typeB chargeB
massB
1       DUM     1  MYR     CAI     1    0.000  12.0110     CR1     0.000
12.0110
2       DUM     1  MYR     HAI     1    0.000   1.0080      HC     0.000
1.0080
3       DUM     1  MYR     CAN     2    0.000  12.0110       C     0.000
12.0110
4       DUM     1  MYR     OAC     2    0.000  15.9994      OA     0.000
15.9994
5       DUM     1  MYR     HAC     2    0.000   1.0080       H     0.000
1.0080

======md-QQ.mdp======
title           = 5-ns free energy MD run
cpp             = /usr/bin/cpp
include         = -I/home/jalemkul/Amyloid_B/MD/Free_Energy/Topologies/
; Run control parameters
integrator      = sd
; start time and timestep in ps
tinit           = 0
dt              = 0.002
nsteps          = 2500000
; number of steps for center of mass motion removal
nstcomm         = 100
; Output control options
; Output frequency for coords (x), velocities (v) and forces (f)
nstxout         = 500
nstvout         = 500
nstfout         = 0
; Output frequency for energies to log file and energy file
nstlog          = 500
nstenergy       = 500
energygrps      = MYR
; Output frequency and precision for xtc file
nstxtcout       = 0
xtc-precision   = 1000
; Neighborsearching parameters
nstlist         = 5
ns_type         = grid
; Periodic boundary conditions: xyz or none
pbc             = xyz
; nblist cut-off
rlist           = 1.0
domain-decomposition     = no
; Options for electrostatics and vdw
; Method for doing electrostatics
coulombtype     = pme
rcoulomb        = 1.0
; Dielectric constant (DC) for cut-off or DC of reaction field
epsilon-r       = 1
; vdw cut-off length
vdwtype         = cut-off
rvdw            = 1.4
; Apply long range dispersion corrections for Energy and Pressure
DispCorr        = no
; Spacing for the PME/PPPM FFT grid
fourierspacing  = 0.12
; FFT grid size, when a value is 0 fourierspacing will be used
fourier_nx      = 0
fourier_ny      = 0
fourier_nz      = 0
; EWALD/PME/PPPM parameters
pme_order       = 6
ewald_rtol      = 1e-06
epsilon_surface = 0
optimize_fft    = no
; Restraints
dihre           = simple
dihre-fc        = 1
nstdihreout     = 1000
disre           = simple
disre_fc        = 1
; Options for temperature coupling
tc_grps         = system
tau_t           = 0.1
ref_t           = 298
; Options for pressure coupling
Pcoupl          = Berendsen
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
sc_alpha        = 0
sc_power        = 1.0
sc_sigma        = 0
; Generate velocities for the startup run
gen_vel         = no
; OPTIONS FOR BONDS
constraints     = hbonds
; Type of constraint algorithm
constraint-algorithm     = Lincs
; Do not constrain the start configuration
unconstrained-start      = no
; Relative tolerance of shake
shake-tol       = 0.0001
; Highest order in the expansion of the constraint coupling matrix
lincs-order     = 12
; Lincs will write a warning to the stderr if in one step a bond
; rotates over more degrees than
lincs-warnangle = 30

======md-LJ.mdp======
(identical to the above, except):
; Free energy control stuff
free_energy     = yes
init_lambda     = 0.0
delta_lambda    = 0
sc_alpha        = 0.5
sc_power        = 1.0
sc_sigma        = 0.3

========================================

Justin A. Lemkul