Use ListOfLists in gmx_mtop_t and gmx_localtop_t
authorBerk Hess <hess@kth.se>
Fri, 6 Dec 2019 09:11:55 +0000 (10:11 +0100)
committerMark Abraham <mark.j.abraham@gmail.com>
Mon, 16 Dec 2019 15:37:47 +0000 (16:37 +0100)
This change uses ListOfLists for exclusions in gmx_mtop_t and
gmx_localtop_t. Exclusions in t_topolgy as well as other atom
and index lists still use t_blocka.

Change-Id: Idfd532ad2a4b19d68e60b05dbd0a3082650ec20a

33 files changed:
src/gromacs/domdec/domdec_topology.cpp
src/gromacs/fileio/tpxio.cpp
src/gromacs/gmxpreprocess/gpp_nextnb.cpp
src/gromacs/gmxpreprocess/gpp_nextnb.h
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/grompp_impl.h
src/gromacs/gmxpreprocess/topio.cpp
src/gromacs/gmxpreprocess/toppush.cpp
src/gromacs/mdlib/dispersioncorrection.cpp
src/gromacs/mdlib/forcerec.cpp
src/gromacs/mdlib/perf_est.cpp
src/gromacs/mdlib/sim_util.cpp
src/gromacs/mdrun/tpi.cpp
src/gromacs/nbnxm/benchmark/bench_setup.cpp
src/gromacs/nbnxm/benchmark/bench_system.cpp
src/gromacs/nbnxm/benchmark/bench_system.h
src/gromacs/nbnxm/nbnxm.h
src/gromacs/nbnxm/pairlist.cpp
src/gromacs/nbnxm/pairlistset.h
src/gromacs/nbnxm/pairlistsets.h
src/gromacs/selection/nbsearch.cpp
src/gromacs/selection/nbsearch.h
src/gromacs/selection/tests/nbsearch.cpp
src/gromacs/tools/convert_tpr.cpp
src/gromacs/topology/block.cpp
src/gromacs/topology/block.h
src/gromacs/topology/exclusionblocks.cpp
src/gromacs/topology/exclusionblocks.h
src/gromacs/topology/mtop_util.cpp
src/gromacs/topology/tests/exclusionblocks.cpp
src/gromacs/topology/topology.cpp
src/gromacs/topology/topology.h
src/gromacs/trajectoryanalysis/modules/rdf.cpp

index 19acbac64bbf6025c254d2c4fba57c9891312fb8..23346db87cdf6f220e4469510edd57250ae5c537 100644 (file)
@@ -75,6 +75,7 @@
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/logger.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/strconvert.h"
@@ -87,6 +88,8 @@
 #include "domdec_vsite.h"
 #include "dump.h"
 
+using gmx::ListOfLists;
+
 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
 #define NITEM_DD_INIT_LOCAL_STATE 5
 
@@ -111,7 +114,7 @@ struct thread_work_t
     t_idef                    idef;       /**< Partial local topology */
     std::unique_ptr<VsitePbc> vsitePbc;   /**< vsite PBC structure */
     int                       nbonded;    /**< The number of bondeds in this struct */
-    t_blocka                  excl;       /**< List of exclusions */
+    ListOfLists<int>          excl;       /**< List of exclusions */
     int                       excl_count; /**< The total exclusion count for \p excl */
 };
 
