Refactor gmx_reverse_top_t
authorMark Abraham <mark.j.abraham@gmail.com>
Fri, 14 May 2021 07:13:52 +0000 (09:13 +0200)
committerPaul Bauer <paul.bauer.q@gmail.com>
Mon, 17 May 2021 08:17:35 +0000 (08:17 +0000)
Made the impl properly private, but to do so had to separate
the responsibility for checking the local topology, which
should never have been here.

This introduces two new headers in preparation for splitting up
domdec_topology.cpp into several parts that can then have their own
headers. But first gmx_reverse_top_t has to be properly encapsulated.

Refs #3887

admin/lsan-suppressions.txt
src/gromacs/domdec/domdec.cpp
src/gromacs/domdec/domdec.h
src/gromacs/domdec/domdec_internal.h
src/gromacs/domdec/domdec_struct.h
src/gromacs/domdec/domdec_topology.cpp
src/gromacs/domdec/localtopologychecker.h [new file with mode: 0644]
src/gromacs/domdec/partition.cpp
src/gromacs/domdec/reversetopology.h [new file with mode: 0644]

index ce8a33ee2b886b0d6f90e511c5b9e0eb58aebbe4..12ceca18c179f61b3f7636d3f5ad75799f252d5a 100644 (file)
@@ -16,6 +16,7 @@ leak:do_inputrec
 leak:do_single_flood
 leak:get_zone_pulse_groups
 leak:gmx::DomainDecompositionBuilder::Impl::build
+leak:gmx::DomainDecompositionBuilder::Impl::Impl
 leak:gmx_check
 leak:gmx_make_edi
 leak:gmx_pme_receive_f
index 93a5f5eeb6bd123f04e89219e3a616705349c55a..c729d4b2b7c4a0246ce5e24318aeed290d04c8a6 100644 (file)
 #include "gromacs/domdec/domdec_network.h"
 #include "gromacs/domdec/ga2la.h"
 #include "gromacs/domdec/gpuhaloexchange.h"
+#include "gromacs/domdec/localtopologychecker.h"
 #include "gromacs/domdec/options.h"
 #include "gromacs/domdec/partition.h"
+#include "gromacs/domdec/reversetopology.h"
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/gmxlib/nrnb.h"
 #include "gromacs/gpu_utils/device_stream_manager.h"
@@ -2995,6 +2997,8 @@ gmx_domdec_t* DomainDecompositionBuilder::Impl::build(LocalAtomSetManager* atomS
 
     dd->atomSets = atomSets;
 
+    dd->localTopologyChecker = std::make_unique<LocalTopologyChecker>();
+
     return dd;
 }
 
index e36334fdc618aa85d6d40e5170a96204733a5990..d157513be386bc2988f0225a1bc155b109d43069 100644 (file)
@@ -258,6 +258,10 @@ 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);
@@ -287,17 +291,19 @@ void checkNumberOfBondedInteractions(const gmx::MDLogger&           mdlog,
                                      gmx::ArrayRef<const gmx::RVec> x,
                                      const matrix                   box);
 
-/*! \brief Generate the local topology and virtual site data */
-void dd_make_local_top(struct gmx_domdec_t*           dd,
-                       struct gmx_domdec_zones_t*     zones,
-                       int                            npbcdim,
-                       matrix                         box,
-                       rvec                           cellsize_min,
-                       const ivec                     npulse,
-                       t_forcerec*                    fr,
-                       gmx::ArrayRef<const gmx::RVec> coordinates,
-                       const gmx_mtop_t&              top,
-                       gmx_localtop_t*                ltop);
+/*! \brief Generate the local topology and virtual site data
+ *
+ * \returns Total count of bonded interactions in the local topology on this domain */
+int dd_make_local_top(struct gmx_domdec_t*           dd,
+                      struct gmx_domdec_zones_t*     zones,
+                      int                            npbcdim,
+                      matrix                         box,
+                      rvec                           cellsize_min,
+                      const ivec                     npulse,
+                      t_forcerec*                    fr,
+                      gmx::ArrayRef<const gmx::RVec> coordinates,
+                      const gmx_mtop_t&              top,
+                      gmx_localtop_t*                ltop);
 
 /*! \brief Sort ltop->ilist when we are doing free energy. */
 void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop);
index a112f8e2564404ab977e84efa4d30c2d2060e0d3..0c30bf2551d09e6f01b2939868839a3809b57b64 100644 (file)
@@ -799,18 +799,4 @@ static constexpr double DD_PRES_SCALE_MARGIN = 1.02;
 
 /*! \endcond */
 
