Error while using position restarints in MD run

User discussions
,
1

GROMACS version: 2023.1
GROMACS modification: Yes/No
Here post your question

Hello all,
I am trying to simulate a DNA-membrane system where I want to apply position restraints only on DNA in the production run. But while using the command:
gmx mdrun -v -deffnm md_0_1,
it is showing LINCS error and the run is stopped. Here is the error:
Step 876, time 0.876 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 31216132.000000, max 3863014656.000000 (between atoms 88473 and 88474)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
39161 39162 112.1 0.1111 0.2950 0.1111
39164 39165 91.2 4.6748 5.1557 0.1111
39164 39166 103.2 0.3457 0.4882 0.1111
39167 39168 92.2 5.0048 2.9605 0.1111
39732 39733 73.3 0.4217 0.1111 0.1111
39734 39735 110.0 0.1011 0.3246 0.1111
39734 39736 98.5 0.3566 0.1376 0.1111
39734 39737 54.0 0.1172 0.1173 0.1111
85119 85121 98.0 0.1111 0.8011 0.1111
85134 85135 36.9 0.1115 0.1153 0.1111
85134 85136 39.3 0.1116 0.1149 0.1111
85134 85137 100.9 0.1104 0.5852 0.1111
85144 85145 40.8 0.1111 0.1111 0.1111
85178 85180 94.0 0.9634 1.5993 0.1111
85181 85182 91.9 0.1187 3.4375 0.1111
85181 85183 123.8 0.6689 0.1997 0.1111
85187 85188 111.7 0.2891 0.3004 0.1111
85189 85190 95.0 0.2718 1.2723 0.1111
85189 85192 101.0 0.5227 0.5831 0.1111
85193 85194 99.5 0.1111 0.6733 0.1111
85193 85195 97.9 0.1111 86.4671 0.1111
85199 85200 52.0 0.5286 0.1013 0.1111
85199 85201 94.2 0.1503 1.5152 0.1111
87736 87737 92.0 0.1111 3.2267 0.1111
87736 87738 90.9 0.1111 7.2435 0.1111
87739 87740 122.3 0.1111 429180832.0000 0.1111
87739 87741 77.0 0.1111 429180864.0000 0.1111
87773 87774 65.4 0.1129 0.1141 0.1111
87773 87775 83.4 0.6898 0.1087 0.1111
87776 87777 102.1 0.2135 701.3235 0.1111
87776 87778 100.0 0.1114 0.6405 0.1111
87779 87780 91.3 0.1116 4.7411 0.1111
87779 87781 95.1 0.3253 1.2405 0.1111
88324 88325 161.3 0.2327 5031.8369 0.1111
88324 88326 144.7 100.1726 4813.1401 0.1111
88324 88327 77.7 0.3314 5032.0215 0.1111
88464 88466 91.7 0.5818 3.6829 0.1111
88469 88470 35.3 1.0136 8601994.0000 0.1111
88473 88474 82.3 143.8495 429180928.0000 0.1111
88473 88475 63.0 0.9098 429180768.0000 0.1111
88476 88477 114.8 0.1111 0.2646 0.1111
88476 88478 102.1 0.1111 40.5930 0.1111
88482 88483 98.8 5.9134 0.1294 0.1111
89192 89193 115.1 0.1111 0.2624 0.1111
89198 89199 39.3 0.1111 48.9639 0.1111
89198 89200 73.6 0.1111 42.9575 0.1111
89207 89208 90.9 0.3581 6.7295 0.1111
89209 89210 96.1 0.1111 0.1240 0.1111
89209 89211 31.3 0.1111 0.1121 0.1111
89209 89212 100.9 0.1111 0.1337 0.1111
89259 89261 92.0 0.1111 3.1863 0.1111
89262 89263 84.2 0.1111 721879.1875 0.1111
89392 89393 54.2 0.2772 0.1111 0.1111
91854 91855 110.1 1.4639 0.3228 0.1111
91854 91856 115.0 0.2031 0.2625 0.1111
91857 91858 44.2 0.1144 231.9760 0.1111
91859 91860 105.3 0.1111 0.1538 0.1111
91859 91861 119.6 0.1111 0.1853 0.1111
91859 91862 92.0 0.1111 3.1109 0.1111
91863 91864 101.1 1.2794 0.5751 0.1111
91866 91867 93.5 0.1376 1.8082 0.1111
91866 91868 91.9 1.6835 3.2789 0.1111
91869 91870 92.5 0.4026 2.5030 0.1111
91869 91871 99.0 1.3699 0.7099 0.1111
91872 91873 90.8 0.1111 7.8312 0.1111
91874 91875 92.1 0.2672 2.9756 0.1111
91874 91876 99.5 0.2225 0.6703 0.1111
91874 91877 99.9 0.2612 0.6491 0.1111
91878 91880 92.9 0.5104 2.2284 0.1111
91878 91881 93.3 0.3811 1.9263 0.1111
93246 93247 90.7 0.3242 9.6808 0.1111
93246 93248 93.6 0.7115 1.7602 0.1111
93252 93253 91.3 2.6884 0.1137 0.1111
93261 93262 92.2 0.1111 2.8810 0.1111
93261 93263 44.2 0.1111 0.1026 0.1111
93269 93272 93.9 0.1111 1.6159 0.1111
96219 96220 99.1 0.1111 0.7003 0.1111
96219 96221 91.5 0.1111 4.2736 0.1111
96284 96286 98.0 0.1111 5222.8994 0.1111
96284 96287 100.0 0.1111 5223.0332 0.1111
97388 97391 103.1 0.1111 0.4917 0.1111
99889 99890 64.4 3.0742 13.3486 0.1111
99889 99891 98.3 0.4694 0.7683 0.1111
99904 99905 91.4 0.2280 4.4763 0.1111
99904 99906 108.6 0.6482 0.3491 0.1111
99964 99965 62.5 0.1986 0.1111 0.1111
99971 99973 112.9 0.2093 0.2860 0.1111

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

