GROMACS version: 2020
GROMACS modification: Yes/No
Here post your question
Hello everyone!! I know that LINCS WARNINGS is a largely discussed topic but after reading the past discussions I couldn’t find the same case I am about to present.
I’m preparing a protein-protein system in explicit water, whose residues I have conditioned at two different pH using playmolecule>proteinprepare.
On the conditioned pdb I have run a cleaning.py script that restores all the atom naming in compliance with aminoacids.rtp of charmm36.
I have then successfully pdb2gmx the system into complex.gro and proceeded with solvatation, ions addition, index creation and em.tpr generation using grompp.
Here is the minimization.mdp:
integrator = steep ; Algorithm (steep = steepest descent minimization)
emtol = 1000.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
emstep = 0.01 ; Energy step size
nsteps = 50000 ; Maximum number of (minimization) steps to perform
; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
nstlist = 1 ; Frequency to update the neighbor list and long range forces
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; Method to determine neighbor list (simple, grid)
rlist = 1.4 ; Cut-off for making neighbor list (short range forces)
coulombtype = PME ; Treatment of long range electrostatic interactions
rcoulomb = 1.4 ; Short-range electrostatic cut-off
rvdw = 1.4 ; Short-range Van der Waals cut-off
pbc = xyz ; Periodic Boundary Conditions
Then I run the minimization using “mdrun -v -deffnm em” which correctly produces em.gro, em.edr, em.trr, em.log.
The em.log state a final potential energy of -4.6038165e+06 reached in 2478 steps.
I would now like to proceed with NVT equilibration but when I run I get:
Step 4, time 0.008 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000647, max 0.009665 (between atoms 3435 and 3436)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3432 3433 66.3 0.1112 0.1110 0.1111
3435 3437 50.2 0.1113 0.1109 0.1111
Step 5, time 0.01 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.026490, max 0.457917 (between atoms 10263 and 10264)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
10266 10267 90.0 0.1111 0.1594 0.1111
10263 10264 90.0 0.1112 0.1620 0.1111
Step 5, time 0.01 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.101927, max 1.751611 (between atoms 3438 and 3440)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
3435 3437 90.0 0.1109 0.2365 0.1111
3438 3440 90.0 0.1111 0.3057 0.1111
Wrote pdb files with previous and current coordinates
Step 6, time 0.012 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000792, max 0.016193 (between atoms 10263 and 10265)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
10266 10267 66.3 0.1594 0.1107 0.1111
Step 6, time 0.012 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 21.915888, max 447.443390 (between atoms 3438 and 3440)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
10263 10264 72.9 0.1620 0.1107 0.1111
3435 3436 39.2 0.1133 0.1158 0.1111
3435 3437 90.0 0.2365 0.1108 0.1111
3438 3439 90.0 0.1120 2.7457 0.1111
3438 3440 90.0 0.3057 49.8221 0.1111
3444 3445 90.0 0.1040 0.1220 0.1040
Wrote pdb files with previous and current coordinates
From here on, I searched this forum for answers and I realized that (as in many cases) something wrong occurred in the input structure
There is a major difference in the atoms interested by the LINCS WARNING between the original structure (pdb input for pdb2gmx) and the output of the minimization (em.gro).
pdb
ATOM 620 HB1 GLN H 153 -6.298 3.119 -3.741 1.00 0.00 H
ATOM 621 HB1 GLN H 153 -7.858 3.603 -3.736 1.00 0.00 H
gro
153GLN HB110264 3.312 9.733 6.686
153GLN HB210265 3.312 9.733 6.686
The two aliphatic hydrogen of CB of GLN residue which have two different set of coordinates in the original pdb, became the same atom after minimization. How is that possible? Has this ever happened to anyone before?
Things I have tried so far:
- decreasing the emstep from 0.01 to 0.005 (won’t work)
Thanks all for the help!!!
PS: This is the NVT.mdp I have used for the equilibration
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.002 ; 2 fs
; Output control
nstxout = 500 ; save coordinates every 1.0 ps
nstvout = 500 ; save velocities every 1.0 ps
nstenergy = 500 ; save energies every 1.0 ps
nstlog = 500 ; update log file every 1.0 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; bonds involving H are constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; also related to accuracy
; Nonbonded settings
cutoff-scheme = Verlet ; Buffered neighbor searching
ns_type = grid ; search neighboring grid cells
nstlist = 10 ; 20 fs, largely irrelevant with Verlet
rcoulomb = 1.0 ; short-range electrostatic cutoff (in nm)
rvdw = 1.0 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; account for cut-off vdW scheme
; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
pme_order = 4 ; cubic interpolation
fourierspacing = 0.16 ; grid spacing for FFT
; Temperature coupling is on
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein Non-Protein ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 310 310 ; reference temperature, one for each group, in K
; Pressure coupling is off
pcoupl = no ; no pressure coupling in NVT
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 300 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed