From b8afebc4bbfc94c71d934fbe488c9bc501645123 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Wed, 5 May 2021 08:53:36 +0200 Subject: [PATCH] Move the code of what will later implement LocalTopologyChecker Some common implementation functionality used in making the reverse topology also needs to be exposed. Otherwise, this is pure code movement with no functionality changes. Refs #3887 --- src/gromacs/domdec/domdec.h | 42 -- src/gromacs/domdec/domdec_topology.cpp | 344 +-------------- src/gromacs/domdec/localtopologychecker.cpp | 409 ++++++++++++++++++ src/gromacs/domdec/localtopologychecker.h | 49 +++ src/gromacs/domdec/partition.cpp | 1 + src/gromacs/domdec/reversetopology.h | 6 + src/gromacs/mdlib/stat.cpp | 3 +- src/gromacs/mdlib/update_vv.cpp | 1 + src/gromacs/mdrun/md.cpp | 1 + src/gromacs/mdrun/mimic.cpp | 1 + src/gromacs/mdrun/rerun.cpp | 1 + .../computeglobalselement.cpp | 1 + 12 files changed, 474 insertions(+), 385 deletions(-) create mode 100644 src/gromacs/domdec/localtopologychecker.cpp diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index d157513be3..bf159fa6ce 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -241,15 +241,6 @@ gmx::ArrayRef dd_constraints_nlocalatoms(const gmx_domdec_t* dd); /* In domdec_top.c */ -/*! \brief Print error output when interactions are missing */ -[[noreturn]] void dd_print_missing_interactions(const gmx::MDLogger& mdlog, - t_commrec* cr, - int local_count, - const gmx_mtop_t& top_global, - const gmx_localtop_t* top_local, - gmx::ArrayRef x, - const matrix box); - /*! \brief Generate and store the reverse topology */ void dd_make_reverse_top(FILE* fplog, gmx_domdec_t* dd, @@ -258,39 +249,6 @@ void dd_make_reverse_top(FILE* fplog, const t_inputrec& inputrec, gmx::DDBondedChecking ddBondedChecking); -/*! \brief Set that the number of bonded interactions in the local - * topology should be checked via observables reduction. */ -void scheduleCheckOfLocalTopology(gmx_domdec_t* dd, int numBondedInteractionsToReduce); - -/*! \brief Return whether the total bonded interaction count across - * domains should be checked this step. */ -bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd); - -//! Return the number of bonded interactions in this domain. -int numBondedInteractions(const gmx_domdec_t& dd); - -/*! \brief Set total bonded interaction count across domains. */ -void setNumberOfBondedInteractionsOverAllDomains(const gmx_domdec_t& dd, int newValue); - -/*! \brief Check whether bonded interactions are missing from the reverse topology - * produced by domain decomposition. - * - * Must only be called when DD is active. - * - * \param[in] mdlog Logger - * \param[in] cr Communication object - * \param[in] top_global Global topology for the error message - * \param[in] top_local Local topology for the error message - * \param[in] x Position vector for the error message - * \param[in] box Box matrix for the error message - */ -void checkNumberOfBondedInteractions(const gmx::MDLogger& mdlog, - t_commrec* cr, - const gmx_mtop_t& top_global, - const gmx_localtop_t* top_local, - gmx::ArrayRef x, - const matrix box); - /*! \brief Generate the local topology and virtual site data * * \returns Total count of bonded interactions in the local topology on this domain */ diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp index 72a2e07ec3..09839f037e 100644 --- a/src/gromacs/domdec/domdec_topology.cpp +++ b/src/gromacs/domdec/domdec_topology.cpp @@ -163,8 +163,7 @@ struct gmx_reverse_top_t::Impl }; -/*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */ -static int nral_rt(int ftype) +int nral_rt(int ftype) { int nral = NRAL(ftype); if (interaction_function[ftype].flags & IF_VSITE) @@ -178,8 +177,7 @@ static int nral_rt(int ftype) return nral; } -/*! \brief Return whether interactions of type \p ftype need to be assigned exactly once */ -static gmx_bool dd_check_ftype(const int ftype, const ReverseTopOptions rtOptions) +bool dd_check_ftype(const int ftype, const ReverseTopOptions& rtOptions) { return ((((interaction_function[ftype].flags & IF_BOND) != 0U) && ((interaction_function[ftype].flags & IF_VSITE) == 0U) @@ -189,280 +187,6 @@ static gmx_bool dd_check_ftype(const int ftype, const ReverseTopOptions rtOption || (rtOptions.includeSettles && ftype == F_SETTLE)); } -/*! \brief Checks whether interactions have been assigned for one function type - * - * Loops over a list of interactions in the local topology of one function type - * and flags each of the interactions as assigned in the global \p isAssigned list. - * Exits with an inconsistency error when an interaction is assigned more than once. - */ -static void flagInteractionsForType(const int ftype, - const InteractionList& il, - const reverse_ilist_t& ril, - const gmx::Range& atomRange, - const int numAtomsPerMolecule, - gmx::ArrayRef globalAtomIndices, - gmx::ArrayRef isAssigned) -{ - const int nril_mol = ril.index[numAtomsPerMolecule]; - const int nral = NRAL(ftype); - - for (int i = 0; i < il.size(); i += 1 + nral) - { - // ia[0] is the interaction type, ia[1, ...] the atom indices - const int* ia = il.iatoms.data() + i; - // Extract the global atom index of the first atom in this interaction - const int a0 = globalAtomIndices[ia[1]]; - /* Check if this interaction is in - * the currently checked molblock. - */ - if (atomRange.isInRange(a0)) - { - // The molecule index in the list of this molecule type - const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule; - const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule; - const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule; - int j_mol = ril.index[atomOffset]; - bool found = false; - while (j_mol < ril.index[atomOffset + 1] && !found) - { - const int j = moleculeIndex * nril_mol + j_mol; - const int ftype_j = ril.il[j_mol]; - /* Here we need to check if this interaction has - * not already been assigned, since we could have - * multiply defined interactions. - */ - if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0) - { - /* Check the atoms */ - found = true; - for (int a = 0; a < nral; a++) - { - if (globalAtomIndices[ia[1 + a]] - != globalAtomStartInMolecule + ril.il[j_mol + 2 + a]) - { - found = false; - } - } - if (found) - { - isAssigned[j] = 1; - } - } - j_mol += 2 + nral_rt(ftype_j); - } - if (!found) - { - gmx_incons("Some interactions seem to be assigned multiple times"); - } - } - } -} - -/*! \brief Help print error output when interactions are missing in a molblock - * - * \note This function needs to be called on all ranks (contains a global summation) - */ -static std::string printMissingInteractionsMolblock(t_commrec* cr, - const gmx_reverse_top_t& rt, - const char* moltypename, - const reverse_ilist_t& ril, - const gmx::Range& atomRange, - const int numAtomsPerMolecule, - const int numMolecules, - const InteractionDefinitions& idef) -{ - const int nril_mol = ril.index[numAtomsPerMolecule]; - std::vector isAssigned(numMolecules * nril_mol, 0); - gmx::StringOutputStream stream; - gmx::TextWriter log(&stream); - - for (int ftype = 0; ftype < F_NRE; ftype++) - { - if (dd_check_ftype(ftype, rt.options())) - { - flagInteractionsForType( - ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned); - } - } - - gmx_sumi(isAssigned.size(), isAssigned.data(), cr); - - const int numMissingToPrint = 10; - int i = 0; - for (int mol = 0; mol < numMolecules; mol++) - { - int j_mol = 0; - while (j_mol < nril_mol) - { - int ftype = ril.il[j_mol]; - int nral = NRAL(ftype); - int j = mol * nril_mol + j_mol; - if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE)) - { - if (DDMASTER(cr->dd)) - { - if (i == 0) - { - log.writeLineFormatted("Molecule type '%s'", moltypename); - log.writeLineFormatted( - "the first %d missing interactions, except for exclusions:", - numMissingToPrint); - } - log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname); - int a = 0; - for (; a < nral; a++) - { - log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1); - } - while (a < 4) - { - log.writeString(" "); - a++; - } - log.writeString(" global"); - for (int a = 0; a < nral; a++) - { - log.writeStringFormatted("%6d", - atomRange.begin() + mol * numAtomsPerMolecule - + ril.il[j_mol + 2 + a] + 1); - } - log.ensureLineBreak(); - } - i++; - if (i >= numMissingToPrint) - { - break; - } - } - j_mol += 2 + nral_rt(ftype); - } - } - - return stream.toString(); -} - -/*! \brief Help print error output when interactions are missing */ -static void printMissingInteractionsAtoms(const gmx::MDLogger& mdlog, - t_commrec* cr, - const gmx_mtop_t& mtop, - const InteractionDefinitions& idef) -{ - const gmx_reverse_top_t& rt = *cr->dd->reverse_top; - - /* Print the atoms in the missing interactions per molblock */ - int a_end = 0; - for (const gmx_molblock_t& molb : mtop.molblock) - { - const gmx_moltype_t& moltype = mtop.moltype[molb.type]; - const int a_start = a_end; - a_end = a_start + molb.nmol * moltype.atoms.nr; - const gmx::Range atomRange(a_start, a_end); - - auto warning = printMissingInteractionsMolblock(cr, - rt, - *(moltype.name), - rt.interactionListForMoleculeType(molb.type), - atomRange, - moltype.atoms.nr, - molb.nmol, - idef); - - GMX_LOG(mdlog.warning).appendText(warning); - } -} - -void dd_print_missing_interactions(const gmx::MDLogger& mdlog, - t_commrec* cr, - int numBondedInteractionsOverAllDomains, - const gmx_mtop_t& top_global, - const gmx_localtop_t* top_local, - gmx::ArrayRef x, - const matrix box) -{ - int cl[F_NRE]; - gmx_domdec_t* dd = cr->dd; - - GMX_LOG(mdlog.warning) - .appendText( - "Not all bonded interactions have been properly assigned to the domain " - "decomposition cells"); - - const int ndiff_tot = numBondedInteractionsOverAllDomains - - dd->reverse_top->expectedNumGlobalBondedInteractions(); - - for (int ftype = 0; ftype < F_NRE; ftype++) - { - const int nral = NRAL(ftype); - cl[ftype] = top_local->idef.il[ftype].size() / (1 + nral); - } - - gmx_sumi(F_NRE, cl, cr); - - if (DDMASTER(dd)) - { - GMX_LOG(mdlog.warning).appendText("A list of missing interactions:"); - int rest_global = dd->reverse_top->expectedNumGlobalBondedInteractions(); - int rest = numBondedInteractionsOverAllDomains; - for (int ftype = 0; ftype < F_NRE; ftype++) - { - /* In the reverse and local top all constraints are merged - * into F_CONSTR. So in the if statement we skip F_CONSTRNC - * and add these constraints when doing F_CONSTR. - */ - if (dd_check_ftype(ftype, dd->reverse_top->options()) && ftype != F_CONSTRNC) - { - int n = gmx_mtop_ftype_count(top_global, ftype); - if (ftype == F_CONSTR) - { - n += gmx_mtop_ftype_count(top_global, F_CONSTRNC); - } - int ndiff = cl[ftype] - n; - if (ndiff != 0) - { - GMX_LOG(mdlog.warning) - .appendTextFormatted("%20s of %6d missing %6d", - interaction_function[ftype].longname, - n, - -ndiff); - } - rest_global -= n; - rest -= cl[ftype]; - } - } - - int ndiff = rest - rest_global; - if (ndiff != 0) - { - GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff); - } - } - - printMissingInteractionsAtoms(mdlog, cr, top_global, top_local->idef); - write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box); - - std::string errorMessage; - - if (ndiff_tot > 0) - { - errorMessage = - "One or more interactions were assigned to multiple domains of the domain " - "decompostion. Please report this bug."; - } - else - { - errorMessage = gmx::formatString( - "%d of the %d bonded interactions could not be calculated because some atoms " - "involved moved further apart than the multi-body cut-off distance (%g nm) or the " - "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds " - "also see option -ddcheck", - -ndiff_tot, - dd->reverse_top->expectedNumGlobalBondedInteractions(), - dd_cutoff_multibody(dd), - dd_cutoff_twobody(dd)); - } - gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str()); -} - /*! \brief Return global topology molecule information for global atom index \p i_gl */ static void global_atomnr_to_moltype_ind(ArrayRef molblockIndices, int i_gl, @@ -1687,70 +1411,6 @@ static int make_local_bondeds_excls(gmx_domdec_t* dd, return numBondedInteractions; } -void scheduleCheckOfLocalTopology(gmx_domdec_t* dd, const int numBondedInteractionsToReduce) -{ - dd->localTopologyChecker->numBondedInteractionsToReduce = numBondedInteractionsToReduce; - // Note that it's possible for this to still be true from the last - // time it was set, e.g. if repartitioning was triggered before - // global communication that would have acted on the true - // value. This could happen for example when replica exchange took - // place soon after a partition. - dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = true; - // Clear the old global value, which is now invalid - dd->localTopologyChecker->numBondedInteractionsOverAllDomains.reset(); -} - -bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd) -{ - return dd.localTopologyChecker->shouldCheckNumberOfBondedInteractions; -} - -int numBondedInteractions(const gmx_domdec_t& dd) -{ - return dd.localTopologyChecker->numBondedInteractionsToReduce; -} - -void setNumberOfBondedInteractionsOverAllDomains(const gmx_domdec_t& dd, int newValue) -{ - GMX_RELEASE_ASSERT(!dd.localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(), - "Cannot set number of bonded interactions because it is already set"); - dd.localTopologyChecker->numBondedInteractionsOverAllDomains.emplace(newValue); -} - -void checkNumberOfBondedInteractions(const gmx::MDLogger& mdlog, - t_commrec* cr, - const gmx_mtop_t& top_global, - const gmx_localtop_t* top_local, - gmx::ArrayRef x, - const matrix box) -{ - GMX_RELEASE_ASSERT( - DOMAINDECOMP(cr), - "No need to check number of bonded interactions when not using domain decomposition"); - if (cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions) - { - GMX_RELEASE_ASSERT(cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(), - "The check for the total number of bonded interactions requires the " - "value to have been reduced across all domains"); - if (cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value() - != cr->dd->reverse_top->expectedNumGlobalBondedInteractions()) - { - dd_print_missing_interactions( - mdlog, - cr, - cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value(), - top_global, - top_local, - x, - box); // Does not return - } - // Now that the value is set and the check complete, future - // global communication should not compute the value until - // after the next partitioning. - cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = false; - } -} - int dd_make_local_top(gmx_domdec_t* dd, gmx_domdec_zones_t* zones, int npbcdim, diff --git a/src/gromacs/domdec/localtopologychecker.cpp b/src/gromacs/domdec/localtopologychecker.cpp new file mode 100644 index 0000000000..15f7b4900a --- /dev/null +++ b/src/gromacs/domdec/localtopologychecker.cpp @@ -0,0 +1,409 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2021, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ + +/*! \internal \file + * + * \brief This file defines functions used by the domdec module + * while managing the construction, use and error checking for + * topologies local to a DD rank. + * + * \ingroup module_domdec + */ + +#include "gmxpre.h" + +#include +#include + +#include "gromacs/domdec/localtopologychecker.h" + +#include "gromacs/domdec/domdec_internal.h" +#include "gromacs/domdec/reversetopology.h" +#include "gromacs/gmxlib/network.h" +#include "gromacs/mdtypes/commrec.h" +#include "gromacs/topology/idef.h" +#include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/fatalerror.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/stringstream.h" +#include "gromacs/utility/textwriter.h" + +#include "dump.h" + +using gmx::ArrayRef; +using gmx::RVec; + +/*! \brief Checks whether interactions have been assigned for one function type + * + * Loops over a list of interactions in the local topology of one function type + * and flags each of the interactions as assigned in the global \p isAssigned list. + * Exits with an inconsistency error when an interaction is assigned more than once. + */ +static void flagInteractionsForType(const int ftype, + const InteractionList& il, + const reverse_ilist_t& ril, + const gmx::Range& atomRange, + const int numAtomsPerMolecule, + ArrayRef globalAtomIndices, + ArrayRef isAssigned) +{ + const int nril_mol = ril.index[numAtomsPerMolecule]; + const int nral = NRAL(ftype); + + for (int i = 0; i < il.size(); i += 1 + nral) + { + // ia[0] is the interaction type, ia[1, ...] the atom indices + const int* ia = il.iatoms.data() + i; + // Extract the global atom index of the first atom in this interaction + const int a0 = globalAtomIndices[ia[1]]; + /* Check if this interaction is in + * the currently checked molblock. + */ + if (atomRange.isInRange(a0)) + { + // The molecule index in the list of this molecule type + const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule; + const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule; + const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule; + int j_mol = ril.index[atomOffset]; + bool found = false; + while (j_mol < ril.index[atomOffset + 1] && !found) + { + const int j = moleculeIndex * nril_mol + j_mol; + const int ftype_j = ril.il[j_mol]; + /* Here we need to check if this interaction has + * not already been assigned, since we could have + * multiply defined interactions. + */ + if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0) + { + /* Check the atoms */ + found = true; + for (int a = 0; a < nral; a++) + { + if (globalAtomIndices[ia[1 + a]] + != globalAtomStartInMolecule + ril.il[j_mol + 2 + a]) + { + found = false; + } + } + if (found) + { + isAssigned[j] = 1; + } + } + j_mol += 2 + nral_rt(ftype_j); + } + if (!found) + { + gmx_incons("Some interactions seem to be assigned multiple times"); + } + } + } +} + +/*! \brief Help print error output when interactions are missing in a molblock + * + * \note This function needs to be called on all ranks (contains a global summation) + */ +static std::string printMissingInteractionsMolblock(t_commrec* cr, + const gmx_reverse_top_t& rt, + const char* moltypename, + const reverse_ilist_t& ril, + const gmx::Range& atomRange, + const int numAtomsPerMolecule, + const int numMolecules, + const InteractionDefinitions& idef) +{ + const int nril_mol = ril.index[numAtomsPerMolecule]; + std::vector isAssigned(numMolecules * nril_mol, 0); + gmx::StringOutputStream stream; + gmx::TextWriter log(&stream); + + for (int ftype = 0; ftype < F_NRE; ftype++) + { + if (dd_check_ftype(ftype, rt.options())) + { + flagInteractionsForType( + ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned); + } + } + + gmx_sumi(isAssigned.size(), isAssigned.data(), cr); + + const int numMissingToPrint = 10; + int i = 0; + for (int mol = 0; mol < numMolecules; mol++) + { + int j_mol = 0; + while (j_mol < nril_mol) + { + int ftype = ril.il[j_mol]; + int nral = NRAL(ftype); + int j = mol * nril_mol + j_mol; + if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE)) + { + if (DDMASTER(cr->dd)) + { + if (i == 0) + { + log.writeLineFormatted("Molecule type '%s'", moltypename); + log.writeLineFormatted( + "the first %d missing interactions, except for exclusions:", + numMissingToPrint); + } + log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname); + int a = 0; + for (; a < nral; a++) + { + log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1); + } + while (a < 4) + { + log.writeString(" "); + a++; + } + log.writeString(" global"); + for (int a = 0; a < nral; a++) + { + log.writeStringFormatted("%6d", + atomRange.begin() + mol * numAtomsPerMolecule + + ril.il[j_mol + 2 + a] + 1); + } + log.ensureLineBreak(); + } + i++; + if (i >= numMissingToPrint) + { + break; + } + } + j_mol += 2 + nral_rt(ftype); + } + } + + return stream.toString(); +} + +/*! \brief Help print error output when interactions are missing */ +static void printMissingInteractionsAtoms(const gmx::MDLogger& mdlog, + t_commrec* cr, + const gmx_mtop_t& mtop, + const InteractionDefinitions& idef) +{ + const gmx_reverse_top_t& rt = *cr->dd->reverse_top; + + /* Print the atoms in the missing interactions per molblock */ + int a_end = 0; + for (const gmx_molblock_t& molb : mtop.molblock) + { + const gmx_moltype_t& moltype = mtop.moltype[molb.type]; + const int a_start = a_end; + a_end = a_start + molb.nmol * moltype.atoms.nr; + const gmx::Range atomRange(a_start, a_end); + + auto warning = printMissingInteractionsMolblock( + cr, + rt, + *(moltype.name), + cr->dd->reverse_top->interactionListForMoleculeType(molb.type), + atomRange, + moltype.atoms.nr, + molb.nmol, + idef); + + GMX_LOG(mdlog.warning).appendText(warning); + } +} + +/*! \brief Print error output when interactions are missing */ +[[noreturn]] static void dd_print_missing_interactions(const gmx::MDLogger& mdlog, + t_commrec* cr, + int numBondedInteractionsOverAllDomains, + const gmx_mtop_t& top_global, + const gmx_localtop_t* top_local, + ArrayRef x, + const matrix box) +{ + int cl[F_NRE]; + gmx_domdec_t* dd = cr->dd; + + GMX_LOG(mdlog.warning) + .appendText( + "Not all bonded interactions have been properly assigned to the domain " + "decomposition cells"); + + const int ndiff_tot = numBondedInteractionsOverAllDomains + - dd->reverse_top->expectedNumGlobalBondedInteractions(); + + for (int ftype = 0; ftype < F_NRE; ftype++) + { + const int nral = NRAL(ftype); + cl[ftype] = top_local->idef.il[ftype].size() / (1 + nral); + } + + gmx_sumi(F_NRE, cl, cr); + + if (DDMASTER(dd)) + { + GMX_LOG(mdlog.warning).appendText("A list of missing interactions:"); + int rest_global = dd->reverse_top->expectedNumGlobalBondedInteractions(); + int rest = numBondedInteractionsOverAllDomains; + for (int ftype = 0; ftype < F_NRE; ftype++) + { + /* In the reverse and local top all constraints are merged + * into F_CONSTR. So in the if statement we skip F_CONSTRNC + * and add these constraints when doing F_CONSTR. + */ + if (dd_check_ftype(ftype, dd->reverse_top->options()) && ftype != F_CONSTRNC) + { + int n = gmx_mtop_ftype_count(top_global, ftype); + if (ftype == F_CONSTR) + { + n += gmx_mtop_ftype_count(top_global, F_CONSTRNC); + } + int ndiff = cl[ftype] - n; + if (ndiff != 0) + { + GMX_LOG(mdlog.warning) + .appendTextFormatted("%20s of %6d missing %6d", + interaction_function[ftype].longname, + n, + -ndiff); + } + rest_global -= n; + rest -= cl[ftype]; + } + } + + int ndiff = rest - rest_global; + if (ndiff != 0) + { + GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff); + } + } + + printMissingInteractionsAtoms(mdlog, cr, top_global, top_local->idef); + write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box); + + std::string errorMessage; + + if (ndiff_tot > 0) + { + errorMessage = + "One or more interactions were assigned to multiple domains of the domain " + "decompostion. Please report this bug."; + } + else + { + errorMessage = gmx::formatString( + "%d of the %d bonded interactions could not be calculated because some atoms " + "involved moved further apart than the multi-body cut-off distance (%g nm) or the " + "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds " + "also see option -ddcheck", + -ndiff_tot, + dd->reverse_top->expectedNumGlobalBondedInteractions(), + dd_cutoff_multibody(dd), + dd_cutoff_twobody(dd)); + } + gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str()); +} + + +void scheduleCheckOfLocalTopology(gmx_domdec_t* dd, const int numBondedInteractionsToReduce) +{ + dd->localTopologyChecker->numBondedInteractionsToReduce = numBondedInteractionsToReduce; + // Note that it's possible for this to still be true from the last + // time it was set, e.g. if repartitioning was triggered before + // global communication that would have acted on the true + // value. This could happen for example when replica exchange took + // place soon after a partition. + dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = true; + // Clear the old global value, which is now invalid + dd->localTopologyChecker->numBondedInteractionsOverAllDomains.reset(); +} + +bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd) +{ + return dd.localTopologyChecker->shouldCheckNumberOfBondedInteractions; +} + +int numBondedInteractions(const gmx_domdec_t& dd) +{ + return dd.localTopologyChecker->numBondedInteractionsToReduce; +} + +void setNumberOfBondedInteractionsOverAllDomains(gmx_domdec_t* dd, int newValue) +{ + GMX_RELEASE_ASSERT(!dd->localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(), + "Cannot set number of bonded interactions because it is already set"); + dd->localTopologyChecker->numBondedInteractionsOverAllDomains.emplace(newValue); +} + +void checkNumberOfBondedInteractions(const gmx::MDLogger& mdlog, + t_commrec* cr, + const gmx_mtop_t& top_global, + const gmx_localtop_t* top_local, + ArrayRef x, + const matrix box) +{ + GMX_RELEASE_ASSERT( + DOMAINDECOMP(cr), + "No need to check number of bonded interactions when not using domain decomposition"); + if (cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions) + { + GMX_RELEASE_ASSERT(cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(), + "The check for the total number of bonded interactions requires the " + "value to have been reduced across all domains"); + if (cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value() + != cr->dd->reverse_top->expectedNumGlobalBondedInteractions()) + { + dd_print_missing_interactions( + mdlog, + cr, + cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value(), + top_global, + top_local, + x, + box); // Does not return + } + // Now that the value is set and the check complete, future + // global communication should not compute the value until + // after the next partitioning. + cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = false; + } +} diff --git a/src/gromacs/domdec/localtopologychecker.h b/src/gromacs/domdec/localtopologychecker.h index 54517dc570..d202622223 100644 --- a/src/gromacs/domdec/localtopologychecker.h +++ b/src/gromacs/domdec/localtopologychecker.h @@ -46,11 +46,28 @@ #include +#include "gromacs/math/vectypes.h" + +struct gmx_domdec_t; +struct gmx_localtop_t; +struct gmx_mtop_t; +struct t_commrec; +struct t_inputrec; + +namespace gmx +{ +class MDLogger; +template +class ArrayRef; +enum class DDBondedChecking : bool; +} // namespace gmx + namespace gmx { struct LocalTopologyChecker { +public: /*! \brief Data to help check local topology construction * * Partitioning could incorrectly miss a bonded interaction. @@ -81,4 +98,36 @@ struct LocalTopologyChecker } // namespace gmx +//! Set that the local topology should be checked via observables reduction +void scheduleCheckOfLocalTopology(gmx_domdec_t* dd, int numBondedInteractionsToReduce); + +/*! \brief Return whether the total bonded interaction count across + * domains should be checked in observables reduction. */ +bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd); + +//! Return the number of bonded interactions in this domain. +int numBondedInteractions(const gmx_domdec_t& dd); + +/*! \brief Set total bonded interaction count across domains. */ +void setNumberOfBondedInteractionsOverAllDomains(gmx_domdec_t* dd, int newValue); + +/*! \brief Check whether bonded interactions are missing from the reverse topology + * produced by domain decomposition. + * + * Must only be called when DD is active. + * + * \param[in] mdlog Logger + * \param[in] cr Communication object + * \param[in] top_global Global topology for the error message + * \param[in] top_local Local topology for the error message + * \param[in] x Position vector for the error message + * \param[in] box Box matrix for the error message + */ +void checkNumberOfBondedInteractions(const gmx::MDLogger& mdlog, + t_commrec* cr, + const gmx_mtop_t& top_global, + const gmx_localtop_t* top_local, + gmx::ArrayRef x, + const matrix box); + #endif diff --git a/src/gromacs/domdec/partition.cpp b/src/gromacs/domdec/partition.cpp index 84713d0aa6..7e4e9f4984 100644 --- a/src/gromacs/domdec/partition.cpp +++ b/src/gromacs/domdec/partition.cpp @@ -60,6 +60,7 @@ #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/ga2la.h" #include "gromacs/domdec/localatomsetmanager.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/domdec/nsgrid.h" #include "gromacs/ewald/pme_pp.h" diff --git a/src/gromacs/domdec/reversetopology.h b/src/gromacs/domdec/reversetopology.h index 0b9359b4b1..f4bc4a631d 100644 --- a/src/gromacs/domdec/reversetopology.h +++ b/src/gromacs/domdec/reversetopology.h @@ -164,4 +164,10 @@ public: std::unique_ptr impl_; }; +/*! \brief Returns the number of atom entries for il in gmx_reverse_top_t */ +int nral_rt(int ftype); + +/*! \brief Return whether interactions of type \p ftype need to be assigned exactly once */ +bool dd_check_ftype(int ftype, const ReverseTopOptions& rtOptions); + #endif diff --git a/src/gromacs/mdlib/stat.cpp b/src/gromacs/mdlib/stat.cpp index 7bfdc6eebe..83ff1665da 100644 --- a/src/gromacs/mdlib/stat.cpp +++ b/src/gromacs/mdlib/stat.cpp @@ -44,6 +44,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/fileio/checkpoint.h" #include "gromacs/fileio/xtcio.h" #include "gromacs/gmxlib/network.h" @@ -373,7 +374,7 @@ void global_stat(const gmx_global_stat& gs, GMX_RELEASE_ASSERT(DOMAINDECOMP(cr), "No need to check number of bonded interactions when not using domain " "decomposition"); - setNumberOfBondedInteractionsOverAllDomains(*cr->dd, gmx::roundToInt(nb)); + setNumberOfBondedInteractionsOverAllDomains(cr->dd, gmx::roundToInt(nb)); } if (!sig.empty()) diff --git a/src/gromacs/mdlib/update_vv.cpp b/src/gromacs/mdlib/update_vv.cpp index 10feea2d1e..e1738b5285 100644 --- a/src/gromacs/mdlib/update_vv.cpp +++ b/src/gromacs/mdlib/update_vv.cpp @@ -44,6 +44,7 @@ #include #include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/math/vec.h" #include "gromacs/mdlib/constr.h" diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 810606ed22..a6a3263c54 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -61,6 +61,7 @@ #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/domdec_struct.h" #include "gromacs/domdec/gpuhaloexchange.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/domdec/partition.h" #include "gromacs/essentialdynamics/edsam.h" diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 24661da9ac..b55f3cf550 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -57,6 +57,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/domdec/partition.h" #include "gromacs/essentialdynamics/edsam.h" diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index fd71258149..d0783eb775 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -56,6 +56,7 @@ #include "gromacs/domdec/domdec.h" #include "gromacs/domdec/domdec_network.h" #include "gromacs/domdec/domdec_struct.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/domdec/mdsetup.h" #include "gromacs/domdec/partition.h" #include "gromacs/essentialdynamics/edsam.h" diff --git a/src/gromacs/modularsimulator/computeglobalselement.cpp b/src/gromacs/modularsimulator/computeglobalselement.cpp index 57b3163fb6..abfa5b3770 100644 --- a/src/gromacs/modularsimulator/computeglobalselement.cpp +++ b/src/gromacs/modularsimulator/computeglobalselement.cpp @@ -44,6 +44,7 @@ #include "computeglobalselement.h" #include "gromacs/domdec/domdec.h" +#include "gromacs/domdec/localtopologychecker.h" #include "gromacs/gmxlib/network.h" #include "gromacs/gmxlib/nrnb.h" #include "gromacs/math/vec.h" -- 2.22.0