Back Off! I just backed up step876c.pdb to ./#step876c.pdb.1#
^CWrote pdb files with previous and current coordinates
^C^CAborted (core dumped)

I am also attaching the topol.top and md.mdp file which I am using.

topol.top
File ‘topol.top’ was generated
; By user: nabanita (1001)
; On host: kashi
; At date: Mon Apr 15 13:47:07 2024
;
; This is a standalone topology file
;
; Created by:
; :-) GROMACS - gmx pdb2gmx, 2023.3 (-:
;
; Executable: /usr/local/gromacs/bin/gmx
; Data prefix: /usr/local/gromacs
; Working dir: /data/nabanita/13th_demo/hj_test/15thapr_demo
; Command line:
; gmx pdb2gmx -f new.pdb -o dna_mem.gro -ter
; Force field was read from current directory or a relative path - path added.
;

; Include forcefield parameters
#include “./charmm36-jul2022.ff/forcefield.itp”

; Include chain topologies
#include “topol_DNA_chain_D.itp”
#include “topol_DNA_chain_D2.itp”
#include “topol_DNA_chain_D3.itp”
#include “topol_DNA_chain_D4.itp”
#include “topol_Other_chain_P.itp”
#include “topol_Other_chain_P2.itp”
#include “topol_Other_chain_P3.itp”
#include “topol_Other_chain_P4.itp”
#include “topol_Other_chain_P5.itp”
#include “topol_Other_chain_P6.itp”

#ifdef POSRES_DNAP
#include “posre_dnaP.itp”
#endif

; Include water topology
#include “./charmm36-jul2022.ff/tip3p.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

; Include topology for ions
#include “./charmm36-jul2022.ff/ions.itp”

[ system ]
; Name
Protein in water

[ molecules ]
; Compound #mols
DNA_chain_D 1
DNA_chain_D2 1
DNA_chain_D3 1
DNA_chain_D4 1
Other_chain_P 1
Other_chain_P2 1
Other_chain_P3 1
Other_chain_P4 1
Other_chain_P5 1
Other_chain_P6 1
SOL 87531
CL 47
MG 47
CL 47
NA 155
CL 63

md.mdp
title = KALP15-DPPC Production MD
define = -DPOSRES -DPOSRES_DNAP ; position restarint the protein
; Run parameters
integrator = sd ; leap-frog integrator
nsteps = 1000000 ; 1 * 1000000 = 1000 ps (1 ns)
dt = 0.001 ; 1 fs
; Output control
nstxout = 100 ; save coordinates every 2 ps
nstvout = 100 ; save velocities every 2 ps
nstxtcout = 100 ; xtc compressed trajectory output every 2 ps
nstenergy = 100 ; save energies every 2 ps
nstlog = 100 ; update log file every 2 ps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 2 ; accuracy of LINCS
lincs_order = 4 ; 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 ; More accurate thermostat
tc-grps = DNA_DFPE Water_and_ions ; two coupling groups - more accurate
tau_t = 0.5 0.5 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm = 100
comm-mode = Linear
comm-grps = DNA_DFPE Water_and_ions
; Scale COM of reference coordinates
refcoord_scaling = com

2

It seems like you need to run longer equilibration with more restraints before turning them off.

There are also a few things that spring to my mind when I see your parameters, but they are not related to the crash.

  • Why are you using the sd integrator for this simulation?
    • If you want to use the sd integrator, why are you also specifying the v-rescale thermostat?
    • If you want to use the sd integrator I would suggest tau-t of 2.
  • Why are you writing your output files so often, every 100 fs? Do you really need 10,000 frames? Do you really need uncompressed output at all? nstenergy = 10000 should be enough.
3

Thank you @MagnusL for your valuable suggestion.
Regarding your questions:

  1. There are atoms at both ends of an angle, connected by constraints and
    with masses that differ by more than a factor of 13. This means that
    there are likely dynamic modes that are only very weakly coupled. To
    ensure good equipartitioning, you need to either not use constraints on
    all bonds (but, if possible, only on bonds involving hydrogens) or use
    integrator = sd or decrease one or more tolerances:
    verlet-buffer-tolerance <= 0.0001, LINCS iterations >= 2, LINCS order >=
    4 or SHAKE tolerance <= 1e-05

  2. The bond in molecule-type Protein_chain_A between atoms 173 O1 and 174 HO1 has an estimated oscillational period of 9.0e-03 ps, which is less than 10 times the time step of 1.0e-03 ps. Maybe you forgot to change the constraints mdp option”. For that reason I am using 1 fs timestep.

4

Regarding the second question:
The bond in molecule-type Protein_chain_A between atoms 79 C2 and 80 O2 has an estimated oscillational period of 1.9.0e-03 ps, which is less than 10 times the time step of 2.0e-03 ps. Maybe you forgot to change the constraints mdp option”. For that reason I am using 1 fs timestep. For that reason I am using 1 fs timestep.

5

Thanks for the explanation. I think you misunderstood my questions a little. But I’ll follow-up on your answers.

  1. You are currently only constraining bonds to hydrogens - and that is usually the best choice. But still, if you want to use the sd integrator it is confusing to set parameters for the v-rescale thermostat. A tau-t of 2 is better for sd.

  2. I didn’t question your choice of timestep. I’m wondering about the output intervals. Why are nstxout, nstvout etc so low? And do you really need uncompressed data at all?

6

Regrading the answer 2
I don’t need output for such small intervals. I kept it because I assumed that the crash is happening for keeping large nstxout, nstvout etc.

7

Regarding the timestep, if I keep 2 fs instead of 1 fs will it pose a problem alongwith the restraints.

8

And also if we use sd, should we not set any parameters for v-rescale? Should we off the pressure coupling

Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)

