Too many LINCS Warnings in MARTINI Lipid Bilayer simulations

User discussions
1

GROMACS version:2022
GROMACS modification: No

Dear GROMACS users,

I have a coarse grained lipid bilayer(MARTINI 2.2), containing DPPC, DOPC and CHOL. I have managed to minimize the system in two steps, first to a tolerance of 1000 kJ/nm/mol followed by another step with a tolerance of 500 kJ/nm/mol. Both minimizations were done using the steepest descent algorithm with a step size of 0.015 and 0.005 respectively. Now while performing equilibriation runs, I am getting the following errors.

Fatal error:
Too many LINCS warnings (1214)
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

The angle deviations shown in the log file are 179.8-180 degrees, but the maximum threshold seems to be 30 degrees.

Can any one tell me how I can resolve this? Thank you for patience.

2

I assume that you are already applying restraints during EM and the crashing equilibration stage. You could see if it helps running a first (short) equilibration simulation with a shorter integration time step. Try dt = 0.001 (yes even with Martini).

3

Dear MagnusL,

Thank you for your input. I actually have multiple equilibration steps, starting with an integration time step of dt = 0.002, and gradually increasing it all the way to dt = 0.02
These are some example mdp files

;define                   = -DBILAYER_LIPIDHEAD_FC=200
integrator               = md
tinit                    = 0.0
dt                       = 0.002
nsteps                   = 500000

nstlog                   = 1000
nstenergy                = 1000
nstxout-compressed       = 1000
compressed-x-precision   = 100

cutoff-scheme            = Verlet
nstlist                  = 20
ns_type                  = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005

epsilon_r                = 15
coulombtype              = reaction-field
rcoulomb                 = 2
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 2

tcoupl                   = v-rescale
tc-grps                  = membrane W
tau_t                    = 1.0  1.0
ref_t                    = 303.15 303.15

; Pressure coupling:
Pcoupl                   = berendsen
Pcoupltype               = semiisotropic
tau_p                    = 5.0
compressibility          = 3e-4 3e-4
ref_p                    = 1.0  1.0

; GENERATE VELOCITIES FOR STARTUP RUN:
gen_vel                  = yes
gen_temp                 = 303.15
gen_seed                 = 2009092072

refcoord_scaling         = all


define                   = -DBILAYER_LIPIDHEAD_FC=10
integrator               = md
tinit                    = 0.0
dt                       = 0.020
nsteps                   = 50000

nstlog                   = 1000
nstenergy                = 1000
nstxout-compressed       = 1000
compressed-x-precision   = 100

cutoff-scheme            = Verlet
nstlist                  = 20
ns_type                  = grid
pbc                      = xyz
verlet-buffer-tolerance  = 0.005

epsilon_r                = 15
coulombtype              = reaction-field
rcoulomb                 = 2
vdw_type                 = cutoff
vdw-modifier             = Potential-shift-verlet
rvdw                     = 2

tcoupl                   = v-rescale
tc-grps                  = membrane W
tau_t                    = 1.0  1.0
ref_t                    = 303.15 303.15

; Pressure coupling:
Pcoupl                   = berendsen
Pcoupltype               = semiisotropic
tau_p                    = 5.0
compressibility          = 3e-4 3e-4
ref_p                    = 1.0  1.0

; GENERATE VELOCITIES FOR STARTUP RUN:
gen_vel                  = no

refcoord_scaling         = all

P.S, I have commented the restrains on the lipid heads as my system is generated using PACKMOL. Some of the packing defects were not closing in earlier attempts, but closed up after removing restrains on the lipid heads

4

P.S,
Further minimization runs were also performed and the system was minimized to a tolerance of 200 kJ/mol/nm

5

So, did the simulations crash in the first or the second of those two example equilibration runs? Have you tried restraining all of the lipids?

