Licns warning

User discussions
1

GROMACS version: 2022.4
GROMACS modification: Yes/No
I am trying to determine the interaction between a protein and a ligand. Both the protein and ligand are using the gromos54a7 force field.

When I simulated the protein and ligand together in water, after energy minimization, equilibration(nvt, npt) steps, I get LICNS WARNING in MD step. And I get residual files like step2935029b.pdb, step2935029c.pdb … and the simulation does not terminate completely.

When I simulate the protein and ligand separately in water with the same .mdp option at all steps(ions, em, nvt, npt,…), I get the same error only for the ligand.

This is the error message I saw when I simulated the protein and ligand together. :
starting mdrun ‘Protein in water’
5530000 steps, 11060.0 ps (continuing from step 530000, 1060.0 ps).

Step 2935010, time 5870.02 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.002271, max 0.034268 (between atoms 1110 and 1111)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1103 1109 32.7 0.1390 0.1414 0.1390
1109 1110 35.0 0.1400 0.1438 0.1400

Step 2935020, time 5870.04 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000301, max 0.007636 (between atoms 1111 and 1112)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1111 1112 32.4 0.1097 0.1098 0.1090

Step 2935021, time 5870.04 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000224, max 0.004744 (between atoms 1103 and 1109)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1111 1112 32.0 0.1098 0.1091 0.1090

Step 2935023, time 5870.05 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.001704, max 0.055889 (between atoms 1111 and 1112)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1111 1112 90.0 0.1087 0.1151 0.1090
1119 1120 37.3 0.1085 0.1105 0.1090

Step 2935024, time 5870.05 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000448, max 0.011007 (between atoms 1119 and 1120)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1111 1112 65.2 0.1151 0.1083 0.1090
1119 1120 39.7 0.1105 0.1102 0.1090

Step 2935025, time 5870.05 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000402, max 0.009328 (between atoms 1111 and 1112)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
1111 1112 61.3 0.1083 0.1080 0.1090
1119 1120 43.2 0.1102 0.1093 0.1090

And this is .mpd file that I used :

; Run parameters
integrator = md ; leap-frog integrator
dt = 0.002 ; 2 fs
nsteps = 5000000 ; 2*5000000 = 10 ns
init-step = 530000 ; the final step of the previous running

; Output control
nstlog = 5000 ; update log file every 10 ps
nstenergy = 0 ; save energies every 10 ps → no save (5000->0)
nstxout-compressed = 5000 ; save .xtc file

; Neighborsearching
cutoff-scheme = Verlet ; Don’t set rlist, verlet-buffer-tolerance
nstlist = 20 ; will be automatically increased by mdrun, usually to 20 or more
ns_type = grid ; search neighboring grid cells

; Periodic boundary conditions
pbc = xyz ; 3-D PBC

; Electrostatics
coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics
rcoulomb = 1.2 ; short-range electrostatic cutoff (in nm)

; Van der Waals
vdw-type = Cut-off
rvdw = 1.2 ; short-range van der Waals cutoff (in nm)
DispCorr = EnerPres ; apply long range dispersion corrections for Energy and Pressure

; Temperature coupling
tcoupl = V-rescale ; modified Berendsen thermostat
tc-grps = Protein_G19K Water_and_ions ; two coupling groups - more accurate
tau_t = 0.1 0.1 ; time constant for coupling in ps
ref_t = 300 300 ; reference temperature for coupling, one for each group, in K

; Pressure coupling
pcoupl = Parrinello-Rahman ; NPT for production MD
pcoupltype = isotropic ; uniform scaling of box for soluble protein simulation, not for membrane protein
tau_p = 2.0 ; coupling constant in ps
compressibility = 4.5e-5 ; isothermal compressibility of water, bar^-1
ref_p = 1.0 ; reference pressure in bar
refcoord_scaling = com

; Velocity generation
gen_vel = no ; continuous running
ld-seed = 54321 ; seed for random number generation for langevan dynamics

; Bond parameters
constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained
constraint_algorithm = lincs ; holonomic constraints
continuation = yes ; continuous running

Thanks for reading. I’m not sure how to ask as I’m a student who hasn’t had much exposure to MD simulation yet… Please let me know if I’m falling short in what I have to offer.

2

From your simulation of the protein and ligand separately, it sounds like the ligand might be the problem. Have you tried energy minimizing the ligand in vacuo (as suggested here https://manual.gromacs.org/current/user-guide/terminology.html#system-diagnosis) to see if the topology is capable of achieving a stable configuration? If it can’t then the problem could be the ligand topology itself.

Could you describe how you got the topology info for the ligand? Also, were there any notes/warnings that were output into the terminal during any of the previous steps before MD?

3

When I simulated the ligand in vacuum, I still got the same error.

I obtained the topology of the ligand from the ATB (https://atb.uq.edu.au/index.py) site. For other molecules obtained from the above site, until now, the simulation completes without error.

Below are the notes/warnings messages I received.

WARNING 1 [file topol.top, line 3574]:
The GROMOS force fields have been parametrized with a physically
incorrect multiple-time-stepping scheme for a twin-range cut-off. When
used with a single-range cut-off (or a correct Trotter
multiple-time-stepping scheme), physical properties, such as the density,
might differ from the intended values. Since there are researchers
actively working on validating GROMOS with modern integrators we have not
yet removed the GROMOS force fields, but you should be aware of these
issues and check if molecules in your system are affected before
proceeding. Further information is available at
GROMACS can not reproduce properties with the GROMOS force fields - Redmine #2884 (#2884) · Issues · GROMACS / GROMACS · GitLab, and a longer
explanation of our decision to remove physically incorrect algorithms can
be found at http s://doi.org/10.2 6434/chem rxiv.11474583.v1 .

NOTE 1 [file nvt.mdp]:
Removing center of mass motion in the presence of position restraints
might cause artifacts. When you are using position restraints to
equilibrate a macro-molecule, the artifacts are usually negligible.

NOTE 1 [file npt.mdp]:
Removing center of mass motion in the presence of position restraints
might cause artifacts. When you are using position restraints to
equilibrate a macro-molecule, the artifacts are usually negligible.

4

Have you tried using a different forcefield? I’m not sure if this issue is actually to do with how GROMOS is parametrized as raised in Warning 1, but it might be useful to see if having the Ligand parameters generated in a different forcefield could solve the issue

It might also help to visualize the structure to see what might be going wrong

5

As it turns out, I had immobilized the ligand in system and proceeded with the EM step. When I unfastened it and simulated it, the problem was solved.

Thank you.

1 Like