Merge release-2021 into master
[alexxy/gromacs.git] / src / gromacs / mdlib / vsite.cpp
index 9be178c243ddc82bbde88639c1c54b1df3f651ad..c3f9edb3b27d65b522f598d0cfdac1f8fa4de1f7 100644 (file)
@@ -60,7 +60,6 @@
 #include "gromacs/math/vec.h"
 #include "gromacs/mdlib/gmx_omp_nthreads.h"
 #include "gromacs/mdtypes/commrec.h"
-#include "gromacs/mdtypes/mdatom.h"
 #include "gromacs/pbcutil/ishift.h"
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/timing/wallcycle.h"
@@ -214,7 +213,9 @@ public:
     //! Set VSites and distribute VSite work over threads, should be called after DD partitioning
     void setVirtualSites(ArrayRef<const InteractionList> ilist,
                          ArrayRef<const t_iparams>       iparams,
-                         const t_mdatoms&                mdatoms,
+                         int                             numAtoms,
+                         int                             homenr,
+                         ArrayRef<const ParticleType>    ptype,
                          bool                            useDomdec);
 
 private:
@@ -232,13 +233,19 @@ class VirtualSitesHandler::Impl
 {
 public:
     //! Constructor, domdec should be nullptr without DD
-    Impl(const gmx_mtop_t& mtop, gmx_domdec_t* domdec, PbcType pbcType);
+    Impl(const gmx_mtop_t&                 mtop,
+         gmx_domdec_t*                     domdec,
+         PbcType                           pbcType,
+         ArrayRef<const RangePartitioning> updateGroupingPerMoleculeType);
 
     //! Returns the number of virtual sites acting over multiple update groups
     int numInterUpdategroupVirtualSites() const { return numInterUpdategroupVirtualSites_; }
 
     //! Set VSites and distribute VSite work over threads, should be called after DD partitioning
-    void setVirtualSites(ArrayRef<const InteractionList> ilist, const t_mdatoms& mdatoms);
+    void setVirtualSites(ArrayRef<const InteractionList> ilist,
+                         int                             numAtoms,
+                         int                             homenr,
+                         ArrayRef<const ParticleType>    ptype);
 
     /*! \brief Create positions of vsite atoms based for the local system
      *
@@ -2530,7 +2537,7 @@ void VirtualSitesHandler::spreadForces(ArrayRef<const RVec> x,
 }
 
 int countInterUpdategroupVsites(const gmx_mtop_t&                           mtop,
-                                gmx::ArrayRef<const gmx::RangePartitioning> updateGroupingPerMoleculetype)
+                                gmx::ArrayRef<const gmx::RangePartitioning> updateGroupingsPerMoleculeType)
 {
     int n_intercg_vsite = 0;
     for (const gmx_molblock_t& molb : mtop.molblock)
@@ -2538,9 +2545,9 @@ int countInterUpdategroupVsites(const gmx_mtop_t&                           mtop
         const gmx_moltype_t& molt = mtop.moltype[molb.type];
 
         std::vector<int> atomToGroup;
-        if (!updateGroupingPerMoleculetype.empty())
+        if (!updateGroupingsPerMoleculeType.empty())
         {
-            atomToGroup = makeAtomToGroupMapping(updateGroupingPerMoleculetype[molb.type]);
+            atomToGroup = makeAtomToGroupMapping(updateGroupingsPerMoleculeType[molb.type]);
         }
         for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
         {
@@ -2574,7 +2581,8 @@ int countInterUpdategroupVsites(const gmx_mtop_t&                           mtop
 
 std::unique_ptr<VirtualSitesHandler> makeVirtualSitesHandler(const gmx_mtop_t& mtop,
                                                              const t_commrec*  cr,
-                                                             PbcType           pbcType)
+                                                             PbcType           pbcType,
+                                                             ArrayRef<const RangePartitioning> updateGroupingPerMoleculeType)
 {
     GMX_RELEASE_ASSERT(cr != nullptr, "We need a valid commrec");
 
@@ -2603,10 +2611,10 @@ std::unique_ptr<VirtualSitesHandler> makeVirtualSitesHandler(const gmx_mtop_t& m
         return vsite;
     }
 
-    return std::make_unique<VirtualSitesHandler>(mtop, cr->dd, pbcType);
+    return std::make_unique<VirtualSitesHandler>(mtop, cr->dd, pbcType, updateGroupingPerMoleculeType);
 }
 
-ThreadingInfo::ThreadingInfo() : numThreads_(gmx_omp_nthreads_get(emntVSITE))
+ThreadingInfo::ThreadingInfo() : numThreads_(gmx_omp_nthreads_get(ModuleMultiThread::VirtualSite))
 {
     if (numThreads_ > 1)
     {
@@ -2632,27 +2640,21 @@ ThreadingInfo::ThreadingInfo() : numThreads_(gmx_omp_nthreads_get(emntVSITE))
     }
 }
 
-//! Returns the number of inter update-group vsites
-static int getNumInterUpdategroupVsites(const gmx_mtop_t& mtop, const gmx_domdec_t* domdec)
-{
-    gmx::ArrayRef<const gmx::RangePartitioning> updateGroupingPerMoleculetype;
-    if (domdec)
-    {
-        updateGroupingPerMoleculetype = getUpdateGroupingPerMoleculetype(*domdec);
-    }
-
-    return countInterUpdategroupVsites(mtop, updateGroupingPerMoleculetype);
-}
-
-VirtualSitesHandler::Impl::Impl(const gmx_mtop_t& mtop, gmx_domdec_t* domdec, const PbcType pbcType) :
-    numInterUpdategroupVirtualSites_(getNumInterUpdategroupVsites(mtop, domdec)),
+VirtualSitesHandler::Impl::Impl(const gmx_mtop_t&                       mtop,
+                                gmx_domdec_t*                           domdec,
+                                const PbcType                           pbcType,
+                                const ArrayRef<const RangePartitioning> updateGroupingPerMoleculeType) :
+    numInterUpdategroupVirtualSites_(countInterUpdategroupVsites(mtop, updateGroupingPerMoleculeType)),
     domainInfo_({ pbcType, pbcType != PbcType::No && numInterUpdategroupVirtualSites_ > 0, domdec }),
     iparams_(mtop.ffparams.iparams)
 {
 }
 
-VirtualSitesHandler::VirtualSitesHandler(const gmx_mtop_t& mtop, gmx_domdec_t* domdec, const PbcType pbcType) :
-    impl_(new Impl(mtop, domdec, pbcType))
+VirtualSitesHandler::VirtualSitesHandler(const gmx_mtop_t&                       mtop,
+                                         gmx_domdec_t*                           domdec,
+                                         const PbcType                           pbcType,
+                                         const ArrayRef<const RangePartitioning> updateGroupingPerMoleculeType) :
+    impl_(new Impl(mtop, domdec, pbcType, updateGroupingPerMoleculeType))
 {
 }
 
@@ -2688,7 +2690,7 @@ static void assignVsitesToThread(VsiteThread*                    tData,
                                  gmx::ArrayRef<int>              taskIndex,
                                  ArrayRef<const InteractionList> ilist,
                                  ArrayRef<const t_iparams>       ip,
-                                 const ParticleType*             ptype)
+                                 ArrayRef<const ParticleType>    ptype)
 {
     for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
     {
@@ -2758,7 +2760,16 @@ static void assignVsitesToThread(VsiteThread*                    tData,
                                    "a constructing atom that does not belong to our task, such "
                                    "vsites should be assigned to the single 'master' task");
 
-                        task = nthread + thread;
+                        if (tData->useInterdependentTask)
+                        {
+                            // Assign to the interdependent task
+                            task = nthread + thread;
+                        }
+                        else
+                        {
+                            // Assign to the separate, non-parallel task
+                            task = 2 * nthread;
+                        }
                     }
                 }
             }
@@ -2847,7 +2858,9 @@ static void assignVsitesToSingleTask(VsiteThread*                    tData,
 
 void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
                                     ArrayRef<const t_iparams>       iparams,
-                                    const t_mdatoms&                mdatoms,
+                                    const int                       numAtoms,
+                                    const int                       homenr,
+                                    ArrayRef<const ParticleType>    ptype,
                                     const bool                      useDomdec)
 {
     if (numThreads_ <= 1)
@@ -2920,15 +2933,15 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
          * When assigning vsites to threads, we should take care that the last
          * threads also covers the non-local range.
          */
-        vsite_atom_range = mdatoms.nr;
-        natperthread     = (mdatoms.homenr + numThreads_ - 1) / numThreads_;
+        vsite_atom_range = numAtoms;
+        natperthread     = (homenr + numThreads_ - 1) / numThreads_;
     }
 
     if (debug)
     {
         fprintf(debug,
                 "virtual site thread dist: natoms %d, range %d, natperthread %d\n",
-                mdatoms.nr,
+                numAtoms,
                 vsite_atom_range,
                 natperthread);
     }
