[gmx-users] links errors at 700 or 800 K

[Justin Lemkul]

On 5/21/15 10:24 AM, Christopher Neale wrote:
> Dear Users:
>
> I have a system with an 8-mer peptide in a box of water. After a 2-ns npt equilibration, I run 100 ns of nvt at various temperatures (300-800 K) in preparation for a REMD simulation. I am using velocity langevin for temperature control.
>
> Everything is fine for <=600 K. Above that, I sometimes run into LINCS errors:
>
> ######
> Step 19126507, time 38253 (ps)  LINCS WARNING

I would say a 2-fs time step is destined to fail at such a high temperature. 
Almost certainly 1-fs (still in conjunction with constraints) would be more stable.

> relative constraint deviation after LINCS:
> rms 0.069230, max 0.677838 (between atoms 49 and 50)
> bonds that rotated more than 30 degrees:
>   atom 1 atom 2  angle  previous, current, constraint length
>       49     50   90.0    0.1111   0.1864      0.1111
>       46     48   90.0    0.1111   0.1722      0.1111
> ######
>
> with eventual crash:
>
> ######
> WARNING: Listed nonbonded interaction between particles 57 and 59
> at distance 58.376 which is larger than the table limit 2.223 nm.
>
> Fatal error:
> 1 particles communicated to PME rank 0 are more than 2/3 times the cut-off out of the domain decomposition cell of their charge group in dimension y.
> This usually means that your system is not well equilibrated.
> ######
>
> I have verified this with three different force fields in gmx4.6.7 (charmm27, amber99sb-ildn, and opls). I have also seen this in gmx5 (only charmm27 force field tested so far in gmx5 -- the error snippits above are from that gmx 5 run).
>
> The gromacs 5 mdp file (with verlet cutoff scheme) is as follows:
>
> ######
> nsteps = 50000000
> constraints = all-bonds
> lincs-iter =  1
> lincs-order =  6
> constraint_algorithm =  lincs
> integrator = sd
> dt = 0.002
> tinit = 0
> nstcomm = 1
> nstxout = 0
> nstvout = 0
> nstfout = 0
> nstxtcout = 100000
> nstenergy = 100000
> nstlist = 10
> nstlog=0
> ns_type = grid
> cutoff-scheme = verlet
> vdwtype = cut-off
> rlist = 1.2
> rvdw = 1.2
> rvdw-switch = 1.0
> rcoulomb = 1.2
> vdw-modifier = force-switch
> coulombtype = PME
> ewald-rtol = 1e-5
> optimize_fft = yes
> fourierspacing = 0.12
> fourier_nx = 0
> fourier_ny = 0
> fourier_nz = 0
> pme_order = 4
> tc_grps             =  System
> tau_t               =  1.0
> ld_seed             =  -1
> ref_t = 800
> gen_temp = 800
> gen_vel = yes
> unconstrained_start = no
> gen_seed = -1
> Pcoupl = no
> ######
>
> I'm guessing that at high temperature I need to increase lincs-iter or lincs-order ? I'm asking because these errors are rare (1-3 out of 20 100-ns simulations will have such an error) so simply testing it is not that quick.
>

Possibly, but dt itself is the more likely culprit.  Simulating at 600-800 K is 
already taking a sledgehammer to the force field, anyway :)

> Very preliminary results suggest that gromacs5.0.5 with verlet cutoff scheme might have fewer such errors than gromacs 4.6.7 with the group cutoff scheme, possibly due to atoms moving out/in of the cutoff within the 20 fs between neighbourlist updates in the group scheme (3/20 crash with gmx4/group and 1/20 crashes with gmx5/verlet, but the variation between different force fields is on this order so I am not yet sure if the difference is significant).
>
> I will also note that I can get rid of these crashes by using constraints = h-bonds instead of constraints = all-bonds, though it has always been a mystery to me why constraints = h-bonds is stable (even at 300 K) with a 2 fs timestep so I am cautious about using that combination.
>

I suppose that means you're getting distortion of bonds between heavy atoms, 
which is what causes the LINCS failures.  Any unphysical geometries in these 
"constraints = h-bonds" runs?

-Justin

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

Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow

Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201

jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul

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