9

The things is that when you are using the sd integrator that takes over as a thermostat, and uses tau-t as the inverse friction constant. When running gmx grompp you also get a note that says that v-rescale is inactivated (or something along those lines) when you are using the sd integrator. So, it’s less confusing to remove v-rescale from the mdp file. I’d suggest something like:

tcoupl = no ; handled by the sd integrator
tc-grps = DNA_DFPE Water_and_ions ; two coupling groups - more accurate
tau_t = 2 2 ; inverse friction for the sd integrator
ref_t = 300 300 ; reference temperature, one for each group, in K

or even

tcoupl = no ; handled by the sd integrator
tc-grps = system
tau_t = 2 ; inverse friction for the sd integrator
ref_t = 300 ; reference temperature, one for each group, in K

The pressure coupling is not affected by the sd integrator. So you can keep your current settings.

10

Thank you so much for your suggestions. I have one more question, can I keep a time of 2 fs despite the error;

“The bond in molecule-type Protein_chain_A between atoms 79 C2 and 80 O2 has an estimated oscillational period of 1.9.0e-03 ps, which is less than 10 times the time step of 2.0e-03 ps. Maybe you forgot to change the constraints mdp option”. For that reason I am using 1 fs timestep. For that reason I am using 1 fs timestep.”

11

