Partial conversion of gmx_domdec_t to C++
authorBerk Hess <hess@kth.se>
Mon, 14 May 2018 13:55:36 +0000 (15:55 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Wed, 13 Jun 2018 12:40:27 +0000 (14:40 +0200)
Most of the simple buffers in gmx_domdec_t now use std::vector.
Changed AtomDistribution to use std::unique_ptr.
Now new is used for gmx_domdec_t with some default zero values.

Renamed charge groups to atom groups for all modified code
in preparation for the introduction of update groups.

Change-Id: Ie54b2e6fc354cd2fbc874023922369261b3cd092

17 files changed:
src/gromacs/domdec/collect.cpp
src/gromacs/domdec/distribute.cpp
src/gromacs/domdec/domdec.cpp
src/gromacs/domdec/domdec_constraints.cpp
src/gromacs/domdec/domdec_specatomcomm.cpp
src/gromacs/domdec/domdec_struct.h
src/gromacs/domdec/domdec_topology.cpp
src/gromacs/domdec/redistribute.cpp
src/gromacs/mdlib/constr.cpp
src/gromacs/mdlib/coupling.cpp
src/gromacs/mdlib/force.cpp
src/gromacs/mdlib/mdsetup.cpp
src/gromacs/mdlib/shellfc.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdlib/update.cpp
src/gromacs/mdlib/wnblist.cpp
src/gromacs/mdrun/minimize.cpp

index f3a90cce4f148fb86a29f6d84d1bde5ba4e7bffa..ecce553e3ae78eb6c30ae268bd11a806d5b575ee 100644 (file)
@@ -71,7 +71,7 @@ static void dd_collect_cg(gmx_domdec_t  *dd,
     if (state_local->ddp_count == dd->ddp_count)
     {
         /* The local state and DD are in sync, use the DD indices */
-        atomGroups = gmx::constArrayRefFromArray(dd->index_gl, dd->ncg_home);
+        atomGroups = gmx::constArrayRefFromArray(dd->globalAtomGroupIndices.data(), dd->ncg_home);
         nat_home   = dd->nat_home;
     }
     else if (state_local->ddp_count_cg_gl == state_local->ddp_count)
@@ -93,7 +93,7 @@ static void dd_collect_cg(gmx_domdec_t  *dd,
         gmx_incons("Attempted to collect a vector for a state for which the charge group distribution is unknown");
     }
 
-    AtomDistribution *ma = dd->ma;
+    AtomDistribution *ma = dd->ma.get();
 
     /* Collect the charge group and atom counts on the master */
     int localBuffer[2] = { static_cast<int>(atomGroups.size()), nat_home };
@@ -150,8 +150,6 @@ static void dd_collect_vec_sendrecv(gmx_domdec_t                  *dd,
                                     gmx::ArrayRef<const gmx::RVec> lv,
                                     gmx::ArrayRef<gmx::RVec>       v)
 {
-    AtomDistribution *ma = dd->ma;
-
     if (!DDMASTER(dd))
     {
 #if GMX_MPI
@@ -161,12 +159,14 @@ static void dd_collect_vec_sendrecv(gmx_domdec_t                  *dd,
     }
     else
     {
+        AtomDistribution &ma = *dd->ma.get();
+
         /* Copy the master coordinates to the global array */
         const t_block &cgs_gl    = dd->comm->cgs_gl;
 
         int            rank      = dd->masterrank;
         int            localAtom = 0;
-        for (const int &i : ma->domainGroups[rank].atomGroups)
+        for (const int &i : ma.domainGroups[rank].atomGroups)
         {
             for (int globalAtom = cgs_gl.index[i]; globalAtom < cgs_gl.index[i + 1]; globalAtom++)
             {
@@ -178,20 +178,20 @@ static void dd_collect_vec_sendrecv(gmx_domdec_t                  *dd,
         {
             if (rank != dd->rank)
             {
-                const auto &domainGroups = ma->domainGroups[rank];
+                const auto &domainGroups = ma.domainGroups[rank];
 
-                GMX_RELEASE_ASSERT(v.data() != ma->rvecBuffer.data(), "We need different communication and return buffers");
+                GMX_RELEASE_ASSERT(v.data() != ma.rvecBuffer.data(), "We need different communication and return buffers");
 
                 /* When we send/recv instead of scatter/gather, we might need
                  * to increase the communication buffer size here.
                  */
-                if (static_cast<size_t>(domainGroups.numAtoms) > ma->rvecBuffer.size())
+                if (static_cast<size_t>(domainGroups.numAtoms) > ma.rvecBuffer.size())
                 {
-                    ma->rvecBuffer.resize(domainGroups.numAtoms);
+                    ma.rvecBuffer.resize(domainGroups.numAtoms);
                 }
 
 #if GMX_MPI
-                MPI_Recv(ma->rvecBuffer.data(), domainGroups.numAtoms*sizeof(rvec), MPI_BYTE, rank,
+                MPI_Recv(ma.rvecBuffer.data(), domainGroups.numAtoms*sizeof(rvec), MPI_BYTE, rank,
                          rank, dd->mpi_comm_all, MPI_STATUS_IGNORE);
 #endif
                 int localAtom = 0;
@@ -199,7 +199,7 @@ static void dd_collect_vec_sendrecv(gmx_domdec_t                  *dd,
                 {
                     for (int globalAtom = cgs_gl.index[cg]; globalAtom < cgs_gl.index[cg + 1]; globalAtom++)
                     {
-                        copy_rvec(ma->rvecBuffer[localAtom++], v[globalAtom]);
+                        copy_rvec(ma.rvecBuffer[localAtom++], v[globalAtom]);
                     }
                 }
             }
@@ -216,7 +216,7 @@ static void dd_collect_vec_gatherv(gmx_domdec_t                  *dd,
 
     if (DDMASTER(dd))
     {
-        get_commbuffer_counts(dd->ma, &recvCounts, &displacements);
+        get_commbuffer_counts(dd->ma.get(), &recvCounts, &displacements);
     }
 
     dd_gatherv(dd, dd->nat_home*sizeof(rvec), lv.data(), recvCounts, displacements,
index 8149f01755967cffa03a304cdbac40c135ff9d09..fd374eef215c5d249178de2d5a5293bcfe20db25 100644 (file)
 #include "domdec_internal.h"
 #include "utility.h"
 
-static void dd_distribute_vec_sendrecv(gmx_domdec_t *dd, const t_block *cgs,
-                                       const rvec *v, rvec *lv)
+static void distributeVecSendrecv(gmx_domdec_t                   *dd,
+                                  const t_block                  *cgs,
+                                  gmx::ArrayRef<const gmx::RVec>  globalVec,
+                                  gmx::ArrayRef<gmx::RVec>        localVec)
 {
-    AtomDistribution *ma;
-    int               nalloc = 0;
-    rvec             *buf    = nullptr;
-
     if (DDMASTER(dd))
     {
-        ma  = dd->ma;
+        std::vector<gmx::RVec> buffer;
 
         for (int rank = 0; rank < dd->nnodes; rank++)
         {
             if (rank != dd->rank)
             {
-                const auto &domainGroups = ma->domainGroups[rank];
+                const auto &domainGroups = dd->ma->domainGroups[rank];
 
-                if (domainGroups.numAtoms > nalloc)
-                {
-                    nalloc = over_alloc_dd(domainGroups.numAtoms);
-                    srenew(buf, nalloc);
-                }
-                /* Use lv as a temporary buffer */
-                int localAtom = 0;
+                buffer.resize(domainGroups.numAtoms);
+
+                int   localAtom = 0;
                 for (const int &cg : domainGroups.atomGroups)
                 {
                     for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
                     {
-                        copy_rvec(v[globalAtom], buf[localAtom++]);
+                        buffer[localAtom++] = globalVec[globalAtom];
                     }
                 }
                 GMX_RELEASE_ASSERT(localAtom == domainGroups.numAtoms, "The index count and number of indices should match");
 
 #if GMX_MPI
-                MPI_Send(buf, domainGroups.numAtoms*sizeof(rvec), MPI_BYTE,
+                MPI_Send(buffer.data(), domainGroups.numAtoms*sizeof(gmx::RVec), MPI_BYTE,
                          rank, rank, dd->mpi_comm_all);
 #endif
             }
         }
-        sfree(buf);
 
-        const auto &domainGroups = ma->domainGroups[dd->masterrank];
+        const auto &domainGroups = dd->ma->domainGroups[dd->masterrank];
         int         localAtom    = 0;
         for (const int &cg : domainGroups.atomGroups)
         {
             for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
             {
-                copy_rvec(v[globalAtom], lv[localAtom++]);
+                localVec[localAtom++] = globalVec[globalAtom];
             }
         }
     }
     else
     {
 #if GMX_MPI
-        MPI_Recv(lv, dd->nat_home*sizeof(rvec), MPI_BYTE, dd->masterrank,
+        MPI_Recv(localVec.data(), dd->nat_home*sizeof(gmx::RVec), MPI_BYTE, dd->masterrank,
                  MPI_ANY_TAG, dd->mpi_comm_all, MPI_STATUS_IGNORE);
 #endif
     }
 }
 
-static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, const t_block *cgs,
-                                       const rvec *v, rvec *lv)
+static void distributeVecScatterv(gmx_domdec_t                   *dd,
+                                  const t_block                  *cgs,
+                                  gmx::ArrayRef<const gmx::RVec>  globalVec,
+                                  gmx::ArrayRef<gmx::RVec>        localVec)
 {
     int *sendCounts    = nullptr;
     int *displacements = nullptr;
 
     if (DDMASTER(dd))
     {
-        AtomDistribution *ma = dd->ma;
+        AtomDistribution &ma = *dd->ma.get();
 
-        get_commbuffer_counts(ma, &sendCounts, &displacements);
+        get_commbuffer_counts(&ma, &sendCounts, &displacements);
 
-        gmx::ArrayRef<gmx::RVec> buffer = ma->rvecBuffer;
+        gmx::ArrayRef<gmx::RVec> buffer = ma.rvecBuffer;
         int localAtom                   = 0;
         for (int rank = 0; rank < dd->nnodes; rank++)
         {
-            const auto &domainGroups = ma->domainGroups[rank];
+            const auto &domainGroups = ma.domainGroups[rank];
             for (const int &cg : domainGroups.atomGroups)
             {
                 for (int globalAtom = cgs->index[cg]; globalAtom < cgs->index[cg + 1]; globalAtom++)
                 {
-                    copy_rvec(v[globalAtom], buffer[localAtom++]);
+                    buffer[localAtom++] = globalVec[globalAtom];
                 }
             }
         }
@@ -150,19 +145,21 @@ static void dd_distribute_vec_scatterv(gmx_domdec_t *dd, const t_block *cgs,
 
     dd_scatterv(dd, sendCounts, displacements,
                 DDMASTER(dd) ? dd->ma->rvecBuffer.data() : nullptr,
-                dd->nat_home*sizeof(rvec), lv);
+                dd->nat_home*sizeof(gmx::RVec), localVec.data());
 }
 
-static void dd_distribute_vec(gmx_domdec_t *dd, const t_block *cgs,
-                              const rvec *v, rvec *lv)
+static void distributeVec(gmx_domdec_t                   *dd,
+                          const t_block                  *cgs,
+                          gmx::ArrayRef<const gmx::RVec>  globalVec,
+                          gmx::ArrayRef<gmx::RVec>        localVec)
 {
     if (dd->nnodes <= c_maxNumRanksUseSendRecvForScatterAndGather)
     {
-        dd_distribute_vec_sendrecv(dd, cgs, v, lv);
+        distributeVecSendrecv(dd, cgs, globalVec, localVec);
     }
     else
     {
-        dd_distribute_vec_scatterv(dd, cgs, v, lv);
+        distributeVecScatterv(dd, cgs, globalVec, localVec);
     }
 }
 
@@ -266,18 +263,15 @@ static void dd_distribute_state(gmx_domdec_t *dd, const t_block *cgs,
 
     if (state_local->flags & (1 << estX))
     {
-        const rvec *xGlobal = (DDMASTER(dd) ? as_rvec_array(state->x.data()) : nullptr);
-        dd_distribute_vec(dd, cgs, xGlobal, as_rvec_array(state_local->x.data()));
+        distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->x) : gmx::EmptyArrayRef(), state_local->x);
     }
     if (state_local->flags & (1 << estV))
     {
-        const rvec *vGlobal = (DDMASTER(dd) ? as_rvec_array(state->v.data()) : nullptr);
-        dd_distribute_vec(dd, cgs, vGlobal, as_rvec_array(state_local->v.data()));
+        distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->v) : gmx::EmptyArrayRef(), state_local->v);
     }
     if (state_local->flags & (1 << estCGP))
     {
-        const rvec *cgpGlobal = (DDMASTER(dd) ? as_rvec_array(state->cg_p.data()) : nullptr);
-        dd_distribute_vec(dd, cgs, cgpGlobal, as_rvec_array(state_local->cg_p.data()));
+        distributeVec(dd, cgs, DDMASTER(dd) ? makeArrayRef(state->cg_p) : gmx::EmptyArrayRef(), state_local->cg_p);
     }
 }
 
@@ -287,12 +281,12 @@ getAtomGroupDistribution(FILE *fplog,
                          const t_block *cgs, rvec pos[],
                          gmx_domdec_t *dd)
 {
-    AtomDistribution *ma = dd->ma;
+    AtomDistribution &ma = *dd->ma.get();
 
     /* Clear the count */
     for (int rank = 0; rank < dd->nnodes; rank++)
     {
-        ma->domainGroups[rank].numAtoms = 0;
+        ma.domainGroups[rank].numAtoms = 0;
     }
 
     matrix tcm;
@@ -396,7 +390,7 @@ getAtomGroupDistribution(FILE *fplog,
         }
         int domainIndex = dd_index(dd->nc, ind);
         indices[domainIndex].push_back(icg);
-        ma->domainGroups[domainIndex].numAtoms += cgindex[icg + 1] - cgindex[icg];
+        ma.domainGroups[domainIndex].numAtoms += cgindex[icg + 1] - cgindex[icg];
     }
 
     if (fplog)
@@ -408,11 +402,11 @@ getAtomGroupDistribution(FILE *fplog,
 
         nat_sum  = 0;
         nat2_sum = 0;
-        nat_min  = ma->domainGroups[0].numAtoms;
-        nat_max  = ma->domainGroups[0].numAtoms;
+        nat_min  = ma.domainGroups[0].numAtoms;
+        nat_max  = ma.domainGroups[0].numAtoms;
         for (int rank = 0; rank < dd->nnodes; rank++)
         {
-            int numAtoms  = ma->domainGroups[rank].numAtoms;
+            int numAtoms  = ma.domainGroups[rank].numAtoms;
             nat_sum      += numAtoms;
             // cast to double to avoid integer overflows when squaring
             nat2_sum     += gmx::square(static_cast<double>(numAtoms));
@@ -437,7 +431,7 @@ static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
                                  const matrix box, const gmx_ddbox_t *ddbox,
                                  rvec pos[])
 {
-    AtomDistribution *ma = nullptr;
+    AtomDistribution *ma = dd->ma.get();
     int              *ibuf, buf2[2] = { 0, 0 };
     gmx_bool          bMaster = DDMASTER(dd);
 
@@ -445,8 +439,6 @@ static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
 
     if (bMaster)
     {
-        ma = dd->ma;
-
         if (dd->bScrewPBC)
         {
             check_screw_box(box);
@@ -469,14 +461,8 @@ static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
 
     dd->ncg_home = buf2[0];
     dd->nat_home = buf2[1];
-    dd->ncg_tot  = dd->ncg_home;
-    dd->nat_tot  = dd->nat_home;
-    if (dd->ncg_home > dd->cg_nalloc || dd->cg_nalloc == 0)
-    {
-        dd->cg_nalloc = over_alloc_dd(dd->ncg_home);
-        srenew(dd->index_gl, dd->cg_nalloc);
-        srenew(dd->cgindex, dd->cg_nalloc+1);
-    }
+    dd->globalAtomGroupIndices.resize(dd->ncg_home);
+    dd->globalAtomIndices.resize(dd->nat_home);
 
     if (bMaster)
     {
@@ -501,15 +487,16 @@ static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
                 bMaster ? ma->intBuffer.data() : nullptr,
                 bMaster ? ma->intBuffer.data() + dd->nnodes : nullptr,
                 bMaster ? ma->atomGroups.data() : nullptr,
-                dd->ncg_home*sizeof(int), dd->index_gl);
+                dd->ncg_home*sizeof(int), dd->globalAtomGroupIndices.data());
 
     /* Determine the home charge group sizes */
-    dd->cgindex[0] = 0;
+    dd->atomGroups_.index.resize(dd->ncg_home + 1);
+    dd->atomGroups_.index[0] = 0;
     for (int i = 0; i < dd->ncg_home; i++)
     {
-        int cg_gl        = dd->index_gl[i];
-        dd->cgindex[i+1] =
-            dd->cgindex[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
+        int cg_gl        = dd->globalAtomGroupIndices[i];
+        dd->atomGroups_.index[i+1] =
+            dd->atomGroups_.index[i] + cgs->index[cg_gl+1] - cgs->index[cg_gl];
     }
 
     if (debug)
@@ -517,7 +504,7 @@ static void distributeAtomGroups(FILE *fplog, gmx_domdec_t *dd,
         fprintf(debug, "Home charge groups:\n");
         for (int i = 0; i < dd->ncg_home; i++)
         {
-            fprintf(debug, " %d", dd->index_gl[i]);
+            fprintf(debug, " %d", dd->globalAtomGroupIndices[i]);
             if (i % 10 == 9)
             {
                 fprintf(debug, "\n");
index 01dcc9e86013d1a26b8d5d276f4cd0878b16412a..c8dbba14839a20d324a85a58ba3522fd726afc98 100644 (file)
@@ -49,6 +49,7 @@
 
 #include <algorithm>
 
+#include "gromacs/compat/make_unique.h"
 #include "gromacs/domdec/collect.h"
 #include "gromacs/domdec/dlbtiming.h"
 #include "gromacs/domdec/domdec_network.h"
@@ -235,7 +236,7 @@ int ddglatnr(const gmx_domdec_t *dd, int i)
         {
             gmx_fatal(FARGS, "glatnr called with %d, which is larger than the local number of atoms (%d)", i, dd->comm->nat[ddnatNR-1]);
         }
-        atnr = dd->gatindex[i] + 1;
+        atnr = dd->globalAtomIndices[i] + 1;
     }
 
     return atnr;
@@ -264,7 +265,7 @@ void dd_store_state(gmx_domdec_t *dd, t_state *state)
     state->cg_gl.resize(dd->ncg_home);
     for (i = 0; i < dd->ncg_home; i++)
     {
-        state->cg_gl[i] = dd->index_gl[i];
+        state->cg_gl[i] = dd->globalAtomGroupIndices[i];
     }
 
     state->ddp_count_cg_gl = dd->ddp_count;
@@ -343,23 +344,21 @@ void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[], gmx_wallcycle *wcycle)
 {
     wallcycle_start(wcycle, ewcMOVEX);
 
-    int                    nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
-    int                   *index, *cgindex;
+    int                    nzone, nat_tot;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_comm_dim_t *cd;
-    gmx_domdec_ind_t      *ind;
     rvec                   shift = {0, 0, 0}, *buf, *rbuf;
     gmx_bool               bPBC, bScrew;
 
     comm = dd->comm;
 
-    cgindex = dd->cgindex;
+    const BlockRanges &atomGroups = dd->atomGroups();
 
     buf = comm->vbuf.v;
 
     nzone   = 1;
     nat_tot = dd->nat_home;
-    for (d = 0; d < dd->ndim; d++)
+    for (int d = 0; d < dd->ndim; d++)
     {
         bPBC   = (dd->ci[dd->dim[d]] == 0);
         bScrew = (bPBC && dd->bScrewPBC && dd->dim[d] == XX);
@@ -368,18 +367,18 @@ void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[], gmx_wallcycle *wcycle)
             copy_rvec(box[dd->dim[d]], shift);
         }
         cd = &comm->cd[d];
-        for (p = 0; p < cd->np; p++)
+        for (int p = 0; p < cd->np; p++)
         {
-            ind   = &cd->ind[p];
-            index = ind->index;
-            n     = 0;
+            const gmx_domdec_ind_t *ind   = &cd->ind[p];
+            const int              *index = ind->index;
+            int                     n     = 0;
             if (!bPBC)
             {
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i] + 1];
+                    for (int j = at0; j < at1; j++)
                     {
                         copy_rvec(x[j], buf[n]);
                         n++;
@@ -388,11 +387,11 @@ void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[], gmx_wallcycle *wcycle)
             }
             else if (!bScrew)
             {
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i] + 1];
+                    for (int j = at0; j < at1; j++)
                     {
                         /* We need to shift the coordinates */
                         rvec_add(x[j], shift, buf[n]);
@@ -402,11 +401,11 @@ void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[], gmx_wallcycle *wcycle)
             }
             else
             {
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i]+1];
+                    for (int j = at0; j < at1; j++)
                     {
                         /* Shift x */
                         buf[n][XX] = x[j][XX] + shift[XX];
@@ -435,10 +434,10 @@ void dd_move_x(gmx_domdec_t *dd, matrix box, rvec x[], gmx_wallcycle *wcycle)
                              rbuf, ind->nrecv[nzone+1]);
             if (!cd->bInPlace)
             {
-                j = 0;
-                for (zone = 0; zone < nzone; zone++)
+                int j = 0;
+                for (int zone = 0; zone < nzone; zone++)
                 {
-                    for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
+                    for (int i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
                     {
                         copy_rvec(rbuf[j], x[i]);
                         j++;
@@ -457,11 +456,9 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
 {
     wallcycle_start(wcycle, ewcMOVEF);
 
-    int                    nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
-    int                   *index, *cgindex;
+    int                    nzone, nat_tot;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_comm_dim_t *cd;
-    gmx_domdec_ind_t      *ind;
     rvec                  *buf, *sbuf;
     ivec                   vis;
     int                    is;
@@ -469,13 +466,13 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
 
     comm = dd->comm;
 
-    cgindex = dd->cgindex;
+    const BlockRanges &atomGroups = dd->atomGroups();
 
     buf = comm->vbuf.v;
 
     nzone   = comm->zones.n/2;
-    nat_tot = dd->nat_tot;
-    for (d = dd->ndim-1; d >= 0; d--)
+    nat_tot = comm->nat[ddnatZONE];
+    for (int d = dd->ndim-1; d >= 0; d--)
     {
         /* Only forces in domains near the PBC boundaries need to
            consider PBC in the treatment of fshift */
@@ -491,21 +488,21 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
         is              = IVEC2IS(vis);
 
         cd = &comm->cd[d];
-        for (p = cd->np-1; p >= 0; p--)
+        for (int p = cd->np-1; p >= 0; p--)
         {
-            ind      = &cd->ind[p];
-            nat_tot -= ind->nrecv[nzone+1];
+            const gmx_domdec_ind_t *ind  = &cd->ind[p];
+            nat_tot                     -= ind->nrecv[nzone+1];
             if (cd->bInPlace)
             {
                 sbuf = f + nat_tot;
             }
             else
             {
-                sbuf = comm->vbuf2.v;
-                j    = 0;
-                for (zone = 0; zone < nzone; zone++)
+                sbuf  = comm->vbuf2.v;
+                int j = 0;
+                for (int zone = 0; zone < nzone; zone++)
                 {
-                    for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
+                    for (int i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
                     {
                         copy_rvec(f[i], sbuf[j]);
                         j++;
@@ -516,16 +513,16 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
             dd_sendrecv_rvec(dd, d, dddirForward,
                              sbuf, ind->nrecv[nzone+1],
                              buf,  ind->nsend[nzone+1]);
-            index = ind->index;
+            const int *index = ind->index;
             /* Add the received forces */
-            n = 0;
+            int        n = 0;
             if (!bShiftForcesNeedPbc)
             {
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i] + 1];
+                    for (int j = at0; j < at1; j++)
                     {
                         rvec_inc(f[j], buf[n]);
                         n++;
@@ -537,11 +534,11 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
                 /* fshift should always be defined if this function is
                  * called when bShiftForcesNeedPbc is true */
                 assert(NULL != fshift);
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i] + 1];
+                    for (int j = at0; j < at1; j++)
                     {
                         rvec_inc(f[j], buf[n]);
                         /* Add this force to the shift force */
@@ -552,11 +549,11 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
             }
             else
             {
-                for (i = 0; i < ind->nsend[nzone]; i++)
+                for (int i = 0; i < ind->nsend[nzone]; i++)
                 {
-                    at0 = cgindex[index[i]];
-                    at1 = cgindex[index[i]+1];
-                    for (j = at0; j < at1; j++)
+                    int at0 = atomGroups.index[index[i]];
+                    int at1 = atomGroups.index[index[i] + 1];
+                    for (int j = at0; j < at1; j++)
                     {
                         /* Rotate the force */
                         f[j][XX] += buf[n][XX];
@@ -579,34 +576,32 @@ void dd_move_f(gmx_domdec_t *dd, rvec f[], rvec *fshift, gmx_wallcycle *wcycle)
 
 void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
 {
-    int                    nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
-    int                   *index, *cgindex;
+    int                    nzone, nat_tot;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_comm_dim_t *cd;
-    gmx_domdec_ind_t      *ind;
     real                  *buf, *rbuf;
 
     comm = dd->comm;
 
-    cgindex = dd->cgindex;
+    const BlockRanges &atomGroups = dd->atomGroups();
 
     buf = &comm->vbuf.v[0][0];
 
     nzone   = 1;
     nat_tot = dd->nat_home;
-    for (d = 0; d < dd->ndim; d++)
+    for (int d = 0; d < dd->ndim; d++)
     {
         cd = &comm->cd[d];
-        for (p = 0; p < cd->np; p++)
+        for (int p = 0; p < cd->np; p++)
         {
-            ind   = &cd->ind[p];
-            index = ind->index;
-            n     = 0;
-            for (i = 0; i < ind->nsend[nzone]; i++)
+            const gmx_domdec_ind_t *ind   = &cd->ind[p];
+            const int              *index = ind->index;
+            int                     n     = 0;
+            for (int i = 0; i < ind->nsend[nzone]; i++)
             {
-                at0 = cgindex[index[i]];
-                at1 = cgindex[index[i]+1];
-                for (j = at0; j < at1; j++)
+                int at0 = atomGroups.index[index[i]];
+                int at1 = atomGroups.index[index[i] + 1];
+                for (int j = at0; j < at1; j++)
                 {
                     buf[n] = v[j];
                     n++;
@@ -627,10 +622,10 @@ void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
                              rbuf, ind->nrecv[nzone+1]);
             if (!cd->bInPlace)
             {
-                j = 0;
-                for (zone = 0; zone < nzone; zone++)
+                int j = 0;
+                for (int zone = 0; zone < nzone; zone++)
                 {
-                    for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
+                    for (int i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
                     {
                         v[i] = rbuf[j];
                         j++;
@@ -645,39 +640,37 @@ void dd_atom_spread_real(gmx_domdec_t *dd, real v[])
 
 void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
 {
-    int                    nzone, nat_tot, n, d, p, i, j, at0, at1, zone;
-    int                   *index, *cgindex;
+    int                    nzone, nat_tot;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_comm_dim_t *cd;
-    gmx_domdec_ind_t      *ind;
     real                  *buf, *sbuf;
 
     comm = dd->comm;
 
-    cgindex = dd->cgindex;
+    const gmx::BlockRanges &atomGroups = dd->atomGroups();
 
     buf = &comm->vbuf.v[0][0];
 
     nzone   = comm->zones.n/2;
-    nat_tot = dd->nat_tot;
-    for (d = dd->ndim-1; d >= 0; d--)
+    nat_tot = comm->nat[ddnatZONE];
+    for (int d = dd->ndim-1; d >= 0; d--)
     {
         cd = &comm->cd[d];
-        for (p = cd->np-1; p >= 0; p--)
+        for (int p = cd->np-1; p >= 0; p--)
         {
-            ind      = &cd->ind[p];
-            nat_tot -= ind->nrecv[nzone+1];
+            const gmx_domdec_ind_t *ind  = &cd->ind[p];
+            nat_tot                     -= ind->nrecv[nzone+1];
             if (cd->bInPlace)
             {
                 sbuf = v + nat_tot;
             }
             else
             {
-                sbuf = &comm->vbuf2.v[0][0];
-                j    = 0;
-                for (zone = 0; zone < nzone; zone++)
+                sbuf  = &comm->vbuf2.v[0][0];
+                int j = 0;
+                for (int zone = 0; zone < nzone; zone++)
                 {
-                    for (i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
+                    for (int i = ind->cell2at0[zone]; i < ind->cell2at1[zone]; i++)
                     {
                         sbuf[j] = v[i];
                         j++;
@@ -688,14 +681,14 @@ void dd_atom_sum_real(gmx_domdec_t *dd, real v[])
             dd_sendrecv_real(dd, d, dddirForward,
                              sbuf, ind->nrecv[nzone+1],
                              buf,  ind->nsend[nzone+1]);
-            index = ind->index;
+            const int *index = ind->index;
             /* Add the received forces */
-            n = 0;
-            for (i = 0; i < ind->nsend[nzone]; i++)
+            int        n = 0;
+            for (int i = 0; i < ind->nsend[nzone]; i++)
             {
-                at0 = cgindex[index[i]];
-                at1 = cgindex[index[i]+1];
-                for (j = at0; j < at1; j++)
+                int at0 = atomGroups.index[index[i]];
+                int at1 = atomGroups.index[index[i] + 1];
+                for (int j = at0; j < at1; j++)
                 {
                     v[j] += buf[n];
                     n++;
@@ -1112,9 +1105,8 @@ void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
 {
     char          fname[STRLEN], buf[22];
     FILE         *out;
-    int           i, ii, resnr, c;
+    int           resnr;
     const char   *atomname, *resname;
-    real          b;
     gmx_domdec_t *dd;
 
     dd = cr->dd;
@@ -1130,14 +1122,16 @@ void write_dd_pdb(const char *fn, gmx_int64_t step, const char *title,
     fprintf(out, "TITLE     %s\n", title);
     gmx_write_pdb_box(out, dd->bScrewPBC ? epbcSCREW : epbcXYZ, box);
     int molb = 0;
-    for (i = 0; i < natoms; i++)
+    for (int i = 0; i < natoms; i++)
     {
-        ii = dd->gatindex[i];
+        int  ii = dd->globalAtomIndices[i];
         mtopGetAtomAndResidueName(mtop, ii, &molb, &atomname, &resnr, &resname, nullptr);
+        int  c;
+        real b;
         if (i < dd->comm->nat[ddnatZONE])
         {
             c = 0;
-            while (i >= dd->cgindex[dd->comm->zones.cg_range[c+1]])
+            while (i >= dd->atomGroups().index[dd->comm->zones.cg_range[c+1]])
             {
                 c++;
             }
@@ -1499,33 +1493,37 @@ static void set_zones_ncg_home(gmx_domdec_t *dd)
     dd->comm->zone_ncg1[0] = dd->ncg_home;
 }
 
-static void rebuild_cgindex(gmx_domdec_t *dd,
-                            const int *gcgs_index, const t_state *state)
+static void restoreAtomGroups(gmx_domdec_t *dd,
+                              const int *gcgs_index, const t_state *state)
 {
-    int * gmx_restrict dd_cg_gl = dd->index_gl;
-    int * gmx_restrict cgindex  = dd->cgindex;
-    int                nat      = 0;
+    gmx::ArrayRef<const int>  atomGroupsState        = state->cg_gl;
+
+    std::vector<int>         &globalAtomGroupIndices = dd->globalAtomGroupIndices;
+    gmx::BlockRanges         &atomGroups             = dd->atomGroups_;
+
+    globalAtomGroupIndices.resize(atomGroupsState.size());
+    atomGroups.index.resize(atomGroupsState.size() + 1);
 
     /* Copy back the global charge group indices from state
      * and rebuild the local charge group to atom index.
      */
-    cgindex[0] = nat;
-    for (unsigned int i = 0; i < state->cg_gl.size(); i++)
+    int atomIndex = 0;
+    for (unsigned int i = 0; i < atomGroupsState.size(); i++)
     {
-        cgindex[i]  = nat;
-        int cg_gl   = state->cg_gl[i];
-        dd_cg_gl[i] = cg_gl;
-        nat        += gcgs_index[cg_gl+1] - gcgs_index[cg_gl];
+        atomGroups.index[i]        = atomIndex;
+        const int atomGroupGlobal  = atomGroupsState[i];
+        globalAtomGroupIndices[i]  = atomGroupGlobal;
+        atomIndex                 += gcgs_index[atomGroupGlobal + 1] - gcgs_index[atomGroupGlobal];
     }
-    cgindex[state->cg_gl.size()] = nat;
+    atomGroups.index[atomGroupsState.size()] = atomIndex;
 
-    dd->ncg_home = state->cg_gl.size();
-    dd->nat_home = nat;
+    dd->ncg_home = atomGroupsState.size();
+    dd->nat_home = atomIndex;
 
     set_zones_ncg_home(dd);
 }
 
-static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
+static void dd_set_cginfo(gmx::ArrayRef<const int> index_gl, int cg0, int cg1,
                           t_forcerec *fr, char *bLocalCG)
 {
     cginfo_mb_t *cginfo_mb;
@@ -1555,22 +1553,13 @@ static void dd_set_cginfo(int *index_gl, int cg0, int cg1,
 static void make_dd_indices(gmx_domdec_t *dd,
                             const int *gcgs_index, int cg_start)
 {
-    int          nzone, zone, zone1, cg0, cg1, cg1_p1, cg, cg_gl, a, a_gl;
-    int         *zone2cg, *zone_ncg1, *index_gl, *gatindex;
-    gmx_bool     bCGs;
-
-    if (dd->nat_tot > dd->gatindex_nalloc)
-    {
-        dd->gatindex_nalloc = over_alloc_dd(dd->nat_tot);
-        srenew(dd->gatindex, dd->gatindex_nalloc);
-    }
+    const int                 numZones               = dd->comm->zones.n;
+    const int                *zone2cg                = dd->comm->zones.cg_range;
+    const int                *zone_ncg1              = dd->comm->zone_ncg1;
+    gmx::ArrayRef<const int>  globalAtomGroupIndices = dd->globalAtomGroupIndices;
+    const gmx_bool            bCGs                   = dd->comm->bCGs;
 
-    nzone      = dd->comm->zones.n;
-    zone2cg    = dd->comm->zones.cg_range;
-    zone_ncg1  = dd->comm->zone_ncg1;
-    index_gl   = dd->index_gl;
-    gatindex   = dd->gatindex;
-    bCGs       = dd->comm->bCGs;
+    std::vector<int>         &globalAtomIndices      = dd->globalAtomIndices;
 
     if (zone2cg[1] != dd->ncg_home)
     {
@@ -1578,9 +1567,11 @@ static void make_dd_indices(gmx_domdec_t *dd,
     }
 
     /* Make the local to global and global to local atom index */
-    a = dd->cgindex[cg_start];
-    for (zone = 0; zone < nzone; zone++)
+    int a = dd->atomGroups().index[cg_start];
+    globalAtomIndices.resize(a);
+    for (int zone = 0; zone < numZones; zone++)
     {
+        int cg0;
         if (zone == 0)
         {
             cg0 = cg_start;
@@ -1589,30 +1580,30 @@ static void make_dd_indices(gmx_domdec_t *dd,
         {
             cg0 = zone2cg[zone];
         }
-        cg1    = zone2cg[zone+1];
-        cg1_p1 = cg0 + zone_ncg1[zone];
+        int cg1    = zone2cg[zone+1];
+        int cg1_p1 = cg0 + zone_ncg1[zone];
 
-        for (cg = cg0; cg < cg1; cg++)
+        for (int cg = cg0; cg < cg1; cg++)
         {
-            zone1 = zone;
+            int zone1 = zone;
             if (cg >= cg1_p1)
             {
                 /* Signal that this cg is from more than one pulse away */
-                zone1 += nzone;
+                zone1 += numZones;
             }
-            cg_gl = index_gl[cg];
+            int cg_gl = globalAtomGroupIndices[cg];
             if (bCGs)
             {
-                for (a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
+                for (int a_gl = gcgs_index[cg_gl]; a_gl < gcgs_index[cg_gl+1]; a_gl++)
                 {
-                    gatindex[a] = a_gl;
+                    globalAtomIndices.push_back(a_gl);
                     ga2la_set(dd->ga2la, a_gl, a, zone1);
                     a++;
                 }
             }
             else
             {
-                gatindex[a] = cg_gl;
+                globalAtomIndices.push_back(cg_gl);
                 ga2la_set(dd->ga2la, cg_gl, a, zone1);
                 a++;
             }
@@ -1623,33 +1614,31 @@ static void make_dd_indices(gmx_domdec_t *dd,
 static int check_bLocalCG(gmx_domdec_t *dd, int ncg_sys, const char *bLocalCG,
                           const char *where)
 {
-    int i, ngl, nerr;
-
-    nerr = 0;
+    int nerr = 0;
     if (bLocalCG == nullptr)
     {
         return nerr;
     }
-    for (i = 0; i < dd->ncg_tot; i++)
+    for (size_t i = 0; i < dd->globalAtomGroupIndices.size(); i++)
     {
-        if (!bLocalCG[dd->index_gl[i]])
+        if (!bLocalCG[dd->globalAtomGroupIndices[i]])
         {
             fprintf(stderr,
-                    "DD rank %d, %s: cg %d, global cg %d is not marked in bLocalCG (ncg_home %d)\n", dd->rank, where, i+1, dd->index_gl[i]+1, dd->ncg_home);
+                    "DD rank %d, %s: atom group %zu, global atom group %d is not marked in bLocalCG (ncg_home %d)\n", dd->rank, where, i + 1, dd->globalAtomGroupIndices[i] + 1, dd->ncg_home);
             nerr++;
         }
     }
-    ngl = 0;
-    for (i = 0; i < ncg_sys; i++)
+    size_t ngl = 0;
+    for (int i = 0; i < ncg_sys; i++)
     {
         if (bLocalCG[i])
         {
             ngl++;
         }
     }
-    if (ngl != dd->ncg_tot)
+    if (ngl != dd->globalAtomGroupIndices.size())
     {
-        fprintf(stderr, "DD rank %d, %s: In bLocalCG %d cgs are marked as local, whereas there are %d\n", dd->rank, where, ngl, dd->ncg_tot);
+        fprintf(stderr, "DD rank %d, %s: In bLocalCG %zu atom groups are marked as local, whereas there are %zu\n", dd->rank, where, ngl, dd->globalAtomGroupIndices.size());
         nerr++;
     }
 
@@ -1660,102 +1649,102 @@ static void check_index_consistency(gmx_domdec_t *dd,
                                     int natoms_sys, int ncg_sys,
                                     const char *where)
 {
-    int   nerr, ngl, i, a, cell;
-    int  *have;
+    int       nerr = 0;
 
-    nerr = 0;
+    const int numAtomsInZones = dd->comm->nat[ddnatZONE];
 
     if (dd->comm->DD_debug > 1)
     {
-        snew(have, natoms_sys);
-        for (a = 0; a < dd->nat_tot; a++)
+        std::vector<int> have(natoms_sys);
+        for (int a = 0; a < numAtomsInZones; a++)
         {
-            if (have[dd->gatindex[a]] > 0)
+            int globalAtomIndex = dd->globalAtomIndices[a];
+            if (have[globalAtomIndex] > 0)
             {
-                fprintf(stderr, "DD rank %d: global atom %d occurs twice: index %d and %d\n", dd->rank, dd->gatindex[a]+1, have[dd->gatindex[a]], a+1);
+                fprintf(stderr, "DD rank %d: global atom %d occurs twice: index %d and %d\n", dd->rank, globalAtomIndex + 1, have[globalAtomIndex], a+1);
             }
             else
             {
-                have[dd->gatindex[a]] = a + 1;
+                have[globalAtomIndex] = a + 1;
             }
         }
-        sfree(have);
     }
 
-    snew(have, dd->nat_tot);
+    std::vector<int> have(numAtomsInZones);
 
-    ngl  = 0;
-    for (i = 0; i < natoms_sys; i++)
+    int              ngl = 0;
+    for (int i = 0; i < natoms_sys; i++)
     {
+        int a;
+        int cell;
         if (ga2la_get(dd->ga2la, i, &a, &cell))
         {
-            if (a >= dd->nat_tot)
+            if (a >= numAtomsInZones)
             {
-                fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which is larger than nat_tot (%d)\n", dd->rank, i+1, a+1, dd->nat_tot);
+                fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which is larger than nat_tot (%d)\n", dd->rank, i+1, a+1, numAtomsInZones);
                 nerr++;
             }
             else
             {
                 have[a] = 1;
-                if (dd->gatindex[a] != i)
+                if (dd->globalAtomIndices[a] != i)
                 {
-                    fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which has global atom index %d\n", dd->rank, i+1, a+1, dd->gatindex[a]+1);
+                    fprintf(stderr, "DD rank %d: global atom %d marked as local atom %d, which has global atom index %d\n", dd->rank, i+1, a+1, dd->globalAtomIndices[a]+1);
                     nerr++;
                 }
             }
             ngl++;
         }
     }
-    if (ngl != dd->nat_tot)
+    if (ngl != numAtomsInZones)
     {
         fprintf(stderr,
                 "DD rank %d, %s: %d global atom indices, %d local atoms\n",
-                dd->rank, where, ngl, dd->nat_tot);
+                dd->rank, where, ngl, numAtomsInZones);
     }
-    for (a = 0; a < dd->nat_tot; a++)
+    for (int a = 0; a < numAtomsInZones; a++)
     {
         if (have[a] == 0)
         {
             fprintf(stderr,
                     "DD rank %d, %s: local atom %d, global %d has no global index\n",
-                    dd->rank, where, a+1, dd->gatindex[a]+1);
+                    dd->rank, where, a + 1, dd->globalAtomIndices[a] + 1);
         }
     }
-    sfree(have);
 
     nerr += check_bLocalCG(dd, ncg_sys, dd->comm->bLocalCG, where);
 
     if (nerr > 0)
     {
-        gmx_fatal(FARGS, "DD rank %d, %s: %d atom/cg index inconsistencies",
+        gmx_fatal(FARGS, "DD rank %d, %s: %d atom(group) index inconsistencies",
                   dd->rank, where, nerr);
     }
 }
 
-static void clear_dd_indices(gmx_domdec_t *dd, int cg_start, int a_start)
+/* Clear all DD global state indices, starting from \p atomGroupStart and \p atomStart */
+static void clearDDStateIndices(gmx_domdec_t *dd,
+                                int           atomGroupStart,
+                                int           atomStart)
 {
-    int   i;
-    char *bLocalCG;
-
-    if (a_start == 0)
+    if (atomStart == 0)
     {
         /* Clear the whole list without searching */
         ga2la_clear(dd->ga2la);
     }
     else
     {
-        for (i = a_start; i < dd->nat_tot; i++)
+        for (int i = 0; i < dd->comm->nat[ddnatZONE]; i++)
         {
-            ga2la_del(dd->ga2la, dd->gatindex[i]);
+            ga2la_del(dd->ga2la, dd->globalAtomIndices[i]);
         }
     }
 
-    bLocalCG = dd->comm->bLocalCG;
+    char *bLocalCG = dd->comm->bLocalCG;
     if (bLocalCG)
     {
-        for (i = cg_start; i < dd->ncg_tot; i++)
+        for (size_t atomGroup = atomGroupStart; atomGroup < dd->globalAtomGroupIndices.size(); atomGroup++)
         {
-            bLocalCG[dd->index_gl[i]] = FALSE;
+            bLocalCG[dd->globalAtomGroupIndices[atomGroup]] = FALSE;
         }
     }
 
@@ -3202,9 +3191,9 @@ static void make_dd_communicators(FILE *fplog, t_commrec *cr,
 
     if (DDMASTER(dd))
     {
-        dd->ma = new AtomDistribution(dd->nc,
-                                      comm->cgs_gl.nr,
-                                      comm->cgs_gl.index[comm->cgs_gl.nr]);
+        dd->ma = gmx::compat::make_unique<AtomDistribution>(dd->nc,
+                                                            comm->cgs_gl.nr,
+                                                            comm->cgs_gl.index[comm->cgs_gl.nr]);
     }
 }
 
@@ -3459,8 +3448,6 @@ static int determineInitialDlbState(FILE *fplog, t_commrec *cr,
 
 static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
 {
-    int dim;
-
     dd->ndim = 0;
     if (getenv("GMX_DD_ORDER_ZYX") != nullptr)
     {
@@ -3469,7 +3456,7 @@ static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
         {
             fprintf(fplog, "Using domain decomposition order z, y, x\n");
         }
-        for (dim = DIM-1; dim >= 0; dim--)
+        for (int dim = DIM-1; dim >= 0; dim--)
         {
             if (dd->nc[dim] > 1)
             {
@@ -3480,7 +3467,7 @@ static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
     else
     {
         /* Decomposition order x,y,z */
-        for (dim = 0; dim < DIM; dim++)
+        for (int dim = 0; dim < DIM; dim++)
         {
             if (dd->nc[dim] > 1)
             {
@@ -3488,6 +3475,12 @@ static void set_dd_dim(FILE *fplog, gmx_domdec_t *dd)
             }
         }
     }
+
+    if (dd->ndim == 0)
+    {
+        /* Set dim[0] to avoid extra checks on ndim in several places */
+        dd->dim[0] = XX;
+    }
 }
 
 static gmx_domdec_comm_t *init_dd_comm()
@@ -4437,10 +4430,14 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
                 "\nInitializing Domain Decomposition on %d ranks\n", cr->nnodes);
     }
 
-    snew(dd, 1);
+    dd = new gmx_domdec_t;
 
     dd->comm = init_dd_comm();
 
+    /* Initialize DD paritioning counters */
+    dd->comm->partition_step = INT_MIN;
+    dd->ddp_count            = 0;
+
     set_dd_envvar_options(fplog, dd, cr->nodeid);
 
     gmx_ddbox_t ddbox = {0};
@@ -4461,10 +4458,6 @@ gmx_domdec_t *init_domain_decomposition(FILE *fplog, t_commrec *cr,
     /* Set overallocation to avoid frequent reallocation of arrays */
     set_over_alloc_dd(TRUE);
 
-    /* Initialize DD paritioning counters */
-    dd->comm->partition_step = INT_MIN;
-    dd->ddp_count            = 0;
-
     /* We currently don't know the number of threads yet, we set this later */
     dd->comm->nth = 0;
 
@@ -4680,9 +4673,10 @@ void dd_dlb_unlock(gmx_domdec_t *dd)
 static void merge_cg_buffers(int ncell,
                              gmx_domdec_comm_dim_t *cd, int pulse,
                              int  *ncg_cell,
-                             int  *index_gl, int  *recv_i,
+                             gmx::ArrayRef<int> index_gl,
+                             int  *recv_i,
                              rvec *cg_cm,    rvec *recv_vr,
-                             int *cgindex,
+                             gmx::ArrayRef<int> cgindex,
                              cginfo_mb_t *cginfo_mb, int *cginfo)
 {
     gmx_domdec_ind_t *ind, *ind_p;
@@ -4762,7 +4756,9 @@ static void merge_cg_buffers(int ncell,
 }
 
 static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
-                               int nzone, int cg0, const int *cgindex)
+                               int                    nzone,
+                               int                    cg0,
+                               const BlockRanges     &atomGroups)
 {
     int cg, zone, p;
 
@@ -4773,9 +4769,9 @@ static void make_cell2at_index(gmx_domdec_comm_dim_t *cd,
     {
         for (p = 0; p < cd->np; p++)
         {
-            cd->ind[p].cell2at0[zone] = cgindex[cg];
+            cd->ind[p].cell2at0[zone] = atomGroups.index[cg];
             cg += cd->ind[p].nrecv[zone];
-            cd->ind[p].cell2at1[zone] = cgindex[cg];
+            cd->ind[p].cell2at1[zone] = atomGroups.index[cg];
         }
     }
 }
@@ -4910,13 +4906,13 @@ static void
 get_zone_pulse_cgs(gmx_domdec_t *dd,
                    int zonei, int zone,
                    int cg0, int cg1,
-                   const int *index_gl,
-                   const int *cgindex,
+                   gmx::ArrayRef<const int> globalAtomGroupIndices,
+                   const gmx::BlockRanges &atomGroups,
                    int dim, int dim_ind,
                    int dim0, int dim1, int dim2,
                    real r_comm2, real r_bcomm2,
                    matrix box,
-                   ivec tric_dist,
+                   bool distanceIsTriclinic,
                    rvec *normal,
                    real skew_fac2_d, real skew_fac_01,
                    rvec *v_d, rvec *v_0, rvec *v_1,
@@ -4958,7 +4954,7 @@ get_zone_pulse_cgs(gmx_domdec_t *dd,
     {
         r2  = 0;
         rb2 = 0;
-        if (tric_dist[dim_ind] == 0)
+        if (!distanceIsTriclinic)
         {
             /* Rectangular direction, easy */
             r = cg_cm[cg][dim] - c->c[dim_ind][zone];
@@ -5129,7 +5125,7 @@ get_zone_pulse_cgs(gmx_domdec_t *dd,
               (bDist2B && r2  < r_bcomm2)) &&
              (!bBondComm ||
               (GET_CGINFO_BOND_INTER(cginfo[cg]) &&
-               missing_link(comm->cglink, index_gl[cg],
+               missing_link(comm->cglink, globalAtomGroupIndices[cg],
                             comm->bLocalCG)))))
         {
             /* Make an index to the local charge groups */
@@ -5144,7 +5140,7 @@ get_zone_pulse_cgs(gmx_domdec_t *dd,
                 srenew(*ibuf, *ibuf_nalloc);
             }
             ind->index[nsend] = cg;
-            (*ibuf)[nsend]    = index_gl[cg];
+            (*ibuf)[nsend]    = globalAtomGroupIndices[cg];
             nsend_z++;
             vec_rvec_check_alloc(vbuf, nsend+1);
 
@@ -5163,7 +5159,7 @@ get_zone_pulse_cgs(gmx_domdec_t *dd,
                 copy_rvec(cg_cm[cg], vbuf->v[nsend]);
             }
             nsend++;
-            nat += cgindex[cg+1] - cgindex[cg];
+            nat += atomGroups.index[cg+1] - atomGroups.index[cg];
         }
     }
 
@@ -5180,7 +5176,7 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     int                    dim_ind, dim, dim0, dim1, dim2, dimd, p, nat_tot;
     int                    nzone, nzone_send, zone, zonei, cg0, cg1;
     int                    c, i, cg, cg_gl, nrcg;
-    int                   *zone_cg_range, pos_cg, *index_gl, *cgindex, *recv_i;
+    int                   *zone_cg_range, pos_cg;
     gmx_domdec_comm_t     *comm;
     gmx_domdec_zones_t    *zones;
     gmx_domdec_comm_dim_t *cd;
@@ -5189,7 +5185,6 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     gmx_bool               bBondComm, bDist2B, bDistMB, bDistBonded;
     real                   r_comm2, r_bcomm2;
     dd_corners_t           corners;
-    ivec                   tric_dist;
     rvec                  *cg_cm, *normal, *v_d, *v_0 = nullptr, *v_1 = nullptr, *recv_vr;
     real                   skew_fac2_d, skew_fac_01;
     rvec                   sf2_round;
@@ -5229,19 +5224,6 @@ static void setup_dd_communication(gmx_domdec_t *dd,
             cg_cm = nullptr;
     }
 
-    for (dim_ind = 0; dim_ind < dd->ndim; dim_ind++)
-    {
-        /* Check if we need to use triclinic distances */
-        tric_dist[dim_ind] = 0;
-        for (i = 0; i <= dim_ind; i++)
-        {
-            if (ddbox->tric_dir[dd->dim[i]])
-            {
-                tric_dist[dim_ind] = 1;
-            }
-        }
-    }
-
     bBondComm = comm->bBondComm;
 
     /* Do we need to determine extra distances for multi-body bondeds? */
@@ -5292,8 +5274,6 @@ static void setup_dd_communication(gmx_domdec_t *dd,
     }
 
     zone_cg_range = zones->cg_range;
-    index_gl      = dd->index_gl;
-    cgindex       = dd->cgindex;
     cginfo_mb     = fr->cginfo_mb;
 
     zone_cg_range[0]   = 0;
@@ -5308,6 +5288,16 @@ static void setup_dd_communication(gmx_domdec_t *dd,
         dim = dd->dim[dim_ind];
         cd  = &comm->cd[dim_ind];
 
+        /* Check if we need to compute triclinic distances along this dim */
+        bool distanceIsTriclinic = false;
+        for (i = 0; i <= dim_ind; i++)
+        {
+            if (ddbox->tric_dir[dim])
+            {
+                distanceIsTriclinic = true;
+            }
+        }
+
         if (dim >= ddbox->npbcdim && dd->ci[dim] == 0)
         {
             /* No pbc in this dimension, the first node should not comm. */
@@ -5334,7 +5324,7 @@ static void setup_dd_communication(gmx_domdec_t *dd,
             nat   = 0;
             for (zone = 0; zone < nzone_send; zone++)
             {
-                if (tric_dist[dim_ind] && dim_ind > 0)
+                if (dim_ind > 0 && distanceIsTriclinic)
                 {
                     /* Determine slightly more optimized skew_fac's
                      * for rounding.
@@ -5433,10 +5423,11 @@ static void setup_dd_communication(gmx_domdec_t *dd,
 
                         /* Get the cg's for this pulse in this zone */
                         get_zone_pulse_cgs(dd, zonei, zone, cg0_th, cg1_th,
-                                           index_gl, cgindex,
+                                           dd->globalAtomGroupIndices,
+                                           dd->atomGroups(),
                                            dim, dim_ind, dim0, dim1, dim2,
                                            r_comm2, r_bcomm2,
-                                           box, tric_dist,
+                                           box, distanceIsTriclinic,
                                            normal, skew_fac2_d, skew_fac_01,
                                            v_d, v_0, v_1, &corners, sf2_round,
                                            bDistBonded, bBondComm,
@@ -5532,17 +5523,14 @@ static void setup_dd_communication(gmx_domdec_t *dd,
             }
 
             /* Make space for the global cg indices */
-            if (pos_cg + ind->nrecv[nzone] > dd->cg_nalloc
-                || dd->cg_nalloc == 0)
-            {
-                dd->cg_nalloc = over_alloc_dd(pos_cg + ind->nrecv[nzone]);
-                srenew(index_gl, dd->cg_nalloc);
-                srenew(cgindex, dd->cg_nalloc+1);
-            }
+            int numAtomGroupsNew = pos_cg + ind->nrecv[nzone];
+            dd->globalAtomGroupIndices.resize(numAtomGroupsNew);
+            dd->atomGroups_.index.resize(numAtomGroupsNew + 1);
             /* Communicate the global cg indices */
+            int *recv_i;
             if (cd->bInPlace)
             {
-                recv_i = index_gl + pos_cg;
+                recv_i =  dd->globalAtomGroupIndices.data() + pos_cg;
             }
             else
             {
@@ -5583,10 +5571,10 @@ static void setup_dd_communication(gmx_domdec_t *dd,
                 {
                     for (cg = 0; cg < ind->nrecv[zone]; cg++)
                     {
-                        cg_gl              = index_gl[pos_cg];
-                        fr->cginfo[pos_cg] = ddcginfo(cginfo_mb, cg_gl);
-                        nrcg               = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
-                        cgindex[pos_cg+1]  = cgindex[pos_cg] + nrcg;
+                        cg_gl                              = dd->globalAtomGroupIndices[pos_cg];
+                        fr->cginfo[pos_cg]                 = ddcginfo(cginfo_mb, cg_gl);
+                        nrcg                               = GET_CGINFO_NATOMS(fr->cginfo[pos_cg]);
+                        dd->atomGroups_.index[pos_cg + 1]  = dd->atomGroups_.index[pos_cg] + nrcg;
                         if (bBondComm)
                         {
                             /* Update the charge group presence,
@@ -5608,8 +5596,9 @@ static void setup_dd_communication(gmx_domdec_t *dd,
             {
                 /* This part of the code is never executed with bBondComm. */
                 merge_cg_buffers(nzone, cd, p, zone_cg_range,
-                                 index_gl, recv_i, cg_cm, recv_vr,
-                                 cgindex, fr->cginfo_mb, fr->cginfo);
+                                 dd->globalAtomGroupIndices, recv_i, cg_cm, recv_vr,
+                                 dd->atomGroups_.index,
+                                 fr->cginfo_mb, fr->cginfo);
                 pos_cg += ind->nrecv[nzone];
             }
             nat_tot += ind->nrecv[nzone+1];
@@ -5617,19 +5606,15 @@ static void setup_dd_communication(gmx_domdec_t *dd,
         if (!cd->bInPlace)
         {
             /* Store the atom block for easy copying of communication buffers */
-            make_cell2at_index(cd, nzone, zone_cg_range[nzone], cgindex);
+            make_cell2at_index(cd, nzone, zone_cg_range[nzone], dd->atomGroups());
         }
         nzone += nzone;
     }
-    dd->index_gl = index_gl;
-    dd->cgindex  = cgindex;
 
-    dd->ncg_tot          = zone_cg_range[zones->n];
-    dd->nat_tot          = nat_tot;
     comm->nat[ddnatHOME] = dd->nat_home;
     for (i = ddnatZONE; i < ddnatNR; i++)
     {
-        comm->nat[i] = dd->nat_tot;
+        comm->nat[i] = nat_tot;
     }
 
     if (!bBondComm)
@@ -5637,7 +5622,8 @@ static void setup_dd_communication(gmx_domdec_t *dd,
         /* We don't need to update cginfo, since that was alrady done above.
          * So we pass NULL for the forcerec.
          */
-        dd_set_cginfo(dd->index_gl, dd->ncg_home, dd->ncg_tot,
+        dd_set_cginfo(dd->globalAtomGroupIndices,
+                      dd->ncg_home, dd->globalAtomGroupIndices.size(),
                       nullptr, comm->bLocalCG);
     }
 
@@ -5976,12 +5962,13 @@ static void order_vec_cg(int n, const gmx_cgsort_t *sort,
     }
 }
 
-static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort,
+static void order_vec_atom(int ncg, const gmx::BlockRanges *atomGroups,
+                           const gmx_cgsort_t *sort,
                            rvec *v, rvec *buf)
 {
     int a, atot, cg, cg0, cg1, i;
 
-    if (cgindex == nullptr)
+    if (atomGroups == nullptr)
     {
         /* Avoid the useless loop of the atoms within a cg */
         order_vec_cg(ncg, sort, v, buf);
@@ -5993,8 +5980,8 @@ static void order_vec_atom(int ncg, const int *cgindex, const gmx_cgsort_t *sort
     a = 0;
     for (cg = 0; cg < ncg; cg++)
     {
-        cg0 = cgindex[sort[cg].ind];
-        cg1 = cgindex[sort[cg].ind+1];
+        cg0 = atomGroups->index[sort[cg].ind];
+        cg1 = atomGroups->index[sort[cg].ind+1];
         for (i = cg0; i < cg1; i++)
         {
             copy_rvec(v[i], buf[a]);
@@ -6093,7 +6080,7 @@ static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
                  * but we set it here anyhow to avoid a conditional.
                  */
                 sort_i->nsc    = a[i];
-                sort_i->ind_gl = dd->index_gl[i];
+                sort_i->ind_gl = dd->globalAtomGroupIndices[i];
                 sort_i->ind    = i;
                 ncg_new++;
             }
@@ -6117,7 +6104,7 @@ static int dd_sort_order(gmx_domdec_t *dd, t_forcerec *fr, int ncg_home_old)
              * and the global topology index
              */
             cgsort[i].nsc    = a[i];
-            cgsort[i].ind_gl = dd->index_gl[i];
+            cgsort[i].ind_gl = dd->globalAtomGroupIndices[i];
             cgsort[i].ind    = i;
             if (cgsort[i].nsc < moved)
             {
@@ -6163,7 +6150,6 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
 {
     gmx_domdec_sort_t *sort;
     gmx_cgsort_t      *cgsort;
-    int               *cgindex;
     int                ncg_new, i, *ibuf, cgsize;
     rvec              *vbuf;
 
@@ -6190,18 +6176,13 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
             ncg_new = 0;
     }
 
+    const gmx::BlockRanges &atomGroups = dd->atomGroups();
+
     /* We alloc with the old size, since cgindex is still old */
-    vec_rvec_check_alloc(&dd->comm->vbuf, dd->cgindex[dd->ncg_home]);
+    vec_rvec_check_alloc(&dd->comm->vbuf, atomGroups.index[dd->ncg_home]);
     vbuf = dd->comm->vbuf.v;
 
-    if (dd->comm->bCGs)
-    {
-        cgindex = dd->cgindex;
-    }
-    else
-    {
-        cgindex = nullptr;
-    }
+    const gmx::BlockRanges *atomGroupsPtr = (dd->comm->bCGs ? &atomGroups : nullptr);
 
     /* Remove the charge groups which are no longer at home here */
     dd->ncg_home = ncg_new;
@@ -6214,15 +6195,15 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
     /* Reorder the state */
     if (state->flags & (1 << estX))
     {
-        order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->x.data()), vbuf);
+        order_vec_atom(dd->ncg_home, atomGroupsPtr, cgsort, as_rvec_array(state->x.data()), vbuf);
     }
     if (state->flags & (1 << estV))
     {
-        order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->v.data()), vbuf);
+        order_vec_atom(dd->ncg_home, atomGroupsPtr, cgsort, as_rvec_array(state->v.data()), vbuf);
     }
     if (state->flags & (1 << estCGP))
     {
-        order_vec_atom(dd->ncg_home, cgindex, cgsort, as_rvec_array(state->cg_p.data()), vbuf);
+        order_vec_atom(dd->ncg_home, atomGroupsPtr, cgsort, as_rvec_array(state->cg_p.data()), vbuf);
     }
 
     if (fr->cutoff_scheme == ecutsGROUP)
@@ -6238,7 +6219,7 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
     }
     ibuf = sort->ibuf;
     /* Reorder the global cg index */
-    order_int_cg(dd->ncg_home, cgsort, dd->index_gl, ibuf);
+    order_int_cg(dd->ncg_home, cgsort, dd->globalAtomGroupIndices.data(), ibuf);
     /* Reorder the cginfo */
     order_int_cg(dd->ncg_home, cgsort, fr->cginfo, ibuf);
     /* Rebuild the local cg index */
@@ -6247,23 +6228,23 @@ static void dd_sort_state(gmx_domdec_t *dd, rvec *cgcm, t_forcerec *fr, t_state
         ibuf[0] = 0;
         for (i = 0; i < dd->ncg_home; i++)
         {
-            cgsize    = dd->cgindex[cgsort[i].ind+1] - dd->cgindex[cgsort[i].ind];
+            cgsize    = atomGroups.index[cgsort[i].ind+1] - atomGroups.index[cgsort[i].ind];
             ibuf[i+1] = ibuf[i] + cgsize;
         }
         for (i = 0; i < dd->ncg_home+1; i++)
         {
-            dd->cgindex[i] = ibuf[i];
+            dd->atomGroups_.index[i] = ibuf[i];
         }
     }
     else
     {
         for (i = 0; i < dd->ncg_home+1; i++)
         {
-            dd->cgindex[i] = i;
+            dd->atomGroups_.index[i] = i;
         }
     }
     /* Set the home atom number */
-    dd->nat_home = dd->cgindex[dd->ncg_home];
+    dd->nat_home = dd->atomGroups().index[dd->ncg_home];
 
     if (fr->cutoff_scheme == ecutsVERLET)
     {
@@ -6611,7 +6592,7 @@ void dd_partition_system(FILE                *fplog,
     if (bMasterState)
     {
         /* Clear the old state */
-        clear_dd_indices(dd, 0, 0);
+        clearDDStateIndices(dd, 0, 0);
         ncgindex_set = 0;
 
         rvec *xGlobal = (DDMASTER(dd) ? as_rvec_array(state_global->x.data()) : nullptr);
@@ -6636,7 +6617,7 @@ void dd_partition_system(FILE                *fplog,
 
         inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
 
-        dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
+        dd_set_cginfo(dd->globalAtomGroupIndices, 0, dd->ncg_home, fr, comm->bLocalCG);
     }
     else if (state_local->ddp_count != dd->ddp_count)
     {
@@ -6651,10 +6632,10 @@ void dd_partition_system(FILE                *fplog,
         }
 
         /* Clear the old state */
-        clear_dd_indices(dd, 0, 0);
+        clearDDStateIndices(dd, 0, 0);
 
-        /* Build the new indices */
-        rebuild_cgindex(dd, cgs_gl->index, state_local);
+        /* Restore the atom group indices from state_local */
+        restoreAtomGroups(dd, cgs_gl->index, state_local);
         make_dd_indices(dd, cgs_gl->index, 0);
         ncgindex_set = dd->ncg_home;
 
@@ -6667,7 +6648,7 @@ void dd_partition_system(FILE                *fplog,
 
         inc_nrnb(nrnb, eNR_CGCM, dd->nat_home);
 
-        dd_set_cginfo(dd->index_gl, 0, dd->ncg_home, fr, comm->bLocalCG);
+        dd_set_cginfo(dd->globalAtomGroupIndices, 0, dd->ncg_home, fr, comm->bLocalCG);
 
         set_ddbox(dd, bMasterState, ir, state_local->box,
                   TRUE, &top_local->cgs, as_rvec_array(state_local->x.data()), &ddbox);
@@ -6679,7 +6660,7 @@ void dd_partition_system(FILE                *fplog,
         /* We have the full state, only redistribute the cgs */
 
         /* Clear the non-home indices */
-        clear_dd_indices(dd, dd->ncg_home, dd->nat_home);
+        clearDDStateIndices(dd, dd->ncg_home, dd->nat_home);
         ncgindex_set = 0;
 
         /* Avoid global communication for dim's without pbc and -gcom */
@@ -6920,7 +6901,7 @@ void dd_partition_system(FILE                *fplog,
         {
             if (EEL_FULL(ir->coulombtype) && dd->n_intercg_excl > 0)
             {
-                nat_f_novirsum = dd->nat_tot;
+                nat_f_novirsum = comm->nat[ddnatZONE];
             }
             else
             {
@@ -6939,8 +6920,8 @@ void dd_partition_system(FILE                *fplog,
      * allocation, zeroing and copying, but this is probably not worth
      * the complications and checking.
      */
-    forcerec_set_ranges(fr, dd->ncg_home, dd->ncg_tot,
-                        dd->nat_tot, comm->nat[ddnatCON], nat_f_novirsum);
+    forcerec_set_ranges(fr, dd->ncg_home, dd->globalAtomGroupIndices.size(),
+                        comm->nat[ddnatZONE], comm->nat[ddnatCON], nat_f_novirsum);
 
     /* Update atom data for mdatoms and several algorithms */
     mdAlgorithmsSetupAtomData(cr, ir, top_global, top_local, fr,
index 0636bef10f32998ef55be5b53007050f148708e4..0fe94bfc951f5da24216edb4e871bbca89d8edb3 100644 (file)
@@ -256,9 +256,9 @@ static void atoms_to_settles(gmx_domdec_t *dd,
     {
         if (GET_CGINFO_SETTLE(cginfo[cg]))
         {
-            for (int a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
+            for (int a = dd->atomGroups().index[cg]; a < dd->atomGroups().index[cg+1]; a++)
             {
-                int a_gl = dd->gatindex[a];
+                int a_gl = dd->globalAtomIndices[a];
                 int a_mol;
                 mtopGetMolblockIndex(mtop, a_gl, &mb, nullptr, &a_mol);
 
@@ -350,9 +350,9 @@ static void atoms_to_constraints(gmx_domdec_t *dd,
     {
         if (GET_CGINFO_CONSTR(cginfo[cg]))
         {
-            for (int a = dd->cgindex[cg]; a < dd->cgindex[cg+1]; a++)
+            for (int a = dd->atomGroups().index[cg]; a < dd->atomGroups().index[cg+1]; a++)
             {
-                int a_gl = dd->gatindex[a];
+                int a_gl = dd->globalAtomIndices[a];
                 int molnr, a_mol;
                 mtopGetMolblockIndex(mtop, a_gl, &mb, &molnr, &a_mol);
 
index 0b91385c258cf0bb1c05390fa12ae8f915e4d7f5..510ca68483b1985ef3ffb593f8c653cb8588433a 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2006,2007,2008,2009,2010,2012,2013,2014,2015,2017, by the GROMACS development team, led by
+ * Copyright (c) 2006,2007,2008,2009,2010,2012,2013,2014,2015,2017,2018, 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.
@@ -531,16 +531,11 @@ int setup_specat_communication(gmx_domdec_t               *dd,
             }
             nrecv_local += buf[0];
             spas->nrecv  = buf[1];
-            if (nat_tot_specat + spas->nrecv > dd->gatindex_nalloc)
-            {
-                dd->gatindex_nalloc =
-                    over_alloc_dd(nat_tot_specat + spas->nrecv);
-                srenew(dd->gatindex, dd->gatindex_nalloc);
-            }
+            dd->globalAtomIndices.resize(nat_tot_specat + spas->nrecv);
             /* Send and receive the indices */
             dd_sendrecv_int(dd, d, dir == 0 ? dddirBackward : dddirForward,
                             spac->ibuf, spas->nsend,
-                            dd->gatindex+nat_tot_specat, spas->nrecv);
+                            dd->globalAtomIndices.data() + nat_tot_specat, spas->nrecv);
             nat_tot_specat += spas->nrecv;
         }
 
@@ -566,7 +561,7 @@ int setup_specat_communication(gmx_domdec_t               *dd,
         /* Make a global to local index for the communication atoms */
         for (i = nat_tot_prev; i < nat_tot_specat; i++)
         {
-            gmx_hash_change_or_set(ga2la_specat, dd->gatindex[i], i);
+            gmx_hash_change_or_set(ga2la_specat, dd->globalAtomIndices[i], i);
         }
     }
 
index de89abf52ee61ec6a73e17338b560a170d45cd97..31898035575cfded68589521f952fdb7a6967eb6 100644 (file)
 
 #include <cstddef>
 
+#include <memory>
+
 #include "gromacs/math/vectypes.h"
+#include "gromacs/topology/block.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/gmxmpi.h"
 #include "gromacs/utility/real.h"
@@ -134,8 +137,8 @@ struct gmx_domdec_t {
     /* Communication with the PME only nodes */
     int                    pme_nodeid;
     gmx_bool               pme_receive_vir_ener;
-    gmx_pme_comm_n_box_t  *cnb;
-    int                    nreq_pme;
+    gmx_pme_comm_n_box_t  *cnb      = nullptr;
+    int                    nreq_pme = 0;
     MPI_Request            req_pme[8];
 
 
@@ -154,7 +157,7 @@ struct gmx_domdec_t {
     int  neighbor[DIM][2];
 
     /* Only available on the master node */
-    AtomDistribution *ma;
+    std::unique_ptr<AtomDistribution> ma;
 
     /* Are there inter charge group constraints */
     gmx_bool bInterCGcons;
@@ -169,33 +172,33 @@ struct gmx_domdec_t {
     int  n_intercg_excl;
 
     /* Vsite stuff */
-    gmx_hash_t                *ga2la_vsite;
-    gmx_domdec_specat_comm_t  *vsite_comm;
+    gmx_hash_t                *ga2la_vsite = nullptr;
+    gmx_domdec_specat_comm_t  *vsite_comm  = nullptr;
 
     /* Constraint stuff */
-    gmx_domdec_constraints_t *constraints;
-    gmx_domdec_specat_comm_t *constraint_comm;
-
-    /* The local to gobal charge group index and local cg to local atom index */
-    int   ncg_home;
-    int   ncg_tot;
-    int  *index_gl;
-    int  *cgindex;
-    int   cg_nalloc;
-    /* Local atom to local cg index, only for special cases */
-    int  *la2lc;
-    int   la2lc_nalloc;
+    gmx_domdec_constraints_t *constraints     = nullptr;
+    gmx_domdec_specat_comm_t *constraint_comm = nullptr;
+
+    /* The number of home atom groups */
+    int                     ncg_home = 0;
+    /* Global atom group indices for the home and all non-home groups */
+    std::vector<int>        globalAtomGroupIndices;
+    /* The atom groups for the home and all non-home groups, todo: make private */
+    gmx::BlockRanges        atomGroups_;
+    const gmx::BlockRanges &atomGroups() const
+    {
+        return atomGroups_;
+    }
+    /* Local atom to local atom-group index, only used for checking bondeds */
+    std::vector<int> localAtomGroupFromAtom;
 
     /* The number of home atoms */
-    int   nat_home;
-    /* The total number of atoms: home and received zones */
-    int   nat_tot;
-    /* Index from the local atoms to the global atoms */
-    int  *gatindex;
-    int   gatindex_nalloc;
+    int              nat_home = 0;
+    /* Index from the local atoms to the global atoms, covers home and received zones */
+    std::vector<int> globalAtomIndices;
 
     /* Global atom number to local atom number list */
-    gmx_ga2la_t  *ga2la;
+    gmx_ga2la_t  *ga2la = nullptr;
 
     /* Communication stuff */
     gmx_domdec_comm_t *comm;
@@ -205,9 +208,8 @@ struct gmx_domdec_t {
 
 
     /* gmx_pme_recv_f buffer */
-    int   pme_recv_f_alloc;
-    rvec *pme_recv_f_buf;
-
+    int   pme_recv_f_alloc = 0;
+    rvec *pme_recv_f_buf   = nullptr;
 };
 
 #endif
index e6fe90fc8cae6213d916a566dc76920301bf8f42..129b87b784c9482ed0debba2c145d861ea245fbe 100644 (file)
@@ -190,7 +190,7 @@ static void print_missing_interactions_mb(FILE *fplog, t_commrec *cr,
     int  nril_mol = ril->index[nat_mol];
     snew(assigned, nmol*nril_mol);
 
-    int *gatindex = cr->dd->gatindex;
+    gmx::ArrayRef<const int> gatindex = cr->dd->globalAtomIndices;
     for (int ftype = 0; ftype < F_NRE; ftype++)
     {
         if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle))
@@ -1093,26 +1093,18 @@ static void add_vsite(gmx_ga2la_t *ga2la, const int *index, const int *rtil,
     }
 }
 
-/*! \brief Update the local atom to local charge group index */
-static void make_la2lc(gmx_domdec_t *dd)
+/*! \brief Build the index that maps each local atom to its local atom group */
+static void makeLocalAtomGroupFromAtom(gmx_domdec_t *dd)
 {
-    int *cgindex, *la2lc, cg, a;
+    const gmx::BlockRanges &atomGroups = dd->atomGroups();
 
-    cgindex = dd->cgindex;
+    dd->localAtomGroupFromAtom.clear();
 
-    if (dd->nat_tot > dd->la2lc_nalloc)
+    for (size_t g = 0; g < dd->globalAtomGroupIndices.size(); g++)
     {
-        dd->la2lc_nalloc = over_alloc_dd(dd->nat_tot);
-        snew(dd->la2lc, dd->la2lc_nalloc);
-    }
-    la2lc = dd->la2lc;
-
-    /* Make the local atom to local cg index */
-    for (cg = 0; cg < dd->ncg_tot; cg++)
-    {
-        for (a = cgindex[cg]; a < cgindex[cg+1]; a++)
+        for (int a = atomGroups.index[g]; a < atomGroups.index[g + 1]; a++)
         {
-            la2lc[a] = cg;
+            dd->localAtomGroupFromAtom.push_back(g);
         }
     }
 }
@@ -1538,7 +1530,7 @@ static int make_bondeds_zone(gmx_domdec_t *dd,
                              int izone,
                              int at_start, int at_end)
 {
-    int                i, i_gl, mb, mt, mol, i_mol;
+    int                mb, mt, mol, i_mol;
     int               *index, *rtil;
     gmx_bool           bBCheck;
     gmx_reverse_top_t *rt;
@@ -1550,10 +1542,10 @@ static int make_bondeds_zone(gmx_domdec_t *dd,
 
     nbonded_local = 0;
 
-    for (i = at_start; i < at_end; i++)
+    for (int i = at_start; i < at_end; i++)
     {
         /* Get the global atom number */
-        i_gl = dd->gatindex[i];
+        const int i_gl = dd->globalAtomIndices[i];
         global_atomnr_to_moltype_ind(rt, i_gl, &mb, &mt, &mol, &i_mol);
         /* Check all intramolecular interactions assigned to this atom */
         index = rt->ril_mt[mt].index;
@@ -1607,12 +1599,10 @@ static int make_bondeds_zone(gmx_domdec_t *dd,
 static void set_no_exclusions_zone(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                    int iz, t_blocka *lexcls)
 {
-    int  a0, a1, a;
+    int a0 = dd->atomGroups().index[zones->cg_range[iz]];
+    int a1 = dd->atomGroups().index[zones->cg_range[iz + 1]];
 
-    a0 = dd->cgindex[zones->cg_range[iz]];
-    a1 = dd->cgindex[zones->cg_range[iz+1]];
-
-    for (a = a0+1; a < a1+1; a++)
+    for (int a = a0 + 1; a < a1 + 1; a++)
     {
         lexcls->index[a] = lexcls->nra;
     }
@@ -1633,48 +1623,49 @@ static int make_exclusions_zone_cg(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                    int iz,
                                    int cg_start, int cg_end)
 {
-    int               n_excl_at_max, n, count, jla0, jla1, jla;
-    int               cg, la0, la1, la, a_gl, mb, mt, mol, a_mol, j, aj_mol;
+    int               n_excl_at_max;
+    int               mb, mt, mol;
     const t_blocka   *excls;
     gmx_ga2la_t      *ga2la;
     int               cell;
 
     ga2la = dd->ga2la;
 
-    jla0 = dd->cgindex[zones->izone[iz].jcg0];
-    jla1 = dd->cgindex[zones->izone[iz].jcg1];
+    int jla0 = dd->atomGroups().index[zones->izone[iz].jcg0];
+    int jla1 = dd->atomGroups().index[zones->izone[iz].jcg1];
 
     n_excl_at_max = dd->reverse_top->n_excl_at_max;
 
     /* We set the end index, but note that we might not start at zero here */
-    lexcls->nr = dd->cgindex[cg_end];
+    lexcls->nr = dd->atomGroups().index[cg_end];
 
-    n     = lexcls->nra;
-    count = 0;
-    for (cg = cg_start; cg < cg_end; cg++)
+    int n     = lexcls->nra;
+    int count = 0;
+    for (int cg = cg_start; cg < cg_end; cg++)
     {
         if (n + (cg_end - cg_start)*n_excl_at_max > lexcls->nalloc_a)
         {
             lexcls->nalloc_a = over_alloc_large(n + (cg_end - cg_start)*n_excl_at_max);
             srenew(lexcls->a, lexcls->nalloc_a);
         }
-        la0 = dd->cgindex[cg];
-        la1 = dd->cgindex[cg+1];
+        int la0 = dd->atomGroups().index[cg];
+        int la1 = dd->atomGroups().index[cg + 1];
         if (GET_CGINFO_EXCL_INTER(cginfo[cg]) ||
             !GET_CGINFO_EXCL_INTRA(cginfo[cg]))
         {
             /* Copy the exclusions from the global top */
-            for (la = la0; la < la1; la++)
+            for (int la = la0; la < la1; la++)
             {
                 lexcls->index[la] = n;
-                a_gl              = dd->gatindex[la];
+                int a_gl          = dd->globalAtomIndices[la];
+                int a_mol;
                 global_atomnr_to_moltype_ind(dd->reverse_top, a_gl, &mb, &mt, &mol, &a_mol);
                 excls = &moltype[mt].excls;
-                for (j = excls->index[a_mol]; j < excls->index[a_mol+1]; j++)
+                for (int j = excls->index[a_mol]; j < excls->index[a_mol+1]; j++)
                 {
-                    aj_mol = excls->a[j];
+                    int aj_mol = excls->a[j];
                     /* This computation of jla is only correct intra-cg */
-                    jla = la + aj_mol - a_mol;
+                    int jla = la + aj_mol - a_mol;
                     if (jla >= la0 && jla < la1)
                     {
                         /* This is an intra-cg exclusion. We can skip
@@ -1734,10 +1725,10 @@ static int make_exclusions_zone_cg(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
              */
             if (iz == 0)
             {
-                for (la = la0; la < la1; la++)
+                for (int la = la0; la < la1; la++)
                 {
                     lexcls->index[la] = n;
-                    for (j = la0; j < la1; j++)
+                    for (int j = la0; j < la1; j++)
                     {
                         lexcls->a[n++] = j;
                     }
@@ -1747,7 +1738,7 @@ static int make_exclusions_zone_cg(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
             else
             {
                 /* We don't need exclusions for this cg */
-                for (la = la0; la < la1; la++)
+                for (int la = la0; la < la1; la++)
                 {
                     lexcls->index[la] = n;
                 }
@@ -1771,13 +1762,12 @@ static void make_exclusions_zone(gmx_domdec_t *dd,
                                  int at_start, int at_end)
 {
     gmx_ga2la_t *ga2la;
-    int          jla0, jla1;
     int          n_excl_at_max, n, at;
 
     ga2la = dd->ga2la;
 
-    jla0 = dd->cgindex[zones->izone[iz].jcg0];
-    jla1 = dd->cgindex[zones->izone[iz].jcg1];
+    int jla0 = dd->atomGroups().index[zones->izone[iz].jcg0];
+    int jla1 = dd->atomGroups().index[zones->izone[iz].jcg1];
 
     n_excl_at_max = dd->reverse_top->n_excl_at_max;
 
@@ -1805,7 +1795,7 @@ static void make_exclusions_zone(gmx_domdec_t *dd,
 
             /* Copy the exclusions from the global top */
             lexcls->index[at] = n;
-            a_gl              = dd->gatindex[at];
+            a_gl              = dd->globalAtomIndices[at];
             global_atomnr_to_moltype_ind(dd->reverse_top, a_gl,
                                          &mb, &mt, &mol, &a_mol);
             excls = &moltype[mt].excls;
@@ -1854,14 +1844,11 @@ static void check_alloc_index(t_blocka *ba, int nindex_max)
 static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                    t_blocka *lexcls)
 {
-    int nr;
-    int thread;
-
-    nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
+    int nr = dd->atomGroups().index[zones->izone[zones->nizone - 1].cg1];
 
     check_alloc_index(lexcls, nr);
 
-    for (thread = 1; thread < dd->reverse_top->nthread; thread++)
+    for (int thread = 1; thread < dd->reverse_top->nthread; thread++)
     {
         check_alloc_index(&dd->reverse_top->th_work[thread].excl, nr);
     }
@@ -1871,17 +1858,15 @@ static void check_exclusions_alloc(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
 static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
                                     t_blocka *lexcls)
 {
-    int la0, la;
-
-    lexcls->nr = dd->cgindex[zones->izone[zones->nizone-1].cg1];
+    lexcls->nr = dd->atomGroups().index[zones->izone[zones->nizone - 1].cg1];
 
     if (dd->n_intercg_excl == 0)
     {
         /* There are no exclusions involving non-home charge groups,
          * but we need to set the indices for neighborsearching.
          */
-        la0 = dd->cgindex[zones->izone[0].cg1];
-        for (la = la0; la < lexcls->nr; la++)
+        int la0 = dd->atomGroups().index[zones->izone[0].cg1];
+        for (int la = la0; la < lexcls->nr; la++)
         {
             lexcls->index[la] = lexcls->nra;
         }
@@ -1889,7 +1874,7 @@ static void finish_local_exclusions(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
         /* nr is only used to loop over the exclusions for Ewald and RF,
          * so we can set it to the number of home atoms for efficiency.
          */
-        lexcls->nr = dd->cgindex[zones->izone[0].cg1];
+        lexcls->nr = dd->atomGroups().index[zones->izone[0].cg1];
     }
 }
 
@@ -2016,7 +2001,8 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd,
                                       idef_t,
                                       vsite_pbc, vsite_pbc_nalloc,
                                       izone,
-                                      dd->cgindex[cg0t], dd->cgindex[cg1t]);
+                                      dd->atomGroups().index[cg0t],
+                                      dd->atomGroups().index[cg1t]);
 
                 if (izone < nzone_excl)
                 {
@@ -2031,7 +2017,8 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd,
                         excl_t->nra = 0;
                     }
 
-                    if (dd->cgindex[dd->ncg_tot] == dd->ncg_tot &&
+                    int numAtomGroups = dd->globalAtomGroupIndices.size();
+                    if (dd->atomGroups().index[numAtomGroups] == numAtomGroups &&
                         !rt->bExclRequired)
                     {
                         /* No charge groups and no distance check required */
@@ -2100,8 +2087,8 @@ static int make_local_bondeds_excls(gmx_domdec_t *dd,
 
 void dd_make_local_cgs(gmx_domdec_t *dd, t_block *lcgs)
 {
-    lcgs->nr    = dd->ncg_tot;
-    lcgs->index = dd->cgindex;
+    lcgs->nr    = dd->globalAtomGroupIndices.size();
+    lcgs->index = dd->atomGroups_.index.data();
 }
 
 void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
@@ -2167,7 +2154,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
         }
         if (bRCheckMB || bRCheck2B)
         {
-            make_la2lc(dd);
+            makeLocalAtomGroupFromAtom(dd);
             if (fr->bMolPBC)
             {
                 pbc_null = set_pbc_dd(&pbc, fr->ePBC, dd->nc, TRUE, box);
@@ -2182,7 +2169,7 @@ void dd_make_local_top(gmx_domdec_t *dd, gmx_domdec_zones_t *zones,
     dd->nbonded_local =
         make_local_bondeds_excls(dd, zones, mtop, fr->cginfo,
                                  bRCheckMB, rcheck, bRCheck2B, rc,
-                                 dd->la2lc,
+                                 dd->localAtomGroupFromAtom.data(),
                                  pbc_null, cgcm_or_x,
                                  &ltop->idef, vsite,
                                  &ltop->excls, &nexcl);
index df8bae94f6716b364483930e6f8d1158fa18e0d8..98f3382415b066b40edd624f9309f47e21655c75 100644 (file)
@@ -71,7 +71,7 @@
 #define DD_FLAG_BW(d) (1<<(16+(d)*2+1))
 
 static int compact_and_copy_vec_at(int ncg, int *move,
-                                   int *cgindex,
+                                   const gmx::BlockRanges &atomGroups,
                                    int nvec, int vec,
                                    rvec *src, gmx_domdec_comm_t *comm,
                                    gmx_bool bCompact)
@@ -90,7 +90,7 @@ static int compact_and_copy_vec_at(int ncg, int *move,
     i0 = 0;
     for (icg = 0; icg < ncg; icg++)
     {
-        i1 = cgindex[icg+1];
+        i1 = atomGroups.index[icg+1];
         m  = move[icg];
         if (m == -1)
         {
@@ -125,7 +125,7 @@ static int compact_and_copy_vec_at(int ncg, int *move,
 }
 
 static int compact_and_copy_vec_cg(int ncg, int *move,
-                                   int *cgindex,
+                                   const gmx::BlockRanges &atomGroups,
                                    int nvec, rvec *src, gmx_domdec_comm_t *comm,
                                    gmx_bool bCompact)
 {
@@ -143,7 +143,7 @@ static int compact_and_copy_vec_cg(int ncg, int *move,
     i0 = 0;
     for (icg = 0; icg < ncg; icg++)
     {
-        i1 = cgindex[icg+1];
+        i1 = atomGroups.index[icg + 1];
         m  = move[icg];
         if (m == -1)
         {
@@ -170,79 +170,81 @@ static int compact_and_copy_vec_cg(int ncg, int *move,
     return home_pos;
 }
 
-static int compact_ind(int ncg, int *move,
-                       int *index_gl, int *cgindex,
-                       int *gatindex,
-                       gmx_ga2la_t *ga2la, char *bLocalCG,
-                       int *cginfo)
+static int compact_ind(int                 ncg,
+                       const int          *move,
+                       gmx::ArrayRef<int>  globalAtomGroupIndices,
+                       gmx::BlockRanges   *atomGroups,
+                       gmx::ArrayRef<int>  globalAtomIndices,
+                       gmx_ga2la_t        *ga2la,
+                       char               *bLocalCG,
+                       int                *cginfo)
 {
-    int cg, nat, a0, a1, a, a_gl;
-    int home_pos;
-
-    home_pos = 0;
-    nat      = 0;
-    for (cg = 0; cg < ncg; cg++)
+    int home_pos = 0;
+    int nat      = 0;
+    for (int cg = 0; cg < ncg; cg++)
     {
-        a0 = cgindex[cg];
-        a1 = cgindex[cg+1];
+        int a0 = atomGroups->index[cg];
+        int a1 = atomGroups->index[cg+1];
         if (move[cg] == -1)
         {
             /* Compact the home arrays in place.
              * Anything that can be done here avoids access to global arrays.
              */
-            cgindex[home_pos] = nat;
-            for (a = a0; a < a1; a++)
+            atomGroups->index[home_pos] = nat;
+            for (int a = a0; a < a1; a++)
             {
-                a_gl          = gatindex[a];
-                gatindex[nat] = a_gl;
+                const int a_gl         = globalAtomIndices[a];
+                globalAtomIndices[nat] = a_gl;
                 /* The cell number stays 0, so we don't need to set it */
                 ga2la_change_la(ga2la, a_gl, nat);
                 nat++;
             }
-            index_gl[home_pos] = index_gl[cg];
-            cginfo[home_pos]   = cginfo[cg];
+            globalAtomGroupIndices[home_pos] = globalAtomGroupIndices[cg];
+            cginfo[home_pos]                 = cginfo[cg];
             /* The charge group remains local, so bLocalCG does not change */
             home_pos++;
         }
         else
         {
             /* Clear the global indices */
-            for (a = a0; a < a1; a++)
+            for (int a = a0; a < a1; a++)
             {
-                ga2la_del(ga2la, gatindex[a]);
+                ga2la_del(ga2la, globalAtomIndices[a]);
             }
             if (bLocalCG)
             {
-                bLocalCG[index_gl[cg]] = FALSE;
+                bLocalCG[globalAtomGroupIndices[cg]] = FALSE;
             }
         }
     }
-    cgindex[home_pos] = nat;
+    atomGroups->index[home_pos] = nat;
 
     return home_pos;
 }
 
-static void clear_and_mark_ind(int ncg, int *move,
-                               int *index_gl, int *cgindex, int *gatindex,
-                               gmx_ga2la_t *ga2la, char *bLocalCG,
-                               int *cell_index)
+static void clear_and_mark_ind(int                       ncg,
+                               const int                *move,
+                               gmx::ArrayRef<const int>  globalAtomGroupIndices,
+                               const gmx::BlockRanges   &atomGroups,
+                               gmx::ArrayRef<const int>  globalAtomIndices,
+                               gmx_ga2la_t              *ga2la,
+                               char                     *bLocalCG,
+                               int                      *cell_index)
 {
-    int cg, a0, a1, a;
-
-    for (cg = 0; cg < ncg; cg++)
+    for (int cg = 0; cg < ncg; cg++)
     {
         if (move[cg] >= 0)
         {
-            a0 = cgindex[cg];
-            a1 = cgindex[cg+1];
+            int a0 = atomGroups.index[cg];
+            int a1 = atomGroups.index[cg + 1];
             /* Clear the global indices */
-            for (a = a0; a < a1; a++)
+            for (int a = a0; a < a1; a++)
             {
-                ga2la_del(ga2la, gatindex[a]);
+                ga2la_del(ga2la, globalAtomIndices[a]);
             }
             if (bLocalCG)
             {
-                bLocalCG[index_gl[cg]] = FALSE;
+                bLocalCG[globalAtomGroupIndices[cg]] = FALSE;
             }
             /* Signal that this cg has moved using the ns cell index.
              * Here we set it to -1. fill_grid will change it
@@ -269,14 +271,14 @@ static void print_cg_move(FILE *fplog,
     {
         fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition (%f) in direction %c\n",
                 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
-                ddglatnr(dd, dd->cgindex[cg]), limitd, dim2char(dim));
+                ddglatnr(dd, dd->atomGroups().index[cg]), limitd, dim2char(dim));
     }
     else
     {
         /* We don't have a limiting distance available: don't print it */
         fprintf(fplog, "%s %d moved more than the distance allowed by the domain decomposition in direction %c\n",
                 dd->comm->bCGs ? "The charge group starting at atom" : "Atom",
-                ddglatnr(dd, dd->cgindex[cg]), dim2char(dim));
+                ddglatnr(dd, dd->atomGroups().index[cg]), dim2char(dim));
     }
     fprintf(fplog, "distance out of cell %f\n",
             dir == 1 ? pos_d - comm->cell_x1[dim] : pos_d - comm->cell_x0[dim]);
@@ -352,7 +354,7 @@ static void calc_cg_move(FILE *fplog, gmx_int64_t step,
                          ivec tric_dir, matrix tcm,
                          rvec cell_x0, rvec cell_x1,
                          rvec limitd, rvec limit0, rvec limit1,
-                         const int *cgindex,
+                         const gmx::BlockRanges &atomGroups,
                          int cg_start, int cg_end,
                          rvec *cg_cm,
                          int *move)
@@ -370,8 +372,8 @@ static void calc_cg_move(FILE *fplog, gmx_int64_t step,
 
     for (cg = cg_start; cg < cg_end; cg++)
     {
-        k0   = cgindex[cg];
-        k1   = cgindex[cg+1];
+        k0   = atomGroups.index[cg];
+        k1   = atomGroups.index[cg + 1];
         nrcg = k1 - k0;
         if (nrcg == 1)
         {
@@ -543,7 +545,6 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
     real               pos_d;
     matrix             tcm;
     rvec              *cg_cm = nullptr, cell_x0, cell_x1, limitd, limit0, limit1;
-    int               *cgindex;
     cginfo_mb_t       *cginfo_mb;
     gmx_domdec_comm_t *comm;
     int               *moved;
@@ -564,9 +565,9 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
     bool bV   = state->flags & (1<<estV);
     bool bCGP = state->flags & (1<<estCGP);
 
-    if (dd->ncg_tot > comm->nalloc_int)
+    if (dd->globalAtomGroupIndices.size() > static_cast<size_t>(comm->nalloc_int))
     {
-        comm->nalloc_int = over_alloc_dd(dd->ncg_tot);
+        comm->nalloc_int = over_alloc_dd(dd->globalAtomGroupIndices.size());
         srenew(comm->buf_int, comm->nalloc_int);
     }
     move = comm->buf_int;
@@ -609,7 +610,7 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
 
     make_tric_corr_matrix(npbcdim, state->box, tcm);
 
-    cgindex = dd->cgindex;
+    const gmx::BlockRanges &atomGroups = dd->atomGroups();
 
     nthread = gmx_omp_nthreads_get(emntDomdec);
 
@@ -623,7 +624,7 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
         {
             calc_cg_move(fplog, step, dd, state, tric_dir, tcm,
                          cell_x0, cell_x1, limitd, limit0, limit1,
-                         cgindex,
+                         atomGroups,
                          ( thread   *dd->ncg_home)/nthread,
                          ((thread+1)*dd->ncg_home)/nthread,
                          fr->cutoff_scheme == ecutsGROUP ? cg_cm : as_rvec_array(state->x.data()),
@@ -646,12 +647,12 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
                 comm->cggl_flag_nalloc[mc] = over_alloc_dd(ncg[mc]+1);
                 srenew(comm->cggl_flag[mc], comm->cggl_flag_nalloc[mc]*DD_CGIBS);
             }
-            comm->cggl_flag[mc][ncg[mc]*DD_CGIBS  ] = dd->index_gl[cg];
+            comm->cggl_flag[mc][ncg[mc]*DD_CGIBS  ] = dd->globalAtomGroupIndices[cg];
             /* We store the cg size in the lower 16 bits
              * and the place where the charge group should go
              * in the next 6 bits. This saves some communication volume.
              */
-            nrcg = cgindex[cg+1] - cgindex[cg];
+            nrcg = atomGroups.index[cg+1] - atomGroups.index[cg];
             comm->cggl_flag[mc][ncg[mc]*DD_CGIBS+1] = nrcg | flag;
             ncg[mc] += 1;
             nat[mc] += nrcg;
@@ -695,7 +696,7 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
              * but that could give rise to rounding issues.
              */
             home_pos_cg =
-                compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
+                compact_and_copy_vec_cg(dd->ncg_home, move, dd->atomGroups(),
                                         nvec, cg_cm, comm, bCompact);
             break;
         case ecutsVERLET:
@@ -704,7 +705,7 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
              * many conditionals for both for with and without charge groups.
              */
             home_pos_cg =
-                compact_and_copy_vec_cg(dd->ncg_home, move, cgindex,
+                compact_and_copy_vec_cg(dd->ncg_home, move, dd->atomGroups(),
                                         nvec, as_rvec_array(state->x.data()), comm, FALSE);
             if (bCompact)
             {
@@ -718,18 +719,18 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
 
     vec         = 0;
     home_pos_at =
-        compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
+        compact_and_copy_vec_at(dd->ncg_home, move, dd->atomGroups(),
                                 nvec, vec++, as_rvec_array(state->x.data()),
                                 comm, bCompact);
     if (bV)
     {
-        compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
+        compact_and_copy_vec_at(dd->ncg_home, move, dd->atomGroups(),
                                 nvec, vec++, as_rvec_array(state->v.data()),
                                 comm, bCompact);
     }
     if (bCGP)
     {
-        compact_and_copy_vec_at(dd->ncg_home, move, cgindex,
+        compact_and_copy_vec_at(dd->ncg_home, move, dd->atomGroups(),
                                 nvec, vec++, as_rvec_array(state->cg_p.data()),
                                 comm, bCompact);
     }
@@ -737,7 +738,7 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
     if (bCompact)
     {
         compact_ind(dd->ncg_home, move,
-                    dd->index_gl, dd->cgindex, dd->gatindex,
+                    dd->globalAtomGroupIndices, &dd->atomGroups_, dd->globalAtomIndices,
                     dd->ga2la, comm->bLocalCG,
                     fr->cginfo);
     }
@@ -758,11 +759,15 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
         }
 
         clear_and_mark_ind(dd->ncg_home, move,
-                           dd->index_gl, dd->cgindex, dd->gatindex,
+                           dd->globalAtomGroupIndices, dd->atomGroups(), dd->globalAtomIndices,
                            dd->ga2la, comm->bLocalCG,
                            moved);
     }
 
+    /* Now we can remove the excess global atom-group indices from the list */
+    dd->globalAtomGroupIndices.resize(home_pos_cg);
+    dd->atomGroups_.index.resize(home_pos_cg + 1);
+
     cginfo_mb = fr->cginfo_mb;
 
     *ncg_stay_home = home_pos_cg;
@@ -876,6 +881,9 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
                                         comm->vbuf.v[buf_pos][d3]*tcm[d3][dim2];
                                 }
                             }
+
+                            GMX_ASSERT(dim2 >= 0 && dim2 < DIM, "Keep the static analyzer happy");
+
                             /* Check of we are not at the box edge.
                              * pbc is only handled in the first step above,
                              * but this check could move over pbc while
@@ -916,15 +924,10 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
             nrcg = flag & DD_FLAG_NRCG;
             if (mc == -1)
             {
-                if (home_pos_cg+1 > dd->cg_nalloc)
-                {
-                    dd->cg_nalloc = over_alloc_dd(home_pos_cg+1);
-                    srenew(dd->index_gl, dd->cg_nalloc);
-                    srenew(dd->cgindex, dd->cg_nalloc+1);
-                }
                 /* Set the global charge group index and size */
-                dd->index_gl[home_pos_cg]  = comm->buf_int[cg*DD_CGIBS];
-                dd->cgindex[home_pos_cg+1] = dd->cgindex[home_pos_cg] + nrcg;
+                const int globalAtomGroupIndex = comm->buf_int[cg*DD_CGIBS];
+                dd->globalAtomGroupIndices.push_back(globalAtomGroupIndex);
+                dd->atomGroups_.index.push_back(dd->atomGroups_.index[home_pos_cg] + nrcg);
                 /* Copy the state from the buffer */
                 if (fr->cutoff_scheme == ecutsGROUP)
                 {
@@ -935,10 +938,10 @@ void dd_redistribute_cg(FILE *fplog, gmx_int64_t step,
 
                 /* Set the cginfo */
                 fr->cginfo[home_pos_cg] = ddcginfo(cginfo_mb,
-                                                   dd->index_gl[home_pos_cg]);
+                                                   globalAtomGroupIndex);
                 if (comm->bLocalCG)
                 {
-                    comm->bLocalCG[dd->index_gl[home_pos_cg]] = TRUE;
+                    comm->bLocalCG[globalAtomGroupIndex] = TRUE;
                 }
 
                 for (i = 0; i < nrcg; i++)
index 9448b953b541c65158ab7dfd715dd9309626ed7c..ad2be05bbae50a133e44f5d14398dfdbbfa2c8a2 100644 (file)
@@ -260,7 +260,7 @@ static void write_constr_pdb(const char *fn, const char *title,
             {
                 continue;
             }
-            ii = dd->gatindex[i];
+            ii = dd->globalAtomIndices[i];
         }
         else
         {
index 1da83d6c6f7b17ac8232260b562e506e4cff16df..4a1cc52eb1971ecb36a2afaddae3210dd3c258f4 100644 (file)
@@ -784,7 +784,7 @@ void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt
 void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
                      const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
 {
-    const int                                 *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr);
+    const int                                 *gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
     int                                        i;
     int                                        gc = 0;
     gmx::ThreeFry2x64<0>                       rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
index 3e07b81f1311cd5a2692719bd3ac2d4bb91f89a0..f37a997a808d92ac063cdb7362c5eaa5706b33cd 100644 (file)
@@ -336,7 +336,7 @@ void do_force_lowlevel(t_forcerec           *fr,
                     idef, (const rvec *) x, hist,
                     forceForUseWithShiftForces, forceWithVirial,
                     fr, &pbc, graph, enerd, nrnb, lambda, md, fcd,
-                    DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr,
+                    DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr,
                     flags);
 
 
index b00123e0a6c8898fdedd0747666393d4b5169ee4..98c8a1e561f76dda3737e7c8107434d0ef148d06 100644 (file)
@@ -77,7 +77,7 @@ void mdAlgorithmsSetupAtomData(const t_commrec   *cr,
     if (usingDomDec)
     {
         numAtomIndex = dd_natoms_mdatoms(cr->dd);
-        atomIndex    = cr->dd->gatindex;
+        atomIndex    = cr->dd->globalAtomIndices.data();
         numHomeAtoms = cr->dd->nat_home;
     }
     else
index 5fe6e1994f89c8b253aac5af9b3bc42914b686b4..06220f2f32c1461ccf6f0617ecdae3dbf59c494b 100644 (file)
@@ -629,7 +629,7 @@ void make_local_shells(const t_commrec *cr,
             }
             if (dd)
             {
-                shell[nshell] = shfc->shell_gl[ind[dd->gatindex[i]]];
+                shell[nshell] = shfc->shell_gl[ind[dd->globalAtomIndices[i]]];
             }
             else
             {
index ae10681136c03fea092c3791e5589b42d8b9cc82..3a1aabb419a4139fad75322e5d0ac75cf75e6ad9 100644 (file)
@@ -1115,7 +1115,7 @@ static void do_force_cutsVERLET(FILE *fplog,
 
     if (DOMAINDECOMP(cr))
     {
-        cg1 = cr->dd->ncg_tot;
+        cg1 = cr->dd->globalAtomGroupIndices.size();
     }
     else
     {
@@ -1785,7 +1785,7 @@ static void do_force_cutsGROUP(FILE *fplog,
     cg0 = 0;
     if (DOMAINDECOMP(cr))
     {
-        cg1 = cr->dd->ncg_tot;
+        cg1 = cr->dd->globalAtomGroupIndices.size();
     }
     else
     {
index d6e4da7decc1ff9977c229724fba47129e7c3dc2..0634b5d5c017a2b96c4f6eb8178a3bac9eeac098 100644 (file)
@@ -1579,7 +1579,7 @@ update_sd_second_half(gmx_int64_t                    step,
                     as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()),
                     as_rvec_array(state->v.data()), nullptr,
                     step, inputrec->ld_seed,
-                    DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr);
+                    DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
             }
             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
         }
@@ -1880,7 +1880,7 @@ void update_coords(gmx_int64_t                    step,
                             md->cFREEZE, md->cACC, md->cTC,
                             x_rvec, xp_rvec, v_rvec, f_rvec,
                             step, inputrec->ld_seed,
-                            DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr);
+                            DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
                     }
                     break;
                 case (eiBD):
@@ -1890,7 +1890,7 @@ void update_coords(gmx_int64_t                    step,
                                  x_rvec, xp_rvec, v_rvec, f_rvec,
                                  inputrec->bd_fric,
                                  upd->sd->bd_rf,
-                                 step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->gatindex : nullptr);
+                                 step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
                     break;
                 case (eiVV):
                 case (eiVVAK):
index 2ef86fbc165329769ef13807d53de6c76de4761d..87aa001d26fde28229ca0315e6cc05296d9b242c 100644 (file)
@@ -75,7 +75,7 @@ static void write_nblist(FILE *out, gmx_domdec_t *dd, t_nblist *nblist, int nDNL
 
             for (zi = 0; zi < dd_zones->n; zi++)
             {
-                ca1[zi] = dd->cgindex[dd_zones->cg_range[zi+1]];
+                ca1[zi] = dd->atomGroups().index[dd_zones->cg_range[zi + 1]];
             }
             i = 0;
             for (zi = 0; zi < dd_zones->nizone && zi < dd_zones->n; zi++)
index 928358a6ace36e13139e8572501f01fe8a3610c3..23cb6bf22fabe857f7d7ab1b5461ac0758d2bfd8 100644 (file)
@@ -271,7 +271,7 @@ static void get_f_norm_max(const t_commrec *cr,
 
     if (la_max >= 0 && DOMAINDECOMP(cr))
     {
-        a_max = cr->dd->gatindex[la_max];
+        a_max = cr->dd->globalAtomIndices[la_max];
     }
     else
     {