[gmx-users] implicit solvent LINCS errors

[Joshua L. Phillips]

I've found that I often get LINCS warnings like this when starting from
highly extended conformations when using implicit solvent. The GBSA
surface tension combined with the lack of viscosity (due to the absence
of explicit water) allows the protein to change conformation much faster
than LINCS likes by default.

Normally, I just run a short vacuum simulation (keep the same settings
as Justin suggested but set GBSA = No) to let the system relax a little
more before starting a GBSA run (EM is often just not enough). Usually
only about 50ps. Also, doing this with position restraints can help slow
down the collapse, and usually results in just enough collapse that the
following GBSA run will satisfy the default LINCS settings.

Another option might be to run a short simulation with GBSA turned on,
but using Langevin dynamics to include some additional friction that
should slow down the initial collapse.

Also, you could fudge with the environment variables associated with
LINCS, but this seems a little dangerous compared to the above two
suggestions.
 
I would be interested in peoples' opinions and suggestions on how to
handle this issue as it is quite common when starting from highly
extended structures using GBSA. Some proteins are more susceptible than
others...

I also tend to find that version 4.5.4 is much more stable than even
4.5.3 for implicit solvent simulations.

-- Josh

On Mon, 2011-05-30 at 18:06 -0500, Michael D. Daily wrote:
> Thanks Justin, this is very helpful.  I'll attempt these fixes tomorrow.
> 
> Mike
> 
> On 5/30/2011 5:50 PM, Justin A. Lemkul wrote:
> >
> >
> > Michael D. Daily wrote:
> >>
> >> Hi all,
> >>
> >> I'm trying to run implicit solvent calculations in gromacs 4.5 with 
> >> the charmm forcefield.  I am able to minimize successfully and 
> >> compile for 
> >
> > When troubleshooting, it is always advisable to try the latest version 
> > (4.5.4) to see if the problem is reproducible.  If a pertinent bug has 
> > been fixed, there's no use troubleshooting the broken version.
> >
> >> mdrun, but soon after starting, mdrun complains about excessive 
> >> rotation in LINCS (see the error printed below that).  I also include 
> >> my mdp file at the bottom.  Can anyone advise me as to the possible 
> >> cause of such errors, as it is difficult to diagnose given that 
> >> grompp worked fine.
> >>
> >
> > For reference:
> >
> > http://www.gromacs.org/Documentation/Terminology/Blowing_Up#Diagnosing_an_Unstable_System 
> >
> >
> > If grompp worked, that just means your coordinate and topology matched 
> > and there were no internal conflicts within the .mdp file.  It is no 
> > guarantee that the resulting simulation will actually work, 
> > unfortunately.
> >
> >> --- lincs error ---
> >>
> >> Step 1, time 0.002 (ps)  LINCS WARNING
> >
> > Typically an instant LINCS failure indicates insufficient 
> > minimization.  You said you minimized successfully, but what does this 
> > mean?  What values did you achieve for Fmax and Epot?
> >
> >> relative constraint deviation after LINCS:
> >> rms 0.000780, max 0.020692 (between atoms 880 and 881)
> >> bonds that rotated more than 30 degrees:
> >>  atom 1 atom 2  angle  previous, current, constraint length
> >>     606    607   36.7    1.0527   0.1080      0.1080
> >>     614    615   35.6    0.8972   0.1125      0.1111
> >>     614    616   75.3    0.1054   0.1121      0.1111
> >>     880    881   58.0    0.1068   0.1134      0.1111
> >>     880    882   50.4    0.9168   0.1121      0.1111
> >>     889    890   55.3    0.1066   0.1122      0.1111
> >>     889    891   35.3    0.8588   0.1118      0.1111
> >>
> >> ---- mdp file ------------
> >>
> >> ; title and include files
> >> title                    = 1EX6-S35P_md1
> >> cpp                      = cpp
> >> include                  = 
> >> -I/home/yoo2/myusr/gromacs-3.3.3/share/gromacs/top/ -I./
> >
> > Any reason you're including files from an ancient Gromacs version?
> >
> >> define                   =
> >> ; integrator and input/output setting up
> >> integrator               = md
> >> nsteps                   = 1000000 ; 2 ns
> >> ;nsteps                   = 5000 ; 2 ns
> >> dt                       = 0.002
> >> nstxout                  = 5000
> >> nstvout                  = 5000
> >> nstenergy                = 500
> >> nstxtcout                = 500
> >> nstlog                   = 500
> >> xtc_grps                 = System
> >> energygrps               = System
> >> comm_mode                = Linear
> >>
> >> ;implicit solvent
> >> implicit_solvent = GBSA
> >> gb_algorithm = Still
> >> gb_saltconc = 0.15
> >
> > FYI, gb_saltconc is nonfunctional.  Don't expect it to do anything :)
> >
> >> rgbradii = 1.0
> >>
> >> ; neighbor searching and vdw/pme setting up
> >> nstlist                  = 10
> >> ns_type                  = grid
> >> pbc                      = xyz
> >> ;rlist                    = 1.4
> >> rlist                    = 1.0
> >>
> >> ;coulombtype              = pme
> >> coulombtype              = Cut-Off
> >> fourierspacing           = 0.1
> >> pme_order                = 6
> >> ;rcoulomb                 = 1.4
> >> rcoulomb                 = 1.1
> >>
> >> ;vdwtype                  = switch
> >> vdwtype                  = Cut-Off
> >> rvdw_switch              = 1.0
> >> ;rvdw                     = 1.2
> >> rvdw                     = 1.1
> >>
> >
> > All of these are potentially problematic.  Running implicit 
> > simulations typically requires longer cutoffs than would normally be 
> > needed for explicit solvent simulations.  Try rlist=rvdw=rcoulomb=2.0 nm.
> >
> >> ; cpt control
> >> tcoupl                   = nose-hoover
> >
> > A better choice for initial equilibration would be either V-rescale or 
> > Berendsen.  I know this can be an issue in explicit solvent, when 
> > velocities can oscillate a lot at the outset of a simulation using 
> > Nose-Hoover and the simulation box can explode; I don't know if this 
> > is such a big deal with implicit, but it can't hurt to try.
> >
> >> tc-grps                  = System
> >> tau_t                    = 0.4
> >> ref_t                    = 300.0
> >> Pcoupl                   = parrinello-rahman
> >
> > I don't know how an implicit box will respond to pressure coupling, 
> > but it would be better to try NVT first and see if it's stable, then 
> > try NPT and see if things break down.
> >
> > One option that might be advantageous is to use the all-vs-all kernels 
> > for a speed upgrade.  You can accomplish this with:
> >
> > rlist = 0
> > nstlist = 0
> > rvdw = 0
> > rcoulomb = 0
> > rgbradii = 0
> > pbc = no
> > comm-mode = angular
> >
> > You'd have to run with mdrun -pd (particle decomposition), but the end 
> > result can be quite fast and you avoid potential periodicity effects.
> >
> > -Justin
> >
> 
> 
> -- 
> Michael D. Daily, Ph.D.
> Postdoctoral Fellow
> Qiang Cui group
> Department of Chemistry
> University of Wisconsin-Madison
> 

-- 
Joshua L. Phillips
Ph.D. Candidate - School of Engineering
University of California, Merced
jphillips7 at ucmerced.edu