From: Berk Hess Date: Thu, 16 Aug 2018 11:51:18 +0000 (+0200) Subject: Ensure domains are large enough for atom motion X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=cd41b0fefbc65a682c6f05ebeaa7c139d45206d7;p=alexxy%2Fgromacs.git Ensure domains are large enough for atom motion 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 much more likely that "atom moved to far" errors appeared at DD and PME redistribution. Now minimum DD cell size setting correctly takes into account atom displacement (when there is a reference temperature). Note that this can significantly increase the minimum DD cell size for solvent systems and slighlty for systems with large molecules. Fixes #2614 Change-Id: Ie41131e9eed3ef828928516a6b8ebfb9b5ba2bdb --- diff --git a/docs/release-notes/2018/2018.3.rst b/docs/release-notes/2018/2018.3.rst index c97856100d..58d8bef243 100644 --- a/docs/release-notes/2018/2018.3.rst +++ b/docs/release-notes/2018/2018.3.rst @@ -33,6 +33,20 @@ systems with constraints. :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 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 33745e1be4..913c2488de 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -64,6 +64,7 @@ #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" @@ -6402,15 +6403,37 @@ static void set_dd_limits_and_grid(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd, } 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) { diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index cf886d8079..c1a6fcdcce 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -1177,3 +1177,117 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, *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; +} diff --git a/src/gromacs/mdlib/calc_verletbuf.h b/src/gromacs/mdlib/calc_verletbuf.h index 9bd1ba8149..65cecc3113 100644 --- a/src/gromacs/mdlib/calc_verletbuf.h +++ b/src/gromacs/mdlib/calc_verletbuf.h @@ -1,7 +1,7 @@ /* * 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. @@ -100,6 +100,24 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, 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.