@@ -119,8 +122,6 @@ struct thread_work_t
 struct gmx_reverse_top_t
 {
     //! @cond Doxygen_Suppress
-    //! \brief The maximum number of exclusions one atom can have
-    int n_excl_at_max = 0;
     //! \brief Are there constraints in this revserse top?
     bool bConstr = false;
     //! \brief Are there settles in this revserse top?
@@ -453,14 +454,15 @@ static void global_atomnr_to_moltype_ind(const gmx_reverse_top_t* rt, int i_gl,
 }
 
 /*! \brief Returns the maximum number of exclusions per atom */
-static int getMaxNumExclusionsPerAtom(const t_blocka& excls)
+static int getMaxNumExclusionsPerAtom(const ListOfLists<int>& excls)
 {
     int maxNumExcls = 0;
-    for (int at = 0; at < excls.nr; at++)
+    for (gmx::index at = 0; at < excls.ssize(); at++)
     {
-        const int numExcls = excls.index[at + 1] - excls.index[at];
+        const auto list     = excls[at];
+        const int  numExcls = list.ssize();
 
-        GMX_RELEASE_ASSERT(numExcls != 1 || excls.a[excls.index[at]] == at,
+        GMX_RELEASE_ASSERT(numExcls != 1 || list[0] == at,
                            "With 1 exclusion we expect a self-exclusion");
 
         maxNumExcls = std::max(maxNumExcls, numExcls);
@@ -717,10 +719,7 @@ void dd_make_reverse_top(FILE*              fplog,
             make_reverse_top(mtop, ir->efep != efepNO, !dd->comm->systemInfo.haveSplitConstraints,
                              !dd->comm->systemInfo.haveSplitSettles, bBCheck, &dd->nbonded_global);
 
-    gmx_reverse_top_t* rt = dd->reverse_top;
-
     dd->haveExclusions = false;
-    rt->n_excl_at_max  = 0;
     for (const gmx_molblock_t& molb : mtop->molblock)
     {
         const int maxNumExclusionsPerAtom = getMaxNumExclusionsPerAtom(mtop->moltype[molb.type].excls);
@@ -729,7 +728,6 @@ void dd_make_reverse_top(FILE*              fplog,
         {
             dd->haveExclusions = true;
         }
-        rt->n_excl_at_max = std::max(rt->n_excl_at_max, maxNumExclusionsPerAtom);
     }
 
     if (vsite && vsite->numInterUpdategroupVsites > 0)
@@ -995,37 +993,12 @@ static real dd_dist2(t_pbc* pbc_null, const rvec* x, const int i, int j)
     return norm2(dx);
 }
 
-/*! \brief Append t_blocka block structures 1 to nsrc in src to *dest */
-static void combine_blocka(t_blocka* dest, gmx::ArrayRef<const thread_work_t> src)
+/*! \brief Append ListOfLists exclusion objects 1 to nsrc in \p src to \p *dest */
+static void combineExclusions(ListOfLists<int>* dest, gmx::ArrayRef<const thread_work_t> src)
 {
-    int ni = src.back().excl.nr;
-    int na = 0;
-    for (const thread_work_t& th_work : src)
-    {
-        na += th_work.excl.nra;
-    }
-    if (ni + 1 > dest->nalloc_index)
-    {
-        dest->nalloc_index = over_alloc_large(ni + 1);
-        srenew(dest->index, dest->nalloc_index);
-    }
-    if (dest->nra + na > dest->nalloc_a)
-    {
-        dest->nalloc_a = over_alloc_large(dest->nra + na);
-        srenew(dest->a, dest->nalloc_a);
-    }
     for (gmx::index s = 1; s < src.ssize(); s++)
     {
-        for (int i = dest->nr + 1; i < src[s].excl.nr + 1; i++)
-        {
-            dest->index[i] = dest->nra + src[s].excl.index[i];
-        }
-        for (int i = 0; i < src[s].excl.nra; i++)
-        {
-            dest->a[dest->nra + i] = src[s].excl.a[i];
-        }
-        dest->nr = src[s].excl.nr;
-        dest->nra += src[s].excl.nra;
+        dest->appendListOfLists(src[s].excl);
     }
 }
 
@@ -1377,11 +1350,11 @@ static int make_bondeds_zone(gmx_domdec_t*                      dd,
 }
 
 /*! \brief Set the exclusion data for i-zone \p iz for the case of no exclusions */
-static void set_no_exclusions_zone(const gmx_domdec_zones_t* zones, int iz, t_blocka* lexcls)
+static void set_no_exclusions_zone(const gmx_domdec_zones_t* zones, int iz, ListOfLists<int>* lexcls)
 {
     for (int a = zones->cg_range[iz]; a < zones->cg_range[iz + 1]; a++)
     {
-        lexcls->index[a + 1] = lexcls->nra;
+        lexcls->pushBack({});
     }
 }
 
@@ -1390,52 +1363,33 @@ static void make_exclusions_zone(gmx_domdec_t*                     dd,
                                  gmx_domdec_zones_t*               zones,
                                  const std::vector<gmx_moltype_t>& moltype,
                                  const int*                        cginfo,
-                                 t_blocka*                         lexcls,
+                                 ListOfLists<int>*                 lexcls,
                                  int                               iz,
                                  int                               at_start,
                                  int                               at_end,
                                  const gmx::ArrayRef<const int>    intermolecularExclusionGroup)
 {
-    int n_excl_at_max, n, at;
-
     const gmx_ga2la_t& ga2la = *dd->ga2la;
 
     const auto& jAtomRange = zones->iZones[iz].jAtomRange;
 
-    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 = at_end;
+    const gmx::index oldNumLists = lexcls->ssize();
 
-    n = lexcls->nra;
-    for (at = at_start; at < at_end; at++)
+    std::vector<int> exclusionsForAtom;
+    for (int at = at_start; at < at_end; at++)
     {
-        if (n + 1000 > lexcls->nalloc_a)
-        {
-            lexcls->nalloc_a = over_alloc_large(n + 1000);
-            srenew(lexcls->a, lexcls->nalloc_a);
-        }
+        exclusionsForAtom.clear();
 
         if (GET_CGINFO_EXCL_INTER(cginfo[at]))
         {
-            int             a_gl, mb, mt, mol, a_mol, j;
-            const t_blocka* excls;
-
-            if (n + n_excl_at_max > lexcls->nalloc_a)
-            {
-                lexcls->nalloc_a = over_alloc_large(n + n_excl_at_max);
-                srenew(lexcls->a, lexcls->nalloc_a);
-            }
+            int a_gl, mb, mt, mol, a_mol;
 
             /* Copy the exclusions from the global top */
-            lexcls->index[at] = n;
-            a_gl              = dd->globalAtomIndices[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;
-            for (j = excls->index[a_mol]; j < excls->index[a_mol + 1]; j++)
+            const auto excls = moltype[mt].excls[a_mol];
+            for (const int aj_mol : excls)
             {
-                const int aj_mol = excls->a[j];
-
                 if (const auto* jEntry = ga2la.find(a_gl + aj_mol - a_mol))
                 {
                     /* This check is not necessary, but it can reduce
@@ -1444,16 +1398,11 @@ static void make_exclusions_zone(gmx_domdec_t*                     dd,
                      */
                     if (jAtomRange.isInRange(jEntry->la))
                     {
-                        lexcls->a[n++] = jEntry->la;
+                        exclusionsForAtom.push_back(jEntry->la);
                     }
                 }
             }
         }
-        else
-        {
-            /* We don't need exclusions for this atom */
-            lexcls->index[at] = n;
-        }
 
         bool isExcludedAtom = !intermolecularExclusionGroup.empty()
                               && std::find(intermolecularExclusionGroup.begin(),
@@ -1462,51 +1411,26 @@ static void make_exclusions_zone(gmx_domdec_t*                     dd,
 
         if (isExcludedAtom)
         {
-            if (n + intermolecularExclusionGroup.ssize() > lexcls->nalloc_a)
-            {
-                lexcls->nalloc_a = over_alloc_large(n + intermolecularExclusionGroup.size());
-                srenew(lexcls->a, lexcls->nalloc_a);
-            }
             for (int qmAtomGlobalIndex : intermolecularExclusionGroup)
             {
                 if (const auto* entry = dd->ga2la->find(qmAtomGlobalIndex))
                 {
-                    lexcls->a[n++] = entry->la;
+                    exclusionsForAtom.push_back(entry->la);
                 }
             }
         }
-    }
-
-    lexcls->index[lexcls->nr] = n;
-    lexcls->nra               = n;
-}
 
-
-/*! \brief Ensure we have enough space in \p ba for \p nindex_max indices */
-static void check_alloc_index(t_blocka* ba, int nindex_max)
-{
-    if (nindex_max + 1 > ba->nalloc_index)
-    {
-        ba->nalloc_index = over_alloc_dd(nindex_max + 1);
-        srenew(ba->index, ba->nalloc_index);
+        /* Append the exclusions for this atom to the topology */
+        lexcls->pushBack(exclusionsForAtom);
     }
-}
 
-/*! \brief Ensure that we have enough space for exclusion storate in \p lexcls */
-static void check_exclusions_alloc(const gmx_domdec_t* dd, const gmx_domdec_zones_t* zones, t_blocka* lexcls)
-{
-    const int nr = zones->iZones.back().iAtomRange.end();
-
-    check_alloc_index(lexcls, nr);
-
-    for (size_t thread = 1; thread < dd->reverse_top->th_work.size(); thread++)
-    {
-        check_alloc_index(&dd->reverse_top->th_work[thread].excl, nr);
-    }
+    GMX_RELEASE_ASSERT(
+            lexcls->ssize() - oldNumLists == at_end - at_start,
+            "The number of exclusion list should match the number of atoms in the range");
 }
 
 /*! \brief Set the total count indexes for the local exclusions, needed by several functions */
-static void finish_local_exclusions(gmx_domdec_t* dd, gmx_domdec_zones_t* zones, t_blocka* lexcls)
+static void finish_local_exclusions(gmx_domdec_t* dd, gmx_domdec_zones_t* zones, ListOfLists<int>* lexcls)
 {
     const gmx::Range<int> nonhomeIzonesAtomRange(zones->iZones[0].iAtomRange.end(),
                                                  zones->iZones.back().iAtomRange.end());
@@ -1516,19 +1440,10 @@ static void finish_local_exclusions(gmx_domdec_t* dd, gmx_domdec_zones_t* zones,
         /* There are no exclusions involving non-home charge groups,
          * but we need to set the indices for neighborsearching.
          */
-        for (int la : nonhomeIzonesAtomRange)
+        for (int gmx_unused la : nonhomeIzonesAtomRange)
         {
-            lexcls->index[la] = lexcls->nra;
+            lexcls->pushBack({});
         }
-
-        /* 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 = nonhomeIzonesAtomRange.begin();
-    }
-    else
-    {
-        lexcls->nr = nonhomeIzonesAtomRange.end();
     }
 }
 
@@ -1556,7 +1471,7 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
                                     t_pbc*              pbc_null,
                                     rvec*               cg_cm,
                                     t_idef*             idef,
-                                    t_blocka*           lexcls,
+                                    ListOfLists<int>*   lexcls,
                                     int*                excl_count)
 {
     int                nzone_bondeds, nzone_excl;
@@ -1588,8 +1503,6 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
         nzone_excl = 1;
     }
 
-    check_exclusions_alloc(dd, zones, lexcls);
-
     rt = dd->reverse_top;
 
     rc2 = rc * rc;
@@ -1598,8 +1511,7 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
     clear_idef(idef);
     nbonded_local = 0;
 
-    lexcls->nr  = 0;
-    lexcls->nra = 0;
+    lexcls->clear();
     *excl_count = 0;
 
     for (int izone = 0; izone < nzone_bondeds; izone++)
@@ -1613,9 +1525,8 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
         {
             try
             {
-                int       cg0t, cg1t;
-                t_idef*   idef_t;
-                t_blocka* excl_t;
+                int     cg0t, cg1t;
+                t_idef* idef_t;
 
                 cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
                 cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
@@ -1636,15 +1547,17 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
 
                 if (izone < nzone_excl)
                 {
+                    ListOfLists<int>* excl_t;
                     if (thread == 0)
                     {
+                        // Thread 0 stores exclusions directly in the final storage
                         excl_t = lexcls;
                     }
                     else
                     {
-                        excl_t      = &rt->th_work[thread].excl;
-                        excl_t->nr  = 0;
-                        excl_t->nra = 0;
+                        // Threads > 0 store in temporary storage, starting at list index 0
+                        excl_t = &rt->th_work[thread].excl;
+                        excl_t->clear();
                     }
 
                     /* No charge groups and no distance check required */
@@ -1669,7 +1582,7 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
         {
             if (rt->th_work.size() > 1)
             {
-                combine_blocka(lexcls, rt->th_work);
+                combineExclusions(lexcls, rt->th_work);
             }
 
             for (const thread_work_t& th_work : rt->th_work)
@@ -1690,7 +1603,7 @@ static int make_local_bondeds_excls(gmx_domdec_t*       dd,
     finish_local_exclusions(dd, zones, lexcls);
     if (debug)
     {
-        fprintf(debug, "We have %d exclusions, check count %d\n", lexcls->nra, *excl_count);
+        fprintf(debug, "We have %d exclusions, check count %d\n", lexcls->numElements(), *excl_count);
     }
 
     return nbonded_local;
@@ -2062,12 +1975,11 @@ static void bonded_cg_distance_mol(const gmx_moltype_t* molt,
     }
     if (bExcl)
     {
-        const t_blocka* excls = &molt->excls;
-        for (int ai = 0; ai < excls->nr; ai++)
+        const auto& excls = molt->excls;
+        for (gmx::index ai = 0; ai < excls.ssize(); ai++)
         {
-            for (int j = excls->index[ai]; j < excls->index[ai + 1]; j++)
+            for (const int aj : excls[ai])
             {
-                int aj = excls->a[j];
                 if (ai != aj)
                 {
                     real rij2 = distance2(cg_cm[ai], cg_cm[aj]);
index 86698e5e82bfceee2d8883f17848023644b28aa4..801174e5e55f312f5f50c5bb489a3336350ef257 100644 (file)
@@ -2125,19 +2125,25 @@ static void do_block(gmx::ISerializer* serializer, t_block* block)
     serializer->doIntArray(block->index, block->nr + 1);
 }
 
-static void do_blocka(gmx::ISerializer* serializer, t_blocka* block)
+static void doListOfLists(gmx::ISerializer* serializer, gmx::ListOfLists<int>* listOfLists)
 {
-    serializer->doInt(&block->nr);
-    serializer->doInt(&block->nra);
+    int numLists = listOfLists->ssize();
+    serializer->doInt(&numLists);
+    int numElements = listOfLists->elementsView().ssize();
+    serializer->doInt(&numElements);
     if (serializer->reading())
     {
-        block->nalloc_index = block->nr + 1;
-        snew(block->index, block->nalloc_index);
-        block->nalloc_a = block->nra;
-        snew(block->a, block->nalloc_a);
+        std::vector<int> listRanges(numLists + 1);
+        serializer->doIntArray(listRanges.data(), numLists + 1);
+        std::vector<int> elements(numElements);
+        serializer->doIntArray(elements.data(), numElements);
+        *listOfLists = gmx::ListOfLists<int>(std::move(listRanges), std::move(elements));
+    }
+    else
+    {
+        serializer->doIntArray(const_cast<int*>(listOfLists->listRangesView().data()), numLists + 1);
+        serializer->doIntArray(const_cast<int*>(listOfLists->elementsView().data()), numElements);
     }
-    serializer->doIntArray(block->index, block->nr + 1);
-    serializer->doIntArray(block->a, block->nra);
 }
 
 /* This is a primitive routine to make it possible to translate atomic numbers
@@ -2450,7 +2456,7 @@ static void do_moltype(gmx::ISerializer* serializer, gmx_moltype_t* molt, t_symt
     sfree(cgs.index);
 
     /* This used to be in the atoms struct */
-    do_blocka(serializer, &molt->excls);
+    doListOfLists(serializer, &molt->excls);
 }
 
 static void do_molblock(gmx::ISerializer* serializer, gmx_molblock_t* molb, int numAtomsPerMolecule)
index 19aa460d27637b7f6bf294e815b7d27a392f8fb9..b411d6a602b2cba8438b8f561059e2b3fa4e9af2 100644 (file)
@@ -172,14 +172,15 @@ void __print_nnb(t_nextnb* nnb, char* s)
 }
 #endif
 
-static void nnb2excl(t_nextnb* nnb, t_blocka* excl)
+static void nnb2excl(t_nextnb* nnb, gmx::ListOfLists<int>* excls)
 {
     int       i, j, j_index;
     int       nre, nrx, nrs, nr_of_sortables;
     sortable* s;
 
-    srenew(excl->index, nnb->nr + 1);
-    excl->index[0] = 0;
+    excls->clear();
+
+    std::vector<int> exclusionsForAtom;
     for (i = 0; (i < nnb->nr); i++)
     {
         /* calculate the total number of exclusions for atom i */
@@ -230,16 +231,13 @@ static void nnb2excl(t_nextnb* nnb, t_blocka* excl)
         nr_of_sortables = j_index;
         prints("after rm-double", j_index, s);
 
-        /* make space for arrays */
-        srenew(excl->a, excl->nra + nr_of_sortables);
-
         /* put the sorted exclusions in the target list */
+        exclusionsForAtom.clear();
         for (nrs = 0; (nrs < nr_of_sortables); nrs++)
         {
-            excl->a[excl->nra + nrs] = s[nrs].aj;
+            exclusionsForAtom.push_back(s[nrs].aj);
         }
-        excl->nra += nr_of_sortables;
-        excl->index[i + 1] = excl->nra;
+        excls->pushBack(exclusionsForAtom);
 
         /* cleanup temporary space */
         sfree(s);
@@ -432,7 +430,7 @@ static void sort_and_purge_nnb(t_nextnb* nnb)
 }
 
 
-void generate_excl(int nrexcl, int nratoms, gmx::ArrayRef<InteractionsOfType> plist, t_blocka* excl)
+void generate_excl(int nrexcl, int nratoms, gmx::ArrayRef<InteractionsOfType> plist, gmx::ListOfLists<int>* excls)
 {
     t_nextnb nnb;
     if (nrexcl < 0)
@@ -441,8 +439,7 @@ void generate_excl(int nrexcl, int nratoms, gmx::ArrayRef<InteractionsOfType> pl
     }
     init_nnb(&nnb, nratoms, nrexcl);
     gen_nnb(&nnb, plist);
-    excl->nr = nratoms;
     sort_and_purge_nnb(&nnb);
-    nnb2excl(&nnb, excl);
+    nnb2excl(&nnb, excls);
     done_nnb(&nnb);
 }
index 4d014cdf97a6afb61f478eb8b7cd902251330212..8e640c1bab5060438ad4ac42e9ef9bc3242508c3 100644 (file)
 
 #include "gromacs/utility/arrayref.h"
 
-struct t_blocka;
 struct InteractionsOfType;
 
+namespace gmx
+{
+template<typename>
+class ListOfLists;
+}
+
 struct t_nextnb
 {
     int nr;   /* nr atoms (0 <= i < nr) (atoms->nr)            */
@@ -75,7 +80,7 @@ void gen_nnb(t_nextnb* nnb, gmx::ArrayRef<InteractionsOfType> plist);
  * initiated using init_nnb.
  */
 
-void generate_excl(int nrexcl, int nratoms, gmx::ArrayRef<InteractionsOfType> plist, t_blocka* excl);
+void generate_excl(int nrexcl, int nratoms, gmx::ArrayRef<InteractionsOfType> plist, gmx::ListOfLists<int>* excls);
 /* Generate an exclusion block from bonds and constraints in
  * plist.
  */
index 3922eeb01c1e046fad881d005632817aa2e477c6..735b5c37d67350b6bc66e4970ed13fe00d2276ff 100644 (file)
@@ -234,7 +234,7 @@ void InteractionOfType::setForceParameter(int pos, real value)
 void MoleculeInformation::initMolInfo()
 {
     init_block(&mols);
-    init_blocka(&excls);
+    excls.clear();
     init_t_atoms(&atoms, 0, FALSE);
 }
 
index 1e06a7d6da5103cef538af118254c62822ce52ca..7eaa8f1c233a00b04df356fe431adaf72e6d0f35 100644 (file)
@@ -47,6 +47,7 @@
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/basedefinitions.h"
 #include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/real.h"
 
 /*! \libinternal \brief
@@ -176,7 +177,7 @@ struct MoleculeInformation
     //! Molecules separated in datastructure.
     t_block mols;
     //! Exclusions in the molecule.
-    t_blocka excls;
+    gmx::ListOfLists<int> excls;
     //! Interactions of a defined type.
     std::array<InteractionsOfType, F_NRE> interactions;
 
index 7058eb623b62caf89134349657154f9ac6c6a96a..041cb951369c199c3d1edc0ab149e35e61be157b 100644 (file)
@@ -1376,7 +1376,7 @@ void generate_qmexcl(gmx_mtop_t* sys, t_inputrec* ir, warninp* wi, GmxQmmmMode q
                     /* Copy the exclusions to a new array, since this is the only
                      * thing that needs to be modified for QMMM.
                      */
-                    copy_blocka(&sys->moltype[molb->type].excls, &sys->moltype.back().excls);
+                    sys->moltype.back().excls = sys->moltype[molb->type].excls;
                     /* Set the molecule type for the QMMM molblock */
                     molb->type = sys->moltype.size() - 1;
                 }
index af199cdc799e50b0184338e1ce7ccc60d8217dc9..f78a72af0282f114c0864a2d2c884f357f44360d 100644 (file)
@@ -2615,10 +2615,8 @@ static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionsOfType> interactio
 
 static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, InteractionsOfType* nbp, warninp* wi)
 {
-    int       n, ntype;
-    t_atom*   atom;
-    t_blocka* excl;
-    bool      bExcl;
+    int     n, ntype;
+    t_atom* atom;
 
     n    = mol->atoms.nr;
     atom = mol->atoms.atom;
@@ -2628,20 +2626,20 @@ static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, Interact
                "Number of pairs of generated non-bonded parameters should be a perfect square");
 
     /* Add a pair interaction for all non-excluded atom pairs */
-    excl = &mol->excls;
+    const auto& excls = mol->excls;
     for (int i = 0; i < n; i++)
     {
         for (int j = i + 1; j < n; j++)
         {
-            bExcl = FALSE;
-            for (int k = excl->index[i]; k < excl->index[i + 1]; k++)
+            bool pairIsExcluded = false;
+            for (const int atomK : excls[i])
             {
-                if (excl->a[k] == j)
+                if (atomK == j)
                 {
-                    bExcl = TRUE;
+                    pairIsExcluded = true;
                 }
             }
-            if (!bExcl)
+            if (!pairIsExcluded)
             {
                 if (nb_funct != F_LJ)
                 {
@@ -2661,24 +2659,20 @@ static void generate_LJCpairsNB(MoleculeInformation* mol, int nb_funct, Interact
     }
 }
 
-static void set_excl_all(t_blocka* excl)
+static void set_excl_all(gmx::ListOfLists<int>* excl)
 {
-    int nat, i, j, k;
-
     /* Get rid of the current exclusions and exclude all atom pairs */
-    nat       = excl->nr;
-    excl->nra = nat * nat;
-    srenew(excl->a, excl->nra);
-    k = 0;
-    for (i = 0; i < nat; i++)
+    const int        numAtoms = excl->ssize();
+    std::vector<int> exclusionsForAtom(numAtoms);
+    for (int i = 0; i < numAtoms; i++)
     {
-        excl->index[i] = k;
-        for (j = 0; j < nat; j++)
-        {
-            excl->a[k++] = j;
-        }
+        exclusionsForAtom[i] = i;
+    }
+    excl->clear();
+    for (int i = 0; i < numAtoms; i++)
+    {
+        excl->pushBack(exclusionsForAtom);
     }
-    excl->index[nat] = k;
 }
 
 static void decouple_atoms(t_atoms*    atoms,
index 7cbd9e41c7c76cc1cf368a1033193d4840fc4b8c..e436bc9d42614f2ad2bed5543efc4d561f37a17c 100644 (file)
@@ -186,17 +186,14 @@ DispersionCorrection::TopologyParams::TopologyParams(const gmx_mtop_t&         m
              */
             for (const gmx_molblock_t& molb : mtop.molblock)
             {
-                const int       nmol  = molb.nmol;
-                const t_atoms*  atoms = &mtop.moltype[molb.type].atoms;
-                const t_blocka* excl  = &mtop.moltype[molb.type].excls;
+                const int      nmol  = molb.nmol;
+                const t_atoms* atoms = &mtop.moltype[molb.type].atoms;
+                const auto&    excl  = mtop.moltype[molb.type].excls;
                 for (int i = 0; (i < atoms->nr); i++)
                 {
                     const int tpi = atomtypeAOrB(atoms->atom[i], q);
-                    const int j1  = excl->index[i];
-                    const int j2  = excl->index[i + 1];
-                    for (int j = j1; j < j2; j++)
+                    for (const int k : excl[i])
                     {
-                        const int k = excl->a[j];
                         if (k > i)
                         {
                             const int tpj = atomtypeAOrB(atoms->atom[k], q);
index 66c4c33b5647d48e047a7d9ff49dc7f9224b7f78..9440cf8f088e99c42c1f1737eb7df1d9705914ef 100644 (file)
@@ -212,7 +212,7 @@ static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_f
     {
         const gmx_molblock_t& molb = mtop->molblock[mb];
         const gmx_moltype_t&  molt = mtop->moltype[molb.type];
-        const t_blocka&       excl = molt.excls;
+        const auto&           excl = molt.excls;
 
         /* Check if the cginfo is identical for all molecules in this block.
          * If so, we only need an array of the size of one molecule.
@@ -285,9 +285,9 @@ static std::vector<cginfo_mb_t> init_cginfo_mb(const gmx_mtop_t* mtop, const t_f
 
                 bool haveExclusions = false;
                 /* Loop over all the exclusions of atom ai */
-                for (int j = excl.index[a]; j < excl.index[a + 1]; j++)
+                for (const int j : excl[a])
                 {
-                    if (excl.a[j] != a)
+                    if (j != a)
                     {
                         haveExclusions = true;
                         break;
index c93de7f1051e1d0f321699dce161a9e4a25f73a5..b2b061061fa0bd0af50130b08992fc919bdeb5cf 100644 (file)
@@ -233,7 +233,7 @@ void count_bonded_distances(const gmx_mtop_t& mtop, const t_inputrec& ir, double
         }
         if (bExcl)
         {
-            ndtot_c += molb.nmol * (molt->excls.nra - molt->atoms.nr) / 2.;
+            ndtot_c += molb.nmol * (molt->excls.numElements() - molt->atoms.nr) / 2.;
         }
     }
 
index 692c2570b66f9dfa613e29dae78d4f183685027f..8e0541bc8634e80e9a2c7999c16b286e60a88948 100644 (file)
@@ -1143,7 +1143,7 @@ void do_force(FILE*                               fplog,
         wallcycle_start_nocount(wcycle, ewcNS);
         wallcycle_sub_start(wcycle, ewcsNBS_SEARCH_LOCAL);
         /* Note that with a GPU the launch overhead of the list transfer is not timed separately */
-        nbv->constructPairlist(InteractionLocality::Local, &top->excls, step, nrnb);
+        nbv->constructPairlist(InteractionLocality::Local, top->excls, step, nrnb);
 
         nbv->setupGpuShortRangeWork(fr->gpuBonded, InteractionLocality::Local);
 
@@ -1242,7 +1242,7 @@ void do_force(FILE*                               fplog,
             wallcycle_start_nocount(wcycle, ewcNS);
             wallcycle_sub_start(wcycle, ewcsNBS_SEARCH_NONLOCAL);
             /* Note that with a GPU the launch overhead of the list transfer is not timed separately */
-            nbv->constructPairlist(InteractionLocality::NonLocal, &top->excls, step, nrnb);
+            nbv->constructPairlist(InteractionLocality::NonLocal, top->excls, step, nrnb);
 
             nbv->setupGpuShortRangeWork(fr->gpuBonded, InteractionLocality::NonLocal);
             wallcycle_sub_stop(wcycle, ewcsNBS_SEARCH_NONLOCAL);
index 2ed8bd4586eb53126c8881264dbea9f9eb3ce17c..1c9066c7be0929d416054a47021176c1e0e3e177 100644 (file)
@@ -666,7 +666,7 @@ void LegacySimulator::do_tpi()
                 /* TODO: Avoid updating all atoms at every bNS step */
                 fr->nbv->setAtomProperties(*mdatoms, fr->cginfo);
 
-                fr->nbv->constructPairlist(InteractionLocality::Local, &top.excls, step, nrnb);
+                fr->nbv->constructPairlist(InteractionLocality::Local, top.excls, step, nrnb);
 
                 bNS = FALSE;
             }
index f3d6b27526a0eee43ed7f1786b5ed1772d0f812c..c2c2328f66508dda20d36da78c8cfa96777c1e63 100644 (file)
@@ -223,7 +223,7 @@ static std::unique_ptr<nonbonded_verlet_t> setupNbnxmForBenchInstance(const Kern
                       { 0, int(system.coordinates.size()) }, atomDensity, atomInfo,
                       system.coordinates, 0, nullptr);
 
-    nbv->constructPairlist(gmx::InteractionLocality::Local, &system.excls, 0, &nrnb);
+    nbv->constructPairlist(gmx::InteractionLocality::Local, system.excls, 0, &nrnb);
 
     t_mdatoms mdatoms;
     // We only use (read) the atom type and charge from mdatoms
index fba71d2b9a0ca899ff30e3257d137485e354da7e..8ed07ece85b992c02c70ef9e88cd298422978ecc 100644 (file)
@@ -167,10 +167,8 @@ BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor)
     charges.resize(numAtoms);
     atomInfoAllVdw.resize(numAtoms);
     atomInfoOxygenVdw.resize(numAtoms);
-    snew(excls.index, numAtoms + 1);
-    snew(excls.a, numAtoms * numAtomsInMolecule);
-    excls.index[0] = 0;
 
+    std::vector<int> exclusionsForAtom;
     for (int a = 0; a < numAtoms; a++)
     {
         if (a % numAtomsInMolecule == 0)
@@ -191,12 +189,13 @@ BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor)
         SET_CGINFO_HAS_Q(atomInfoAllVdw[a]);
         SET_CGINFO_HAS_Q(atomInfoOxygenVdw[a]);
 
+        exclusionsForAtom.clear();
         const int firstAtomInMolecule = a - (a % numAtomsInMolecule);
         for (int aj = 0; aj < numAtomsInMolecule; aj++)
         {
-            excls.a[a * numAtomsInMolecule + aj] = firstAtomInMolecule + aj;
+            exclusionsForAtom.push_back(firstAtomInMolecule + aj);
         }
-        excls.index[a + 1] = (a + 1) * numAtomsInMolecule;
+        excls.pushBack(exclusionsForAtom);
     }
 
     forceRec.ntype = numAtomTypes;
index 512368ddda507fe8937ab58d4e797a08d813553a..adcc85d4ffacc443b36b48dc8040e09ffa2fb1aa 100644 (file)
@@ -48,7 +48,7 @@
 
 #include "gromacs/math/vectypes.h"
 #include "gromacs/mdtypes/forcerec.h"
-#include "gromacs/topology/block.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 
 namespace gmx
@@ -80,7 +80,7 @@ struct BenchmarkSystem
     //! Atom info where only oxygen atoms are marked to have Van der Waals interactions
     std::vector<int> atomInfoOxygenVdw;
     //! Information about exclusions.
-    t_blocka excls;
+    ListOfLists<int> excls;
     //! Storage for atom positions.
     std::vector<gmx::RVec> coordinates;
     //! System simulation box.
index 1548f3704ff9c74b83d01b82d29a29ce1e2a1cd4..cb5da28eede392e4da2e658a7c4fbab6eeba48ff 100644 (file)
@@ -133,7 +133,6 @@ struct interaction_const_t;
 struct nonbonded_verlet_t;
 class PairSearch;
 class PairlistSets;
-struct t_blocka;
 struct t_commrec;
 struct t_lambda;
 struct t_mdatoms;
@@ -153,6 +152,8 @@ class GpuEventSynchronizer;
 namespace gmx
 {
 class ForceWithShiftForces;
+template<typename>
+class ListOfLists;
 class MDLogger;
 class UpdateGroupsCog;
 } // namespace gmx
@@ -251,7 +252,10 @@ public:
     gmx::ArrayRef<const int> getGridIndices() const;
 
     //! Constructs the pairlist for the given locality
-    void constructPairlist(gmx::InteractionLocality iLocality, const t_blocka* excl, int64_t step, t_nrnb* nrnb);
+    void constructPairlist(gmx::InteractionLocality     iLocality,
+                           const gmx::ListOfLists<int>& exclusions,
+                           int64_t                      step,
+                           t_nrnb*                      nrnb);
 
     //! Updates all the atom properties in Nbnxm
     void setAtomProperties(const t_mdatoms& mdatoms, gmx::ArrayRef<const int> atomInfo);
index 6174905d578fffe2f68497beee7ea8e84cf12139..f5efbef63864fe3af81188c9637cc4fd93d4b1e4 100644 (file)
 #include "gromacs/pbcutil/pbc.h"
 #include "gromacs/simd/simd.h"
 #include "gromacs/simd/vector_operations.h"
-#include "gromacs/topology/block.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 
 #include "atomdata.h"
@@ -1363,12 +1363,12 @@ static nbnxn_sci_t* getOpenIEntry(NbnxnPairlistGpu* nbl)
  * Set all atom-pair exclusions from the topology stored in exclusions
  * as masks in the pair-list for simple list entry iEntry.
  */
-static void setExclusionsForIEntry(const Nbnxm::GridSet& gridSet,
-                                   NbnxnPairlistCpu*     nbl,
-                                   gmx_bool              diagRemoved,
-                                   int                   na_cj_2log,
-                                   const nbnxn_ci_t&     iEntry,
-                                   const t_blocka&       exclusions)
+static void setExclusionsForIEntry(const Nbnxm::GridSet&   gridSet,
+                                   NbnxnPairlistCpu*       nbl,
+                                   gmx_bool                diagRemoved,
+                                   int                     na_cj_2log,
+                                   const nbnxn_ci_t&       iEntry,
+                                   const ListOfLists<int>& exclusions)
 {
     if (iEntry.cj_ind_end == iEntry.cj_ind_start)
     {
@@ -1391,11 +1391,8 @@ static void setExclusionsForIEntry(const Nbnxm::GridSet& gridSet,
         if (iAtom >= 0)
         {
             /* Loop over the topology-based exclusions for this i-atom */
-            for (int exclIndex = exclusions.index[iAtom]; exclIndex < exclusions.index[iAtom + 1];
-                 exclIndex++)
+            for (const int jAtom : exclusions[iAtom])
             {
-                const int jAtom = exclusions.a[exclIndex];
-
                 if (jAtom == iAtom)
                 {
                     /* The self exclusion are already set, save some time */
@@ -1876,9 +1873,9 @@ static void make_fep_list(gmx::ArrayRef<const int> atomIndices,
 static void setExclusionsForIEntry(const Nbnxm::GridSet& gridSet,
                                    NbnxnPairlistGpu*     nbl,
                                    gmx_bool              diagRemoved,
-                                   int gmx_unused     na_cj_2log,
-                                   const nbnxn_sci_t& iEntry,
-                                   const t_blocka&    exclusions)
+                                   int gmx_unused          na_cj_2log,
+                                   const nbnxn_sci_t&      iEntry,
+                                   const ListOfLists<int>& exclusions)
 {
     if (iEntry.numJClusterGroups() == 0)
     {
@@ -1912,11 +1909,8 @@ static void setExclusionsForIEntry(const Nbnxm::GridSet& gridSet,
             const int iCluster = i / c_clusterSize;
 
             /* Loop over the topology-based exclusions for this i-atom */
-            for (int exclIndex = exclusions.index[iAtom]; exclIndex < exclusions.index[iAtom + 1];
-                 exclIndex++)
+            for (const int jAtom : exclusions[iAtom])
             {
-                const int jAtom = exclusions.a[exclIndex];
-
                 if (jAtom == iAtom)
                 {
                     /* The self exclusions are already set, save some time */
@@ -3083,7 +3077,7 @@ static void nbnxn_make_pairlist_part(const Nbnxm::GridSet&   gridSet,
                                      const Grid&             jGrid,
                                      PairsearchWork*         work,
                                      const nbnxn_atomdata_t* nbat,
-                                     const t_blocka&         exclusions,
+                                     const ListOfLists<int>& exclusions,
                                      real                    rlist,
                                      const PairlistType      pairlistType,
                                      int                     ci_block,
@@ -3924,7 +3918,7 @@ static void prepareListsForDynamicPruning(gmx::ArrayRef<NbnxnPairlistCpu> lists)
 void PairlistSet::constructPairlists(const Nbnxm::GridSet&         gridSet,
                                      gmx::ArrayRef<PairsearchWork> searchWork,
                                      nbnxn_atomdata_t*             nbat,
-                                     const t_blocka*               excl,
+                                     const ListOfLists<int>&       exclusions,
                                      const int                     minimumIlistCountForGpuBalancing,
                                      t_nrnb*                       nrnb,
                                      SearchCycleCounting*          searchCycleCounting)
@@ -4037,14 +4031,14 @@ void PairlistSet::constructPairlists(const Nbnxm::GridSet&         gridSet,
                     /* Divide the i cells equally over the pairlists */
                     if (isCpuType_)
                     {
-                        nbnxn_make_pairlist_part(gridSet, iGrid, jGrid, &work, nbat, *excl, rlist,
+                        nbnxn_make_pairlist_part(gridSet, iGrid, jGrid, &work, nbat, exclusions, rlist,
                                                  params_.pairlistType, ci_block, nbat->bUseBufferFlags,
                                                  nsubpair_target, progBal, nsubpair_tot_est, th,
                                                  numLists, &cpuLists_[th], fepListPtr);
                     }
                     else
                     {
-                        nbnxn_make_pairlist_part(gridSet, iGrid, jGrid, &work, nbat, *excl, rlist,
+                        nbnxn_make_pairlist_part(gridSet, iGrid, jGrid, &work, nbat, exclusions, rlist,
                                                  params_.pairlistType, ci_block, nbat->bUseBufferFlags,
                                                  nsubpair_target, progBal, nsubpair_tot_est, th,
                                                  numLists, &gpuLists_[th], fepListPtr);
@@ -4202,12 +4196,12 @@ void PairlistSet::constructPairlists(const Nbnxm::GridSet&         gridSet,
 void PairlistSets::construct(const InteractionLocality iLocality,
                              PairSearch*               pairSearch,
                              nbnxn_atomdata_t*         nbat,
-                             const t_blocka*           excl,
+                             const ListOfLists<int>&   exclusions,
                              const int64_t             step,
                              t_nrnb*                   nrnb)
 {
-    pairlistSet(iLocality).constructPairlists(pairSearch->gridSet(), pairSearch->work(), nbat, excl,
-                                              minimumIlistCountForGpuBalancing_, nrnb,
+    pairlistSet(iLocality).constructPairlists(pairSearch->gridSet(), pairSearch->work(), nbat,
+                                              exclusions, minimumIlistCountForGpuBalancing_, nrnb,
                                               &pairSearch->cycleCounting_);
 
     if (iLocality == InteractionLocality::Local)
@@ -4234,11 +4228,11 @@ void PairlistSets::construct(const InteractionLocality iLocality,
 }
 
 void nonbonded_verlet_t::constructPairlist(const InteractionLocality iLocality,
-                                           const t_blocka*           excl,
+                                           const ListOfLists<int>&   exclusions,
                                            int64_t                   step,
                                            t_nrnb*                   nrnb)
 {
-    pairlistSets_->construct(iLocality, pairSearch_.get(), nbat.get(), excl, step, nrnb);
+    pairlistSets_->construct(iLocality, pairSearch_.get(), nbat.get(), exclusions, step, nrnb);
 
     if (useGpu())
     {
index d0e930330c54202269dabe7b4878d6dc9332f38c..1133846a6bccdcf774b632c5af7be9742cae762c 100644 (file)
@@ -62,9 +62,14 @@ struct nbnxn_atomdata_t;
 struct PairlistParams;
 struct PairsearchWork;
 struct SearchCycleCounting;
-struct t_blocka;
 struct t_nrnb;
 
+namespace gmx
+{
+template<typename>
+class ListOfLists;
+}
+
 namespace Nbnxm
 {
 class GridSet;
@@ -85,7 +90,7 @@ public:
     void constructPairlists(const Nbnxm::GridSet&         gridSet,
                             gmx::ArrayRef<PairsearchWork> searchWork,
                             nbnxn_atomdata_t*             nbat,
-                            const t_blocka*               excl,
+                            const gmx::ListOfLists<int>&  exclusions,
                             int                           minimumIlistCountForGpuBalancing,
                             t_nrnb*                       nrnb,
                             SearchCycleCounting*          searchCycleCounting);
index 2349159a1dd328df5537ebe8166a395f0490eff9..960f5f96791dbf4e2a529672a28b1b779f6524a9 100644 (file)
@@ -57,9 +57,13 @@ struct nbnxn_atomdata_t;
 class PairlistSet;
 enum class PairlistType;
 class PairSearch;
-struct t_blocka;
 struct t_nrnb;
 
+namespace gmx
+{
+template<typename>
+class ListOfLists;
+}
 
 class PairlistSets
 {
@@ -69,12 +73,12 @@ public:
                  int                   minimumIlistCountForGpuBalancing);
 
     //! Construct the pairlist set for the given locality
-    void construct(gmx::InteractionLocality iLocality,
-                   PairSearch*              pairSearch,
-                   nbnxn_atomdata_t*        nbat,
-                   const t_blocka*          excl,
-                   int64_t                  step,
-                   t_nrnb*                  nrnb);
+    void construct(gmx::InteractionLocality     iLocality,
+                   PairSearch*                  pairSearch,
+                   nbnxn_atomdata_t*            nbat,
+                   const gmx::ListOfLists<int>& exclusions,
+                   int64_t                      step,
+                   t_nrnb*                      nrnb);
 
     //! Dispatches the dynamic pruning kernel for the given locality
     void dispatchPruneKernel(gmx::InteractionLocality iLocality,
index ddc8bd67d341b1dee67fb8009a6f21076bf6ed13..08d8f6cc7b0d028933598d4c3ec8577aba4ca374 100644 (file)
 #include "gromacs/math/functions.h"
 #include "gromacs/math/vec.h"
 #include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/block.h"
 #include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/exceptions.h"
 #include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/mutex.h"
 #include "gromacs/utility/stringutil.h"
 
@@ -137,7 +137,7 @@ public:
      */
     void                  init(AnalysisNeighborhood::SearchMode     mode,
                                bool                                 bXY,
-                               const t_blocka*                      excls,
+                               const ListOfLists<int>*              excls,
                                const t_pbc*                         pbc,
                                const AnalysisNeighborhoodPositions& positions);
     PairSearchImplPointer getPairSearch();
@@ -280,7 +280,7 @@ private:
     //! Reference position indices (NULL if no indices).
     const int* refIndices_;
     //! Exclusions.
-    const t_blocka* excls_;
+    const ListOfLists<int>* excls_;
     //! PBC data.
     t_pbc pbc_;
 
@@ -337,8 +337,6 @@ public:
         testPositions_    = nullptr;
         testExclusionIds_ = nullptr;
         testIndices_      = nullptr;
-        nexcl_            = 0;
-        excl_             = nullptr;
         clear_rvec(xtest_);
         clear_rvec(testcell_);
         clear_ivec(currCell_);
@@ -376,10 +374,8 @@ private:
     const int* testExclusionIds_;
     //! Reference to the test position indices.
     const int* testIndices_;
-    //! Number of excluded reference positions for current test particle.
-    int nexcl_;
     //! Exclusions for current test particle.
-    const int* excl_;
+    ArrayRef<const int> excl_;
     //! Index of the currently active test position in \p testPositions_.
     int testIndex_;
     //! Stores test position during a pair loop.
@@ -862,7 +858,7 @@ int AnalysisNeighborhoodSearchImpl::shiftCell(const ivec cell, rvec shift) const
 
 void AnalysisNeighborhoodSearchImpl::init(AnalysisNeighborhood::SearchMode     mode,
                                           bool                                 bXY,
-                                          const t_blocka*                      excls,
+                                          const ListOfLists<int>*              excls,
                                           const t_pbc*                         pbc,
                                           const AnalysisNeighborhoodPositions& positions)
 {
@@ -988,16 +984,13 @@ void AnalysisNeighborhoodPairSearchImpl::reset(int testIndex)
         if (search_.excls_ != nullptr)
         {
             const int exclIndex = testExclusionIds_[index];
-            if (exclIndex < search_.excls_->nr)
+            if (exclIndex < search_.excls_->ssize())
             {
-                const int startIndex = search_.excls_->index[exclIndex];
-                nexcl_               = search_.excls_->index[exclIndex + 1] - startIndex;
-                excl_                = &search_.excls_->a[startIndex];
+                excl_ = (*search_.excls_)[exclIndex];
             }
             else
             {
-                nexcl_ = 0;
-                excl_  = nullptr;
+                excl_ = ArrayRef<const int>();
             }
         }
     }
@@ -1014,15 +1007,16 @@ void AnalysisNeighborhoodPairSearchImpl::nextTestPosition()
 
 bool AnalysisNeighborhoodPairSearchImpl::isExcluded(int j)
 {
-    if (exclind_ < nexcl_)
+    const int nexcl = excl_.ssize();
+    if (exclind_ < nexcl)
     {
         const int index = (search_.refIndices_ != nullptr ? search_.refIndices_[j] : j);
         const int refId = search_.refExclusionIds_[index];
-        while (exclind_ < nexcl_ && excl_[exclind_] < refId)
+        while (exclind_ < nexcl && excl_[exclind_] < refId)
         {
             ++exclind_;
         }
-        if (exclind_ < nexcl_ && refId == excl_[exclind_])
+        if (exclind_ < nexcl && refId == excl_[exclind_])
         {
             ++exclind_;
             return true;
@@ -1258,12 +1252,12 @@ public:
 
     SearchImplPointer getSearch();
 
-    Mutex           createSearchMutex_;
-    SearchList      searchList_;
-    real            cutoff_;
-    const t_blocka* excls_;
-    SearchMode      mode_;
-    bool            bXY_;
+    Mutex                   createSearchMutex_;
+    SearchList              searchList_;
+    real                    cutoff_;
+    const ListOfLists<int>* excls_;
+    SearchMode              mode_;
+    bool                    bXY_;
 };
 
 AnalysisNeighborhood::Impl::SearchImplPointer AnalysisNeighborhood::Impl::getSearch()
@@ -1304,7 +1298,7 @@ void AnalysisNeighborhood::setXYMode(bool bXY)
     impl_->bXY_ = bXY;
 }
 
-void AnalysisNeighborhood::setTopologyExclusions(const t_blocka* excls)
+void AnalysisNeighborhood::setTopologyExclusions(const ListOfLists<int>* excls)
 {
     GMX_RELEASE_ASSERT(impl_->searchList_.empty(),
                        "Changing the exclusions after initSearch() not currently supported");
index ffe93854b644f4ddba04e4a2a41c4ac90f37f7ba..7066cafd12d14bedad6061bdf7e26e821a58c8ff 100644 (file)
 #include "gromacs/utility/gmxassert.h"
 #include "gromacs/utility/real.h"
 
-struct t_blocka;
 struct t_pbc;
 
 namespace gmx
 {
+template<typename>
+class ListOfLists;
 
 namespace internal
 {
@@ -282,7 +283,7 @@ public:
      *
      * \see AnalysisNeighborhoodPositions::exclusionIds()
      */
-    void setTopologyExclusions(const t_blocka* excls);
+    void setTopologyExclusions(const ListOfLists<int>* excls);
     /*! \brief
      * Sets the algorithm to use for searching.
      *
index 7e85995a1dda77e71e6e2e62fd10a1b9812eef4b..ec81ab55889f14a2fc8a1cbb706aa87d32ce577a 100644 (file)
@@ -64,6 +64,7 @@
 #include "gromacs/random/threefry.h"
 #include "gromacs/random/uniformrealdistribution.h"
 #include "gromacs/topology/block.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/stringutil.h"
 
@@ -309,13 +310,13 @@ void NeighborhoodSearchTestData::computeReferencesInternal(t_pbc* pbc, bool bXY)
 class ExclusionsHelper
 {
 public:
-    static void markExcludedPairs(RefPairList* refPairs, int testIndex, const t_blocka* excls);
+    static void markExcludedPairs(RefPairList* refPairs, int testIndex, const gmx::ListOfLists<int>* excls);
 
     ExclusionsHelper(int refPosCount, int testPosCount);
 
     void generateExclusions();
 
-    const t_blocka* exclusions() const { return &excls_; }
+    const gmx::ListOfLists<int>* exclusions() const { return &excls_; }
 
     gmx::ArrayRef<const int> refPosIds() const
     {
@@ -327,21 +328,18 @@ public:
     }
 
 private:
-    int              refPosCount_;
-    int              testPosCount_;
-    std::vector<int> exclusionIds_;
-    std::vector<int> exclsIndex_;
-    std::vector<int> exclsAtoms_;
-    t_blocka         excls_;
+    int                   refPosCount_;
+    int                   testPosCount_;
+    std::vector<int>      exclusionIds_;
+    gmx::ListOfLists<int> excls_;
 };
 
 // static
-void ExclusionsHelper::markExcludedPairs(RefPairList* refPairs, int testIndex, const t_blocka* excls)
+void ExclusionsHelper::markExcludedPairs(RefPairList* refPairs, int testIndex, const gmx::ListOfLists<int>* excls)
 {
     int count = 0;
-    for (int i = excls->index[testIndex]; i < excls->index[testIndex + 1]; ++i)
+    for (const int excludedIndex : (*excls)[testIndex])
     {
-        const int                           excludedIndex = excls->a[i];
         NeighborhoodSearchTestData::RefPair searchPair(excludedIndex, 0.0);
         RefPairList::iterator               excludedRefPair =
                 std::lower_bound(refPairs->begin(), refPairs->end(), searchPair);
@@ -364,13 +362,6 @@ ExclusionsHelper::ExclusionsHelper(int refPosCount, int testPosCount) :
     exclusionIds_.resize(std::max(refPosCount, testPosCount), 1);
     exclusionIds_[0] = 0;
     std::partial_sum(exclusionIds_.begin(), exclusionIds_.end(), exclusionIds_.begin());
-
-    excls_.nr           = 0;
-    excls_.index        = nullptr;
-    excls_.nra          = 0;
-    excls_.a            = nullptr;
-    excls_.nalloc_index = 0;
-    excls_.nalloc_a     = 0;
 }
 
 void ExclusionsHelper::generateExclusions()
@@ -379,21 +370,16 @@ void ExclusionsHelper::generateExclusions()
     // particles would be higher, or where the exclusions would not be random,
     // to make a higher percentage of the exclusions to actually be within the
     // cutoff.
-    exclsIndex_.reserve(testPosCount_ + 1);
-    exclsAtoms_.reserve(testPosCount_ * 20);
-    exclsIndex_.push_back(0);
+    std::vector<int> exclusionsForAtom;
     for (int i = 0; i < testPosCount_; ++i)
     {
+        exclusionsForAtom.clear();
         for (int j = 0; j < 20; ++j)
         {
-            exclsAtoms_.push_back(i + j * 3);
+            exclusionsForAtom.push_back(i + j * 3);
         }
-        exclsIndex_.push_back(exclsAtoms_.size());
+        excls_.pushBack(exclusionsForAtom);
     }
-    excls_.nr    = exclsIndex_.size();
-    excls_.index = exclsIndex_.data();
-    excls_.nra   = exclsAtoms_.size();
-    excls_.a     = exclsAtoms_.data();
 }
 
 /********************************************************************
@@ -414,7 +400,7 @@ public:
     void testPairSearchFull(gmx::AnalysisNeighborhoodSearch*          search,
                             const NeighborhoodSearchTestData&         data,
                             const gmx::AnalysisNeighborhoodPositions& pos,
-                            const t_blocka*                           excls,
+                            const gmx::ListOfLists<int>*              excls,
                             const gmx::ArrayRef<const int>&           refIndices,
                             const gmx::ArrayRef<const int>&           testIndices,
                             bool                                      selfPairs);
@@ -522,7 +508,7 @@ void NeighborhoodSearchTest::testPairSearchIndexed(gmx::AnalysisNeighborhood*
 void NeighborhoodSearchTest::testPairSearchFull(gmx::AnalysisNeighborhoodSearch*          search,
                                                 const NeighborhoodSearchTestData&         data,
                                                 const gmx::AnalysisNeighborhoodPositions& pos,
-                                                const t_blocka*                           excls,
+                                                const gmx::ListOfLists<int>*              excls,
                                                 const gmx::ArrayRef<const int>& refIndices,
                                                 const gmx::ArrayRef<const int>& testIndices,
                                                 bool                            selfPairs)
index 09adddc79534ef0660d78c4fd3eda462db9ff7cb..d4d4685a5878de5d60f684147358f9ae04935d73 100644 (file)
@@ -129,39 +129,33 @@ static void reduce_block(const gmx_bool bKeep[], t_block* block, const char* nam
     block->nr    = newi;
 }
 
-static void reduce_blocka(const int invindex[], const gmx_bool bKeep[], t_blocka* block, const char* name)
+static gmx::ListOfLists<int>
+reduce_blocka(const int invindex[], const gmx_bool bKeep[], const t_blocka& block, const char* name)
 {
-    int *index, *a;
-    int  i, j, k, newi, newj;
+    gmx::ListOfLists<int> lists;
 
-    snew(index, block->nr);
-    snew(a, block->nra);
-
-    newi = newj = 0;
-    for (i = 0; (i < block->nr); i++)
+    std::vector<int> exclusionsForAtom;
+    for (int i = 0; i < block.nr; i++)
     {
-        for (j = block->index[i]; (j < block->index[i + 1]); j++)
+        if (bKeep[i])
         {
-            k = block->a[j];
-            if (bKeep[k])
+            exclusionsForAtom.clear();
+            for (int j = block.index[i]; j < block.index[i + 1]; j++)
             {
-                a[newj] = invindex[k];
-                newj++;
+                const int k = block.a[j];
+                if (bKeep[k])
+                {
+                    exclusionsForAtom.push_back(invindex[k]);
+                }
             }
-        }
-        if (newj > index[newi])
-        {
-            newi++;
-            index[newi] = newj;
+            lists.pushBack(exclusionsForAtom);
         }
     }
 
-    fprintf(stderr, "Reduced block %8s from %6d to %6d index-, %6d to %6d a-entries\n", name,
-            block->nr, newi, block->nra, newj);
-    block->index = index;
-    block->a     = a;
-    block->nr    = newi;
-    block->nra   = newj;
+    fprintf(stderr, "Reduced block %8s from %6d to %6zu index-, %6d to %6d a-entries\n", name,
+            block.nr, lists.size(), block.nra, lists.numElements());
+
+    return lists;
 }
 
 static void reduce_rvec(int gnx, const int index[], rvec vv[])
@@ -271,7 +265,7 @@ static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[],
     invindex = invind(gnx, top.atoms.nr, index);
 
     reduce_block(bKeep, &(top.mols), "mols");
-    reduce_blocka(invindex, bKeep, &(top.excls), "excls");
+    gmx::ListOfLists<int> excls = reduce_blocka(invindex, bKeep, top.excls, "excls");
     reduce_rvec(gnx, index, x);
     reduce_rvec(gnx, index, v);
     reduce_atom(gnx, index, top.atoms.atom, top.atoms.atomname, &(top.atoms.nres), top.atoms.resinfo);
@@ -297,7 +291,7 @@ static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[],
         }
     }
     mtop->moltype[0].atoms = top.atoms;
-    mtop->moltype[0].excls = top.excls;
+    mtop->moltype[0].excls = excls;
 
     mtop->molblock.resize(1);
     mtop->molblock[0].type = 0;
index 7c1cc84357f905775b0f3f38749d84ea473567a1..f0a91414fe2bc066042e865bd57792ca96ecae93 100644 (file)
@@ -42,6 +42,7 @@
 
 #include <algorithm>
 
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/txtdump.h"
 
@@ -218,6 +219,19 @@ static int pr_blocka_title(FILE* fp, int indent, const char* title, const t_bloc
     return indent;
 }
 
+static int pr_listoflists_title(FILE* fp, int indent, const char* title, const gmx::ListOfLists<int>* lists)
+{
+    if (available(fp, lists, indent, title))
+    {
+        indent = pr_title(fp, indent, title);
+        pr_indent(fp, indent);
+        fprintf(fp, "numLists=%zu\n", lists->size());
+        pr_indent(fp, indent);
+        fprintf(fp, "numElements=%d\n", lists->numElements());
+    }
+    return indent;
+}
+
 static void low_pr_blocka(FILE* fp, int indent, const char* title, const t_blocka* block, gmx_bool bShowNumbers)
 {
     int i;
@@ -325,6 +339,43 @@ void pr_blocka(FILE* fp, int indent, const char* title, const t_blocka* block, g
     }
 }
 
+void pr_listoflists(FILE* fp, int indent, const char* title, const gmx::ListOfLists<int>* lists, gmx_bool bShowNumbers)
+{
+    if (available(fp, lists, indent, title))
+    {
+        indent = pr_listoflists_title(fp, indent, title, lists);
+        for (gmx::index i = 0; i < lists->ssize(); i++)
+        {
+            int                      size = pr_indent(fp, indent);
+            gmx::ArrayRef<const int> list = (*lists)[i];
+            if (list.empty())
+            {
+                size += fprintf(fp, "%s[%d]={", title, int(i));
+            }
+            else
+            {
+                size += fprintf(fp, "%s[%d][num=%zu]={", title, bShowNumbers ? int(i) : -1, list.size());
+            }
+            bool isFirst = true;
+            for (const int j : list)
+            {
+                if (!isFirst)
+                {
+                    size += fprintf(fp, ", ");
+                }
+                if ((size) > (USE_WIDTH))
+                {
+                    fprintf(fp, "\n");
+                    size = pr_indent(fp, indent + INDENT);
+                }
+                size += fprintf(fp, "%d", j);
+                isFirst = false;
+            }
+            fprintf(fp, "}\n");
+        }
+    }
+}
+
 void copy_block(const t_block* src, t_block* dst)
 {
     dst->nr = src->nr;
index 8db0164b90c8e6ef9ca008c3666b98b62c1ac3df..fe0a30a56a15661f993bb13d56e711c90b85925a 100644 (file)
@@ -48,6 +48,9 @@
 namespace gmx
 {
 
+template<typename>
+class ListOfLists;
+
 /*! \brief Division of a range of indices into consecutive blocks
  *
  * A range of consecutive indices 0 to full.range.end() is divided
@@ -217,5 +220,6 @@ void stupid_fill_blocka(t_blocka* grp, int natom);
 
 void pr_block(FILE* fp, int indent, const char* title, const t_block* block, gmx_bool bShowNumbers);
 void pr_blocka(FILE* fp, int indent, const char* title, const t_blocka* block, gmx_bool bShowNumbers);
+void pr_listoflists(FILE* fp, int indent, const char* title, const gmx::ListOfLists<int>* block, gmx_bool bShowNumbers);
 
 #endif
index 9472c43781e67e9c831f26ca6521ca807ae0a00f..b73b68097e676ee50a04dc2ac816c123d8a1ca0f 100644 (file)
 
 #include "gromacs/topology/block.h"
 #include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/stringutil.h"
 
 namespace gmx
 {
 
+namespace
+{
+
+//! Converts ListOfLists to a list of ExclusionBlocks
+void listOfListsToExclusionBlocks(const ListOfLists<int>& b, gmx::ArrayRef<ExclusionBlock> b2)
+{
+    for (gmx::index i = 0; i < b.ssize(); i++)
+    {
+        for (int jAtom : b[i])
+        {
+            b2[i].atomNumber.push_back(jAtom);
+        }
+    }
+}
+
+//! Converts a list of ExclusionBlocks to ListOfLists
+void exclusionBlocksToListOfLists(gmx::ArrayRef<const ExclusionBlock> b2, ListOfLists<int>* b)
+{
+    b->clear();
+
+    for (const auto& block : b2)
+    {
+        b->pushBack(block.atomNumber);
+    }
+}
+
+} // namespace
+
 void blockaToExclusionBlocks(const t_blocka* b, gmx::ArrayRef<ExclusionBlock> b2)
 {
     for (int i = 0; (i < b->nr); i++)
@@ -79,21 +108,12 @@ void exclusionBlocksToBlocka(gmx::ArrayRef<const ExclusionBlock> b2, t_blocka* b
     b->index[i] = nra;
 }
 
-void mergeExclusions(t_blocka* excl, gmx::ArrayRef<ExclusionBlock> b2)
+namespace
 {
-    if (b2.empty())
-    {
-        return;
-    }
-    GMX_RELEASE_ASSERT(b2.ssize() == excl->nr,
-                       "Cannot merge exclusions for "
-                       "blocks that do not describe the same number "
-                       "of particles");
 
-    /* Convert the t_blocka entries to ExclusionBlock form */
-    blockaToExclusionBlocks(excl, b2);
-
-    /* Count and sort the exclusions */
+//! Counts and sorts the exclusions
+int countAndSortExclusions(gmx::ArrayRef<ExclusionBlock> b2)
+{
     int nra = 0;
     for (auto& block : b2)
     {
@@ -118,10 +138,29 @@ void mergeExclusions(t_blocka* excl, gmx::ArrayRef<ExclusionBlock> b2)
             nra += block.nra();
         }
     }
-    excl->nra = nra;
-    srenew(excl->a, excl->nra);
 
-    exclusionBlocksToBlocka(b2, excl);
+    return nra;
+}
+
+} // namespace
+
+void mergeExclusions(ListOfLists<int>* excl, gmx::ArrayRef<ExclusionBlock> b2)
+{
+    if (b2.empty())
+    {
+        return;
+    }
+    GMX_RELEASE_ASSERT(b2.ssize() == excl->ssize(),
+                       "Cannot merge exclusions for "
+                       "blocks that do not describe the same number "
+                       "of particles");
+
+    /* Convert the t_blocka entries to ExclusionBlock form */
+    listOfListsToExclusionBlocks(*excl, b2);
+
+    countAndSortExclusions(b2);
+
+    exclusionBlocksToListOfLists(b2, excl);
 }
 
 } // namespace gmx
index b27db36093d71f0c40a945ea8562a80079d86a93..69c9c54c21047ce8bc8add7e2403c8df4c1a37d9 100644 (file)
@@ -44,6 +44,8 @@ struct t_blocka;
 
 namespace gmx
 {
+template<typename>
+class ListOfLists;
 
 /*! \libinternal \brief
  * Describes exclusions for a single atom.
@@ -61,7 +63,7 @@ struct ExclusionBlock
  * Requires that \c b2 and \c excl describe the same number of
  * particles, if \c b2 describes a non-zero number.
  */
-void mergeExclusions(t_blocka* excl, gmx::ArrayRef<ExclusionBlock> b2);
+void mergeExclusions(ListOfLists<int>* excl, gmx::ArrayRef<ExclusionBlock> b2);
 
 /*! \brief
  * Convert the exclusions.
index 4bc2e2946e229397becd7045b3bf00bcf5b9cd6b..541788aa602fa85b26c2ca7bfef057f0f1b10545 100644 (file)
@@ -650,39 +650,41 @@ t_atoms gmx_mtop_global_atoms(const gmx_mtop_t* mtop)
  * The cat routines below are old code from src/kernel/topcat.c
  */
 
-static void blockacat(t_blocka* dest, const t_blocka* src, int copies, int dnum, int snum)
+static void blockacat(t_blocka* dest, const gmx::ListOfLists<int>& src, int copies, int dnum, int snum)
 {
-    int i, j, l, size;
+    int j, l;
     int destnr  = dest->nr;
     int destnra = dest->nra;
 
-    if (src->nr)
+    if (!src.empty())
     {
-        size = (dest->nr + copies * src->nr + 1);
+        size_t size = (dest->nr + copies * src.size() + 1);
         srenew(dest->index, size);
     }
-    if (src->nra)
+    gmx::ArrayRef<const int> srcElements = src.elementsView();
+    if (!srcElements.empty())
     {
-        size = (dest->nra + copies * src->nra);
+        size_t size = (dest->nra + copies * srcElements.size());
         srenew(dest->a, size);
     }
 
+    gmx::ArrayRef<const int> srcListRanges = src.listRangesView();
     for (l = destnr, j = 0; (j < copies); j++)
     {
-        for (i = 0; (i < src->nr); i++)
+        for (gmx::index i = 0; i < src.ssize(); i++)
         {
-            dest->index[l++] = dest->nra + src->index[i];
+            dest->index[l++] = dest->nra + srcListRanges[i];
         }
-        dest->nra += src->nra;
+        dest->nra += src.ssize();
     }
     for (l = destnra, j = 0; (j < copies); j++)
     {
-        for (i = 0; (i < src->nra); i++)
+        for (const int srcElement : srcElements)
         {
-            dest->a[l++] = dnum + src->a[i];
+            dest->a[l++] = dnum + srcElement;
         }
         dnum += snum;
-        dest->nr += src->nr;
+        dest->nr += src.ssize();
     }
     dest->index[dest->nr] = dest->nra;
     dest->nalloc_index    = dest->nr;
@@ -941,6 +943,31 @@ static void copyAtomtypesFromMtop(const gmx_mtop_t& mtop, t_atomtypes* atomtypes
  * \param[in] mtop  Reference to input mtop.
  * \param[in] excls Pointer to final excls data structure.
  */
+static void copyExclsFromMtop(const gmx_mtop_t& mtop, gmx::ListOfLists<int>* excls)
+{
+    excls->clear();
+    int atomIndex = 0;
+    for (const gmx_molblock_t& molb : mtop.molblock)
+    {
+        const gmx_moltype_t& molt = mtop.moltype[molb.type];
+
+        for (int mol = 0; mol < molb.nmol; mol++)
+        {
+            excls->appendListOfLists(molt.excls, atomIndex);
+
+            atomIndex += molt.atoms.nr;
+        }
+    }
+}
+
+/*! \brief Copy excls from mtop.
+ *
+ * Makes a deep copy of excls(t_blocka) from gmx_mtop_t.
+ * Used to initialize legacy t_topology.
+ *
+ * \param[in] mtop  Reference to input mtop.
+ * \param[in] excls Pointer to final excls data structure.
+ */
 static void copyExclsFromMtop(const gmx_mtop_t& mtop, t_blocka* excls)
 {
     init_blocka(excls);
@@ -952,7 +979,7 @@ static void copyExclsFromMtop(const gmx_mtop_t& mtop, t_blocka* excls)
         int srcnr  = molt.atoms.nr;
         int destnr = natoms;
 
-        blockacat(excls, &molt.excls, molb.nmol, destnr, srcnr);
+        blockacat(excls, molt.excls, molb.nmol, destnr, srcnr);
 
         natoms += molb.nmol * srcnr;
     }
@@ -966,21 +993,21 @@ static void copyExclsFromMtop(const gmx_mtop_t& mtop, t_blocka* excls)
  * \param[inout]    excls   existing exclusions in local topology
  * \param[in]       ids     list of global IDs of atoms
  */
-static void addMimicExclusions(t_blocka* excls, const gmx::ArrayRef<const int> ids)
+static void addMimicExclusions(gmx::ListOfLists<int>* excls, const gmx::ArrayRef<const int> ids)
 {
     t_blocka inter_excl{};
     init_blocka(&inter_excl);
     size_t n_q = ids.size();
 
-    inter_excl.nr  = excls->nr;
+    inter_excl.nr  = excls->ssize();
     inter_excl.nra = n_q * n_q;
 
     size_t total_nra = n_q * n_q;
 
-    snew(inter_excl.index, excls->nr + 1);
+    snew(inter_excl.index, excls->ssize() + 1);
     snew(inter_excl.a, total_nra);
 
-    for (int i = 0; i < excls->nr; ++i)
+    for (int i = 0; i < inter_excl.nr; ++i)
     {
         inter_excl.index[i] = 0;
     }
@@ -1012,7 +1039,7 @@ static void addMimicExclusions(t_blocka* excls, const gmx::ArrayRef<const int> i
 
     inter_excl.index[inter_excl.nr] = n_q * n_q;
 
-    std::vector<gmx::ExclusionBlock> qmexcl2(excls->nr);
+    std::vector<gmx::ExclusionBlock> qmexcl2(excls->size());
     gmx::blockaToExclusionBlocks(&inter_excl, qmexcl2);
 
     // Merge the created exclusion list with the existing one
index ab9f1ad1227d98909f09b37739599105213112c2..53a2a05fbcb04605c034338b9e1eabcb09f3db09 100644 (file)
@@ -47,6 +47,7 @@
 
 #include "gromacs/topology/block.h"
 #include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/smalloc.h"
 
 #include "testutils/cmdlinetest.h"
@@ -101,6 +102,21 @@ void makeTestBlockAData(t_blocka* ba)
     addGroupToBlocka(ba, indices);
 }
 
+//! Return ListOfLists filled with some datastructures
+ListOfLists<int> makeTestListOfLists()
+{
+    ListOfLists<int> list;
+
+    std::vector<int> indices = { 12, 11, 9, 6, 2 };
+    list.pushBack(indices);
+    indices = { 10, 8, 5, 1 };
+    list.pushBack(indices);
+    indices = { 7, 4, 0 };
+    list.pushBack(indices);
+
+    return list;
+}
+
 class ExclusionBlockTest : public ::testing::Test
 {
 public:
@@ -108,6 +124,7 @@ public:
     {
         const int natom = 3;
         makeTestBlockAData(&ba_);
+        list_ = makeTestListOfLists();
         b_.resize(natom);
     }
     ~ExclusionBlockTest() override { done_blocka(&ba_); }
@@ -126,8 +143,23 @@ public:
         }
     }
 
+    void compareBlocksAndList()
+    {
+        GMX_RELEASE_ASSERT(ssize(b_) == list_.ssize(), "The list counts should match");
+        for (index i = 0; i < ssize(b_); i++)
+        {
+            gmx::ArrayRef<const int> jList = list_[i];
+            EXPECT_EQ(b_[i].nra(), jList.ssize()) << "Block size mismatch at " << i << ".";
+            for (int j = 0; j < b_[i].nra(); j++)
+            {
+                EXPECT_EQ(b_[i].atomNumber[j], jList[j]) << "Block mismatch at " << i << " , " << j << ".";
+            }
+        }
+    }
+
 protected:
     t_blocka                    ba_;
+    ListOfLists<int>            list_;
     std::vector<ExclusionBlock> b_;
 };
 
@@ -148,8 +180,8 @@ TEST_F(ExclusionBlockTest, ConvertExclusionBlockToBlocka)
 
 TEST_F(ExclusionBlockTest, MergeExclusions)
 {
-    mergeExclusions(&ba_, b_);
-    compareBlocks();
+    mergeExclusions(&list_, b_);
+    compareBlocksAndList();
 }
 
 } // namespace
index e01993300bb2ab27429e88d045912376cb8f1816..3a1b82eb502faccf6bb25a456a31a7b2c8f33ec1 100644 (file)
@@ -76,7 +76,7 @@ void init_top(t_topology* top)
 }
 
 
-gmx_moltype_t::gmx_moltype_t() : name(nullptr), excls()
+gmx_moltype_t::gmx_moltype_t() : name(nullptr)
 {
     init_t_atoms(&atoms, 0, FALSE);
 }
@@ -84,7 +84,6 @@ gmx_moltype_t::gmx_moltype_t() : name(nullptr), excls()
 gmx_moltype_t::~gmx_moltype_t()
 {
     done_atom(&atoms);
-    done_blocka(&excls);
 }
 
 gmx_mtop_t::gmx_mtop_t()
@@ -136,7 +135,6 @@ void done_top_mtop(t_topology* top, gmx_mtop_t* mtop)
 
 gmx_localtop_t::gmx_localtop_t()
 {
-    init_blocka_null(&excls);
     init_idef(&idef);
     init_atomtypes(&atomtypes);
 }
@@ -146,7 +144,6 @@ gmx_localtop_t::~gmx_localtop_t()
     if (!useInDomainDecomp_)
     {
         done_idef(&idef);
-        done_blocka(&excls);
         done_atomtypes(&atomtypes);
     }
 }
@@ -287,7 +284,7 @@ static void pr_moltype(FILE*                 fp,
     pr_indent(fp, indent);
     fprintf(fp, "name=\"%s\"\n", *(molt->name));
     pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
-    pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers);
+    pr_listoflists(fp, indent, "excls", &molt->excls, bShowNumbers);
     for (j = 0; (j < F_NRE); j++)
     {
         pr_ilist(fp, indent, interaction_function[j].longname, ffparams->functype.data(),
@@ -457,15 +454,18 @@ static void cmp_cmap(FILE* fp, const gmx_cmap_t* cmap1, const gmx_cmap_t* cmap2,
     }
 }
 
-static void cmp_blocka(FILE* fp, const t_blocka* b1, const t_blocka* b2, const char* s)
+static void cmp_listoflists(FILE*                        fp,
+                            const gmx::ListOfLists<int>& list1,
+                            const gmx::ListOfLists<int>& list2,
+                            const char*                  s)
 {
     char buf[32];
 
     fprintf(fp, "comparing blocka %s\n", s);
-    sprintf(buf, "%s.nr", s);
-    cmp_int(fp, buf, -1, b1->nr, b2->nr);
-    sprintf(buf, "%s.nra", s);
-    cmp_int(fp, buf, -1, b1->nra, b2->nra);
+    sprintf(buf, "%s.numLists", s);
+    cmp_int(fp, buf, -1, list1.ssize(), list2.ssize());
+    sprintf(buf, "%s.numElements", s);
+    cmp_int(fp, buf, -1, list1.numElements(), list2.numElements());
 }
 
 static void compareFfparams(FILE*                 fp,
@@ -534,7 +534,7 @@ static void compareMoltypes(FILE*                              fp,
         compareAtoms(fp, &mt1[i].atoms, &mt2[i].atoms, relativeTolerance, absoluteTolerance);
         compareInteractionLists(fp, &mt1[i].ilist, &mt2[i].ilist);
         std::string buf = gmx::formatString("excls[%d]", i);
-        cmp_blocka(fp, &mt1[i].excls, &mt2[i].excls, buf.c_str());
+        cmp_listoflists(fp, mt1[i].excls, mt2[i].excls, buf.c_str());
     }
 }
 
@@ -674,8 +674,8 @@ int getGroupType(const SimulationGroups& group, SimulationAtomGroupType type, in
 
 void copy_moltype(const gmx_moltype_t* src, gmx_moltype_t* dst)
 {
-    dst->name = src->name;
-    copy_blocka(&src->excls, &dst->excls);
+    dst->name          = src->name;
+    dst->excls         = src->excls;
     t_atoms* atomsCopy = copy_t_atoms(&src->atoms);
     dst->atoms         = *atomsCopy;
     sfree(atomsCopy);
index 710bbbd2d8ecbf6e1123d95521198a5e07d99cc8..0e5c0d9e17e099b4b137f9c270e8450ae7dc3a8e 100644 (file)
@@ -48,6 +48,7 @@
 #include "gromacs/topology/idef.h"
 #include "gromacs/topology/symtab.h"
 #include "gromacs/utility/enumerationhelpers.h"
+#include "gromacs/utility/listoflists.h"
 #include "gromacs/utility/unique_cptr.h"
 
 enum class SimulationAtomGroupType : int
@@ -83,10 +84,10 @@ struct gmx_moltype_t
     /*! \brief Default copy constructor */
     gmx_moltype_t(const gmx_moltype_t&) = default;
 
-    char**           name;  /**< Name of the molecule type            */
-    t_atoms          atoms; /**< The atoms in this molecule           */
-    InteractionLists ilist; /**< Interaction list with local indices  */
-    t_blocka         excls; /**< The exclusions                       */
+    char**                name;  /**< Name of the molecule type            */
+    t_atoms               atoms; /**< The atoms in this molecule           */
+    InteractionLists      ilist; /**< Interaction list with local indices  */
+    gmx::ListOfLists<int> excls; /**< The exclusions                       */
 };
 
 /*! \brief Block of molecules of the same type, used in gmx_mtop_t */
@@ -214,7 +215,7 @@ struct gmx_localtop_t
     //! Atomtype properties
     t_atomtypes atomtypes;
     //! The exclusions
-    t_blocka excls;
+    gmx::ListOfLists<int> excls;
     //! Flag for domain decomposition so we don't free already freed memory.
     bool useInDomainDecomp_ = false;
 };
index cc45dadebcee5dcbef639e6e4dce08c45b02bd10..f9a58d2f2f0e3c78881cfd37d685bd29bc236350 100644 (file)
@@ -386,7 +386,7 @@ void Rdf::initAnalysis(const TrajectoryAnalysisSettings& settings, const Topolog
             }
         }
         localTop_ = top.expandedTopology();
-        if (localTop_->excls.nr == 0)
+        if (localTop_->excls.empty())
         {
             GMX_THROW(InconsistentInputError(
                     "-excl is set, but the file provided to -s does not define exclusions"));