[gmx-users] implicit solvent LINCS errors

[Justin A. Lemkul]

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

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

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

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