[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======
; Production MD run for free energy calculation
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
Graduate Research Assistant
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul at vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/
========================================
More information about the gromacs.org_gmx-users
mailing list