LINCS Error during Equilibration Step

User discussions
1

GROMACS version: Latest as of 9/14/2022
GROMACS modification: Yes/No
Here post your question

Both the mdp file and the error are listed below. I am trying to run a simulation using OPLS. I am using the standard OPLS files for mdp files. What is going on? The errors I receive are generally about the PME electrostatics which are about -.0005 which is insignificant. Not sure what other issues there are with this.

Here is my NPT file

title = NVT equilibration for KALP15-DPPC
define = -DPOSRES ; position restrain the protein
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 50000 ; 2 * 50000 = 100 ps
dt = 0.001 ; 2 fs
; Output control
nstxout = 100 ; save coordinates every 0.2 ps
nstvout = 100 ; save velocities every 0.2 ps
nstenergy = 100 ; save energies every 0.2 ps
nstlog = 100 ; update log file every 0.2 ps
; Bond parameters
continuation = no ; first dynamics run
constraint_algorithm = lincs ; holonomic constraints
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 2 ; accuracy of LINCS
lincs_order = 8 ; also related to accuracy
; Neighborsearching
ns_type = grid ; search neighboring grid cels
nstlist = 5 ; 10 fs
rlist = 1.2 ; short-range neighborlist cutoff (in nm)
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
; 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_GTP Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant, in ps
ref_t = 323 323 ; 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
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = yes ; assign velocities from Maxwell distribution
gen_temp = 323 ; temperature for Maxwell distribution
gen_seed = -1 ; generate a random seed
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm = 1
comm-mode = Linear
comm-grps = Protein_GTP Water_and_ions

Here is the error I am receiving

Step 316, time 0.316 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.031434, max 1.650226 (between atoms 2720 and 2690)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2719 2689 45.5 0.0947 0.0946 0.0945
2720 2690 112.2 0.2258 0.2504 0.0945
2721 2694 54.5 0.0945 0.0946 0.0945
2722 2698 36.6 0.0945 0.0946 0.0945

Back Off! I just backed up step316b.pdb to ./#step316b.pdb.1#

Back Off! I just backed up step316c.pdb to ./#step316c.pdb.1#
Wrote pdb files with previous and current coordinates

Step 317, time 0.317 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.025719, max 1.350201 (between atoms 2720 and 2690)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2719 2689 47.7 0.0946 0.0947 0.0945
2720 2690 115.2 0.2504 0.2221 0.0945
2721 2694 59.0 0.0946 0.0943 0.0945
2722 2698 37.2 0.0946 0.0945 0.0945

Back Off! I just backed up step317b.pdb to ./#step317b.pdb.1#

Back Off! I just backed up step317c.pdb to ./#step317c.pdb.1#
Wrote pdb files with previous and current coordinates

Step 318, time 0.318 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.030322, max 1.591839 (between atoms 2720 and 2690)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2719 2689 45.0 0.0947 0.0941 0.0945
2720 2690 112.7 0.2221 0.2449 0.0945
2721 2694 62.1 0.0943 0.0945 0.0945
2722 2698 39.0 0.0945 0.0944 0.0945

Back Off! I just backed up step318b.pdb to ./#step318b.pdb.1#

Back Off! I just backed up step318c.pdb to ./#step318c.pdb.1#
Wrote pdb files with previous and current coordinates

Step 319, time 0.319 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.028260, max 1.483552 (between atoms 2720 and 2690)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2719 2689 44.8 0.0941 0.0946 0.0945
2720 2690 113.8 0.2449 0.2347 0.0945
2721 2694 61.3 0.0945 0.0948 0.0945
2722 2698 39.9 0.0944 0.0945 0.0945

Back Off! I just backed up step319b.pdb to ./#step319b.pdb.1#

Back Off! I just backed up step319c.pdb to ./#step319c.pdb.1#
Wrote pdb files with previous and current coordinates

Step 320, time 0.32 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.028147, max 1.477654 (between atoms 2720 and 2690)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
2719 2689 45.5 0.0946 0.0942 0.0945
2720 2690 113.8 0.2347 0.2341 0.0945
2721 2694 61.3 0.0948 0.0942 0.0945
2722 2698 39.2 0.0945 0.0945 0.0945

Program: gmx mdrun, version 2020.1-Ubuntu-2020.1-1
Source file: src/gromacs/mdlib/constr.cpp (line 224)

Fatal error:
Too many LINCS warnings (1003)
If you know what you are doing you can adjust the lincs warning threshold in
your mdp file
or set the environment variable GMX_MAXCONSTRWARN to -1,
but normally it is better to fix the problem

2

Don’t be so sure. If you got your ligand topology from LigParGen, this is a well-known issue. The topology can, in fact, sum to a non-integer value that should be corrected.

Beware using my GROMOS membrane inputs for non-membrane OPLS simulations. The cutoffs, temperature, constraints, etc. are wrong for the kind of simulation you’re doing.

If the simulation is blowing up, there are lots of troubleshooting tips here (even a specific example of a protein-ligand complex!): https://manual.gromacs.org/current/user-guide/terminology.html#system-diagnosis

3

How do I correct the error for ligpargen ligand?

4

Grompp will tell you what your overall system charge is. Justin will probably disagree with me on this one, but if it’s +/-0.01 charge units then I am of a mind to neglect the issue. However, I doubt the overall system charge has anything to do with your system crashing, which is also why Justin suggested you look at systems blowing up (with a link you should check out). A system with a net charge of even +/-N for any value of N may not be physically realistic but it should be able to run without crashing even with PME (counterintuitively, I suppose). I am not sure what you mean about errors being with PME, all I see are LINCS warning, which implies general system blowing up, usually.

My general advice when people have systems blowing up is (a) check if you made a huge mistake, and then (b) re-run EM with define=-DFLEXIBLE. I’ve never understood why, but this can help to get your system over barriers that seem to get stuck on EM with rigid waters. If it makes you feel better, you can always EM again with settle on before your MD, but I’ve never seen the point in doing so.

5

Yes, I disagree very strongly :) Errors of that magnitude are always suggestive of an underlying problem. I suspect the OP’s issue is not related to the -0.0005 e charge, but it should still be addressed since LigParGen is known to add up partial charges incorrectly in some cases.

Here, I will point out that if a large force is acting on a water molecule (or between a water and something else), oftentimes the O-H bond is the easiest thing to perturb and it will extend a long way. If you then try to run dynamics with rigid water, the constraint algorithm immediately fails. If you want to run with rigid water, it is advisable to re-minimize with constraints on to resolve the issue before dynamics fails. I’ve had this happen to me a couple of times, to the point that I never attempt flexible water during EM because it seems to be more harm than good, and a better solution may exist. But mileage clearly varies on this approach.

Without seeing the topology, it’s hard to say. One easy way to do this is simply to redistribute +0.0001 e on five different atoms, as long as the resulting charges are still sensible and reasonably close to similar species in the force field. You may also be able to just stick +0.0005 e onto one atom and call it done. The magnitude of the issue is quite small, so you won’t really go wrong either way.