Decouple update group handling from domain decomposition module
authorMark Abraham <mark.j.abraham@gmail.com>
Thu, 13 May 2021 07:36:56 +0000 (07:36 +0000)
committerArtem Zhmurov <zhmurov@gmail.com>
Thu, 13 May 2021 07:36:56 +0000 (07:36 +0000)
src/gromacs/domdec/builder.h
src/gromacs/domdec/domdec.cpp
src/gromacs/domdec/domdec.h
src/gromacs/domdec/domdec_internal.h
src/gromacs/mdlib/constr.cpp
src/gromacs/mdlib/constr.h
src/gromacs/mdlib/tests/updategroups.cpp
src/gromacs/mdlib/updategroups.cpp
src/gromacs/mdlib/updategroups.h
src/gromacs/mdrun/runner.cpp
src/testutils/include/testutils/loggertest.h

index 0c6e764664a2e6ac11be4a297412cfeac8699fb5..1ca7408dd3c5b9560a0355d7e453d05ce1cd5401 100644 (file)
@@ -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<const RVec> 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<const RangePartitioning> updateGroupingPerMoleculeType,
+                               bool                              useUpdateGroups,
+                               real                              maxUpdateGroupRadius,
+                               ArrayRef<const RVec>              xGlobal);
     //! Destructor
     ~DomainDecompositionBuilder();
     //! Build the resulting DD manager
index a9304bbab60126f27fe00f449b4c2b320a73c21a..93a5f5eeb6bd123f04e89219e3a616705349c55a 100644 (file)
 #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<const gmx::RangePartitioning> 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<const RangePartitioning> 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<double>(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<double>(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<const gmx::RVec> 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<const RangePartitioning> updateGroupingPerMoleculeType,
+                                  const bool                        useUpdateGroups,
+                                  const real                        maxUpdateGroupRadius,
+                                  gmx::ArrayRef<const gmx::RVec>    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<const RVec> 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<const RangePartitioning> updateGroupingPerMoleculeType,
+         bool                              useUpdateGroups,
+         real                              maxUpdateGroupRadius,
+         ArrayRef<const RVec>              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<const RVec> 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<const RangePartitioning> updateGroupingPerMoleculeType,
+                                       const bool                        useUpdateGroups,
+                                       const real                        maxUpdateGroupRadius,
+                                       ArrayRef<const RVec>              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<const RangePartitioning> updateGroupingPerMoleculeType,
+                                                       const bool           useUpdateGroups,
+                                                       const real           maxUpdateGroupRadius,
                                                        ArrayRef<const RVec> 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))
 {
 }
 
index f60c429b6ad1a09507663246f81fd3e5bc2d08ee..e36334fdc618aa85d6d40e5170a96204733a5990 100644 (file)
@@ -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<const gmx::RangePartitioning> 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);
 
index ddfa6b847ffa2b027bcf03c0207639774599efb1..a112f8e2564404ab977e84efa4d30c2d2060e0d3 100644 (file)
@@ -423,7 +423,7 @@ struct DDSystemInfo
     //! True when update groups are used
     bool useUpdateGroups = false;
     //! Update atom grouping for each molecule type
-    std::vector<gmx::RangePartitioning> updateGroupingsPerMoleculeType;
+    gmx::ArrayRef<const gmx::RangePartitioning> updateGroupingsPerMoleculeType;
     //! The maximum radius over all update groups
     real maxUpdateGroupRadius;
 
index 0b8d94234c739dd2a2e79a5bbab263b2a21b22d0..ceee50f658bccc616013b2ef56316417dcb12639 100644 (file)
@@ -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(
index f5bc84ea2ce427aa9f4f51ef734cf0e69d77ac3d..14cd95193a8b9763aa027df328286460264b3509 100644 (file)
@@ -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,
index 4b75b16a99971991a73a105afdf7c6c4e10de37e..8934d0b74bcc6ff352639348ddf74aedf3195603 100644 (file)
@@ -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
index fe040b6a98f4639f420393fe462ceb0f23d653e4..64a0a724dd5c9cc80a2f163fc0a738dd585220ae 100644 (file)
 #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<RangePartitioning>&& 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<RangePartitioning>&& 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
index 6daa2e22a805fa960eb6fd0d583fa8bb323984fa..150d2e527b21595d6835baa7f33f94d9acb30e61 100644 (file)
 
 #include <vector>
 
+#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<RangePartitioning> 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<const RangePartitioning> 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<RangePartitioning>&& updateGroupingPerMoleculeType, real maxUpdateGroupRadius);
+
+    bool                              useUpdateGroups() const { return useUpdateGroups_; }
+    real                              maxUpdateGroupRadius() const { return maxUpdateGroupRadius_; }
+    ArrayRef<const RangePartitioning> 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<RangePartitioning> 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<const RangePartitioning> updateGroupingsPerMoleculeType,
-                                 real                                   temperature);
+UpdateGroups makeUpdateGroups(const gmx::MDLogger&             mdlog,
+                              std::vector<RangePartitioning>&& updateGroupingPerMoleculeType,
+                              real                             maxUpdateGroupRadius,
+                              bool                             useDomainDecomposition,
+                              bool                             systemHasConstraintsOrVsites,
+                              real                             cutoffMargin);
 
 } // namespace gmx
 
index 8046ef7cf35329b95e4915209dd98746a3d4687e..7c4bcaf2dc43691100c749894a69e5b368d17973 100644 (file)
@@ -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<const gmx::RangePartitioning> 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(),
index a597b574e9d6d9fd91c3437479a6cd6af00b8efd..3418febd253dcfc3879c6e4995f358c6f9c73baf 100644 (file)
@@ -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);