Add InteractionDefinitions
[alexxy/gromacs.git] / src / gromacs / mdlib / vsite.cpp
index 74c82f64598d9a298624094441c119e2a6965bc5..1973d2088c2da7c09e2b066a4bc6d59179b82407 100644 (file)
  * Any remaining vsites are assigned to a separate master thread task.
  */
 
+using gmx::ArrayRef;
 using gmx::RVec;
 
-static void init_ilist(t_ilist* ilist)
-{
-    for (int i = 0; i < F_NRE; i++)
-    {
-        ilist[i].nr     = 0;
-        ilist[i].nalloc = 0;
-        ilist[i].iatoms = nullptr;
-    }
-}
-
 /*! \brief List of atom indices belonging to a task */
 struct AtomIndex
 {
@@ -109,7 +100,7 @@ struct AtomIndex
 struct InterdependentTask
 {
     //! The interaction lists, only vsite entries are used
-    t_ilist ilist[F_NRE];
+    InteractionLists ilist;
     //! Thread/task-local force buffer
     std::vector<RVec> force;
     //! The atom indices of the vsites of our task
@@ -117,19 +108,13 @@ struct InterdependentTask
     //! Flags if elements in force are spread to or not
     std::vector<bool> use;
     //! The number of entries set to true in use
-    int nuse;
+    int nuse = 0;
     //! Array of atoms indices, size nthreads, covering all nuse set elements in use
     std::vector<AtomIndex> atomIndex;
     //! List of tasks (force blocks) this task spread forces to
     std::vector<int> spreadTask;
     //! List of tasks that write to this tasks force block range
     std::vector<int> reduceTask;
-
-    InterdependentTask()
-    {
-        init_ilist(ilist);
-        nuse = 0;
-    }
 };
 
 /*! \brief Vsite thread task data structure */
@@ -140,7 +125,7 @@ struct VsiteThread
     //! End of atom range of this task
     int rangeEnd;
     //! The interaction lists, only vsite entries are used
-    t_ilist ilist[F_NRE];
+    std::array<InteractionList, F_NRE> ilist;
     //! Local fshift accumulation buffer
     rvec fshift[SHIFTS];
     //! Local virial dx*df accumulation buffer
@@ -155,7 +140,6 @@ struct VsiteThread
     {
         rangeStart = -1;
         rangeEnd   = -1;
-        init_ilist(ilist);
         clear_rvecs(SHIFTS, fshift);
         clear_mat(dxdf);
         useInterdependentTask = false;
@@ -166,8 +150,7 @@ struct VsiteThread
  *
  * \param[in] ilist  The interaction list
  */
-template<typename T>
-static int vsiteIlistNrCount(const T* ilist)
+static int vsiteIlistNrCount(ArrayRef<const InteractionList> ilist)
 {
     int nr = 0;
     for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
@@ -417,7 +400,7 @@ static void constr_vsite4FDN(const rvec   xi,
 }
 
 
-static int constr_vsiten(const t_iatom* ia, const t_iparams ip[], rvec* x, const t_pbc* pbc)
+static int constr_vsiten(const t_iatom* ia, ArrayRef<const t_iparams> ip, rvec* x, const t_pbc* pbc)
 {
     rvec x1, dx;
     dvec dsum;
@@ -477,12 +460,12 @@ static PbcMode getPbcMode(const t_pbc* pbcPtr)
     }
 }
 
-static void construct_vsites_thread(rvec            x[],
-                                    real            dt,
-                                    rvec*           v,
-                                    const t_iparams ip[],
-                                    const t_ilist   ilist[],
-                                    const t_pbc*    pbc_null)
+static void construct_vsites_thread(rvec                            x[],
+                                    real                            dt,
+                                    rvec*                           v,
+                                    ArrayRef<const t_iparams>       ip,
+                                    ArrayRef<const InteractionList> ilist,
+                                    const t_pbc*                    pbc_null)
 {
     real inv_dt;
     if (v != nullptr)
@@ -500,7 +483,7 @@ static void construct_vsites_thread(rvec            x[],
 
     for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
     {
-        if (ilist[ftype].nr == 0)
+        if (ilist[ftype].empty())
         {
             continue;
         }
@@ -508,9 +491,9 @@ static void construct_vsites_thread(rvec            x[],
         { // TODO remove me
             int nra = interaction_function[ftype].nratoms;
             int inc = 1 + nra;
-            int nr  = ilist[ftype].nr;
+            int nr  = ilist[ftype].size();
 
-            const t_iatom* ia = ilist[ftype].iatoms;
+            const t_iatom* ia = ilist[ftype].iatoms.data();
 
             for (int i = 0; i < nr;)
             {
@@ -609,16 +592,16 @@ static void construct_vsites_thread(rvec            x[],
     }
 }
 
-void construct_vsites(const gmx_vsite_t* vsite,
-                      rvec               x[],
-                      real               dt,
-                      rvec*              v,
-                      const t_iparams    ip[],
-                      const t_ilist      ilist[],
-                      PbcType            pbcType,
-                      gmx_bool           bMolPBC,
-                      const t_commrec*   cr,
-                      const matrix       box)
+void construct_vsites(const gmx_vsite_t*              vsite,
+                      rvec                            x[],
+                      real                            dt,
+                      rvec*                           v,
+                      ArrayRef<const t_iparams>       ip,
+                      ArrayRef<const InteractionList> ilist,
+                      PbcType                         pbcType,
+                      gmx_bool                        bMolPBC,
+                      const t_commrec*                cr,
+                      const matrix                    box)
 {
     const bool useDomdec = (vsite != nullptr && vsite->useDomdec);
     GMX_ASSERT(!useDomdec || (cr != nullptr && DOMAINDECOMP(cr)),
@@ -749,20 +732,13 @@ void constructVsitesGlobal(const gmx_mtop_t& mtop, gmx::ArrayRef<gmx::RVec> x)
     {
         const gmx_molblock_t& molb = mtop.molblock[mb];
         const gmx_moltype_t&  molt = mtop.moltype[molb.type];
-        if (vsiteIlistNrCount(molt.ilist.data()) > 0)
+        if (vsiteIlistNrCount(molt.ilist) > 0)
         {
             int atomOffset = mtop.moleculeBlockIndices[mb].globalAtomStart;
             for (int mol = 0; mol < molb.nmol; mol++)
             {
-                t_ilist ilist[F_NRE];
-                for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
-                {
-                    ilist[ftype].nr     = molt.ilist[ftype].size();
-                    ilist[ftype].iatoms = const_cast<t_iatom*>(molt.ilist[ftype].iatoms.data());
-                }
-
                 construct_vsites(nullptr, as_rvec_array(x.data()) + atomOffset, 0.0, nullptr,
-                                 mtop.ffparams.iparams.data(), ilist, PbcType::No, TRUE, nullptr, nullptr);
+                                 mtop.ffparams.iparams, molt.ilist, PbcType::No, TRUE, nullptr, nullptr);
                 atomOffset += molt.atoms.nr;
             }
         }
@@ -1553,13 +1529,13 @@ static void spread_vsite4FDN(const t_iatom  ia[],
 }
 
 
-static int spread_vsiten(const t_iatom   ia[],
-                         const t_iparams ip[],
-                         const rvec      x[],
-                         rvec            f[],
-                         rvec            fshift[],
-                         const t_pbc*    pbc,
-                         const t_graph*  g)
+static int spread_vsiten(const t_iatom             ia[],
+                         ArrayRef<const t_iparams> ip,
+                         const rvec                x[],
+                         rvec                      f[],
+                         rvec                      fshift[],
+                         const t_pbc*              pbc,
+                         const t_graph*            g)
 {
     rvec xv, dx, fi;
     int  n3, av, i, ai;
@@ -1602,27 +1578,27 @@ static int spread_vsiten(const t_iatom   ia[],
 }
 
 
-static int vsite_count(const t_ilist* ilist, int ftype)
+static int vsite_count(ArrayRef<const InteractionList> ilist, int ftype)
 {
     if (ftype == F_VSITEN)
     {
-        return ilist[ftype].nr / 3;
+        return ilist[ftype].size() / 3;
     }
     else
     {
-        return ilist[ftype].nr / (1 + interaction_function[ftype].nratoms);
+        return ilist[ftype].size() / (1 + interaction_function[ftype].nratoms);
     }
 }
 
-static void spread_vsite_f_thread(const rvec     x[],
-                                  rvec           f[],
-                                  rvec*          fshift,
-                                  gmx_bool       VirCorr,
-                                  matrix         dxdf,
-                                  t_iparams      ip[],
-                                  const t_ilist  ilist[],
-                                  const t_graph* g,
-                                  const t_pbc*   pbc_null)
+static void spread_vsite_f_thread(const rvec                      x[],
+                                  rvec                            f[],
+                                  rvec*                           fshift,
+                                  gmx_bool                        VirCorr,
+                                  matrix                          dxdf,
+                                  ArrayRef<const t_iparams>       ip,
+                                  ArrayRef<const InteractionList> ilist,
+                                  const t_graph*                  g,
+                                  const t_pbc*                    pbc_null)
 {
     const PbcMode pbcMode = getPbcMode(pbc_null);
     /* We need another pbc pointer, as with charge groups we switch per vsite */
@@ -1633,7 +1609,7 @@ static void spread_vsite_f_thread(const rvec     x[],
      * higher type vsites from lower types         */
     for (int ftype = c_ftypeVsiteEnd - 1; ftype >= c_ftypeVsiteStart; ftype--)
     {
-        if (ilist[ftype].nr == 0)
+        if (ilist[ftype].empty())
         {
             continue;
         }
@@ -1641,9 +1617,9 @@ static void spread_vsite_f_thread(const rvec     x[],
         { // TODO remove me
             int nra = interaction_function[ftype].nratoms;
             int inc = 1 + nra;
-            int nr  = ilist[ftype].nr;
+            int nr  = ilist[ftype].size();
 
-            const t_iatom* ia = ilist[ftype].iatoms;
+            const t_iatom* ia = ilist[ftype].iatoms.data();
 
             if (pbcMode == PbcMode::all)
             {
@@ -1724,17 +1700,17 @@ static void clearTaskForceBufferUsedElements(InterdependentTask* idTask)
 void spread_vsite_f(const gmx_vsite_t* vsite,
                     const rvec* gmx_restrict x,
                     rvec* gmx_restrict f,
-                    rvec* gmx_restrict fshift,
-                    gmx_bool           VirCorr,
-                    matrix             vir,
-                    t_nrnb*            nrnb,
-                    const t_idef*      idef,
-                    PbcType            pbcType,
-                    gmx_bool           bMolPBC,
-                    const t_graph*     g,
-                    const matrix       box,
-                    const t_commrec*   cr,
-                    gmx_wallcycle*     wcycle)
+                    rvec* gmx_restrict            fshift,
+                    gmx_bool                      VirCorr,
+                    matrix                        vir,
+                    t_nrnb*                       nrnb,
+                    const InteractionDefinitions& idef,
+                    PbcType                       pbcType,
+                    gmx_bool                      bMolPBC,
+                    const t_graph*                g,
+                    const matrix                  box,
+                    const t_commrec*              cr,
+                    gmx_wallcycle*                wcycle)
 {
     wallcycle_start(wcycle, ewcVSITESPREAD);
     const bool useDomdec = vsite->useDomdec;
@@ -1768,7 +1744,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite,
         {
             clear_mat(dxdf);
         }
-        spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef->iparams, idef->il, g, pbc_null);
+        spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef.iparams, idef.il, g, pbc_null);
 
         if (VirCorr)
         {
@@ -1789,7 +1765,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite,
             clear_mat(vsite->tData[vsite->nthreads]->dxdf);
         }
         spread_vsite_f_thread(x, f, fshift, VirCorr, vsite->tData[vsite->nthreads]->dxdf,
-                              idef->iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null);
+                              idef.iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null);
 
 #pragma omp parallel num_threads(vsite->nthreads)
         {
@@ -1836,7 +1812,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite,
                         copy_rvec(f[idTask->vsite[i]], idTask->force[idTask->vsite[i]]);
                     }
                     spread_vsite_f_thread(x, as_rvec_array(idTask->force.data()), fshift_t, VirCorr,
-                                          tData.dxdf, idef->iparams, tData.idTask.ilist, g, pbc_null);
+                                          tData.dxdf, idef.iparams, tData.idTask.ilist, g, pbc_null);
 
                     /* We need a barrier before reducing forces below
                      * that have been produced by a different thread above.
@@ -1874,7 +1850,7 @@ void spread_vsite_f(const gmx_vsite_t* vsite,
                 }
 
                 /* Spread the vsites that spread locally only */
-                spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef->iparams,
+                spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef.iparams,
                                       tData.ilist, g, pbc_null);
             }
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
@@ -1914,15 +1890,15 @@ void spread_vsite_f(const gmx_vsite_t* vsite,
         dd_move_f_vsites(cr->dd, f, fshift);
     }
 
-    inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef->il, F_VSITE2));
-    inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef->il, F_VSITE2FD));
-    inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef->il, F_VSITE3));
-    inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef->il, F_VSITE3FD));
-    inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef->il, F_VSITE3FAD));
-    inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef->il, F_VSITE3OUT));
-    inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef->il, F_VSITE4FD));
-    inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef->il, F_VSITE4FDN));
-    inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef->il, F_VSITEN));
+    inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef.il, F_VSITE2));
+    inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef.il, F_VSITE2FD));
+    inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef.il, F_VSITE3));
+    inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef.il, F_VSITE3FD));
+    inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef.il, F_VSITE3FAD));
+    inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef.il, F_VSITE3OUT));
+    inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef.il, F_VSITE4FD));
+    inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef.il, F_VSITE4FDN));
+    inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef.il, F_VSITEN));
 
     wallcycle_stop(wcycle, ewcVSITESPREAD);
 }
