[gmx-users] atom moved too far

[Justin Lemkul]

On 12/30/19 5:34 AM, Christos Deligkaris wrote:
> dear all,
>
> While running my production run of a small molecule noncovalently bound to
> DNA in water, I got the error:
>
> Step 8921600:
>
> Atom 20338 moved more than the distance allowed by the domain decomposition
> (0.457738) in direction X
>
> distance out of cell 0.495788
>
> New coordinates:    3.304    4.938    2.969
>
> Old cell boundaries in direction X:    0.000    0.715
>
> New cell boundaries in direction X:    0.000    0.709
>
>
> -------------------------------------------------------
>
> Program:     gmx mdrun, version 2018.1
>
> Source file: src/gromacs/domdec/domdec.cpp (line 4076)
>
> MPI rank:    0 (out of 12)
>
>
> Fatal error:
>
> An atom moved too far between two domain decomposition steps
>
> This usually means that your system is not well equilibrated
>
> The error occured at about time=18 ns. I erased the output files,
> resubmitted the calculation again from the beginning time=0 and in that
> attempt the error appeared at about 16 ns. In both cases it was an atom of
> a water molecule that moved too far.
>
> I read the "gromacs blow up" webpage and thought maybe the harmonic
> intermolecular constraints I have on the cap base pairs are too strong for
> the time step of 0.002 so I decreased the time step to 0.001ps and that
> calculation is now at about 60 ns (still running). I also tried using only
> 6 cores with a step of 0.002 (I used 12 cores in the two failed
> calculations) and that calculation is now at time=50ns (still running).
>
> The energy minimization reached the cutoff for maximum force 1000 in 346
> steps. My NVT equilibration was 50,000 steps of dt = 0.002 and I used 8 NPT
> equilibration calculations, each with 31,250 steps of dt=0.002 (total NPT
> equilibration time 500ps) where I slowly decreased the position restraints
> on DNA and the small molecule, as well as the harmonic restraint between
> the two.
>
> I implemented the harmonic restraint between the small molecule-DNA and the
> cap base pairs using the following in the topology (just showing one
> example)
>
> [ intermolecular_interactions ]
>
> [ bonds ]
>
> 455     483      6       0.18564  200

It is easier to use the pull code to perform the same operation, and is 
more likely to be compatible with domain decomposition.

> For the production run, I use
>
> comm-mode = Linear
>
> nstcomm = 100
>
> comm-grps = System
>
>
> cutoff-scheme = Verlet
>
> nstlist = 10
>
> coulombtype = PME
>
> rcoulomb = 1.2
>
> vdwtype = Cut-off
>
> vdw-modifier = Potential-shift
>
> rvdw = 1.2
>
> rvdw-switch = 0
>
> DispCorr = no
>
> tcoupl = nose-hoover
>
> nh-chain-length = 1
>
> nsttcouple = -1
>
> tc-grps = DNA_NNK Water_and_ions
>
> tau-t = 0.4 0.4
>
> ref-t = 300 300
>
> pcoupl = Parrinello-Rahman
>
> pcoupltype = isotropic
>
> tau-p = 1.0
>
> ref-p = 1.0
>
> refcoord-scaling = com
>
> constraints = h-bonds
>
> constraint-algorithm = lincs
>
> I believe SETTLES is included in the TIP3P topology by default so I am
> using that also.
>
> I have run three 50 ns production runs of the small molecule in water
> without any problems (NPT equilibration was done in 1 calculation for these
> three).
>
> I was surprised that the two failed calculations did not fail at exactly
> the same time step (start with same positions, velocities and solve
> determinist equations of motion).

Domain decomposition and other optimizations introduce small 
differences, so you should not expect different runs to behave 
identically, even with the same .tpr file.

> The error message suggests that the system was not well equilibrated.
> Should I try doing my NPT equilibrations over 1 ns instead of 500ps?

What has been your preparation protocol so far?

> Is it possible that using as a temperature coupling group the
> macromolecule-small molecule results in the blow up if the small molecule
> ends up leaving the macromolecule? (I did see the small molecule becoming
> fully solvated at some point in the trajectory of the failed calculations)
>
> Using a force constant of 200 kJ/mol/nm is perhaps too large for a time
> step of 0.002 ps? (I think I tried removing the harmonic restraints between
> the two bases on the cap base pairs and that calculation also failed). And
> I have seen in the literature a time step of 0.002 ps is reasonable for DNA
> and DNA-ligand calculations.
>
> It is not clear to me why using less cores allows the calculation to
> proceed (unless it got to 50ns by chance and it will blow up eventually).

Probably because with fewer cores you are not using domain 
decomposition, which points to some kind of incompatibility with your 
use of [intermolecular_interactions] and domain decomposition. If you 
really need some kind of restraint, try the pull code or merge the DNA 
chains into a single topology and apply a simple distance restraint.

-Justin

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

Justin A. Lemkul, Ph.D.
Assistant Professor
Office: 301 Fralin Hall
Lab: 303 Engel Hall

Virginia Tech Department of Biochemistry
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.thelemkullab.com

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