-//! \internal \brief Reverse topology class
-struct gmx_reverse_top_t
-{
-    //! Constructor
-    gmx_reverse_top_t(const gmx_mtop_t& mtop, bool useFreeEnergy, const ReverseTopOptions& reverseTopOptions);
-    //! Destructor
-    ~gmx_reverse_top_t();
-
-    //! Private implementation definition
-    struct Impl;
-    //! Private implementation declaration
-    std::unique_ptr<Impl> impl_;
-};
-
 #endif
index e4e2c784cece6581e1d5fb129dd3c6e6bd1c4e98..002dd66f0530c09d2a5effc8ffba2a0baef2c146 100644 (file)
@@ -69,7 +69,7 @@ struct gmx_domdec_specat_comm_t;
 class gmx_ga2la_t;
 struct gmx_pme_comm_n_box_t;
 struct t_inputrec;
-struct gmx_reverse_top_t;
+class gmx_reverse_top_t;
 struct gmx_mtop_t;
 struct ReverseTopOptions;
 
@@ -78,6 +78,7 @@ namespace gmx
 template<typename T>
 class HashedMap;
 class LocalAtomSetManager;
+struct LocalTopologyChecker;
 class GpuHaloExchange;
 } // namespace gmx
 
@@ -233,6 +234,9 @@ struct gmx_domdec_t
     /* The managed atom sets that are updated in domain decomposition */
     gmx::LocalAtomSetManager* atomSets = nullptr;
 
+    //! The handler for checking whether the local topology is missing interactions
+    std::unique_ptr<gmx::LocalTopologyChecker> localTopologyChecker;
+
     /* gmx_pme_recv_f buffer */
     std::vector<gmx::RVec> pmeForceReceiveBuffer;
 
index 3a9f715cf40f11011599982dfeb40facdc0392e4..72a2e07ec33be09743ef24d8b9c88d9067d6bb55 100644 (file)
@@ -58,7 +58,9 @@
 #include "gromacs/domdec/domdec.h"
 #include "gromacs/domdec/domdec_network.h"
 #include "gromacs/domdec/ga2la.h"
+#include "gromacs/domdec/localtopologychecker.h"
 #include "gromacs/domdec/options.h"
+#include "gromacs/domdec/reversetopology.h"
 #include "gromacs/gmxlib/network.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/forcerec.h"
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/pbcutil/mshift.h"
 #include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/ifunc.h"
 #include "gromacs/topology/mtop_util.h"
 #include "gromacs/topology/topsort.h"
 #include "gromacs/utility/cstringutil.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
-#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/logger.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/strconvert.h"
@@ -100,57 +100,6 @@ using gmx::RVec;
 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
 #define NITEM_DD_INIT_LOCAL_STATE 5
 