It seems like PACKMOL does not give very stable configurations since many people try to work-around using different methods. So, if you apply restraints, make sure that you do not always use the PACKMOL coordinates as reference. It is better to use the output from the previous simulation instead. But that might also mean that you need more simulations, e.g.:

  1. EM, restrain whole lipids
  2. Equilibration NVT, dt = 0.001, restrain whole lipids
  3. Equilibration NVT, dt = 0.002, restrain whole lipids
  4. Equilibraiton NVT, dt = 0.004, restrain whole lipids
  5. Equilibration NVT, dt = 0.005, restrain lipid heads
  6. Equilibration NPT, dt = 0.010, restrain lipid heads
  7. Equilibration NPT, dt = 0.020, restrain lipid heads
  8. Equilibration NPT, dt = 0.020, no restraints

If it turns out that the simulations crash with dt = 0.020 you might have to settle for 0.015 or even 0.010.

6

Dear MagnusL,

The LINCS error comes up in the very first equilibriation run using a time step of dt = 0.002. The example mdp files that I had posted in my previous reply was for the first and last equilibration steps respectively. I have 6 such steps, not including minimization, with the following parameters

  1. Minimization 1, steepest descent with a tolerance of 1000 kJ/mol/nm
  2. Minimization 2, steepest descent with a tolerance of 500 kJ/mol/nm
  3. Minimization 3, steepest descent with a tolerance of 300 kJ/mol/nm
  4. Minimization 1, steepest descent with a tolerance of 200 kJ/mol/nm
  5. Equilibration 1, NPT, dt = 0.002, no restrains
  6. Equilibration 2, NPT, dt = 0.005, no restrains
  7. Equilibration 3, NPT, dt = 0.010, no restrains
  8. Equilibration 4, NPT, dt = 0.015, no restrains
  9. Equilibration 5, NPT, dt = 0.020, no restrains
  10. Production Run, NPT, dt = 0.020, no restrains

I have managed to minimize the system to below a tolerance of 200 kJ/mol/nm, but in the very first equilibration run, I am getting LINCS errors. I will try out your procedure as well. Thank you for your inputs

7

I have attempted to run the equilibration runs as specified by you, but with not luck. I am still receiving LINCS errors in the very first step. I also tried reducing the value of dt to dt = 0.0005 and dt = 0.0001 with no luck. Do you think running further steps of energy minimization could help?

8

I think one stage of EM with the default tolerance is almost always enough (5000 to 50000 steps depending on the system). How many EM steps are performed before your EM simulations finish? I remember some cases when EM finished prematurely when run on GPU. Are your simulations running on GPU? If so, does it make a difference if you run with -nb cpu?

9

Dear MagnusL,

These are the minimization outputs for the different systems I have. The minimization runs were done on a CPU machine and the equilibration runs are being carried out on GPU.

Screenshot from 2024-09-17 15-05-05
Screenshot from 2024-09-17 15-05-13
Screenshot from 2024-09-17 15-05-17
Screenshot from 2024-09-17 15-05-25

10

I have received LINCS errors for all systems, bar the 5_LPE one which has already ran for about 880ns with dt = 0.02.

11

I must admit I don’t know what could be the cause of these problems. Four rounds of EM should definitely be (more than?) enough. It might be worth switching to one temperature coupling group, but I don’t think that will make any significant difference.

If NVT equilibration with restraints and dt = 0.0001 does not work after EM there must be something strange going on.

12

Does it help running the equilibrations on CPU only?

13

Or does it work in GROMACS 2024.3?

14

I tried doing this as well. It still throws the same LINCS error. I am unsure of what to try at this point. During my minimization runs, the first minimization was done as a soft-core minimization, and all the following steps were not. Do you think this could be causing a problem?

15

Can you generate the starting configurations without using PACKMOL?

16

No, I cannot. I have an initial phase-separated system which requires packing of different sets of lipids in different regions of the membrane.

17

Double_Raft
Double_Raft1850×913 280 KB

I have very different compositions in the middle as compared to the rest of the membrane. To my knowledge, I cannot build such a system on any other platform than PACKMOL

18

Then I’m afraid I can’t come up with any other suggestions right now.

19

Or one more thing, what radius are you using for the Martini beads in PACKMOL?

20

I have not set any radii… I have taken the already martinized PDBs from the CHARMM-GUI MARTINI Maker and used the APL values given there to construct the system on PACKMOL