[gmx-users] LINCS WARNING after good minimization and equilibration (NPT and NVT)

[Justin A. Lemkul]

zeppelin zeppelin wrote:
> Dear GROMACS users,
> I only begin to work with gromacs and now have some problems with tetra 
> protein.
> I used AMBER force field and made several changes (in ffamber.rtp, 
> specbond.dat, ffamber99bon.itp) to create an isopeptide bonds between 
> monomers.
> Minimization and equilibration (NPT and NVT) were OK.
> But in MD (1.2 ns) I met a problem:
> __________________________________________________________________
> Step 241687, time 483.374 (ps)  LINCS WARNING
> relative constraint deviation after LINCS:
> rms nan, max inf (between atoms 2431 and 2432)
> bonds that rotated more than 30 degrees:
>  atom 1 atom 2  angle  previous, current, constraint length
>    2431   2432   90.0    0.1244      inf      0.1229
>    2428   2430  113.1    0.1010 64418817769472.0000      0.1010
>    2428   2429  105.7    0.1010 65196395593728.0000      0.1010
>   .........................................
> t = 483.374 ps: Water molecule starting at atom 57549 can not be settled.
> Check for bad contacts and/or reduce the timestep.
> Wrote pdb files with previous and current coordinates
> Segmentation fault
> ________________________________________________________________________
> When I looked on step241686c.pdb in pymol, I found, that GLY that form 
> one of isopeptide bonds exploded (it's atoms were far from the water box).
> I saw that a lot of people have a close problem, but usually they have 
> it on the first steps. Also in my case minimization and equilibration 
> were OK. I note some letters with the close problem, but I don't 
> understand what to start with in my case. I kindly ask you to give me an 
> idea or a link to a letter with the best solution.
> 

LINCS errors have many causes, and as such, there is no universal solution.  The 
normal advice for solving this problem is as follows:

http://www.gromacs.org/Documentation/Terminology/Blowing_Up#Diagnosing_an_Unstable_System

I suspect with nstxout = 1000 you should have plenty of frames to watch at the 
point where the system becomes unstable.  Looking at a single .pdb file from the 
crash will not show you how it happened; it just shows you what the result of 
the unstable system was.  It sounds to me like whatever parameters you're using 
for your isopeptide bond are not stable, but you should decide that for yourself 
based on what you see in the trajectory.

> Simulation is based on tutorial: "Lysozyme in Water. Justin Lemkul"
> 

Note that the .mdp files I provide are for OPLS-AA, and may not be universally 
applicable.  I make this caveat clear in the summary of the tutorial.  You may 
have to change some of the parameters (particularly cutoffs, nstlist, time step, 
etc) to be compatible with AMBER.

-Justin

> Thanks,
> Yulian
> 
> 
> HISTORY
> 
>    /gromacs-4.0.5/bin/pdb2gmx -f 3ALB_AMBER.pdb -o processed.gro -water 
> tip3p -missing -merge
>    /gromacs-4.0.5/bin/editconf -f processed.gro -o newbox.gro -c -d 1.0 
> -bt cubic
>    /gromacs-4.0.5/bin/genbox -cp newbox.gro -cs spc216.gro -o solv.gro 
> -p topol.top
>    /gromacs-4.0.5/bin/grompp -f ions.mdp -c solv.gro -p topol.top -o 
> ions.tpr
>    /gromacs-4.0.5/bin/genion -s ions.tpr -o solv.gro -p topol.top -pname 
> Na+ -np 24
>    /gromacs-4.0.5/bin/grompp -f minim.mdp -c solv.gro -p topol.top -o 
> em.tpr
>    /gromacs-4.0.5/bin/mdrun -v -deffnm em
>    vmd em.gro
>    /gromacs-4.0.5/bin/g_energy -f em.edr -o potential.xvg
>    xmgrace potential.xvg
>    /gromacs-4.0.5/bin/grompp -f nvt.mdp -c em.gro -p topol.top -o nvt.tpr
>    /gromacs-4.0.5/bin/mdrun -v -deffnm nvt
>    /gromacs-4.0.5/bin/g_energy -f nvt.edr
>    xmgrace nvt.edr
>    xmgrace energy.xvg
>    /gromacs-4.0.5/bin/grompp -f npt.mdp -c nvt.gro -t nvt.cpt -p 
> topol.top -o npt.tpr
>    /gromacs-4.0.5/bin/mdrun -v -deffnm npt
>    vmd npt.gro
>    /gromacs-4.0.5/bin/g_energy -f npt.edr -o pressure.xvg
>    xmgrace pressure.xvg
>    /gromacs-4.0.5/bin/g_energy -f npt.edr -o density.xvg
>    xmgrace density.xvg
>    /gromacs-4.0.5/bin/grompp -f md.mdp -c npt.gro -t npt.cpt -p 
> topol.top -o md_2ns.tpr
>    /gromacs-4.0.5/bin/mdrun -v -deffnm md_2ns
> ____________________________________________________________________
> 
> md.mdp
> 
> title        = AMBER Ub_tetra48
> ; Run parameters
> integrator    = md        ; leap-frog integrator
> nsteps        = 600000    ; 2 * 500000 = 1000 ps, 1 ns
> dt        = 0.002        ; 2 fs
> ; Output control
> nstxout        = 1000        ; save coordinates every 2 ps
> nstvout        = 1000        ; save velocities every 2 ps
> nstxtcout    = 1000        ; xtc compressed trajectory output every 2 ps
> nstenergy    = 1000        ; save energies every 2 ps
> nstlog        = 1000        ; update log file every 2 ps
> ; Bond parameters
> continuation    = yes        ; Restarting after NPT
> constraint_algorithm = lincs    ; holonomic constraints
> constraints    = all-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 cells
> nstlist        = 5        ; 10 fs
> rlist        = 1.0        ; short-range neighborlist cutoff (in nm)
> rcoulomb    = 1.0        ; short-range electrostatic cutoff (in nm)
> rvdw        = 1.0        ; 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    ; modified Berendsen thermostat
> tc-grps        = Protein Non-Protein    ; two coupling groups - more 
> accurate
> tau_t        = 0.1    0.1    ; time constant, in ps
> ref_t        = 300     300    ; reference temperature, one for each 
> group, in K
> ; Pressure coupling is on
> pcoupl        = Parrinello-Rahman    ; Pressure coupling on in NPT
> pcoupltype    = isotropic    ; uniform scaling of box vectors
> tau_p        = 2.0        ; time constant, in ps
> ref_p        = 1.0        ; reference pressure, in bar
> compressibility = 4.5e-5    ; isothermal compressibility of water, 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
> 

-- 
========================================

Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

========================================