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

[Christopher Neale]

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
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.

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.

Please let me note in advance that although 800K is quite hot, this is the system that I am trying to fix.

Thank you for your assistance,
Chris.