-struct reverse_ilist_t
-{
-    std::vector<int> index;              /* Index for each atom into il          */
-    std::vector<int> il;                 /* ftype|type|a0|...|an|ftype|...       */
-    int              numAtomsInMolecule; /* The number of atoms in this molecule */
-};
-
-struct MolblockIndices
-{
-    int a_start;
-    int a_end;
-    int natoms_mol;
-    int type;
-};
-
-/*! \brief Struct for thread local work data for local topology generation */
-struct thread_work_t
-{
-    /*! \brief Constructor
-     *
-     * \param[in] ffparams  The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
-     */
-    thread_work_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
-
-    InteractionDefinitions         idef;               /**< Partial local topology */
-    std::unique_ptr<gmx::VsitePbc> vsitePbc = nullptr; /**< vsite PBC structure */
-    int              numBondedInteractions  = 0; /**< The number of bonded interactions observed */
-    ListOfLists<int> excl;                       /**< List of exclusions */
-};
-
-/*! \brief Options for setting up gmx_reverse_top_t */
-struct ReverseTopOptions
-{
-    //! Constructor, constraints and settles are not including with a single argument
-    ReverseTopOptions(DDBondedChecking ddBondedChecking,
-                      bool             includeConstraints = false,
-                      bool             includeSettles     = false) :
-        ddBondedChecking(ddBondedChecking),
-        includeConstraints(includeConstraints),
-        includeSettles(includeSettles)
-    {
-    }
-
-    //! \brief For which bonded interactions to check assignments
-    const DDBondedChecking ddBondedChecking;
-    //! \brief Whether constraints are stored in this reverse top
-    const bool includeConstraints;
-    //! \brief Whether settles are stored in this reverse top
-    const bool includeSettles;
-};
-
 /*! \brief Struct for the reverse topology: links bonded interactions to atoms */
 struct gmx_reverse_top_t::Impl
 {
@@ -161,7 +110,7 @@ struct gmx_reverse_top_t::Impl
     //! Options for the setup of this reverse topology
     const ReverseTopOptions options;
     //! Are there interaction of type F_POSRES and/or F_FBPOSRES
-    bool havePositionRestraints;
+    bool hasPositionRestraints;
     //! \brief Are there bondeds/exclusions between atoms?
     bool bInterAtomicInteractions = false;
     //! \brief Reverse ilist for all moltypes
@@ -329,7 +278,7 @@ static std::string printMissingInteractionsMolblock(t_commrec*               cr,
 
     for (int ftype = 0; ftype < F_NRE; ftype++)
     {
-        if (dd_check_ftype(ftype, rt.impl_->options))
+        if (dd_check_ftype(ftype, rt.options()))
         {
             flagInteractionsForType(
                     ftype, idef.il[ftype], ril, atomRange, numAtomsPerMolecule, cr->dd->globalAtomIndices, isAssigned);
@@ -409,8 +358,14 @@ static void printMissingInteractionsAtoms(const gmx::MDLogger&          mdlog,
         a_end                        = a_start + molb.nmol * moltype.atoms.nr;
         const gmx::Range<int> atomRange(a_start, a_end);
 
-        auto warning = printMissingInteractionsMolblock(
-                cr, rt, *(moltype.name), rt.impl_->ril_mt[molb.type], atomRange, moltype.atoms.nr, molb.nmol, idef);
+        auto warning = printMissingInteractionsMolblock(cr,
+                                                        rt,
+                                                        *(moltype.name),
+                                                        rt.interactionListForMoleculeType(molb.type),
+                                                        atomRange,
+                                                        moltype.atoms.nr,
+                                                        molb.nmol,
+                                                        idef);
 
         GMX_LOG(mdlog.warning).appendText(warning);
     }
@@ -433,7 +388,7 @@ void dd_print_missing_interactions(const gmx::MDLogger&  mdlog,
                     "decomposition cells");
 
     const int ndiff_tot = numBondedInteractionsOverAllDomains
-                          - dd->reverse_top->impl_->expectedNumGlobalBondedInteractions;
+                          - dd->reverse_top->expectedNumGlobalBondedInteractions();
 
     for (int ftype = 0; ftype < F_NRE; ftype++)
     {
@@ -446,7 +401,7 @@ void dd_print_missing_interactions(const gmx::MDLogger&  mdlog,
     if (DDMASTER(dd))
     {
         GMX_LOG(mdlog.warning).appendText("A list of missing interactions:");
-        int rest_global = dd->reverse_top->impl_->expectedNumGlobalBondedInteractions;
+        int rest_global = dd->reverse_top->expectedNumGlobalBondedInteractions();
         int rest        = numBondedInteractionsOverAllDomains;
         for (int ftype = 0; ftype < F_NRE; ftype++)
         {
@@ -454,7 +409,7 @@ void dd_print_missing_interactions(const gmx::MDLogger&  mdlog,
              * 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->impl_->options) && ftype != F_CONSTRNC)
+            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)
@@ -501,7 +456,7 @@ void dd_print_missing_interactions(const gmx::MDLogger&  mdlog,
                 "two-body cut-off distance (%g nm), see option -rdd, for pairs and tabulated bonds "
                 "also see option -ddcheck",
                 -ndiff_tot,
-                dd->reverse_top->impl_->expectedNumGlobalBondedInteractions,
+                dd->reverse_top->expectedNumGlobalBondedInteractions(),
                 dd_cutoff_multibody(dd),
                 dd_cutoff_twobody(dd));
     }
@@ -565,13 +520,6 @@ static int getMaxNumExclusionsPerAtom(const ListOfLists<int>& excls)
     return maxNumExcls;
 }
 
-//! Options for linking atoms in make_reverse_ilist
-enum class AtomLinkRule
-{
-    FirstAtom,        //!< Link all interactions to the first atom in the atom list
-    AllAtomsInBondeds //!< Link bonded interactions to all atoms involved, don't link vsites
-};
-
 /*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */
 static int low_make_reverse_ilist(const InteractionLists&  il_mt,
                                   const t_atom*            atom,
@@ -694,13 +642,62 @@ gmx_reverse_top_t::gmx_reverse_top_t(const gmx_mtop_t&        mtop,
 
 gmx_reverse_top_t::~gmx_reverse_top_t() {}
 
+const ReverseTopOptions& gmx_reverse_top_t::options() const
+{
+    return impl_->options;
+}
+
+const reverse_ilist_t& gmx_reverse_top_t::interactionListForMoleculeType(int moleculeType) const
+{
+    return impl_->ril_mt[moleculeType];
+}
+
+int gmx_reverse_top_t::expectedNumGlobalBondedInteractions() const
+{
+    return impl_->expectedNumGlobalBondedInteractions;
+}
+
+ArrayRef<const MolblockIndices> gmx_reverse_top_t::molblockIndices() const
+{
+    return impl_->mbi;
+}
+
+bool gmx_reverse_top_t::hasIntermolecularInteractions() const
+{
+    return impl_->bIntermolecularInteractions;
+}
+
+const reverse_ilist_t& gmx_reverse_top_t::interactionListForIntermolecularInteractions() const
+{
+    return impl_->ril_intermol;
+}
+
+bool gmx_reverse_top_t::hasInterAtomicInteractions() const
+{
+    return impl_->bInterAtomicInteractions;
+}
+
+bool gmx_reverse_top_t::hasPositionRestraints() const
+{
+    return impl_->hasPositionRestraints;
+}
+
+ArrayRef<thread_work_t> gmx_reverse_top_t::threadWorkObjects() const
+{
+    return impl_->th_work;
+}
+
+bool gmx_reverse_top_t::doSorting() const
+{
+    return impl_->ilsort != ilsortNO_FE;
+}
+
 /*! \brief Generate the reverse topology */
 gmx_reverse_top_t::Impl::Impl(const gmx_mtop_t&        mtop,
                               const bool               useFreeEnergy,
                               const ReverseTopOptions& reverseTopOptions) :
     options(reverseTopOptions),
-    havePositionRestraints(
-            gmx_mtop_ftype_count(mtop, F_POSRES) + gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0),
+    hasPositionRestraints(gmx_mtop_ftype_count(mtop, F_POSRES) + gmx_mtop_ftype_count(mtop, F_FBPOSRES) > 0),
     bInterAtomicInteractions(mtop.bIntermolecularInteractions)
 {
     bInterAtomicInteractions = mtop.bIntermolecularInteractions;
@@ -1388,7 +1385,7 @@ static int make_bondeds_zone(gmx_reverse_top_t*                 rt,
     int mol                 = 0;
     int atomIndexInMolecule = 0;
 
-    const auto ddBondedChecking = rt->impl_->options.ddBondedChecking;
+    const auto ddBondedChecking = rt->options().ddBondedChecking;
 
     int numBondedInteractions = 0;
 
@@ -1396,10 +1393,11 @@ static int make_bondeds_zone(gmx_reverse_top_t*                 rt,
     {
         /* Get the global atom number */
         const int globalAtomIndex = globalAtomIndices[i];
-        global_atomnr_to_moltype_ind(rt->impl_->mbi, globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule);
+        global_atomnr_to_moltype_ind(
+                rt->molblockIndices(), globalAtomIndex, &mb, &mt, &mol, &atomIndexInMolecule);
         /* Check all intramolecular interactions assigned to this atom */
-        gmx::ArrayRef<const int>     index = rt->impl_->ril_mt[mt].index;
-        gmx::ArrayRef<const t_iatom> rtil  = rt->impl_->ril_mt[mt].il;
+        gmx::ArrayRef<const int>     index = rt->interactionListForMoleculeType(mt).index;
+        gmx::ArrayRef<const t_iatom> rtil  = rt->interactionListForMoleculeType(mt).il;
 
         numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intramolecular>(
                 i,
@@ -1422,26 +1420,26 @@ static int make_bondeds_zone(gmx_reverse_top_t*                 rt,
                 ddBondedChecking);
 
         // Assign position restraints, when present, for the home zone
-        if (izone == 0 && rt->impl_->havePositionRestraints)
+        if (izone == 0 && rt->hasPositionRestraints())
         {
-            numBondedInteractions +=
-                    assignPositionRestraintsForAtom(i,
-                                                    mol,
-                                                    atomIndexInMolecule,
-                                                    rt->impl_->ril_mt[mt].numAtomsInMolecule,
-                                                    rtil,
-                                                    index[atomIndexInMolecule],
-                                                    index[atomIndexInMolecule + 1],
-                                                    molb[mb],
-                                                    ip_in,
-                                                    idef);
+            numBondedInteractions += assignPositionRestraintsForAtom(
+                    i,
+                    mol,
+                    atomIndexInMolecule,
+                    rt->interactionListForMoleculeType(mt).numAtomsInMolecule,
+                    rtil,
+                    index[atomIndexInMolecule],
+                    index[atomIndexInMolecule + 1],
+                    molb[mb],
+                    ip_in,
+                    idef);
         }
 
-        if (rt->impl_->bIntermolecularInteractions)
+        if (rt->hasIntermolecularInteractions())
         {
             /* Check all intermolecular interactions assigned to this atom */
-            index = rt->impl_->ril_intermol.index;
-            rtil  = rt->impl_->ril_intermol.il;
+            index = rt->interactionListForIntermolecularInteractions().index;
+            rtil  = rt->interactionListForIntermolecularInteractions().il;
 
             numBondedInteractions += assignInteractionsForAtom<InteractionConnectivity::Intermolecular>(
                     i,
@@ -1543,23 +1541,25 @@ static void make_exclusions_zone(ArrayRef<const int>               globalAtomInd
             "The number of exclusion list should match the number of atoms in the range");
 }
 
-/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */
-static void make_local_bondeds_excls(gmx_domdec_t*           dd,
-                                     gmx_domdec_zones_t*     zones,
-                                     const gmx_mtop_t&       mtop,
-                                     const int*              cginfo,
-                                     const bool              checkDistanceMultiBody,
-                                     const ivec              rcheck,
-                                     const gmx_bool          checkDistanceTwoBody,
-                                     const real              cutoff,
-                                     const t_pbc*            pbc_null,
-                                     ArrayRef<const RVec>    coordinates,
-                                     InteractionDefinitions* idef,
-                                     ListOfLists<int>*       lexcls)
+/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls
+ *
+ * \returns Total count of bonded interactions in the local topology on this domain */
+static int make_local_bondeds_excls(gmx_domdec_t*           dd,
+                                    gmx_domdec_zones_t*     zones,
+                                    const gmx_mtop_t&       mtop,
+                                    const int*              cginfo,
+                                    const bool              checkDistanceMultiBody,
+                                    const ivec              rcheck,
+                                    const gmx_bool          checkDistanceTwoBody,
+                                    const real              cutoff,
+                                    const t_pbc*            pbc_null,
+                                    ArrayRef<const RVec>    coordinates,
+                                    InteractionDefinitions* idef,
+                                    ListOfLists<int>*       lexcls)
 {
     int nzone_bondeds = 0;
 
-    if (dd->reverse_top->impl_->bInterAtomicInteractions)
+    if (dd->reverse_top->hasInterAtomicInteractions())
     {
         nzone_bondeds = zones->n;
     }
@@ -1580,7 +1580,7 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
 
     /* Clear the counts */
     idef->clear();
-    dd->reverse_top->impl_->numBondedInteractions = 0;
+    int numBondedInteractions = 0;
 
     lexcls->clear();
 
@@ -1589,7 +1589,8 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
         const int cg0 = zones->cg_range[izone];
         const int cg1 = zones->cg_range[izone + 1];
 
-        const int numThreads = rt->impl_->th_work.size();
+        gmx::ArrayRef<thread_work_t> threadWorkObjects = rt->threadWorkObjects();
+        const int                    numThreads        = threadWorkObjects.size();
 #pragma omp parallel for num_threads(numThreads) schedule(static)
         for (int thread = 0; thread < numThreads; thread++)
         {
@@ -1606,11 +1607,11 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
                 }
                 else
                 {
-                    idef_t = &rt->impl_->th_work[thread].idef;
+                    idef_t = &threadWorkObjects[thread].idef;
                     idef_t->clear();
                 }
 
-                rt->impl_->th_work[thread].numBondedInteractions =
+                threadWorkObjects[thread].numBondedInteractions =
                         make_bondeds_zone(rt,
                                           dd->globalAtomIndices,
                                           *dd->ga2la,
@@ -1638,7 +1639,7 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
                     else
                     {
                         // Threads > 0 store in temporary storage, starting at list index 0
-                        excl_t = &rt->impl_->th_work[thread].excl;
+                        excl_t = &threadWorkObjects[thread].excl;
                         excl_t->clear();
                     }
 
@@ -1646,7 +1647,7 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
                     make_exclusions_zone(dd->globalAtomIndices,
                                          *dd->ga2la,
                                          zones,
-                                         rt->impl_->mbi,
+                                         rt->molblockIndices(),
                                          mtop.moltype,
                                          cginfo,
                                          excl_t,
@@ -1659,55 +1660,61 @@ static void make_local_bondeds_excls(gmx_domdec_t*           dd,
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
         }
 
-        if (rt->impl_->th_work.size() > 1)
+        if (threadWorkObjects.size() > 1)
         {
-            combine_idef(idef, rt->impl_->th_work);
+            combine_idef(idef, threadWorkObjects);
         }
 
-        for (const thread_work_t& th_work : rt->impl_->th_work)
+        for (const thread_work_t& th_work : threadWorkObjects)
         {
-            dd->reverse_top->impl_->numBondedInteractions += th_work.numBondedInteractions;
+            numBondedInteractions += th_work.numBondedInteractions;
         }
 
         if (izone < numIZonesForExclusions)
         {
-            for (std::size_t th = 1; th < rt->impl_->th_work.size(); th++)
+            for (std::size_t th = 1; th < threadWorkObjects.size(); th++)
             {
-                lexcls->appendListOfLists(rt->impl_->th_work[th].excl);
+                lexcls->appendListOfLists(threadWorkObjects[th].excl);
             }
         }
     }
 
+    if (debug)
+    {
+        fprintf(debug, "We have %d exclusions\n", lexcls->numElements());
+    }
+
+    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->reverse_top->impl_->shouldCheckNumberOfBondedInteractions = true;
+    dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = true;
     // Clear the old global value, which is now invalid
-    dd->reverse_top->impl_->numBondedInteractionsOverAllDomains.reset();
-
-    if (debug)
-    {
-        fprintf(debug, "We have %d exclusions\n", lexcls->numElements());
-    }
+    dd->localTopologyChecker->numBondedInteractionsOverAllDomains.reset();
 }
 
 bool shouldCheckNumberOfBondedInteractions(const gmx_domdec_t& dd)
 {
-    return dd.reverse_top->impl_->shouldCheckNumberOfBondedInteractions;
+    return dd.localTopologyChecker->shouldCheckNumberOfBondedInteractions;
 }
 
 int numBondedInteractions(const gmx_domdec_t& dd)
 {
-    return dd.reverse_top->impl_->numBondedInteractions;
+    return dd.localTopologyChecker->numBondedInteractionsToReduce;
 }
 
 void setNumberOfBondedInteractionsOverAllDomains(const gmx_domdec_t& dd, int newValue)
 {
-    GMX_RELEASE_ASSERT(!dd.reverse_top->impl_->numBondedInteractionsOverAllDomains.has_value(),
+    GMX_RELEASE_ASSERT(!dd.localTopologyChecker->numBondedInteractionsOverAllDomains.has_value(),
                        "Cannot set number of bonded interactions because it is already set");
-    dd.reverse_top->impl_->numBondedInteractionsOverAllDomains.emplace(newValue);
+    dd.localTopologyChecker->numBondedInteractionsOverAllDomains.emplace(newValue);
 }
 
 void checkNumberOfBondedInteractions(const gmx::MDLogger&           mdlog,
@@ -1720,18 +1727,18 @@ void checkNumberOfBondedInteractions(const gmx::MDLogger&           mdlog,
     GMX_RELEASE_ASSERT(
             DOMAINDECOMP(cr),
             "No need to check number of bonded interactions when not using domain decomposition");
-    if (cr->dd->reverse_top->impl_->shouldCheckNumberOfBondedInteractions)
+    if (cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions)
     {
-        GMX_RELEASE_ASSERT(cr->dd->reverse_top->impl_->numBondedInteractionsOverAllDomains.has_value(),
+        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->reverse_top->impl_->numBondedInteractionsOverAllDomains.value()
-            != cr->dd->reverse_top->impl_->expectedNumGlobalBondedInteractions)
+        if (cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value()
+            != cr->dd->reverse_top->expectedNumGlobalBondedInteractions())
         {
             dd_print_missing_interactions(
                     mdlog,
                     cr,
-                    cr->dd->reverse_top->impl_->numBondedInteractionsOverAllDomains.value(),
+                    cr->dd->localTopologyChecker->numBondedInteractionsOverAllDomains.value(),
                     top_global,
                     top_local,
                     x,
@@ -1740,20 +1747,20 @@ void checkNumberOfBondedInteractions(const gmx::MDLogger&           mdlog,
         // 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->reverse_top->impl_->shouldCheckNumberOfBondedInteractions = false;
+        cr->dd->localTopologyChecker->shouldCheckNumberOfBondedInteractions = false;
     }
 }
 
-void dd_make_local_top(gmx_domdec_t*        dd,
-                       gmx_domdec_zones_t*  zones,
-                       int                  npbcdim,
-                       matrix               box,
-                       rvec                 cellsize_min,
-                       const ivec           npulse,
-                       t_forcerec*          fr,
-                       ArrayRef<const RVec> coordinates,
-                       const gmx_mtop_t&    mtop,
-                       gmx_localtop_t*      ltop)
+int dd_make_local_top(gmx_domdec_t*        dd,
+                      gmx_domdec_zones_t*  zones,
+                      int                  npbcdim,
+                      matrix               box,
+                      rvec                 cellsize_min,
+                      const ivec           npulse,
+                      t_forcerec*          fr,
+                      ArrayRef<const RVec> coordinates,
+                      const gmx_mtop_t&    mtop,
+                      gmx_localtop_t*      ltop)
 {
     real  rc = -1;
     ivec  rcheck;
@@ -1767,7 +1774,7 @@ void dd_make_local_top(gmx_domdec_t*        dd,
     bool checkDistanceMultiBody = false;
     bool checkDistanceTwoBody   = false;
 
-    if (dd->reverse_top->impl_->bInterAtomicInteractions)
+    if (dd->reverse_top->hasInterAtomicInteractions())
     {
         /* We need to check to which cell bondeds should be assigned */
         rc = dd_cutoff_twobody(dd);
@@ -1821,34 +1828,36 @@ void dd_make_local_top(gmx_domdec_t*        dd,
         }
     }
 
-    make_local_bondeds_excls(dd,
-                             zones,
-                             mtop,
-                             fr->cginfo.data(),
-                             checkDistanceMultiBody,
-                             rcheck,
-                             checkDistanceTwoBody,
-                             rc,
-                             pbc_null,
-                             coordinates,
-                             &ltop->idef,
-                             &ltop->excls);
+    int numBondedInteractionsToReduce = make_local_bondeds_excls(dd,
+                                                                 zones,
+                                                                 mtop,
+                                                                 fr->cginfo.data(),
+                                                                 checkDistanceMultiBody,
+                                                                 rcheck,
+                                                                 checkDistanceTwoBody,
+                                                                 rc,
+                                                                 pbc_null,
+                                                                 coordinates,
+                                                                 &ltop->idef,
+                                                                 &ltop->excls);
 
     /* The ilist is not sorted yet,
      * we can only do this when we have the charge arrays.
      */
     ltop->idef.ilsort = ilsortUNKNOWN;
+
+    return numBondedInteractionsToReduce;
 }
 
 void dd_sort_local_top(const gmx_domdec_t& dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop)
 {
-    if (dd.reverse_top->impl_->ilsort == ilsortNO_FE)
+    if (dd.reverse_top->doSorting())
     {
-        ltop->idef.ilsort = ilsortNO_FE;
+        gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
     }
     else
     {
-        gmx_sort_ilist_fe(&ltop->idef, mdatoms->chargeA, mdatoms->chargeB);
+        ltop->idef.ilsort = ilsortNO_FE;
     }
 }
 
diff --git a/src/gromacs/domdec/localtopologychecker.h b/src/gromacs/domdec/localtopologychecker.h
new file mode 100644 (file)
index 0000000..54517dc
--- /dev/null
@@ -0,0 +1,84 @@
+/*
+ * 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.
+ */
+/*! \libinternal \file
+ *
+ * \brief This file declares functionality for checking whether
+ * local topologies describe all bonded interactions.
+ *
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_LOCALTOPOLOGYCHECKER_H
+#define GMX_DOMDEC_LOCALTOPOLOGYCHECKER_H
+
+#include <optional>
+
+namespace gmx
+{
+
+struct LocalTopologyChecker
+{
+    /*! \brief Data to help check local topology construction
+     *
+     * Partitioning could incorrectly miss a bonded interaction.
+     * However, checking for that requires a global communication
+     * stage, which does not otherwise happen during partitioning. So,
+     * for performance, we do that alongside the first global energy
+     * reduction after a new DD is made. These variables handle
+     * whether the check happens, its input for this domain, output
+     * across all domains, and the expected value it should match. */
+    /*! \{ */
+    /*! \brief Number of bonded interactions found in the local
+     * topology for this domain. */
+    int numBondedInteractionsToReduce = 0;
+    /*! \brief Whether to check at the next global communication
+     * stage the total number of bonded interactions found.
+     *
+     * Cleared after that number is found. */
+    bool shouldCheckNumberOfBondedInteractions = false;
+    /*! \brief The total number of bonded interactions found in
+     * the local topology across all domains.
+     *
+     * Only has a value after reduction across all ranks, which is
+     * removed when it is again time to check after a new
+     * partition. */
+    std::optional<int> numBondedInteractionsOverAllDomains;
+    /*! \} */
+};
+
+} // namespace gmx
+
+#endif
index 9a9409422c21accab930841751f0f0a07e4c4291..84713d0aa64986b4f1b41ba480eafeb8da3351b6 100644 (file)
@@ -3191,16 +3191,17 @@ void dd_partition_system(FILE*                     fplog,
     {
         numPulses[dd->dim[i]] = comm->cd[i].numPulses();
     }
-    dd_make_local_top(dd,
-                      &comm->zones,
-                      dd->unitCellInfo.npbcdim,
-                      state_local->box,
-                      comm->cellsize_min,
-                      numPulses,
-                      fr,
-                      state_local->x,
-                      top_global,
-                      top_local);
+    int numBondedInteractionsToReduce = dd_make_local_top(dd,
+                                                          &comm->zones,
+                                                          dd->unitCellInfo.npbcdim,
+                                                          state_local->box,
+                                                          comm->cellsize_min,
+                                                          numPulses,
+                                                          fr,
+                                                          state_local->x,
+                                                          top_global,
+                                                          top_local);
+    scheduleCheckOfLocalTopology(dd, numBondedInteractionsToReduce);
 
     wallcycle_sub_stop(wcycle, WallCycleSubCounter::DDMakeTop);
 
diff --git a/src/gromacs/domdec/reversetopology.h b/src/gromacs/domdec/reversetopology.h
new file mode 100644 (file)
index 0000000..62568fe
--- /dev/null
@@ -0,0 +1,167 @@
+/*
+ * 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.
+ */
+/*! \libinternal \file
+ *
+ * \brief This file makes declarations used for building
+ * the reverse topology
+ *
+ * \inlibraryapi
+ * \ingroup module_domdec
+ */
+
+#ifndef GMX_DOMDEC_REVERSETOPOLOGY_H
+#define GMX_DOMDEC_REVERSETOPOLOGY_H
+
+#include <cstdio>
+
+#include <memory>
+#include <vector>
+
+#include "gromacs/mdlib/vsite.h"
+#include "gromacs/topology/idef.h"
+#include "gromacs/utility/listoflists.h"
+
+struct gmx_domdec_t;
+struct gmx_ffparams_t;
+struct gmx_mtop_t;
+struct t_atoms;
+struct t_inputrec;
+struct ReverseTopOptions;
+
+namespace gmx
+{
+class VirtualSitesHandler;
+enum class DDBondedChecking : bool;
+} // namespace gmx
+
+//! Options for linking atoms in make_reverse_ilist
+enum class AtomLinkRule
+{
+    FirstAtom,        //!< Link all interactions to the first atom in the atom list
+    AllAtomsInBondeds //!< Link bonded interactions to all atoms involved, don't link vsites
+};
+
+struct MolblockIndices
+{
+    int a_start;
+    int a_end;
+    int natoms_mol;
+    int type;
+};
+
+struct reverse_ilist_t
+{
+    std::vector<int> index;              /* Index for each atom into il          */
+    std::vector<int> il;                 /* ftype|type|a0|...|an|ftype|...       */
+    int              numAtomsInMolecule; /* The number of atoms in this molecule */
+};
+
+/*! \brief Struct for thread local work data for local topology generation */
+struct thread_work_t
+{
+    /*! \brief Constructor
+     *
+     * \param[in] ffparams  The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
+     */
+    thread_work_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
+
+    InteractionDefinitions         idef;               /**< Partial local topology */
+    std::unique_ptr<gmx::VsitePbc> vsitePbc = nullptr; /**< vsite PBC structure */
+    int numBondedInteractions               = 0; /**< The number of bonded interactions observed */
+    gmx::ListOfLists<int> excl;                  /**< List of exclusions */
+    int                   excl_count = 0;        /**< The total exclusion count for \p excl */
+};
+
+/*! \brief Options for setting up gmx_reverse_top_t */
+struct ReverseTopOptions
+{
+    //! Constructor, constraints and settles are not including with a single argument
+    ReverseTopOptions(gmx::DDBondedChecking ddBondedChecking,
+                      bool                  includeConstraints = false,
+                      bool                  includeSettles     = false) :
+        ddBondedChecking(ddBondedChecking),
+        includeConstraints(includeConstraints),
+        includeSettles(includeSettles)
+    {
+    }
+
+    //! \brief For which bonded interactions to check assignments
+    const gmx::DDBondedChecking ddBondedChecking;
+    //! \brief Whether constraints are stored in this reverse top
+    const bool includeConstraints;
+    //! \brief Whether settles are stored in this reverse top
+    const bool includeSettles;
+};
+
+//! \internal \brief Reverse topology class
+class gmx_reverse_top_t
+{
+public:
+    //! Constructor
+    gmx_reverse_top_t(const gmx_mtop_t& mtop, bool useFreeEnergy, const ReverseTopOptions& reverseTopOptions);
+    //! Destructor
+    ~gmx_reverse_top_t();
+
+    //! Gets the options that configured the construction
+    const ReverseTopOptions& options() const;
+
+    //! Gets the interaction list for the given molecule type
+    const reverse_ilist_t& interactionListForMoleculeType(int moleculeType) const;
+
+    //! Returns the total count of bonded interactions, used for checking partitioning
+    int expectedNumGlobalBondedInteractions() const;
+
+    //! Returns the molecule block indices
+    gmx::ArrayRef<const MolblockIndices> molblockIndices() const;
+    //! Returns whether the reverse topology describes intermolecular interactions
+    bool hasIntermolecularInteractions() const;
+    //! Gets the interaction list for any intermolecular interactions
+    const reverse_ilist_t& interactionListForIntermolecularInteractions() const;
+    //! Returns whether the reverse topology describes interatomic interactions
+    bool hasInterAtomicInteractions() const;
+    //! Returns whether there are interactions of type F_POSRES and/or F_FBPOSRES
+    bool hasPositionRestraints() const;
+    //! Returns the per-thread working structures for making the local topology
+    gmx::ArrayRef<thread_work_t> threadWorkObjects() const;
+    //! Returns whether the local topology interactions should be sorted
+    bool doSorting() const;
+
+    //! Private implementation definition
+    struct Impl;
+    //! Private implementation declaration
+    std::unique_ptr<Impl> impl_;
+};
+
+#endif