[gmx-users] LINCS warnings for a single molecule in otherwise stable system

[Cara Kreck]

Hi everyone,

I am trying to run a simulation of a lipid bilayer with a multi-component solution. One of the components is ethylene glycol. Some of the components have a tendency to wriggle their way into the lipids. When an EG molecule does this, it tends to get uncomfortable and trigger some LINCS warnings. It only happens with a single EG molecule at a time. It never happens with any of the other components or the EG molecules in the solution. It never leads to anything else becoming unstable or the system blowing up, it just causes a crash when the warning limit is reached. The strangest thing is when the simulation is continued from the checkpoint, the problem goes away for a while but eventually reoccurs within a few hours with a completely different EG molecule wedged in the membrane.

I have tried increasing both lincs-order and lincs-iter, but it still happens. I have tried increasing lincs-warnangle, but that just lets the angle fluctuations get bigger before the job is killed. I tried using a slightly different topology for EG and I thought that had helped for a while. The first time the angle fluctuations started, they stayed around the 30 degrees mark and went away by themselves after around 600 warnings, so the job was then able to run for the full 24 hours allowed. The last 2 jobs I have submitted have been killed prematurely by LINCS warnings that didn't go away by themselves though. This is very frustrating as I am running on a very busy cluster and it can take 24-48 hours for a job to progress through the queue, only to then be killed within a couple of hours of running. The only thing I can think of now is to increase lincs-warnangle by a small amount so that the smaller fluctuations are given a chance to sort themselves out before being killed.

If anyone else could give me any suggestions, I would really appreciate it. I am using GROMACS 4.6.7 and the GROMOS53A6 FF. Here is a copy of my mdp file:

integrator               = md
dt                       = 0.002 ; 2fs
nsteps                   = 500000000 ; 1000 ns
comm_grps                = DOPC !DOPC
nstxout                  = 1000
nstvout                  = 0
nstlog                   = 1000
nstenergy                = 1000
energygrps               = DOPC !DOPC
nstcalcenergy            = 5
nstlist                  = 5
ns_type                  = grid
pbc                      = xyz
rlist                    = 0.8 
coulombtype              = Reaction-Field
rcoulomb                 = 1.4
epsilon_rf               = 62
vdwtype                  = Cut-off
rvdw                     = 1.4 
tcoupl                   = berendsen 
tc-grps                  = DOPC !DOPC
tau_t                    = 0.1 0.1
ref_t                    = 303 303
Pcoupl                   = berendsen
pcoupltype               = semiisotropic
tau_p                    = 1.0 1.0
compressibility          = 4.6e-5 4.6e-5
ref_p                    = 1.0 1.0
gen_vel                  = no
gen_temp                 = 150
gen_seed                 = 173529
constraints              = all-bonds
constraint_algorithm     = lincs
lincs-order              = 6
lincs-iter               = 2
continuation             = yes

Thanks,

Cara