From d8f52d842d33370a16f390b266e273f6677fe0c4 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Fri, 21 Sep 2018 16:32:14 +0200 Subject: [PATCH] Add update groups to DD atom displacement When using update groups, the DD cell size is limited by the displacement of centers of geometry of update groups raher than the displacement of atoms. This allows for smaller DD cells. Change-Id: I3bf840d09d2ce770bd049c33624e577f7aa3b6a7 --- src/gromacs/domdec/domdec.cpp | 4 +- src/gromacs/mdlib/calc_verletbuf.cpp | 279 ++++++++++++++++++++++----- src/gromacs/mdlib/calc_verletbuf.h | 22 ++- 3 files changed, 247 insertions(+), 58 deletions(-) diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index 1cc507d119..5f556966e4 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -2235,7 +2235,9 @@ static void set_dd_limits_and_grid(const gmx::MDLogger &mdlog, */ constexpr real c_chanceThatAtomMovesBeyondDomain = 1e-12; const real limitForAtomDisplacement = - minCellSizeForAtomDisplacement(*mtop, *ir, c_chanceThatAtomMovesBeyondDomain); + minCellSizeForAtomDisplacement(*mtop, *ir, + comm->updateGroupingPerMoleculetype, + c_chanceThatAtomMovesBeyondDomain); GMX_LOG(mdlog.info).appendTextFormatted( "Minimum cell size due to atom displacement: %.3f nm", limitForAtomDisplacement); diff --git a/src/gromacs/mdlib/calc_verletbuf.cpp b/src/gromacs/mdlib/calc_verletbuf.cpp index cff9084eb9..87b1f6ecb7 100644 --- a/src/gromacs/mdlib/calc_verletbuf.cpp +++ b/src/gromacs/mdlib/calc_verletbuf.cpp @@ -50,9 +50,11 @@ #include "gromacs/mdlib/nbnxn_util.h" #include "gromacs/mdtypes/inputrec.h" #include "gromacs/mdtypes/md_enums.h" +#include "gromacs/topology/block.h" #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/real.h" #include "gromacs/utility/strconvert.h" /* The code in this file estimates a pairlist buffer length @@ -1132,9 +1134,222 @@ 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) +static real +chanceOfAtomCrossingCell(gmx::ArrayRef atomtypes, + real kT_fac, + real cellSize) +{ + /* We assume atoms 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 top 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 (const VerletbufAtomtype &att : atomtypes) + { + const atom_nonbonded_kinetic_prop_t &propAtom = att.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.n*chancePerAtom; + } + + return chance; +} + +/* Struct for storing constraint properties of atoms */ +struct AtomConstraintProps +{ + void addConstraint(real length) + { + numConstraints += 1; + sumLengths += length; + } + + int numConstraints = 0; /* The number of constraints of an atom */ + real sumLengths = 0; /* The sum of constraint length over the constraints */ +}; + +/* Constructs and returns a list of constraint properties per atom */ +static std::vector +getAtomConstraintProps(const gmx_moltype_t &moltype, + const gmx_ffparams_t &ffparams) +{ + const t_atoms &atoms = moltype.atoms; + std::vector props(atoms.nr); + + for (const auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT)) + { + // Settles are handled separately + if (ilist.functionType == F_SETTLE) + { + continue; + } + + for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist)) + { + int type = ilist.iatoms[i]; + int a1 = ilist.iatoms[i + 1]; + int a2 = ilist.iatoms[i + 2]; + real length = ffparams.iparams[type].constr.dA; + props[a1].addConstraint(length); + props[a2].addConstraint(length); + } + } + + return props; +} + +/* Return the chance of at least one update group in a molecule crossing a cell of size cellSize */ +static real +chanceOfUpdateGroupCrossingCell(const gmx_moltype_t &moltype, + const gmx_ffparams_t &ffparams, + const gmx::RangePartitioning &updateGrouping, + real kT_fac, + real cellSize) +{ + const t_atoms &atoms = moltype.atoms; + GMX_ASSERT(updateGrouping.fullRange().end() == atoms.nr, "The update groups should match the molecule type"); + + const pot_derivatives_t boundaryInteraction = { 1/cellSize, 0, 0 }; + + const auto atomConstraintProps = getAtomConstraintProps(moltype, ffparams); + + real chance = 0; + for (int group = 0; group < updateGrouping.numBlocks(); group++) + { + const auto &block = updateGrouping.block(group); + /* Determine the number of atoms with constraints and the mass of the COG */ + int numAtomsWithConstraints = 0; + real massSum = 0; + for (const int &atom : block) + { + if (atomConstraintProps[atom].numConstraints > 0) + { + numAtomsWithConstraints++; + } + massSum += moltype.atoms.atom[atom].m; + } + /* Determine the maximum possible distance between the center of mass + * and the center of geometry of the update group + */ + real maxComCogDistance = 0; + if (numAtomsWithConstraints == 2) + { + for (const int &atom : block) + { + if (atomConstraintProps[atom].numConstraints > 0) + { + GMX_ASSERT(atomConstraintProps[atom].numConstraints == 1, + "Two atoms should be connected by one constraint"); + maxComCogDistance = std::abs(atoms.atom[atom].m/massSum - 0.5)*atomConstraintProps[atom].sumLengths; + break; + } + } + } + else if (numAtomsWithConstraints > 2) + { + for (const int &atom : block) + { + if (atomConstraintProps[atom].numConstraints == numAtomsWithConstraints - 1) + { + real comCogDistance = atomConstraintProps[atom].sumLengths/numAtomsWithConstraints; + maxComCogDistance = std::max(maxComCogDistance, comCogDistance); + } + } + } + else if (block.size() > 1) + { + // All normal atoms must be connected by SETTLE + for (const int &atom : block) + { + const auto &ilist = moltype.ilist[F_SETTLE]; + GMX_RELEASE_ASSERT(ilist.size() > 0, "There should be at least one settle in this moltype"); + for (int i = 0; i < ilist.size(); i += 1 + NRAL(F_SETTLE)) + { + if (atom == ilist.iatoms[i + 1]) + { + const t_iparams &iparams = ffparams.iparams[ilist.iatoms[i]]; + real dOH = iparams.settle.doh; + real dHH = iparams.settle.dhh; + real dOMidH = std::sqrt(dOH*dOH - 0.25_real*dHH*dHH); + maxComCogDistance = std::abs(atoms.atom[atom].m/massSum - 1.0_real/3.0_real)*dOMidH; + } + } + } + } + real s2_3d = kT_fac/massSum; + chance += energyDriftAtomPair(false, false, + s2_3d, 0, 0, + cellSize - 2*maxComCogDistance, + &boundaryInteraction); + } + + return chance; +} + +/* Return the chance of at least one update group in the system crossing a cell of size cellSize */ +static real +chanceOfUpdateGroupCrossingCell(const gmx_mtop_t &mtop, + PartitioningPerMoltype updateGrouping, + real kT_fac, + real cellSize) +{ + GMX_RELEASE_ASSERT(static_cast(updateGrouping.size()) == mtop.moltype.size(), + "The update groups should match the topology"); + + real chance = 0; + for (const gmx_molblock_t &molblock : mtop.molblock) + { + const gmx_moltype_t &moltype = mtop.moltype[molblock.type]; + chance += + molblock.nmol*chanceOfUpdateGroupCrossingCell(moltype, mtop.ffparams, updateGrouping[molblock.type], + kT_fac, cellSize); + } + + return chance; +} + +real +minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop, + const t_inputrec &ir, + PartitioningPerMoltype updateGrouping, + real chanceRequested) { if (!EI_DYNAMICS(ir.eI) || (EI_MD(ir.eI) && ir.etc == etcNO)) { @@ -1164,57 +1379,17 @@ real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop, 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 }; + int ib = (ib0 + ib1)/2; + cellSize = ib*resolution; - real chance = 0; - for (const VerletbufAtomtype &att : atomtypes) + real chance; + if (updateGrouping.empty()) { - const atom_nonbonded_kinetic_prop_t &propAtom = att.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.n*chancePerAtom; + chance = chanceOfAtomCrossingCell(atomtypes, kT_fac, cellSize); + } + else + { + chance = chanceOfUpdateGroupCrossingCell(mtop, updateGrouping, kT_fac, cellSize); } /* Note: chance is for every nstlist steps */ diff --git a/src/gromacs/mdlib/calc_verletbuf.h b/src/gromacs/mdlib/calc_verletbuf.h index 58c0ab860a..a444b2e097 100644 --- a/src/gromacs/mdlib/calc_verletbuf.h +++ b/src/gromacs/mdlib/calc_verletbuf.h @@ -36,12 +36,18 @@ #ifndef GMX_MDLIB_CALC_VERLETBUF_H #define GMX_MDLIB_CALC_VERLETBUF_H +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/real.h" struct gmx_mtop_t; struct t_inputrec; +namespace gmx +{ +class RangePartitioning; +} // namespace gmx + struct VerletbufListSetup { int cluster_size_i; /* Cluster pair-list i-cluster size atom count */ @@ -100,11 +106,15 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, int *n_nonlin_vsite, real *rlist); +/* Convenience type */ +using PartitioningPerMoltype = gmx::ArrayRef; + /* 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. + * an atom or update group crosses to non nearest-neighbor cells + * is <= chanceRequested within ir.nstlist steps. + * Update groups are used when !updateGrouping.empty(). * 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. * @@ -114,9 +124,11 @@ void calc_verlet_buffer_size(const gmx_mtop_t *mtop, real boxvol, * * Note: This size increases (very slowly) with system size. */ -real minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop, - const t_inputrec &ir, - real chanceRequested); +real +minCellSizeForAtomDisplacement(const gmx_mtop_t &mtop, + const t_inputrec &ir, + PartitioningPerMoltype updateGrouping, + real chanceRequested); /* Struct for unique atom type for calculating the energy drift. * The atom displacement depends on mass and constraints. -- 2.22.0