I would not use a 2 fs time step when you get that warning.

But I’m a bit surprised that you get such a long oscillation period for a bond in a protein. What force field is this?

12

I am using Charmm36 force-field and sorry for not changing the name, I am simulating DNA with membrane

13

I’m not used to simulating DNA, but I’m surprised if the Charmm36 force field would not be able to run with a 2 fs time step when you have constrained the bonds to hydrogens. I assume you haven’t modified any of the force field parameters, right?

14

No I have not modified any of the force-field parameters. But from nvt step, this error is coming up

15

Sorry it’s a note which pops up

16

Your settings are very close to the threshold for getting a note (and with reasonable margin from getting a warning. From src/gromacs/gmxpreprocess/grompp.cpp

The stability limit of leap-frog or velocity verlet is 4.44 steps per oscillational period.
But accurate bonds distributions are lost far before that limit.
To allow relatively common schemes (although not common with Gromacs)
of dt=1 fs without constraints and dt=2 fs with only H-bond constraints
we set the note limit to 10.

You are currently at 9 steps per oscillational period. You could test if it works OK with a time step of 2 ns. Run a few nanoseconds in an NVE ensemble and check the energy conservation. If the energy and the system itself seems stable, it should be OK.

17

Hello,
I want to restraint only the DNA in my case. This is the md.mdp file which I am using. But after running for 100 ns I don’t see any bending in the membrane as well. Does it mean my lipids are also restrained.
title = KALP15-DPPC Production MD
define = -DPOSRES -DPOSRES_DNA
; Run parameters
integrator = md ; leap-frog integrator
nsteps = 500000000 ; 1 * 500000000 = 500000 ps (500 ns)
dt = 0.001 ; 1 fs
; Output control
nstxout = 1000000 ; save coordinates every 2 ps
nstvout = 1000000 ; save velocities every 2 ps
nstxtcout = 1000000 ; xtc compressed trajectory output every 2 ps
nstenergy = 1000000 ; save energies every 2 ps
nstlog = 1000000 ; update log file every 2 ps
; Bond parameters
continuation = yes ; Restarting after NPT
constraint_algorithm = lincs ; holonomic constraints
constraints = h-bonds ; all bonds (even heavy atom-H bonds) constrained
lincs_iter = 1 ; accuracy of LINCS
lincs_order = 4 ; 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 = Nose-Hoover ; More accurate thermostat
tc-grps = DNA_DFPE Water_and_ions ; two coupling groups - more accurate
tau_t = 0.5 0.5 ; time constant, in ps
ref_t = 300 300 ; reference temperature, one for each group, in K
; Pressure coupling is on
pcoupl = C-rescale ; Pressure coupling on in NPT
pcoupltype = semiisotropic ; uniform scaling of x-y box vectors, independent z
tau_p = 2.0 ; time constant, in ps
ref_p = 1.0 1.0 ; reference pressure, x-y, z (in bar)
compressibility = 4.5e-5 4.5e-5 ; isothermal compressibility, bar^-1
; Periodic boundary conditions
pbc = xyz ; 3-D PBC
; Dispersion correction
DispCorr = EnerPres ; account for cut-off vdW scheme
; Velocity generation
gen_vel = no ; Velocity generation is off
; COM motion removal
; These options remove motion of the protein/bilayer relative to the solvent/ions
nstcomm = 100
comm-mode = Linear
comm-grps = DNA_DFPE Water_and_ions
; Scale COM of reference coordinates
refcoord_scaling = all

18

If you look for #ifdef POSRES in your topology (and all files that are included from it) you will see which molecules/atoms are restrained by that define statement.

Without knowing your force field or topology, I would presume that -DPOSRES would restrain proteins and possibly lipids, whereas -DPOSRES_DNA would restrain your DNA molecules. If you only want to restrain the DNA you can almost certainly remove -DPOSRES.

19

@MagnusL thank you so much for your suggetions. I am using Charmm 36 force- field and below mentioned is the topol.top file.

File ‘topol.top’ was generated
; By user: nabanita (1001)
; On host: kashi
; At date: Wed Apr 17 10:32:08 2024
;
; This is a standalone topology file
;
; Created by:
; :-) GROMACS - gmx pdb2gmx, 2023.3 (-:
;
; Executable: /usr/local/gromacs/bin/gmx
; Data prefix: /usr/local/gromacs
; Working dir: /data/nabanita/13th_demo/hj_test/17th_apr
; Command line:
; gmx pdb2gmx -f new.pdb -o dna_mem.gro -ter
; Force field was read from current directory or a relative path - path added.
;

; Include forcefield parameters
#include “./charmm36-jul2022.ff/forcefield.itp”

; Include chain topologies
#include “topol_DNA_chain_D.itp”
#include “topol_DNA_chain_D2.itp”
#include “topol_DNA_chain_D3.itp”
#include “topol_DNA_chain_D4.itp”
#include “topol_Other_chain_P.itp”
#include “topol_Other_chain_P2.itp”
#include “topol_Other_chain_P3.itp”
#include “topol_Other_chain_P4.itp”
#include “topol_Other_chain_P5.itp”
#include “topol_Other_chain_P6.itp”

; Include position restrain file
#ifdef POSRES_DNA
#include “posre_dna.itp”
#endif

; Include water topology
#include “./charmm36-jul2022.ff/tip3p.itp”

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
; i funct fcx fcy fcz
1 1 1000 1000 1000
#endif

; Include topology for ions
#include “./charmm36-jul2022.ff/ions.itp”

[ system ]
; Name
Protein in water

[ molecules ]
; Compound #mols
DNA_chain_D 1
DNA_chain_D2 1
DNA_chain_D3 1
DNA_chain_D4 1
Other_chain_P 1
Other_chain_P2 1
Other_chain_P3 1
Other_chain_P4 1
Other_chain_P5 1
Other_chain_P6 1
SOL 87531
CL 47
MG 47
CL 47
NA 155

20

Just based on the contents you copied, -DPOSRES does not have any effect at all. However, I guess there might be sections similar to:

#ifdef POSRES
#include "posre_topol_Other_chain_P.itp"
#endif

in the included .itp files.

Still, if you only want to restrain the DNA you can definitely set: define = -DPOSRES_DNA in the mdp file.