[gmx-users] LINCS Warnings in the middle of an MD run

[K.A.Feenstra]

Christoph Meier wrote:
> 
> Dear all,
> 
> One of my GROMACS MD simulations crashed after running stably for almost 4ns (=several
> weeks on our system). The error message I get is this:
> 
> "...
> LINCS WARNING
> relative constraint deviation after LINCS:
> max 0.000987 (between atoms 3286 and 3289) rms 0.000026
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
> ..."
> 
> I got almost 100 of these LINCS WARNINGs before the system eventually crashed.
> The atoms involved (3286 and 3289) belong to a protonated lysine residue, the topology
> of which was generated in the standard (pdb2gmx) way.
> Does anybody have suggestions as to what I could do to get round this problem ?
> Do you think there is anything wrong with the lysine residue or anything else ?
> 
> I would be grateful for any advice you could give me. Thanks very much !

[.. after reading David's comment, I felt I could add something more...]

This is a 'normal' error, in the sense that a rotatable group like
your protonated NH3 group in the lysine sidechain, can build up a large
rotational velocity. Since this motion couples only weakly with the
rest of the system, once you have the high rot. velocity, it stays there
for some time. If the velocity builds up too high, lincs will get problems
re-setting the bond lengths that get elongeted due to centripetal forces.
There is a finite chance to get this buil-up, which obviously increases
with simulation length...

In short, my view is that this is a more or less physical phenomenon
(i.e. not a simulation artefact), but the problem is that eventually at
very high velocities your timestep is actually too large. The workaround
I usually take is to re-start the simulation using a frame a few 10's of
ps before the first lincs warning, and introduce slight round-off errors
by using only 6 decimal places (in stead of 7 which are the full precision). 
This works as follows:
1) use trjconv to write out one frame from your .trr as a .gro file, and
   set -ndec 6 for the number of decimal places, lile:
   trjconv -f traj.trr -dump 3890 -o frame.gro -ndec 6
2) use trjconv to convert the .gro file back to .trr:
   trjconv -f frame.gro -o frame.trr
   (this is because tpbconv will only read from a .trr file)
3) use tpbconv to generate a .tpr file with the coordinates/velocities
   from the selected frame:
   tpconv -s topol.tpr -f frame.trr -o tpxout.tpr
   you might want to tune the remaining run-time using the -extend option.
4) re-start mdrun using tpxout.tpr

This trick works because the slight round-off you introduce (6 vs the 7
decimal places in a full precision .trr) will cause the simulation to
'veer' off its original path very slightly (not normally noticable), with
a good chance of avoiding the problematic velocity build-up.

-- 
Groetjes,

Anton
 _____________ _____________________________________________________
|             | Anton Feenstra       Dept. of Biophysical Chemistry |
|  _   _  ___,| K.A.Feenstra at chem.rug.nl    University of Groningen |
| / \ / \'| | | Tel +31 50 363 4327                    Nijenborgh 4 |
|(   |   )| | | Fax +31 50 363 4800               9747 AG Groningen |
| \_/ \_/ | | | http://md.chem.rug.nl/~anton        The Netherlands |
|             | "Fly to the Court of England and Unfold" (Macbeth, A|
|_____________|_____________________________________________________|