:issue:`2613`
+Avoid "atom moved to far" errors
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+The introduction of the dual pair list has led to larger nstlist values,
+which leads to larger atom displacements between domain decomposition
+steps. This has made it more likely that the errors
+"An atom moved too far between two domain decomposition steps" and
+"N particles communicated to PME rank M are more than 2/3 times the cut-off
+out of the domain decomposition cell ..." appear for stable systems.
+Now atom displacements are correctly taken into account when determining
+the minimum cell size, so these errors should only appear for unstable systems.
+
+:issue:`2614`
+
grompp now checks that pull groups are not close to half the box size
""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
#include "gromacs/math/functions.h"
#include "gromacs/math/vec.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdlib/calc_verletbuf.h"
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdlib/force.h"
#include "gromacs/mdlib/forcerec.h"
}
comm->cutoff_mbody = 0;
+ /* Determine the minimum cell size limit, affected by many factors */
comm->cellsize_limit = 0;
comm->bBondComm = FALSE;
- /* Atoms should be able to move by up to half the list buffer size (if > 0)
- * within nstlist steps. Since boundaries are allowed to displace by half
- * a cell size, DD cells should be at least the size of the list buffer.
+ /* We do not allow home atoms to move beyond the neighboring domain
+ * between domain decomposition steps, which limits the cell size.
+ * Get an estimate of cell size limit due to atom displacement.
+ * In most cases this is a large overestimate, because it assumes
+ * non-interaction atoms.
+ * We set the chance to 1 in a trillion steps.
*/
+ constexpr real c_chanceThatAtomMovesBeyondDomain = 1e-12;
+ const real limitForAtomDisplacement =
+ minCellSizeForAtomDisplacement(*mtop, *ir, c_chanceThatAtomMovesBeyondDomain);
+ if (fplog)
+ {
+ fprintf(fplog,
+ "Minimum cell size due to atom displacement: %.3f nm\n",
+ limitForAtomDisplacement);
+ }
comm->cellsize_limit = std::max(comm->cellsize_limit,
- ir->rlist - std::max(ir->rvdw, ir->rcoulomb));
+ limitForAtomDisplacement);
+
+ /* TODO: PME decomposition currently requires atoms not to be more than
+ * 2/3 of comm->cutoff, which is >=rlist, outside of their domain.
+ * In nearly all cases, limitForAtomDisplacement will be smaller
+ * than 2/3*rlist, so the PME requirement is satisfied.
+ * But it would be better for both correctness and performance
+ * to use limitForAtomDisplacement instead of 2/3*comm->cutoff.
+ * Note that we would need to improve the pairlist buffer case.
+ */
if (comm->bInterCGBondeds)
{
*rlist = std::max(ir->rvdw, ir->rcoulomb) + ib1*resolution;
}
+
+/* Returns the pairlist buffer size for use as a minimum buffer size
+ *
+ * Note that this is a rather crude estimate. It is ok for a buffer
+ * set for good energy conservation or RF electrostatics. But it is
+ * too small with PME and the buffer set with the default tolerance.
+ */
+static real minCellSizeFromPairlistBuffer(const t_inputrec &ir)
+{
+ return ir.rlist - std::max(ir.rvdw, ir.rcoulomb);
+}
+
+real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ real chanceRequested)
+{
+ if (!EI_DYNAMICS(ir.eI) || (EI_MD(ir.eI) && ir.etc == etcNO))
+ {
+ return minCellSizeFromPairlistBuffer(ir);
+ }
+
+ /* We use the maximum temperature with multiple T-coupl groups.
+ * We could use a per particle temperature, but since particles
+ * interact, this might underestimate the displacements.
+ */
+ const real temperature = maxReferenceTemperature(ir);
+
+ const bool setMassesToOne = (ir.eI == eiBD && ir.bd_fric > 0);
+
+ verletbuf_atomtype_t *att = nullptr;
+ int natt = -1;
+ get_verlet_buffer_atomtypes(&mtop, setMassesToOne, &att, &natt, nullptr);
+
+ const real kT_fac = displacementVariance(ir, temperature,
+ ir.nstlist*ir.delta_t);
+
+ /* Resolution of the cell size */
+ real resolution = 0.001;
+
+ /* Search using bisection, avoid 0 and start at 1 */
+ int ib0 = 0;
+ /* The chance will be neglible at 10 times the max sigma */
+ int ib1 = (int)(10*maxSigma(kT_fac, natt, att)/resolution) + 1;
+ real cellSize = 0;
+ while (ib1 - ib0 > 1)
+ {
+ int ib = (ib0 + ib1)/2;
+ cellSize = ib*resolution;
+
+ /* We assumes atom are distributed uniformly over the cell width.
+ * Once an atom has moved by more than the cellSize (as passed
+ * as the buffer argument to energyDriftAtomPair() below),
+ * the chance of crossing the boundary of the neighbor cell
+ * thus increases as 1/cellSize with the additional displacement
+ * on to of cellSize. We thus create a linear interaction with
+ * derivative = -1/cellSize. Using this in the energyDriftAtomPair
+ * function will return the chance of crossing the next boundary.
+ */
+ const pot_derivatives_t boundaryInteraction = { 1/cellSize, 0, 0 };
+
+ real chance = 0;
+ for (int i = 0; i < natt; i++)
+ {
+ const atom_nonbonded_kinetic_prop_t &propAtom = att[i].prop;
+ real s2_2d;
+ real s2_3d;
+ get_atom_sigma2(kT_fac, &propAtom, &s2_2d, &s2_3d);
+
+ real chancePerAtom = energyDriftAtomPair(propAtom.bConstr, false,
+ s2_2d + s2_3d, s2_2d, 0,
+ cellSize,
+ &boundaryInteraction);
+
+ if (propAtom.bConstr)
+ {
+ /* energyDriftAtomPair() uses an unlimited Gaussian displacement
+ * distribution for constrained atoms, whereas they can
+ * actually not move more than the COM of the two constrained
+ * atoms plus twice the distance from the COM.
+ * Use this maximum, limited displacement when this results in
+ * a smaller chance (note that this is still an overestimate).
+ */
+ real massFraction = propAtom.con_mass/(propAtom.mass + propAtom.con_mass);
+ real comDistance = propAtom.con_len*massFraction;
+
+ real chanceWithMaxDistance =
+ energyDriftAtomPair(false, false,
+ s2_3d, 0, 0,
+ cellSize - 2*comDistance,
+ &boundaryInteraction);
+ chancePerAtom = std::min(chancePerAtom, chanceWithMaxDistance);
+ }
+
+ /* Take into account the line density of the boundary */
+ chancePerAtom /= cellSize;
+
+ chance += att[i].n*chancePerAtom;
+ }
+
+ /* Note: chance is for every nstlist steps */
+ if (chance > chanceRequested*ir.nstlist)
+ {
+ ib0 = ib;
+ }
+ else
+ {
+ ib1 = ib;
+ }
+ }
+
+ sfree(att);
+
+ return cellSize;
+}
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2012,2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
int *n_nonlin_vsite,
real *rlist);
+/* Determines the mininum cell size based on atom displacement
+ *
+ * The value returned is the minimum size for which the chance that
+ * an atom crosses to non nearest-neighbor cells is <= chanceRequested
+ * within ir.nstlist steps.
+ * Without T-coupling, SD or BD, we can not estimate atom displacements
+ * and fall back to the, crude, estimate of using the pairlist buffer size.
+ *
+ * Note: Like the Verlet buffer estimate, this estimate is based on
+ * non-interacting atoms and constrained atom-pairs. Therefore for
+ * any system that is not an ideal gas, this will be an overestimate.
+ *
+ * Note: This size increases (very slowly) with system size.
+ */
+real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop,
+ const t_inputrec &ir,
+ real chanceRequested);
+
/* Struct for unique atom type for calculating the energy drift.
* The atom displacement depends on mass and constraints.
* The energy jump for given distance depend on LJ type and q.
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff