[gmx-users] Problems in putting restraints during free energy perturbation calculations
Searle Duay
searle.duay at uconn.edu
Mon Mar 12 17:41:04 CET 2018
Hi Mark,
Thanks for reminding me about this. Yes, I've read up on this on a tutorial
and I probably interchanged which interactions should be turned on first.
Thank you again. Have a wonderful day!
Searle
On Mon, Mar 12, 2018 at 12:32 PM, Mark Abraham <mark.j.abraham at gmail.com>
wrote:
> Hi,
>
> Don't do that. You have charges on atoms that have no VDW, so there is
> nothing to stop them flying all over the place and imparting crazy forces
> on other atoms that have VDW, etc. Turn the VDW fully on before adding any
> charge. This is a fairly well known issue but something that people
> regularly get wrong, too. (Is there maybe a tutorial you did that missed
> giving a pointer on this issue?)
>
> Mark
>
> On Mon, Mar 12, 2018 at 5:21 PM Searle Duay <searle.duay at uconn.edu> wrote:
>
> > 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
> > --
> > Gromacs Users mailing list
> >
> > * Please search the archive at
> > http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> > posting!
> >
> > * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
> >
> > * For (un)subscribe requests visit
> > https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> > send a mail to gmx-users-request at gromacs.org.
> >
> --
> Gromacs Users mailing list
>
> * Please search the archive at http://www.gromacs.org/
> Support/Mailing_Lists/GMX-Users_List before posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>
--
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