From 86b442ebef2ca6bd30d16254aca7aae86ea522ad Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Thu, 13 May 2021 07:36:56 +0000 Subject: [PATCH] Decouple update group handling from domain decomposition module --- src/gromacs/domdec/builder.h | 20 ++- src/gromacs/domdec/domdec.cpp | 154 ++++++++----------- src/gromacs/domdec/domdec.h | 6 - src/gromacs/domdec/domdec_internal.h | 2 +- src/gromacs/mdlib/constr.cpp | 9 +- src/gromacs/mdlib/constr.h | 2 +- src/gromacs/mdlib/tests/updategroups.cpp | 116 +++++++++++++- src/gromacs/mdlib/updategroups.cpp | 82 +++++++++- src/gromacs/mdlib/updategroups.h | 99 ++++++++++-- src/gromacs/mdrun/runner.cpp | 29 ++-- src/testutils/include/testutils/loggertest.h | 6 + 11 files changed, 388 insertions(+), 137 deletions(-) diff --git a/src/gromacs/domdec/builder.h b/src/gromacs/domdec/builder.h index 0c6e764664..1ca7408dd3 100644 --- a/src/gromacs/domdec/builder.h +++ b/src/gromacs/domdec/builder.h @@ -60,6 +60,7 @@ namespace gmx { class MDLogger; class LocalAtomSetManager; +class RangePartitioning; struct DomdecOptions; struct MdrunOptions; @@ -76,14 +77,17 @@ class DomainDecompositionBuilder { public: //! Constructor - DomainDecompositionBuilder(const MDLogger& mdlog, - t_commrec* cr, - const DomdecOptions& options, - const MdrunOptions& mdrunOptions, - const gmx_mtop_t& mtop, - const t_inputrec& ir, - const matrix box, - ArrayRef xGlobal); + DomainDecompositionBuilder(const MDLogger& mdlog, + t_commrec* cr, + const DomdecOptions& options, + const MdrunOptions& mdrunOptions, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const matrix box, + ArrayRef updateGroupingPerMoleculeType, + bool useUpdateGroups, + real maxUpdateGroupRadius, + ArrayRef xGlobal); //! Destructor ~DomainDecompositionBuilder(); //! Build the resulting DD manager diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index a9304bbab6..93a5f5eeb6 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -116,9 +116,11 @@ #include "utility.h" // TODO remove this when moving domdec into gmx namespace +using gmx::ArrayRef; using gmx::DdRankOrder; using gmx::DlbOption; using gmx::DomdecOptions; +using gmx::RangePartitioning; static const char* enumValueToString(DlbState enumValue) { @@ -206,12 +208,6 @@ int ddglatnr(const gmx_domdec_t* dd, int i) return atnr; } -gmx::ArrayRef getUpdateGroupingsPerMoleculeType(const gmx_domdec_t& dd) -{ - GMX_RELEASE_ASSERT(dd.comm, "Need a valid dd.comm"); - return dd.comm->systemInfo.updateGroupingsPerMoleculeType; -} - void dd_store_state(const gmx_domdec_t& dd, t_state* state) { if (state->ddp_count != dd.ddp_count) @@ -982,11 +978,6 @@ int dd_pme_maxshift_y(const gmx_domdec_t& dd) } } -bool ddMayHaveSplitConstraints(const gmx_domdec_t& dd) -{ - return dd.comm->systemInfo.mayHaveSplitConstraints; -} - bool ddUsesUpdateGroups(const gmx_domdec_t& dd) { return dd.comm->systemInfo.useUpdateGroups; @@ -1916,32 +1907,16 @@ static gmx_domdec_comm_t* init_dd_comm() return comm; } -/* Returns whether mtop contains constraints and/or vsites */ -static bool systemHasConstraintsOrVsites(const gmx_mtop_t& mtop) +static void setupUpdateGroups(const gmx::MDLogger& mdlog, + const gmx_mtop_t& mtop, + ArrayRef updateGroupingsPerMoleculeType, + const bool useUpdateGroups, + const real maxUpdateGroupRadius, + DDSystemInfo* systemInfo) { - return std::any_of(IListRange(mtop).begin(), IListRange(mtop).end(), [](const auto ilist) { - return !extractILists(ilist.list(), IF_CONSTRAINT | IF_VSITE).empty(); - }); -} - -static void setupUpdateGroups(const gmx::MDLogger& mdlog, - const gmx_mtop_t& mtop, - const t_inputrec& inputrec, - const real cutoffMargin, - DDSystemInfo* systemInfo) -{ - /* When we have constraints and/or vsites, it is beneficial to use - * update groups (when possible) to allow independent update of groups. - */ - if (!systemHasConstraintsOrVsites(mtop)) - { - /* No constraints or vsites, atoms can be updated independently */ - return; - } - - systemInfo->updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop); - systemInfo->useUpdateGroups = (!systemInfo->updateGroupingsPerMoleculeType.empty() - && getenv("GMX_NO_UPDATEGROUPS") == nullptr); + systemInfo->updateGroupingsPerMoleculeType = updateGroupingsPerMoleculeType; + systemInfo->useUpdateGroups = useUpdateGroups; + systemInfo->maxUpdateGroupRadius = maxUpdateGroupRadius; if (systemInfo->useUpdateGroups) { @@ -1952,32 +1927,13 @@ static void setupUpdateGroups(const gmx::MDLogger& mdlog, * systemInfo->updateGroupingsPerMoleculeType[molblock.type].numBlocks(); } - systemInfo->maxUpdateGroupRadius = computeMaxUpdateGroupRadius( - mtop, systemInfo->updateGroupingsPerMoleculeType, maxReferenceTemperature(inputrec)); - - /* To use update groups, the large domain-to-domain cutoff distance - * should be compatible with the box size. - */ - systemInfo->useUpdateGroups = (atomToAtomIntoDomainToDomainCutoff(*systemInfo, 0) < cutoffMargin); - - if (systemInfo->useUpdateGroups) - { - GMX_LOG(mdlog.info) - .appendTextFormatted( - "Using update groups, nr %d, average size %.1f atoms, max. radius %.3f " - "nm\n", - numUpdateGroups, - mtop.natoms / static_cast(numUpdateGroups), - systemInfo->maxUpdateGroupRadius); - } - else - { - GMX_LOG(mdlog.info) - .appendTextFormatted( - "The combination of rlist and box size prohibits the use of update " - "groups\n"); - systemInfo->updateGroupingsPerMoleculeType.clear(); - } + GMX_LOG(mdlog.info) + .appendTextFormatted( + "Using update groups, nr %d, average size %.1f atoms, max. radius %.3f " + "nm\n", + numUpdateGroups, + mtop.natoms / static_cast(numUpdateGroups), + systemInfo->maxUpdateGroupRadius); } } @@ -2021,26 +1977,24 @@ static bool moleculesAreAlwaysWhole(const gmx_mtop_t& } /*! \brief Generate the simulation system information */ -static DDSystemInfo getSystemInfo(const gmx::MDLogger& mdlog, - DDRole ddRole, - MPI_Comm communicator, - const DomdecOptions& options, - const gmx_mtop_t& mtop, - const t_inputrec& ir, - const matrix box, - gmx::ArrayRef xGlobal) +static DDSystemInfo getSystemInfo(const gmx::MDLogger& mdlog, + DDRole ddRole, + MPI_Comm communicator, + const DomdecOptions& options, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const matrix box, + ArrayRef updateGroupingPerMoleculeType, + const bool useUpdateGroups, + const real maxUpdateGroupRadius, + gmx::ArrayRef xGlobal) { const real tenPercentMargin = 1.1; DDSystemInfo systemInfo; - /* We need to decide on update groups early, as this affects communication distances */ - systemInfo.useUpdateGroups = false; - if (ir.cutoff_scheme == CutoffScheme::Verlet) - { - real cutoffMargin = std::sqrt(max_cutoff2(ir.pbcType, box)) - ir.rlist; - setupUpdateGroups(mdlog, mtop, ir, cutoffMargin, &systemInfo); - } + setupUpdateGroups( + mdlog, mtop, updateGroupingPerMoleculeType, useUpdateGroups, maxUpdateGroupRadius, &systemInfo); systemInfo.moleculesAreAlwaysWhole = moleculesAreAlwaysWhole( mtop, systemInfo.useUpdateGroups, systemInfo.updateGroupingsPerMoleculeType); @@ -2873,14 +2827,17 @@ class DomainDecompositionBuilder::Impl { public: //! Constructor - Impl(const MDLogger& mdlog, - t_commrec* cr, - const DomdecOptions& options, - const MdrunOptions& mdrunOptions, - const gmx_mtop_t& mtop, - const t_inputrec& ir, - const matrix box, - ArrayRef xGlobal); + Impl(const MDLogger& mdlog, + t_commrec* cr, + const DomdecOptions& options, + const MdrunOptions& mdrunOptions, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const matrix box, + ArrayRef updateGroupingPerMoleculeType, + bool useUpdateGroups, + real maxUpdateGroupRadius, + ArrayRef xGlobal); //! Build the resulting DD manager gmx_domdec_t* build(LocalAtomSetManager* atomSets); @@ -2920,14 +2877,17 @@ public: //! } }; -DomainDecompositionBuilder::Impl::Impl(const MDLogger& mdlog, - t_commrec* cr, - const DomdecOptions& options, - const MdrunOptions& mdrunOptions, - const gmx_mtop_t& mtop, - const t_inputrec& ir, - const matrix box, - ArrayRef xGlobal) : +DomainDecompositionBuilder::Impl::Impl(const MDLogger& mdlog, + t_commrec* cr, + const DomdecOptions& options, + const MdrunOptions& mdrunOptions, + const gmx_mtop_t& mtop, + const t_inputrec& ir, + const matrix box, + ArrayRef updateGroupingPerMoleculeType, + const bool useUpdateGroups, + const real maxUpdateGroupRadius, + ArrayRef xGlobal) : mdlog_(mdlog), cr_(cr), options_(options), mtop_(mtop), ir_(ir) { GMX_LOG(mdlog_.info).appendTextFormatted("\nInitializing Domain Decomposition on %d ranks", cr_->sizeOfDefaultCommunicator); @@ -2947,6 +2907,9 @@ DomainDecompositionBuilder::Impl::Impl(const MDLogger& mdlog, mtop_, ir_, box, + updateGroupingPerMoleculeType, + useUpdateGroups, + maxUpdateGroupRadius, xGlobal); const int numRanksRequested = cr_->sizeOfDefaultCommunicator; @@ -3042,8 +3005,11 @@ DomainDecompositionBuilder::DomainDecompositionBuilder(const MDLogger& mdlo const gmx_mtop_t& mtop, const t_inputrec& ir, const matrix box, + ArrayRef updateGroupingPerMoleculeType, + const bool useUpdateGroups, + const real maxUpdateGroupRadius, ArrayRef xGlobal) : - impl_(new Impl(mdlog, cr, options, mdrunOptions, mtop, ir, box, xGlobal)) + impl_(new Impl(mdlog, cr, options, mdrunOptions, mtop, ir, box, updateGroupingPerMoleculeType, useUpdateGroups, maxUpdateGroupRadius, xGlobal)) { } diff --git a/src/gromacs/domdec/domdec.h b/src/gromacs/domdec/domdec.h index f60c429b6a..e36334fdc6 100644 --- a/src/gromacs/domdec/domdec.h +++ b/src/gromacs/domdec/domdec.h @@ -103,9 +103,6 @@ enum class DDBondedChecking : bool; */ int ddglatnr(const gmx_domdec_t* dd, int i); -/*! \brief Returns a list of update group partitioning for each molecule type or empty when update groups are not used */ -gmx::ArrayRef getUpdateGroupingsPerMoleculeType(const gmx_domdec_t& dd); - /*! \brief Store the global cg indices of the home cgs in state, * * This means it can be reset, even after a new DD partitioning. @@ -149,9 +146,6 @@ int dd_pme_maxshift_x(const gmx_domdec_t& dd); /*! \brief Returns the maximum shift for coordinate communication in PME, dim y */ int dd_pme_maxshift_y(const gmx_domdec_t& dd); -/*! \brief Return whether constraints, not including settles, may cross domain boundaries */ -bool ddMayHaveSplitConstraints(const gmx_domdec_t& dd); - /*! \brief Return whether update groups are used */ bool ddUsesUpdateGroups(const gmx_domdec_t& dd); diff --git a/src/gromacs/domdec/domdec_internal.h b/src/gromacs/domdec/domdec_internal.h index ddfa6b847f..a112f8e256 100644 --- a/src/gromacs/domdec/domdec_internal.h +++ b/src/gromacs/domdec/domdec_internal.h @@ -423,7 +423,7 @@ struct DDSystemInfo //! True when update groups are used bool useUpdateGroups = false; //! Update atom grouping for each molecule type - std::vector updateGroupingsPerMoleculeType; + gmx::ArrayRef updateGroupingsPerMoleculeType; //! The maximum radius over all update groups real maxUpdateGroupRadius; diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp index 0b8d94234c..ceee50f658 100644 --- a/src/gromacs/mdlib/constr.cpp +++ b/src/gromacs/mdlib/constr.cpp @@ -1097,7 +1097,7 @@ Constraints::Impl::Impl(const gmx_mtop_t& mtop_p, pull_t* pull_work, FILE* log_p, const t_commrec* cr_p, - const bool mayHaveSplitConstraints, + const bool useUpdateGroups, const gmx_multisim_t* ms_p, t_nrnb* nrnb_p, gmx_wallcycle* wcycle_p, @@ -1154,8 +1154,11 @@ Constraints::Impl::Impl(const gmx_mtop_t& mtop_p, } } - // When there are multiple PP domains and update groups are not in use, - // the constraints might be split across them. + // When there are multiple PP domains and update groups are + // not in use, the constraints might be split across the + // domains, needing particular handling. + const bool mayHaveSplitConstraints = DOMAINDECOMP(cr) && !useUpdateGroups; + if (ir.eConstrAlg == ConstraintAlgorithm::Lincs) { lincsd = init_lincs( diff --git a/src/gromacs/mdlib/constr.h b/src/gromacs/mdlib/constr.h index f5bc84ea2c..14cd95193a 100644 --- a/src/gromacs/mdlib/constr.h +++ b/src/gromacs/mdlib/constr.h @@ -107,7 +107,7 @@ private: pull_t* pull_work, FILE* log, const t_commrec* cr, - bool haveHaveSplitConstraints, + bool useUpdateGroups, const gmx_multisim_t* ms, t_nrnb* nrnb, gmx_wallcycle* wcycle, diff --git a/src/gromacs/mdlib/tests/updategroups.cpp b/src/gromacs/mdlib/tests/updategroups.cpp index 4b75b16a99..8934d0b74b 100644 --- a/src/gromacs/mdlib/tests/updategroups.cpp +++ b/src/gromacs/mdlib/tests/updategroups.cpp @@ -47,6 +47,7 @@ #include "gromacs/topology/topology.h" +#include "testutils/loggertest.h" #include "testutils/testasserts.h" namespace gmx @@ -57,6 +58,17 @@ namespace /* TODO: Actually initialize moltype.atoms.atom when this is converted to C++ */ +/*! \brief Returns a flexible ethane united-atom molecule */ +gmx_moltype_t flexibleEthaneUA() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 2; + moltype.ilist[F_BONDS].iatoms = { 0, 0, 1 }; + + return moltype; +} + /*! \brief Returns an ethane united-atom molecule */ gmx_moltype_t ethaneUA() { @@ -92,7 +104,7 @@ gmx_moltype_t ethane() return moltype; } -/*! \brief Returns a butane united-atom molecule */ +/*! \brief Returns a butane fully-constrained united-atom molecule */ gmx_moltype_t butaneUA() { gmx_moltype_t moltype = {}; @@ -153,6 +165,8 @@ public: gmx_mtop_t mtop_; //! Default temperature for tests real temperature_ = 298; + //! Logger to use in tests + test::LoggerTestHelper logHelper_; }; TEST_F(UpdateGroupsTest, WithEthaneUA) @@ -171,6 +185,10 @@ TEST_F(UpdateGroupsTest, WithEthaneUA) real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); EXPECT_FLOAT_EQ(maxRadius, 0.3 / 2); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } TEST_F(UpdateGroupsTest, WithMethane) @@ -189,7 +207,12 @@ TEST_F(UpdateGroupsTest, WithMethane) real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); EXPECT_FLOAT_EQ(maxRadius, 0.14); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } + TEST_F(UpdateGroupsTest, WithEthane) { mtop_.moltype.emplace_back(ethane()); @@ -208,6 +231,10 @@ TEST_F(UpdateGroupsTest, WithEthane) real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); EXPECT_FLOAT_EQ(maxRadius, 0.094746813); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } TEST_F(UpdateGroupsTest, CheckRadiusCalculationAtDifferentTemperaturesWithEthane) @@ -239,7 +266,7 @@ TEST_F(UpdateGroupsTest, CheckRadiusCalculationAtDifferentTemperaturesWithEthane EXPECT_FLOAT_EQ(maxRadius, 0.125); } -TEST_F(UpdateGroupsTest, WithButaneUA) +TEST_F(UpdateGroupsTest, WithButaneUALogsThatUnsuitableForUpdateGroups) { mtop_.moltype.emplace_back(butaneUA()); { @@ -251,6 +278,15 @@ TEST_F(UpdateGroupsTest, WithButaneUA) auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_); EXPECT_EQ(updateGroupingsPerMoleculeType.size(), 0); + + real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); + EXPECT_FLOAT_EQ(maxRadius, 0.0); + + logHelper_.expectEntryMatchingRegex( + MDLogger::LogLevel::Info, + "At least one moleculetype does not conform to the requirements"); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } TEST_F(UpdateGroupsTest, WithWaterThreeSite) @@ -269,6 +305,10 @@ TEST_F(UpdateGroupsTest, WithWaterThreeSite) real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); EXPECT_FLOAT_EQ(maxRadius, 0.083887339); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } // Tests update group with virtual site @@ -287,17 +327,38 @@ TEST_F(UpdateGroupsTest, WithWaterFourSite) ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1); EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1); + + real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); + EXPECT_FLOAT_EQ(maxRadius, 0.083887339); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } TEST_F(UpdateGroupsTest, WithFourAtomsWithSettle) { mtop_.moltype.emplace_back(waterThreeSite()); mtop_.moltype.back().atoms.nr = 4; + { + t_iparams iparams[2]; + iparams[0].settle = { 0.1, 0.1633 }; + iparams[1].vsite = { 0.128, 0.128 }; + mtop_.ffparams.iparams.push_back(iparams[0]); + mtop_.ffparams.iparams.push_back(iparams[1]); + } auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_); ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1); EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2); + + real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); + EXPECT_FLOAT_EQ(maxRadius, 0.083887339); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } // Tests groups with two constraints and an angle potential @@ -319,6 +380,10 @@ TEST_F(UpdateGroupsTest, WithWaterFlexAngle) real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); EXPECT_FLOAT_EQ(maxRadius, 0.090824135); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); } TEST_F(UpdateGroupsTest, CheckRadiusCalculationAtDifferentTemperaturesWithWaterFlexAngle) @@ -367,8 +432,55 @@ TEST_F(UpdateGroupsTest, WithTwoMoltypes) ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 2); EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 1); EXPECT_EQ(updateGroupingsPerMoleculeType[1].numBlocks(), 1); + + real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); + EXPECT_FLOAT_EQ(maxRadius, 0.14); + + logHelper_.expectNoEntries(MDLogger::LogLevel::Info); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, true, 1e6_real); +} + +TEST_F(UpdateGroupsTest, LogsWhenSizesAreInvalid) +{ + mtop_.moltype.emplace_back(methane()); + { + t_iparams iparams; + iparams.constr = { 0.1, 0.1 }; + mtop_.ffparams.iparams.push_back(iparams); + } + + auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_); + logHelper_.expectEntryMatchingRegex(MDLogger::LogLevel::Info, + "The combination of rlist and box size prohibits"); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), 1e9_real, true, true, 1e6_real); } +TEST_F(UpdateGroupsTest, LogsWhenUpdateGroupsAreNotUseful) +{ + mtop_.moltype.emplace_back(flexibleEthaneUA()); + { + t_iparams iparams; + iparams.harmonic = { 0.1, 10.0, 0.1, 10.0 }; + mtop_.ffparams.iparams.push_back(iparams); + } + + auto updateGroupingsPerMoleculeType = gmx::makeUpdateGroupingsPerMoleculeType(mtop_); + + ASSERT_EQ(updateGroupingsPerMoleculeType.size(), 1); + EXPECT_EQ(updateGroupingsPerMoleculeType[0].numBlocks(), 2); + + real maxRadius = computeMaxUpdateGroupRadius(mtop_, updateGroupingsPerMoleculeType, temperature_); + EXPECT_FLOAT_EQ(maxRadius, 0); + + logHelper_.expectEntryMatchingRegex(MDLogger::LogLevel::Info, + "No constraints or virtual sites are in use"); + UpdateGroups updateGroups = makeUpdateGroups( + logHelper_.logger(), std::move(updateGroupingsPerMoleculeType), maxRadius, true, false, 1e6_real); +} + + } // namespace } // namespace gmx diff --git a/src/gromacs/mdlib/updategroups.cpp b/src/gromacs/mdlib/updategroups.cpp index fe040b6a98..64a0a724dd 100644 --- a/src/gromacs/mdlib/updategroups.cpp +++ b/src/gromacs/mdlib/updategroups.cpp @@ -52,11 +52,14 @@ #include "gromacs/math/units.h" #include "gromacs/math/utilities.h" #include "gromacs/mdlib/constr.h" -#include "gromacs/mdtypes/inputrec.h" +#include "gromacs/pbcutil/pbc.h" #include "gromacs/topology/idef.h" #include "gromacs/topology/ifunc.h" +#include "gromacs/topology/mtop_util.h" #include "gromacs/topology/topology.h" #include "gromacs/utility/listoflists.h" +#include "gromacs/utility/logger.h" +#include "gromacs/utility/messagestringcollector.h" namespace gmx { @@ -751,4 +754,81 @@ real computeMaxUpdateGroupRadius(const gmx_mtop_t& mtop, return maxRadius; } +real computeCutoffMargin(PbcType pbcType, matrix box, const real rlist) +{ + return std::sqrt(max_cutoff2(pbcType, box)) - rlist; +} + +UpdateGroups::UpdateGroups(std::vector&& updateGroupingPerMoleculeType, + const real maxUpdateGroupRadius) : + useUpdateGroups_(true), + updateGroupingPerMoleculeType_(std::move(updateGroupingPerMoleculeType)), + maxUpdateGroupRadius_(maxUpdateGroupRadius) +{ +} + +bool systemHasConstraintsOrVsites(const gmx_mtop_t& mtop) +{ + IListRange ilistRange(mtop); + return std::any_of(ilistRange.begin(), ilistRange.end(), [](const auto& ilists) { + return !extractILists(ilists.list(), IF_CONSTRAINT | IF_VSITE).empty(); + }); +} + +UpdateGroups makeUpdateGroups(const gmx::MDLogger& mdlog, + std::vector&& updateGroupingPerMoleculeType, + const real maxUpdateGroupRadius, + const bool useDomainDecomposition, + const bool systemHasConstraintsOrVsites, + const real cutoffMargin) +{ + MessageStringCollector messages; + + messages.startContext("When checking whether update groups are usable:"); + + if (!useDomainDecomposition) + { + messages.append( + "Domain decomposition is not active, so there is no need for update groups"); + } + + if (!systemHasConstraintsOrVsites) + { + messages.append( + "No constraints or virtual sites are in use, so it is best not to use update " + "groups"); + } + + if (updateGroupingPerMoleculeType.empty()) + { + messages.append( + "At least one moleculetype does not conform to the requirements for using update " + "groups"); + } + + if (getenv("GMX_NO_UPDATEGROUPS") != nullptr) + { + messages.append( + "Environment variable GMX_NO_UPDATEGROUPS prohibited the use of update groups"); + } + + // To use update groups, the large domain-to-domain cutoff + // distance should be compatible with the box size. + if (2 * maxUpdateGroupRadius >= cutoffMargin) + { + messages.append("The combination of rlist and box size prohibits the use of update groups"); + } + + if (!messages.isEmpty()) + { + // Log why we can't use update groups + GMX_LOG(mdlog.info).appendText(messages.toString()); + return UpdateGroups(); + } + + // Success! + return UpdateGroups(std::move(updateGroupingPerMoleculeType), maxUpdateGroupRadius); +} + + } // namespace gmx diff --git a/src/gromacs/mdlib/updategroups.h b/src/gromacs/mdlib/updategroups.h index 6daa2e22a8..150d2e527b 100644 --- a/src/gromacs/mdlib/updategroups.h +++ b/src/gromacs/mdlib/updategroups.h @@ -45,14 +45,16 @@ #include +#include "gromacs/math/vectypes.h" #include "gromacs/utility/arrayref.h" #include "gromacs/utility/real.h" +enum class PbcType : int; struct gmx_mtop_t; -struct t_inputrec; namespace gmx { +class MDLogger; class RangePartitioning; /*! \brief Returns a vector with update groups for each moleculetype in \p mtop @@ -62,28 +64,101 @@ class RangePartitioning; * the restrictions of update groups, e.g. more than two constraints in a row. * * Currently valid update groups are: - * a single atom which is not a virtual site and does not have constraints; - * or a group of atoms where all virtual sites are constructed from atoms - * within the group and at least one non-vsite atom is constrained to - * all other non-vsite atoms. + * - a single atom which is not a virtual site and does not have constraints; + * - or a group of atoms where all virtual sites are constructed from atoms + * within the group and at least one non-vsite atom is constrained to + * all other non-vsite atoms. + * * To have update groups, all virtual sites should be linear 2 or 3 atom * constructions with coefficients >= 0 and sum of coefficients <= 1. * + * This vector is generally consumed in constructing an UpdateGroups object. + * * \param[in] mtop The system topology */ std::vector makeUpdateGroupingsPerMoleculeType(const gmx_mtop_t& mtop); /*! \brief Returns the maximum update group radius * - * \note When \p updateGroupingsPerMoleculeType is empty, 0 is returned. + * \note When \p updateGroups is empty, 0 is returned. + * + * \param[in] mtop The system topology + * \param[in] updateGroupingPerMoleculeType List of update group, size should match the + * number of moltypes in \p mtop or be 0 + * \param[in] temperature The maximum reference temperature, pass -1 + * when unknown or not applicable + */ +real computeMaxUpdateGroupRadius(const gmx_mtop_t& mtop, + ArrayRef updateGroupingPerMoleculeType, + real temperature); + +/*! \brief Return the margin required for successful domain decomposition + * + * \param[in] pbcType The PBC type in use + * \param[in] box The box in use + * \param[in] rlist The list size in use + */ +real computeCutoffMargin(PbcType pbcType, matrix box, real rlist); + +/*! \brief Returns whether mtop contains any constraints and/or vsites + * + * When we have constraints and/or vsites, it is beneficial to use + * update groups (when possible) to allow independent update of + * groups.*/ +bool systemHasConstraintsOrVsites(const gmx_mtop_t& mtop); + +/*! \libinternal + * \brief Owns the update grouping and related data */ +class UpdateGroups +{ +public: + //! Default constructor + UpdateGroups() = default; + //! Constructor when update groups are active + UpdateGroups(std::vector&& updateGroupingPerMoleculeType, real maxUpdateGroupRadius); + + bool useUpdateGroups() const { return useUpdateGroups_; } + real maxUpdateGroupRadius() const { return maxUpdateGroupRadius_; } + ArrayRef updateGroupingPerMoleculeType() const + { + return updateGroupingPerMoleculeType_; + } + +private: + //! Whether update groups are in use + const bool useUpdateGroups_ = false; + //! The update groupings within each respective molecule type, empty when not in use + const std::vector updateGroupingPerMoleculeType_ = {}; + //! The maximum radius of any update group, 0 when not in use + const real maxUpdateGroupRadius_ = 0.0_real; +}; + +/*! \brief Builder for update groups. + * + * Checks the conditions for using update groups, and logs a message + * if they cannot be used, along with the reason why not. + * + * If PP domain decomposition is not in use, there is no reason to use + * update groups. + * + * All molecule types in the system topology must be conform to the + * requirements, such that makeUpdateGroupingsPerMoleculeType() + * returns a non-empty vector. + * + * When we have constraints and/or vsites, it is beneficial to use + * update groups (when possible) to allow independent update of + * groups. But if there are no constraints or vsites, then there is no + * need to use update groups at all. * - * \param[in] mtop The system topology - * \param[in] updateGroupingsPerMoleculeType List of update group, size should match the number of moltypes in \p mtop or be 0 - * \param[in] temperature The maximum reference temperature, pass -1 when unknown or not applicable + * To use update groups, the large domain-to-domain cutoff distance + * should be compatible with the box size. */ -real computeMaxUpdateGroupRadius(const gmx_mtop_t& mtop, - gmx::ArrayRef updateGroupingsPerMoleculeType, - real temperature); +UpdateGroups makeUpdateGroups(const gmx::MDLogger& mdlog, + std::vector&& updateGroupingPerMoleculeType, + real maxUpdateGroupRadius, + bool useDomainDecomposition, + bool systemHasConstraintsOrVsites, + real cutoffMargin); } // namespace gmx diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp index 8046ef7cf3..7c4bcaf2dc 100644 --- a/src/gromacs/mdrun/runner.cpp +++ b/src/gromacs/mdrun/runner.cpp @@ -1335,6 +1335,19 @@ int Mdrunner::mdrunner() useGpuForNonbonded || (emulateGpuNonbonded == EmulateGpuNonbonded::Yes), *hwinfo_->cpuInfo); + // We need to decide on update groups early, as this affects + // inter-domain communication distances. + auto updateGroupingsPerMoleculeType = makeUpdateGroupingsPerMoleculeType(mtop); + const real maxUpdateGroupRadius = computeMaxUpdateGroupRadius( + mtop, updateGroupingsPerMoleculeType, maxReferenceTemperature(*inputrec)); + const real cutoffMargin = std::sqrt(max_cutoff2(inputrec->pbcType, box)) - inputrec->rlist; + UpdateGroups updateGroups = makeUpdateGroups(mdlog, + std::move(updateGroupingsPerMoleculeType), + maxUpdateGroupRadius, + useDomainDecomposition, + systemHasConstraintsOrVsites(mtop), + cutoffMargin); + // This builder is necessary while we have multi-part construction // of DD. Before DD is constructed, we use the existence of // the builder object to indicate that further construction of DD @@ -1350,6 +1363,9 @@ int Mdrunner::mdrunner() mtop, *inputrec, box, + updateGroups.updateGroupingPerMoleculeType(), + updateGroups.useUpdateGroups(), + updateGroups.maxUpdateGroupRadius(), positionsFromStatePointer(globalState.get())); } else @@ -1427,11 +1443,10 @@ int Mdrunner::mdrunner() // update is done so late. try { - const bool useUpdateGroups = cr->dd ? ddUsesUpdateGroups(*cr->dd) : false; const bool haveFrozenAtoms = inputrecFrozenAtoms(inputrec.get()); useGpuForUpdate = decideWhetherToUseGpuForUpdate(useDomainDecomposition, - useUpdateGroups, + updateGroups.useUpdateGroups(), pmeRunMode, domdecOptions.numPmeRanks > 0, useGpuForNonbonded, @@ -1727,12 +1742,8 @@ int Mdrunner::mdrunner() } /* Initialize the virtual site communication */ - gmx::ArrayRef updateGroupingsPerMoleculeType; - if (DOMAINDECOMP(cr)) - { - updateGroupingsPerMoleculeType = getUpdateGroupingsPerMoleculeType(*cr->dd); - } - vsite = makeVirtualSitesHandler(mtop, cr, fr->pbcType, updateGroupingsPerMoleculeType); + vsite = makeVirtualSitesHandler( + mtop, cr, fr->pbcType, updateGroups.updateGroupingPerMoleculeType()); calc_shifts(box, fr->shift_vec); @@ -1940,7 +1951,7 @@ int Mdrunner::mdrunner() doEssentialDynamics, fplog, cr, - DOMAINDECOMP(cr) && ddMayHaveSplitConstraints(*cr->dd), + updateGroups.useUpdateGroups(), ms, &nrnb, wcycle.get(), diff --git a/src/testutils/include/testutils/loggertest.h b/src/testutils/include/testutils/loggertest.h index a597b574e9..3418febd25 100644 --- a/src/testutils/include/testutils/loggertest.h +++ b/src/testutils/include/testutils/loggertest.h @@ -77,6 +77,9 @@ public: * * If not called for a log level, all entries for that level are * accepted. + * + * Note that this expectation should be set up before the logger is + * used in the test. */ void expectEntryMatchingRegex(gmx::MDLogger::LogLevel level, const char* re); @@ -85,6 +88,9 @@ public: * * If not called for a log level, all entries for that level are * accepted. + * + * Note that this expectation should be set up before the logger is + * used in the test. */ void expectNoEntries(gmx::MDLogger::LogLevel level); -- 2.22.0