[gmx-users] Problems in putting restraints during free energy perturbation calculations (cont)
Searle Duay
searle.duay at uconn.edu
Wed Mar 28 18:30:54 CEST 2018
Hi,
I would like to follow-up on this:
https://www.mail-archive.com/gromacs.org_gmx-users@maillist.sys.kth.se/msg31701.html.
I already fixed the problem with the coulombic and vdW interactions with
the help of Mark (thank you!). However, my problem now is on turning off
the restraints from states 41 to 60. I am not sure if I am using the
parameters for restraint lambdas correctly. From states 0 to 40, the
restraint lambdas are all 1.00, then I started turning it off with
decrements of 0.05. The restraints that I applied are pull restraints to
keep a zinc ion close to two residues in my peptide. However, when I
analyze the data using gmx bar, I expect a change in free energy, but there
is no change in free energy. I feel that there's something wrong with my
parameters because when I visualized the trajectory using VMD at the final
state (where restraints should be turned off completely), the zinc ion is
still close to the two residues.
Here is the copy of my MDP file:
; Run control
integrator = sd ; Langevin dynamics
tinit = 0
dt = 0.002
nsteps = 5000000 ; 10 ns
nstcomm = 100
; Output control
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstenergy = 500
tc-grps = PROT SOL_ION
tau-t = 1.0 1.0
ref-t = 300 300
; Pressure coupling is on for NPT
pcoupl = Parrinello-Rahman
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 = 0
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
coul-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
vdw-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 = 100
; Do not generate velocities
gen_vel = no
; 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 NPT
continuation = yes
; 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 = 20710.8
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 = 20710.8
Thanks!
--
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