2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2021, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
38 * \brief This file defines functions used by the domdec module
39 * while managing the construction, use and error checking for
40 * topologies local to a DD rank.
42 * \ingroup module_domdec
47 #include "gromacs/domdec/localtopologychecker.h"
53 #include "gromacs/domdec/domdec_internal.h"
54 #include "gromacs/domdec/reversetopology.h"
55 #include "gromacs/gmxlib/network.h"
56 #include "gromacs/mdtypes/commrec.h"
57 #include "gromacs/mdtypes/observablesreducer.h"
58 #include "gromacs/mdtypes/state.h"
59 #include "gromacs/topology/idef.h"
60 #include "gromacs/topology/ifunc.h"
61 #include "gromacs/topology/mtop_util.h"
62 #include "gromacs/utility/arrayref.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/logger.h"
65 #include "gromacs/utility/stringstream.h"
66 #include "gromacs/utility/textwriter.h"
73 /*! \brief Checks whether interactions have been assigned for one function type
75 * Loops over a list of interactions in the local topology of one function type
76 * and flags each of the interactions as assigned in the global \p isAssigned list.
77 * Exits with an inconsistency error when an interaction is assigned more than once.
79 static void flagInteractionsForType(const int ftype,
80 const InteractionList& il,
81 const reverse_ilist_t& ril,
82 const Range<int>& atomRange,
83 const int numAtomsPerMolecule,
84 ArrayRef<const int> globalAtomIndices,
85 ArrayRef<int> isAssigned)
87 const int nril_mol = ril.index[numAtomsPerMolecule];
88 const int nral = NRAL(ftype);
90 for (int i = 0; i < il.size(); i += 1 + nral)
92 // ia[0] is the interaction type, ia[1, ...] the atom indices
93 const int* ia = il.iatoms.data() + i;
94 // Extract the global atom index of the first atom in this interaction
95 const int a0 = globalAtomIndices[ia[1]];
96 /* Check if this interaction is in
97 * the currently checked molblock.
99 if (atomRange.isInRange(a0))
101 // The molecule index in the list of this molecule type
102 const int moleculeIndex = (a0 - atomRange.begin()) / numAtomsPerMolecule;
103 const int atomOffset = (a0 - atomRange.begin()) - moleculeIndex * numAtomsPerMolecule;
104 const int globalAtomStartInMolecule = atomRange.begin() + moleculeIndex * numAtomsPerMolecule;
105 int j_mol = ril.index[atomOffset];
107 while (j_mol < ril.index[atomOffset + 1] && !found)
109 const int j = moleculeIndex * nril_mol + j_mol;
110 const int ftype_j = ril.il[j_mol];
111 /* Here we need to check if this interaction has
112 * not already been assigned, since we could have
113 * multiply defined interactions.
115 if (ftype == ftype_j && ia[0] == ril.il[j_mol + 1] && isAssigned[j] == 0)
117 /* Check the atoms */
119 for (int a = 0; a < nral; a++)
121 if (globalAtomIndices[ia[1 + a]]
122 != globalAtomStartInMolecule + ril.il[j_mol + 2 + a])
132 j_mol += 2 + nral_rt(ftype_j);
136 gmx_incons("Some interactions seem to be assigned multiple times");
142 /*! \brief Help print error output when interactions are missing in a molblock
144 * \note This function needs to be called on all ranks (contains a global summation)
146 static std::string printMissingInteractionsMolblock(const t_commrec* cr,
147 const gmx_reverse_top_t& rt,
148 const char* moltypename,
149 const reverse_ilist_t& ril,
150 const Range<int>& atomRange,
151 const int numAtomsPerMolecule,
152 const int numMolecules,
153 const InteractionDefinitions& idef)
155 const int nril_mol = ril.index[numAtomsPerMolecule];
156 std::vector<int> isAssigned(numMolecules * nril_mol, 0);
157 StringOutputStream stream;
158 TextWriter log(&stream);
160 for (int ftype = 0; ftype < F_NRE; ftype++)
162 if (dd_check_ftype(ftype, rt.options()))
164 flagInteractionsForType(
165 ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned);
169 gmx_sumi(isAssigned.size(), isAssigned.data(), cr);
171 const int numMissingToPrint = 10;
173 for (int mol = 0; mol < numMolecules; mol++)
176 while (j_mol < nril_mol)
178 int ftype = ril.il[j_mol];
179 int nral = NRAL(ftype);
180 int j = mol * nril_mol + j_mol;
181 if (isAssigned[j] == 0 && !(interaction_function[ftype].flags & IF_VSITE))
183 if (DDMASTER(cr->dd))
187 log.writeLineFormatted("Molecule type '%s'", moltypename);
188 log.writeLineFormatted(
189 "the first %d missing interactions, except for exclusions:",
192 log.writeStringFormatted("%20s atoms", interaction_function[ftype].longname);
194 for (; a < nral; a++)
196 log.writeStringFormatted("%5d", ril.il[j_mol + 2 + a] + 1);
200 log.writeString(" ");
203 log.writeString(" global");
204 for (int a = 0; a < nral; a++)
206 log.writeStringFormatted("%6d",
207 atomRange.begin() + mol * numAtomsPerMolecule
208 + ril.il[j_mol + 2 + a] + 1);
210 log.ensureLineBreak();
213 if (i >= numMissingToPrint)
218 j_mol += 2 + nral_rt(ftype);
222 return stream.toString();
225 /*! \brief Help print error output when interactions are missing */
226 static void printMissingInteractionsAtoms(const MDLogger& mdlog,
228 const gmx_mtop_t& mtop,
229 const InteractionDefinitions& idef)
231 const gmx_reverse_top_t& rt = *cr->dd->reverse_top;
233 /* Print the atoms in the missing interactions per molblock */
235 for (const gmx_molblock_t& molb : mtop.molblock)
237 const gmx_moltype_t& moltype = mtop.moltype[molb.type];
238 const int a_start = a_end;
239 a_end = a_start + molb.nmol * moltype.atoms.nr;
240 const Range<int> atomRange(a_start, a_end);
242 auto warning = printMissingInteractionsMolblock(
246 cr->dd->reverse_top->interactionListForMoleculeType(molb.type),
252 GMX_LOG(mdlog.warning).appendText(warning);
256 /*! \brief Print error output when interactions are missing */
257 [[noreturn]] static void dd_print_missing_interactions(const MDLogger& mdlog,
259 const int numBondedInteractionsOverAllDomains,
260 const int expectedNumGlobalBondedInteractions,
261 const gmx_mtop_t& top_global,
262 const gmx_localtop_t& top_local,
263 ArrayRef<const RVec> x,
267 gmx_domdec_t* dd = cr->dd;
269 GMX_LOG(mdlog.warning)
271 "Not all bonded interactions have been properly assigned to the domain "
272 "decomposition cells");
274 const int ndiff_tot = numBondedInteractionsOverAllDomains - expectedNumGlobalBondedInteractions;
276 for (int ftype = 0; ftype < F_NRE; ftype++)
278 const int nral = NRAL(ftype);
279 cl[ftype] = top_local.idef.il[ftype].size() / (1 + nral);
282 gmx_sumi(F_NRE, cl, cr);
286 GMX_LOG(mdlog.warning).appendText("A list of missing interactions:");
287 int rest_global = expectedNumGlobalBondedInteractions;
288 int rest = numBondedInteractionsOverAllDomains;
289 for (int ftype = 0; ftype < F_NRE; ftype++)
291 /* In the reverse and local top all constraints are merged
292 * into F_CONSTR. So in the if statement we skip F_CONSTRNC
293 * and add these constraints when doing F_CONSTR.
295 if (dd_check_ftype(ftype, dd->reverse_top->options()) && ftype != F_CONSTRNC)
297 int n = gmx_mtop_ftype_count(top_global, ftype);
298 if (ftype == F_CONSTR)
300 n += gmx_mtop_ftype_count(top_global, F_CONSTRNC);
302 int ndiff = cl[ftype] - n;
305 GMX_LOG(mdlog.warning)
306 .appendTextFormatted("%20s of %6d missing %6d",
307 interaction_function[ftype].longname,
316 int ndiff = rest - rest_global;
319 GMX_LOG(mdlog.warning).appendTextFormatted("%20s of %6d missing %6d", "exclusions", rest_global, -ndiff);
323 printMissingInteractionsAtoms(mdlog, cr, top_global, top_local.idef);
324 write_dd_pdb("dd_dump_err", 0, "dump", top_global, cr, -1, as_rvec_array(x.data()), box);
326 std::string errorMessage;
331 "One or more interactions were assigned to multiple domains of the domain "
332 "decompostion. Please report this bug.";
336 errorMessage = formatString(
337 "%d of the %d bonded interactions could not be calculated because some atoms "
338 "involved moved further apart than the multi-body cut-off distance (%g nm) or the "
339 "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds "
340 "also see option -ddcheck",
342 expectedNumGlobalBondedInteractions,
343 dd_cutoff_multibody(dd),
344 dd_cutoff_twobody(dd));
346 gmx_fatal_collective(FARGS, cr->mpi_comm_mygroup, MASTER(cr), "%s", errorMessage.c_str());
349 /*! \brief Data to help check local topology construction
351 * Partitioning could incorrectly miss a bonded interaction.
352 * However, checking for that requires a global communication
353 * stage, which does not otherwise happen during partitioning. So,
354 * for performance, we do that alongside the first global energy
355 * reduction after a new DD is made. These variables handle
356 * whether the check happens, its input for this domain, output
357 * across all domains, and the expected value it should match. */
358 class LocalTopologyChecker::Impl
362 Impl(const MDLogger& mdlog,
364 const gmx_mtop_t& mtop,
365 const gmx_localtop_t& localTopology,
366 const t_state& localState,
367 bool useUpdateGroups);
368 //! Objects used when reporting that interactions are missing
371 const MDLogger& mdlog_;
372 //! Communication record
373 const t_commrec* cr_;
374 //! Global system topology
375 const gmx_mtop_t& mtop_;
377 const gmx_localtop_t& localTopology_;
379 const t_state& localState_;
382 /*! \brief View used for computing the global number of bonded interactions.
384 * Can be written any time, but that is only useful when followed
385 * by a call of the callbackToRequireReduction. Useful to read
386 * only from the callback that the ObservablesReducer will later
387 * make after reduction. */
388 gmx::ArrayRef<double> reductionBuffer_;
389 /*! \brief Callback used after repartitioning to require reduction
390 * of numBondedInteractionsToReduce so that the total number of
391 * bonded interactions can be checked. */
392 gmx::ObservablesReducerBuilder::CallbackToRequireReduction callbackToRequireReduction_;
393 /*! \brief The expected number of global bonded interactions from the system topology */
394 int expectedNumGlobalBondedInteractions_;
398 /*! \brief Compute the total bonded interaction count
400 * \param[in] mtop The global system topology
401 * \param[in] useUpdateGroups Whether update groups are in use
403 * When using domain decomposition without update groups,
404 * constraint-type interactions can be split across domains, and so we
405 * do not consider them in this correctness check. Otherwise, we
408 static int computeExpectedNumGlobalBondedInteractions(const gmx_mtop_t& mtop, const bool useUpdateGroups)
410 int expectedNumGlobalBondedInteractions = gmx_mtop_interaction_count(mtop, IF_BOND);
413 expectedNumGlobalBondedInteractions += gmx_mtop_interaction_count(mtop, IF_CONSTRAINT);
415 return expectedNumGlobalBondedInteractions;
418 LocalTopologyChecker::Impl::Impl(const MDLogger& mdlog,
420 const gmx_mtop_t& mtop,
421 const gmx_localtop_t& localTopology,
422 const t_state& localState,
423 bool useUpdateGroups) :
427 localTopology_(localTopology),
428 localState_(localState),
429 expectedNumGlobalBondedInteractions_(computeExpectedNumGlobalBondedInteractions(mtop, useUpdateGroups))
433 LocalTopologyChecker::LocalTopologyChecker(const MDLogger& mdlog,
435 const gmx_mtop_t& mtop,
436 const gmx_localtop_t& localTopology,
437 const t_state& localState,
438 const bool useUpdateGroups,
439 ObservablesReducerBuilder* observablesReducerBuilder) :
440 impl_(std::make_unique<Impl>(mdlog, cr, mtop, localTopology, localState, useUpdateGroups))
442 Impl* impl = impl_.get();
443 ObservablesReducerBuilder::CallbackFromBuilder callbackFromBuilder =
444 [impl](ObservablesReducerBuilder::CallbackToRequireReduction c, gmx::ArrayRef<double> v) {
445 impl->callbackToRequireReduction_ = std::move(c);
446 impl->reductionBuffer_ = v;
449 // Make the callback that runs afer reduction.
450 ObservablesReducerBuilder::CallbackAfterReduction callbackAfterReduction = [impl](gmx::Step /*step*/) {
451 // Get the total after reduction
452 int numTotalBondedInteractionsFound = impl->reductionBuffer_[0];
453 if (numTotalBondedInteractionsFound != impl->expectedNumGlobalBondedInteractions_)
455 // Give error and exit
456 dd_print_missing_interactions(impl->mdlog_,
458 numTotalBondedInteractionsFound,
459 impl->expectedNumGlobalBondedInteractions_,
461 impl->localTopology_,
463 impl->localState_.box); // Does not return
467 observablesReducerBuilder->addSubscriber(
468 1, std::move(callbackFromBuilder), std::move(callbackAfterReduction));
471 LocalTopologyChecker::~LocalTopologyChecker() = default;
473 LocalTopologyChecker::LocalTopologyChecker(LocalTopologyChecker&&) noexcept = default;
475 LocalTopologyChecker& LocalTopologyChecker::operator=(LocalTopologyChecker&& other) noexcept
477 impl_ = std::move(other.impl_);
481 void LocalTopologyChecker::scheduleCheckOfLocalTopology(const int numBondedInteractionsToReduce)
483 // When we have a single domain, we don't need to reduce and we algorithmically can not miss
484 // any interactions, so we can assert here.
485 if (!havePPDomainDecomposition(impl_->cr_))
487 GMX_RELEASE_ASSERT(numBondedInteractionsToReduce == impl_->expectedNumGlobalBondedInteractions_,
488 "With a single domain the number of assigned bonded interactions should "
489 "always match the global number");
493 // Fill the reduction buffer with the value from this domain to reduce
494 impl_->reductionBuffer_[0] = double(numBondedInteractionsToReduce);
496 // Pass the post-reduction callback to the ObservablesReducer via
497 // the callback it gave us for the purpose.
499 // Note that it's possible that the callbackAfterReduction is already
500 // outstanding, e.g. if repartitioning was triggered before
501 // observables were reduced. This could happen for example when
502 // replica exchange took place soon after a partition. If so, the
503 // callback will be called again. So long as there is no race
504 // between the calls to this function and the calls to
505 // ObservablesReducer for reduction, this will work correctly. It
506 // could be made safer e.g. with checks against duplicate
507 // callbacks, but there is no problem to solve.
509 // There is no need to check the return value from this callback,
510 // as it is not an error to request reduction at a future step.
511 impl_->callbackToRequireReduction_(ReductionRequirement::Eventually);