@@ -2103,24 +2079,24 @@ static inline void flagAtom(InterdependentTask* idTask, int atom, int thread, in
  * taskIndex[] is set for all vsites in our range, either to our local tasks
  * or to the single last task as taskIndex[]=2*nthreads.
  */
-static void assignVsitesToThread(VsiteThread*          tData,
-                                 int                   thread,
-                                 int                   nthread,
-                                 int                   natperthread,
-                                 gmx::ArrayRef<int>    taskIndex,
-                                 const t_ilist*        ilist,
-                                 const t_iparams*      ip,
-                                 const unsigned short* ptype)
+static void assignVsitesToThread(VsiteThread*                    tData,
+                                 int                             thread,
+                                 int                             nthread,
+                                 int                             natperthread,
+                                 gmx::ArrayRef<int>              taskIndex,
+                                 ArrayRef<const InteractionList> ilist,
+                                 ArrayRef<const t_iparams>       ip,
+                                 const unsigned short*           ptype)
 {
     for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
     {
-        tData->ilist[ftype].nr        = 0;
-        tData->idTask.ilist[ftype].nr = 0;
+        tData->ilist[ftype].clear();
+        tData->idTask.ilist[ftype].clear();
 
-        int      nral1 = 1 + NRAL(ftype);
-        int      inc   = nral1;
-        t_iatom* iat   = ilist[ftype].iatoms;
-        for (int i = 0; i < ilist[ftype].nr;)
+        int        nral1 = 1 + NRAL(ftype);
+        int        inc   = nral1;
+        const int* iat   = ilist[ftype].iatoms.data();
+        for (int i = 0; i < ilist[ftype].size();)
         {
             if (ftype == F_VSITEN)
             {
@@ -2193,7 +2169,7 @@ static void assignVsitesToThread(VsiteThread*          tData,
             if (task == thread || task == nthread + thread)
             {
                 /* Copy this vsite to the thread data struct of thread */
-                t_ilist* il_task;
+                InteractionList* il_task;
                 if (task == thread)
                 {
                     il_task = &tData->ilist[ftype];
@@ -2202,17 +2178,8 @@ static void assignVsitesToThread(VsiteThread*          tData,
                 {
                     il_task = &tData->idTask.ilist[ftype];
                 }
-                /* Ensure we have sufficient memory allocated */
-                if (il_task->nr + inc > il_task->nalloc)
-                {
-                    il_task->nalloc = over_alloc_large(il_task->nr + inc);
-                    srenew(il_task->iatoms, il_task->nalloc);
-                }
                 /* Copy the vsite data to the thread-task local array */
-                for (int j = i; j < i + inc; j++)
-                {
-                    il_task->iatoms[il_task->nr++] = iat[j];
-                }
+                il_task->push_back(iat[i], nral1 - 1, iat + i + 1);
                 if (task == nthread + thread)
                 {
                     /* This vsite write outside our own task force block.
@@ -2243,23 +2210,23 @@ static void assignVsitesToThread(VsiteThread*          tData,
 }
 
 /*! \brief Assign all vsites with taskIndex[]==task to task tData */
-static void assignVsitesToSingleTask(VsiteThread*             tData,
-                                     int                      task,
-                                     gmx::ArrayRef<const int> taskIndex,
-                                     const t_ilist*           ilist,
-                                     const t_iparams*         ip)
+static void assignVsitesToSingleTask(VsiteThread*                    tData,
+                                     int                             task,
+                                     gmx::ArrayRef<const int>        taskIndex,
+                                     ArrayRef<const InteractionList> ilist,
+                                     ArrayRef<const t_iparams>       ip)
 {
     for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
     {
-        tData->ilist[ftype].nr        = 0;
-        tData->idTask.ilist[ftype].nr = 0;
+        tData->ilist[ftype].clear();
+        tData->idTask.ilist[ftype].clear();
 
-        int      nral1   = 1 + NRAL(ftype);
-        int      inc     = nral1;
-        t_iatom* iat     = ilist[ftype].iatoms;
-        t_ilist* il_task = &tData->ilist[ftype];
+        int              nral1   = 1 + NRAL(ftype);
+        int              inc     = nral1;
+        const int*       iat     = ilist[ftype].iatoms.data();
+        InteractionList* il_task = &tData->ilist[ftype];
 
-        for (int i = 0; i < ilist[ftype].nr;)
+        for (int i = 0; i < ilist[ftype].size();)
         {
             if (ftype == F_VSITEN)
             {
@@ -2269,17 +2236,8 @@ static void assignVsitesToSingleTask(VsiteThread*             tData,
             /* Check if the vsite is assigned to our task */
             if (taskIndex[iat[1 + i]] == task)
             {
-                /* Ensure we have sufficient memory allocated */
-                if (il_task->nr + inc > il_task->nalloc)
-                {
-                    il_task->nalloc = over_alloc_large(il_task->nr + inc);
-                    srenew(il_task->iatoms, il_task->nalloc);
-                }
                 /* Copy the vsite data to the thread-task local array */
-                for (int j = i; j < i + inc; j++)
-                {
-                    il_task->iatoms[il_task->nr++] = iat[j];
-                }
+                il_task->push_back(iat[i], inc - 1, iat + i + 1);
             }
 
             i += inc;
@@ -2287,7 +2245,10 @@ static void assignVsitesToSingleTask(VsiteThread*             tData,
     }
 }
 
-void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const t_mdatoms* mdatoms, gmx_vsite_t* vsite)
+void split_vsites_over_threads(ArrayRef<const InteractionList> ilist,
+                               ArrayRef<const t_iparams>       ip,
+                               const t_mdatoms*                mdatoms,
+                               gmx_vsite_t*                    vsite)
 {
     int vsite_atom_range, natperthread;
 
@@ -2315,9 +2276,9 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const
             { // TODO remove me
                 if (ftype != F_VSITEN)
                 {
-                    int            nral1 = 1 + NRAL(ftype);
-                    const t_iatom* iat   = ilist[ftype].iatoms;
-                    for (int i = 0; i < ilist[ftype].nr; i += nral1)
+                    int                 nral1 = 1 + NRAL(ftype);
+                    ArrayRef<const int> iat   = ilist[ftype].iatoms;
+                    for (int i = 0; i < ilist[ftype].size(); i += nral1)
                     {
                         for (int j = i + 1; j < i + nral1; j++)
                         {
@@ -2329,10 +2290,10 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const
                 {
                     int vs_ind_end;
 
-                    const t_iatom* iat = ilist[ftype].iatoms;
+                    ArrayRef<const int> iat = ilist[ftype].iatoms;
 
                     int i = 0;
-                    while (i < ilist[ftype].nr)
+                    while (i < ilist[ftype].size())
                     {
                         /* The 3 below is from 1+NRAL(ftype)=3 */
                         vs_ind_end = i + ip[iat[i]].vsiten.n * 3;
@@ -2518,13 +2479,13 @@ void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const
 
         for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
         {
-            if (ilist[ftype].nr > 0)
+            if (!ilist[ftype].empty())
             {
                 fprintf(debug, "%-20s thread dist:", interaction_function[ftype].longname);
                 for (int th = 0; th < vsite->nthreads + 1; th++)
                 {
-                    fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].nr,
-                            vsite->tData[th]->idTask.ilist[ftype].nr);
+                    fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].size(),
+                            vsite->tData[th]->idTask.ilist[ftype].size());
                 }
                 fprintf(debug, "\n");
             }