[gmx-users] atom moved too far

[Christos Deligkaris]

thank you Justin.

I saw on your umbrella sampling tutorial how to implement the restraints
using the pull code.

The protocol I used is (if I understand correctly your question): 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.
What are the best temperature coupling groups to use when we are not
certain whether the small molecule will spend the entire simulation period
physically bound to the macromolecule or whether it will become fully
solvated at some point? Is Macromolecule and Non-macromolecule the best
option since the small molecule will always (to a small or large extent)
interact with water (versus the Macromolecule_and_small_molecule and
everything else grouping option)?

Best wishes,

Christos Deligkaris, PhD
Assistant Professor of Physics, University of Southern Indiana
SC2220, 8600 University Blvd, Evansville IN, 47712
Office Phone: (812) 228-5056
www.deligkaris.org @DeligkarisGroup


On Mon, Dec 30, 2019 at 7:30 AM Justin Lemkul <jalemkul at vt.edu> wrote:

>
>
> 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
>
> ==================================================
>
> --
> Gromacs Users mailing list
>
> * Please search the archive at
> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before
> posting!
>
> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
>
> * For (un)subscribe requests visit
> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or
> send a mail to gmx-users-request at gromacs.org.
>