[gmx-users] Problems in putting restraints during free energy perturbation calculations
Mark Abraham
mark.j.abraham at gmail.com
Mon Mar 12 17:32:18 CET 2018
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.
>
More information about the gromacs.org_gmx-users
mailing list