How to solve "Too many LINCS warnings"?

GROMACS version: 5.1.2
GROMACS modification: Yes/No
Here post your question
Hi, I am simulating MraY protein which is a protein found in thermophiles. At the step of NVT i am getting error for “Too many LINCS warnings (5691)
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”
Before this, I didn’t get this type of error. How should I solve this error ? I am using POPE lipid and TIP3 water. My system looks fine (noting unusal).
NVT parameters
Input Parameters:
integrator = md
tinit = 0
dt = 0.02
nsteps = 500000
init-step = 0
simulation-part = 1
comm-mode = Linear
nstcomm = 100
bd-fric = 0
ld-seed = 1386368483
emtol = 10
emstep = 0.01
niter = 20
fcstep = 0
nstcgsteep = 1000
nbfgscorr = 10
rtpi = 0.05
nstxout = 500
nstvout = 500
nstfout = 0
nstlog = 500
nstcalcenergy = 100
nstenergy = 500
nstxout-compressed = 0
compressed-x-precision = 1000
cutoff-scheme = Verlet
nstlist = 10
ns-type = Grid
pbc = xyz
periodic-molecules = FALSE
verlet-buffer-tolerance = 0.005
rlist = 1.345
rlistlong = 1.345
nstcalclr = 10
coulombtype = PME
coulomb-modifier = Potential-shift
rcoulomb-switch = 0
rcoulomb = 1.2
epsilon-r = 1
epsilon-rf = inf
vdw-type = Cut-off
vdw-modifier = Potential-shift
rvdw-switch = 0
rvdw = 1.2
DispCorr = EnerPres
table-extension = 1
fourierspacing = 0.16
fourier-nx = 64
fourier-ny = 64
fourier-nz = 72
pme-order = 4
ewald-rtol = 1e-05
ewald-rtol-lj = 0.001
lj-pme-comb-rule = Geometric
ewald-geometry = 0
epsilon-surface = 0
implicit-solvent = No
gb-algorithm = Still
nstgbradii = 1
rgbradii = 1
gb-epsilon-solvent = 80
gb-saltconc = 0
gb-obc-alpha = 1
gb-obc-beta = 0.8
gb-obc-gamma = 4.85
gb-dielectric-offset = 0.009
sa-algorithm = Ace-approximation
sa-surface-tension = 2.05016
tcoupl = V-rescale
nsttcouple = 10
nh-chain-length = 0
print-nose-hoover-chain-variables = FALSE
pcoupl = No
pcoupltype = Isotropic
nstpcouple = -1
tau-p = 1
compressibility (3x3):
compressibility[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
compressibility[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p (3x3):
ref-p[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
ref-p[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
refcoord-scaling = No
posres-com (3):
posres-com[0]= 0.00000e+00
posres-com[1]= 0.00000e+00
posres-com[2]= 0.00000e+00
posres-comB (3):
posres-comB[0]= 0.00000e+00
posres-comB[1]= 0.00000e+00
posres-comB[2]= 0.00000e+00
QMMM = FALSE
QMconstraints = 0
QMMMscheme = 0
MMChargeScaleFactor = 1
qm-opts:
ngQM = 0
constraint-algorithm = Lincs
continuation = FALSE
Shake-SOR = FALSE
shake-tol = 0.0001
lincs-order = 4
lincs-iter = 1
lincs-warnangle = 30
nwall = 0
wall-type = 9-3
wall-r-linpot = -1
wall-atomtype[0] = -1
wall-atomtype[1] = -1
wall-density[0] = 0
wall-density[1] = 0
wall-ewald-zfac = 3
pull = FALSE
rotation = FALSE
interactiveMD = FALSE
disre = No
disre-weighting = Conservative
disre-mixed = FALSE
dr-fc = 1000
dr-tau = 0
nstdisreout = 100
orire-fc = 0
orire-tau = 0
nstorireout = 100
free-energy = no
cos-acceleration = 0
deform (3x3):
deform[ 0]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 1]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
deform[ 2]={ 0.00000e+00, 0.00000e+00, 0.00000e+00}
simulated-tempering = FALSE
E-x:
n = 0
E-xt:
n = 0
E-y:
n = 0
E-yt:
n = 0
E-z:
n = 0
E-zt:
n = 0
swapcoords = no
adress = FALSE
userint1 = 0
userint2 = 0
userint3 = 0
userint4 = 0
userreal1 = 0
userreal2 = 0
userreal3 = 0
userreal4 = 0
grpopts:
nrdf: 81349 118038
ref-t: 323 323
tau-t: 0.1 0.1
annealing: No No
annealing-npoints: 0 0
acc: 0 0 0
nfreeze: N N N
energygrp-flags[ 0]: 0

Using 1 MPI process
Using 1 OpenMP thread

Will do PME sum in reciprocal space for electrostatic interactions.

++++ PLEASE READ AND CITE THE FOLLOWING REFERENCE ++++
U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen
A smooth particle mesh Ewald method
J. Chem. Phys. 103 (1995) pp. 8577-8592
-------- -------- — Thank You — -------- --------

Will do ordinary reciprocal space Ewald sum.
Using a Gaussian width (1/beta) of 0.384195 nm for Ewald
Cut-off’s: NS: 1.345 Coulomb: 1.2 LJ: 1.2
Long Range LJ corr.: 4.2451e-04
System total charge: 0.000
Generated table with 4689 data points for Ewald.
Tabscale = 2000 points/nm
Generated table with 4689 data points for LJ6.
Tabscale = 2000 points/nm
Generated table with 4689 data points for LJ12.
Tabscale = 2000 points/nm
Generated table with 4689 data points for 1-4 COUL.
Tabscale = 2000 points/nm
Generated table with 4689 data points for 1-4 LJ6.
Tabscale = 2000 points/nm
Generated table with 4689 data points for 1-4 LJ12.
Tabscale = 2000 points/nm
Potential shift: LJ r^-12: -1.122e-01 r^-6: -3.349e-01, Ewald -1.000e-05

Thanks in advance !

This time step is about 10x too large to be stable. A dt of this value is only suitable for coarse-grain simulations.

Hi @ Sumedha,

Did you solve your problem?
I have been also facing the same error but still I couldn’t solve it.
I would appreciate any help.

Deniz

Hi @Sumedha, @Deniz and @jalemkul

I am also facing the same error, eventhough my timestep is of 0.001 ps (and is a coarse-grain simulation). The previous energy minimization didn’t achieve the force I set as acceptable but achieved machine precission (I was told to go ahead and try to run the equilibration), but the problems with the LINKS error comes up every time, even once I decreased the timestep.
Could any of you solve it and if so can you give me some advice on how to fix my error? Any help is appreciated.

Thank you very much,

Andreu

.mdp file:
; RUN CONTROL PARAMETERS =
; MARTINI - Most simulations are stable with dt=40 fs,
; some (especially rings) require 20-30 fs.
; The range of time steps used for parametrization
; is 20-40 fs, using smaller time steps is therefore not recommended.

integrator = md
; start time and timestep in ps
tinit = 0.0
dt = 0.001
nsteps = 500000
; number of steps for center of mass motion removal =
nstcomm = 100
comm-grps =

; OUTPUT CONTROL OPTIONS =
; Output frequency for coords (x), velocities (v) and forces (f) =
nstxout = 0
nstvout = 0
nstfout = 0
; Output frequency for energies to log file and energy file =
nstlog = 5000
nstenergy = 5000
; Output frequency and precision for xtc file =
nstxout-compressed = 5000
compressed-x-precision = 100
; This selects the subset of atoms for the xtc file. You can =
; select multiple groups. By default all atoms will be written. =
xtc-grps =
; Selection of energy groups =
energygrps =

; NEIGHBORSEARCHING PARAMETERS =
; MARTINI - no need for more frequent updates
; or larger neighborlist cut-off due
; to the use of shifted potential energy functions.

cutoff-scheme = Verlet
; nblist update frequency =
nstlist = 20
; ns algorithm (simple or grid) =
ns_type = grid
; Periodic boundary conditions: xyz or no =
pbc = xyz
; nblist cut-off =
rlist = 1.1
verlet-buffer-tolerance = 0.01

; OPTIONS FOR ELECTROSTATICS AND VDW =
; MARTINI - vdw and electrostatic interactions are used
; in their shifted forms. Changing to other types of
; electrostatics will affect the general performance of
; the model.

; Method for doing electrostatics =
coulombtype = reaction-field
rcoulomb = 1.1
; Dielectric constant (DC) for cut-off or DC of reaction field =
epsilon_r = 15
epsilon_rf = 0
; Method for doing Van der Waals =
vdw_type = cut-off
; cut-off lengths =
rvdw = 1.1
vdw-modifier = Potential-shift-verlet
; Apply long range dispersion corrections for Energy and Pressure =
DispCorr = No

; OPTIONS FOR WEAK COUPLING ALGORITHMS =
; MARTINI - normal temperature and pressure coupling schemes
; can be used. It is recommended to couple individual groups
; in your system seperately.

; Temperature coupling
tcoupl = v-rescale
tau_t = 1.0 1.0
;tau_t = 1.0
tc-grps = protein non-protein
ref_t = 303 303
;tc-grps = water
;ref_t = 303
Pcoupl = Berendsen
Pcoupltype = isotropic
tau_p = 5.0
compressibility = 3e-4
ref_p = 1.0

; GENERATE VELOCITIES FOR STARTUP RUN =
gen_vel = yes
gen_temp = 303
gen_seed = 1234

; OPTIONS FOR BONDS =
; MARTINI - for ring systems constraints are defined
; which are best handled using Lincs.

constraints = none
; Type of constraint algorithm
constraint_algorithm = Lincs
; Do not constrain the start configuration =
continuation = no
; Highest order in the expansion of the constraint coupling matrix =
lincs_order = 4
; Lincs will write a warning to the stderr if in one step a bond =
; rotates over more degrees than =
lincs_warnangle = 30

.top file:

#include “/home/ubuntu/thesis/martini_v2.2.itp”
#include “/home/ubuntu/thesis/Peptide_assembly_GMX5-2016/3_Done/martini_v2.0_ions.itp”
#include “/home/ubuntu/thesis/ini_files/martini v2.0 solvents.itp”

#include “/home/ubuntu/thesis/Rectangularbox/GFD3E/GFD3E_martini_A.itp”

[ system ]
; name
Martini system from /home/ubuntu/thesis/ini_files/GFD3E_try.pdb

[ molecules ]
; name number
W 10196
EOL 8000
GFD3E_martini_A 30
NA+ 90