[gmx-users] Problems in putting restraints during free energy perturbation calculations
Searle Duay
searle.duay at uconn.edu
Mon Mar 12 17:20:57 CET 2018
Good day!
I am doing some free energy perturbation to calculate the free energy of
binding of an ion to a peptide. For one of my processes, I am slowly
turning on the coulombic interactions by an increment of 0.05 lambda,
followed by turning on the van der Waals interactions by an increment of
0.05 lambda. This is being done while a pull restraint is present to keep
the ion bound to the peptide. Following is a copy of the MDP file of my
equilibration run under NPT ensemble:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.0005
nsteps = 200000 ; 100 ps
nstcomm = 200
; Output control
nstxout = 1000
nstvout = 1000
nstfout = 0
nstlog = 1000
nstenergy = 1000
nstxout-compressed = 0
; Neighborsearching and short-range nonbonded interactions
cutoff-scheme = Verlet
nstlist = 20
ns-type = grid
pbc = xyz
rlist = 1.2
; Electrostatics
coulombtype = pme
rcoulomb = 1.2
; van der Waals
vdwtype = Cut-off
vdw-modifier = Force-switch
rvdw-switch = 1.0
rvdw = 1.2
; Apply long range dispersion corrections for Energy and Pressure
DispCorr = EnerPres
; Spacing for the PME/PPPM FFT grid
fourierspacing = 0.12
; EWALD/PME/PPPM parameters
pme-order = 6
ewald-rtol = 1e-06
epsilon-surface = 0
; Temperature coupling
; tcoupl is implicitly handled by the sd integrator
tcoupl = berendsen
tc-grps = PROT SOL_ION
tau-t = 1.0 1.0
ref-t = 300 300
; Pressure coupling is on for NPT
pcoupl = berendsen
pcoupltype = isotropic
tau-p = 5.0
compressibility = 4.5e-5 4.5e-5
ref-p = 1.0 1.0
; Free energy control stuff
free-energy = yes
init-lambda-state = 1
delta-lambda = 0
calc-lambda-neighbors = 1 ; only immediate neighboring windows
; Vectors of lambda specified here
; Each combination is an index that is retrieved from init_lambda_state for
each simulation
; init_lambda_state 0 1 2 3 4 5 6 7 8 9
10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
25 26 27 28 29 30 31 32 33 34 35 36 37 38 39
40 41 42 43 44 45 46 47 48 49 50 51 52 53 54
55 56 57 58 59 60
vdw-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.05 0.10 0.15
0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90
0.95 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00
coul-lambdas = 0.00 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40
0.45 0.50 0.55 0.60 0.65 0.70 0.75 0.80 0.85 0.90 0.95 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00
; We are not transforming any bonded or restrained interactions
bonded-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00
restraint-lambdas = 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00
1.00 1.00 0.95 0.90 0.85 0.80 0.75 0.70 0.65 0.60 0.55 0.50 0.45 0.40 0.35
0.30 0.25 0.20 0.15 0.10 0.05 0.00
; Masses are not changing (particle identities are the same at lambda = 0
and lambda = 1)
mass-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Not doing simulated temperting here
temperature-lambdas = 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00
; Options for the decoupling
sc-alpha = 0.5
sc-coul = no ; linear interpolation of Coulomb (none
in this case)
sc-power = 1
sc-sigma = 0.3
couple-moltype = HETA ; name of moleculetype to decouple
couple-lambda0 = none ; only van der Waals interactions
couple-lambda1 = vdw-q ; turn off everything, in this case
only vdW
couple-intramol = no
nstdhdl = 10
; Do not generate velocities
gen-vel = yes
gen-temp = 300
gen-seed = -1
; options for bonds
constraints = h-bonds ; we only have C-H bonds here
; Type of constraint algorithm
;constraint-algorithm = LINCS
; Constrain the starting configuration
; since we are continuing from NVT
continuation = no
; Highest order in the expansion of the constraint coupling matrix
lincs-order = 12
;
pull = yes
pull-nstxout = 0
pull-nstfout = 0
pull-ngroups = 3
pull-ncoords = 2
pull-group1-name = ZN
pull-group2-name = HIS_17
pull-group3-name = HIS_21
pull-coord1-type = umbrella
pull-coord1-geometry = distance
pull-coord1-groups = 2 1
pull-coord1-dim = Y Y Y
pull-coord1-init = 0.2095
pull-coord1-start = no
pull-coord1-k = 1000
pull-coord2-type = umbrella
pull-coord2-geometry = distance
pull-coord2-groups = 3 1
pull-coord2-dim = Y Y Y
pull-coord2-init = 0.2095
pull-coord2-start = no
pull-coord2-k = 1000
During NPT equilibriation, I am getting the following error:
Step 6556, time 3.278 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 12685.671875, max 173473.953125 (between atoms 276 and 277)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
276 277 83.6 0.3435 17347.4961 0.1000
Back Off! I just backed up step6556b.pdb to ./#step6556b.pdb.1#
Back Off! I just backed up step6556c.pdb to ./#step6556c.pdb.1#
Wrote pdb files with previous and current coordinates
WARNING: Listed nonbonded interaction between particles 277 and 280
at distance 17347.604 which is larger than the table limit 2.200 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
WARNING: Listed nonbonded interaction between particles 273 and 277
at distance 17348.008 which is larger than the table limit 2.200 nm.
This is likely either a 1,4 interaction, or a listed interaction inside
a smaller molecule you are decoupling during a free energy calculation.
Since interactions at distances beyond the table cannot be computed,
they are skipped until they are inside the table limit again. You will
only see this message once, even if it occurs for several interactions.
IMPORTANT: This should not happen in a stable simulation, so there is
probably something wrong with your system. Only change the table-extension
distance in the mdp file if you are really sure that is the reason.
1. Am I right that I have to give a value of 1.00 for the restraint-lambda
to have it turned on?
2. I've been doing this before without the restraints on and everything
went smooth, but when I try to have the restraints turned on, I am getting
this problem. I am not sure if I am doing this correctly and how to address
the problem.
Thank you!
--
Searle Aichelle S. Duay
Ph.D. Student
Chemistry Department, University of Connecticut
searle.duay at uconn.edu
More information about the gromacs.org_gmx-users
mailing list