@@ -2936,7 +2949,7 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
     /* To simplify the vsite assignment, we make an index which tells us
      * to which task particles, both non-vsites and vsites, are assigned.
      */
-    taskIndex_.resize(mdatoms.nr);
+    taskIndex_.resize(numAtoms);
 
     /* Initialize the task index array. Here we assign the non-vsite
      * particles to task=thread, so we easily figure out if vsites
@@ -2944,9 +2957,9 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
      */
     {
         int thread = 0;
-        for (int i = 0; i < mdatoms.nr; i++)
+        for (int i = 0; i < numAtoms; i++)
         {
-            if (mdatoms.ptype[i] == ParticleType::VSite)
+            if (ptype[i] == ParticleType::VSite)
             {
                 /* vsites are not assigned to a task yet */
                 taskIndex_[i] = -1;
@@ -2996,12 +3009,13 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
             }
 
             /* To avoid large f_buf allocations of #threads*vsite_atom_range
-             * we don't use task2 with more than 200000 atoms. This doesn't
-             * affect performance, since with such a large range relatively few
-             * vsites will end up in the separate task.
-             * Note that useTask2 should be the same for all threads.
+             * we don't use the interdependent tasks with more than 200000 atoms.
+             * This doesn't affect performance, since with such a large range
+             * relatively few vsites will end up in the separate task.
+             * Note that useInterdependentTask should be the same for all threads.
              */
-            tData.useInterdependentTask = (vsite_atom_range <= 200000);
+            const int c_maxNumLocalAtomsForInterdependentTask = 200000;
+            tData.useInterdependentTask = (vsite_atom_range <= c_maxNumLocalAtomsForInterdependentTask);
             if (tData.useInterdependentTask)
             {
                 size_t              natoms_use_in_vsites = vsite_atom_range;
@@ -3025,10 +3039,10 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
             else
             {
                 /* The last thread should cover up to the end of the range */
-                tData.rangeEnd = mdatoms.nr;
+                tData.rangeEnd = numAtoms;
             }
             assignVsitesToThread(
-                    &tData, thread, numThreads_, natperthread, taskIndex_, ilists, iparams, mdatoms.ptype);
+                    &tData, thread, numThreads_, natperthread, taskIndex_, ilists, iparams, ptype);
 
             if (tData.useInterdependentTask)
             {
@@ -3110,16 +3124,21 @@ void ThreadingInfo::setVirtualSites(ArrayRef<const InteractionList> ilists,
 }
 
 void VirtualSitesHandler::Impl::setVirtualSites(ArrayRef<const InteractionList> ilists,
-                                                const t_mdatoms&                mdatoms)
+                                                const int                       numAtoms,
+                                                const int                       homenr,
+                                                ArrayRef<const ParticleType>    ptype)
 {
     ilists_ = ilists;
 
-    threadingInfo_.setVirtualSites(ilists, iparams_, mdatoms, domainInfo_.useDomdec());
+    threadingInfo_.setVirtualSites(ilists, iparams_, numAtoms, homenr, ptype, domainInfo_.useDomdec());
 }
 
-void VirtualSitesHandler::setVirtualSites(ArrayRef<const InteractionList> ilists, const t_mdatoms& mdatoms)
+void VirtualSitesHandler::setVirtualSites(ArrayRef<const InteractionList> ilists,
+                                          const int                       numAtoms,
+                                          const int                       homenr,
+                                          ArrayRef<const ParticleType>    ptype)
 {
-    impl_->setVirtualSites(ilists, mdatoms);
+    impl_->setVirtualSites(ilists, numAtoms, homenr, ptype);
 }
 
 } // namespace gmx