Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / selection / indexutil.cpp
index 5436f970a472c2f5a4e36174d026f2064a2cab24..a38fe3620509088df357a6f44458da298d89c430 100644 (file)
@@ -1,32 +1,36 @@
 /*
+ * This file is part of the GROMACS molecular simulation package.
  *
- *                This source code is part of
+ * Copyright (c) 2009,2010,2011,2012,2013,2014, 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.
  *
- *                 G   R   O   M   A   C   S
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
  *
- *          GROningen MAchine for Chemical Simulations
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
  *
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2009, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
  *
  * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
+ * the research papers on the package. Check out http://www.gromacs.org.
  */
 /*! \internal \file
  * \brief
  * \author Teemu Murtola <teemu.murtola@gmail.com>
  * \ingroup module_selection
  */
-#include "gromacs/legacyheaders/index.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
-#include "gromacs/legacyheaders/smalloc.h"
-#include "gromacs/legacyheaders/string2.h"
-#include "gromacs/legacyheaders/typedefs.h"
+#include "gmxpre.h"
+
+#include "indexutil.h"
 
-#include "gromacs/selection/indexutil.h"
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <string>
+#include <vector>
+
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
 
 /********************************************************************
  * gmx_ana_indexgrps_t functions
  */
 struct gmx_ana_indexgrps_t
 {
+    //! Initializes an empty set of groups.
+    explicit gmx_ana_indexgrps_t(int nr) : nr(nr), g(NULL)
+    {
+        names.reserve(nr);
+        snew(g, nr);
+    }
+    ~gmx_ana_indexgrps_t()
+    {
+        for (int i = 0; i < nr; ++i)
+        {
+            gmx_ana_index_deinit(&g[i]);
+        }
+        sfree(g);
+    }
+
     /** Number of index groups. */
-    int                 nr;
+    int                       nr;
     /** Array of index groups. */
-    gmx_ana_index_t    *g;
+    gmx_ana_index_t          *g;
+    /** Group names. */
+    std::vector<std::string>  names;
 };
 
-/*!
- * \param[out] g     Index group structure.
- * \param[in]  ngrps Number of groups for which memory is allocated.
- */
-void
-gmx_ana_indexgrps_alloc(gmx_ana_indexgrps_t **g, int ngrps)
-{
-    snew(*g, 1);
-    (*g)->nr = ngrps;
-    snew((*g)->g,    ngrps);
-}
-
 /*!
  * \param[out] g     Index group structure.
  * \param[in]  top   Topology structure.
@@ -89,7 +107,6 @@ gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
 {
     t_blocka *block = NULL;
     char    **names = NULL;
-    int       i, j;
 
     if (fnm)
     {
@@ -102,30 +119,45 @@ gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
     }
     else
     {
-        snew(*g, 1);
-        (*g)->nr    = 0;
-        (*g)->g     = NULL;
+        *g = new gmx_ana_indexgrps_t(0);
         return;
     }
 
-    gmx_ana_indexgrps_alloc(g, block->nr);
-    for (i = 0; i < block->nr; ++i)
+    try
     {
-        gmx_ana_index_t *grp = &(*g)->g[i];
+        *g = new gmx_ana_indexgrps_t(block->nr);
+        for (int i = 0; i < block->nr; ++i)
+        {
+            gmx_ana_index_t *grp = &(*g)->g[i];
 
-        grp->isize = block->index[i+1] - block->index[i];
-        snew(grp->index, grp->isize);
-        for (j = 0; j < grp->isize; ++j)
+            grp->isize = block->index[i+1] - block->index[i];
+            snew(grp->index, grp->isize);
+            for (int j = 0; j < grp->isize; ++j)
+            {
+                grp->index[j] = block->a[block->index[i]+j];
+            }
+            grp->nalloc_index = grp->isize;
+            (*g)->names.push_back(names[i]);
+        }
+    }
+    catch (...)
+    {
+        for (int i = 0; i < block->nr; ++i)
         {
-            grp->index[j] = block->a[block->index[i]+j];
+            sfree(names[i]);
         }
-        grp->name         = names[i];
-        grp->nalloc_index = grp->isize;
+        sfree(names);
+        done_blocka(block);
+        sfree(block);
+        throw;
     }
-
+    for (int i = 0; i < block->nr; ++i)
+    {
+        sfree(names[i]);
+    }
+    sfree(names);
     done_blocka(block);
     sfree(block);
-    sfree(names);
 }
 
 /*!
@@ -136,39 +168,7 @@ gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
 void
 gmx_ana_indexgrps_free(gmx_ana_indexgrps_t *g)
 {
-    int  i;
-
-    if (g->nr == 0)
-    {
-        sfree(g);
-        return;
-    }
-    for (i = 0; i < g->nr; ++i)
-    {
-        gmx_ana_index_deinit(&g->g[i]);
-    }
-    sfree(g->g);
-    g->nr    = 0;
-    g->g     = NULL;
-    sfree(g);
-}
-
-/*!
- * \param[out] dest Destination index groups.
- * \param[in]  src  Source index groups.
- *
- * A deep copy is made for all fields, including the group names.
- */
-void
-gmx_ana_indexgrps_clone(gmx_ana_indexgrps_t **dest, gmx_ana_indexgrps_t *src)
-{
-    int g;
-
-    gmx_ana_indexgrps_alloc(dest, src->nr);
-    for (g = 0; g < src->nr; ++g)
-    {
-        gmx_ana_index_copy(&(*dest)->g[g], &src->g[g], true);
-    }
+    delete g;
 }
 
 /*!
@@ -199,53 +199,64 @@ gmx_ana_indexgrps_get_grp(gmx_ana_indexgrps_t *g, int n)
 }
 
 /*!
- * \param[out] dest Output structure.
- * \param[in]  src  Input index groups.
- * \param[in]  n    Number of the group to extract.
+ * \param[out] dest     Output structure.
+ * \param[out] destName Receives the name of the group if found.
+ * \param[in]  src      Input index groups.
+ * \param[in]  n        Number of the group to extract.
  * \returns true if \p n is a valid group in \p src, false otherwise.
  */
 bool
-gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, int n)
+gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, std::string *destName,
+                          gmx_ana_indexgrps_t *src, int n)
 {
+    destName->clear();
     if (n < 0 || n >= src->nr)
     {
         dest->isize = 0;
         return false;
     }
 
+    if (destName != NULL)
+    {
+        *destName = src->names[n];
+    }
     gmx_ana_index_copy(dest, &src->g[n], true);
     return true;
 }
 
 /*!
- * \param[out] dest Output structure.
- * \param[in]  src  Input index groups.
- * \param[in]  name Name (or part of the name) of the group to extract.
+ * \param[out] dest     Output structure.
+ * \param[out] destName Receives the name of the group if found.
+ * \param[in]  src      Input index groups.
+ * \param[in]  name     Name (or part of the name) of the group to extract.
  * \returns true if \p name is a valid group in \p src, false otherwise.
  *
  * Uses the Gromacs routine find_group() to find the actual group;
  * the comparison is case-insensitive.
  */
 bool
-gmx_ana_indexgrps_find(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, char *name)
+gmx_ana_indexgrps_find(gmx_ana_index_t *dest, std::string *destName,
+                       gmx_ana_indexgrps_t *src,
+                       const char *name)
 {
-    int    i;
-    char **names;
+    const char **names;
 
+    destName->clear();
     snew(names, src->nr);
-    for (i = 0; i < src->nr; ++i)
+    for (int i = 0; i < src->nr; ++i)
     {
-        names[i] = src->g[i].name;
+        names[i] = src->names[i].c_str();
     }
-    i = find_group(name, src->nr, names);
+    int n = find_group(const_cast<char *>(name), src->nr,
+                       const_cast<char **>(names));
     sfree(names);
-    if (i == NOTSET)
+    if (n < 0)
     {
         dest->isize = 0;
         return false;
     }
 
-    return gmx_ana_indexgrps_extract(dest, src, i);
+    return gmx_ana_indexgrps_extract(dest, destName, src, n);
 }
 
 /*!
@@ -257,12 +268,10 @@ gmx_ana_indexgrps_find(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, char *na
 void
 gmx_ana_indexgrps_print(FILE *fp, gmx_ana_indexgrps_t *g, int maxn)
 {
-    int  i;
-
-    for (i = 0; i < g->nr; ++i)
+    for (int i = 0; i < g->nr; ++i)
     {
-        fprintf(fp, " %2d: ", i);
-        gmx_ana_index_dump(fp, &g->g[i], i, maxn);
+        fprintf(fp, " Group %2d \"%s\" ", i, g->names[i].c_str());
+        gmx_ana_index_dump(fp, &g->g[i], maxn);
     }
 }
 
@@ -307,7 +316,6 @@ gmx_ana_index_clear(gmx_ana_index_t *g)
 {
     g->isize        = 0;
     g->index        = NULL;
-    g->name         = NULL;
     g->nalloc_index = 0;
 }
 
@@ -315,29 +323,25 @@ gmx_ana_index_clear(gmx_ana_index_t *g)
  * \param[out] g      Output structure.
  * \param[in]  isize  Number of atoms in the new group.
  * \param[in]  index  Array of \p isize atoms (can be NULL if \p isize is 0).
- * \param[in]  name   Name for the new group (can be NULL).
  * \param[in]  nalloc Number of elements allocated for \p index
  *   (if 0, \p index is not freed in gmx_ana_index_deinit())
  *
  * No copy if \p index is made.
  */
 void
-gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, char *name,
-                  int nalloc)
+gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, int nalloc)
 {
     g->isize        = isize;
     g->index        = index;
-    g->name         = name;
     g->nalloc_index = nalloc;
 }
 
 /*!
  * \param[out] g      Output structure.
  * \param[in]  natoms Number of atoms.
- * \param[in]  name   Name for the new group (can be NULL).
  */
 void
-gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms, char *name)
+gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms)
 {
     int  i;
 
@@ -347,7 +351,6 @@ gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms, char *name)
     {
         g->index[i] = i;
     }
-    g->name         = name;
     g->nalloc_index = natoms;
 }
 
@@ -363,7 +366,6 @@ gmx_ana_index_deinit(gmx_ana_index_t *g)
     {
         sfree(g->index);
     }
-    sfree(g->name);
     gmx_ana_index_clear(g);
 }
 
@@ -386,38 +388,21 @@ gmx_ana_index_copy(gmx_ana_index_t *dest, gmx_ana_index_t *src, bool bAlloc)
             snew(dest->index, dest->isize);
             dest->nalloc_index = dest->isize;
         }
-        memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
-    }
-    if (bAlloc && src->name)
-    {
-        dest->name = strdup(src->name);
-    }
-    else if (bAlloc || src->name)
-    {
-        dest->name = src->name;
+        std::memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
     }
 }
 
 /*!
  * \param[in]  fp     Where to print the output.
  * \param[in]  g      Index group to print.
- * \param[in]  i      Group number to use if the name is NULL.
  * \param[in]  maxn   Maximum number of indices to print (-1 = print all).
  */
 void
-gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int i, int maxn)
+gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int maxn)
 {
     int  j, n;
 
-    if (g->name)
-    {
-        fprintf(fp, "\"%s\"", g->name);
-    }
-    else
-    {
-        fprintf(fp, "Group %d", i+1);
-    }
-    fprintf(fp, " (%d atoms)", g->isize);
+    fprintf(fp, "(%d atoms)", g->isize);
     if (maxn != 0)
     {
         fprintf(fp, ":");
@@ -438,32 +423,16 @@ gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int i, int maxn)
     fprintf(fp, "\n");
 }
 
-/*!
- * \param[in]  g      Input index group.
- * \param[in]  natoms Number of atoms to check against.
- *
- * If any atom index in the index group is less than zero or >= \p natoms,
- * gmx_fatal() is called.
- */
-void
-gmx_ana_index_check(gmx_ana_index_t *g, int natoms)
+int
+gmx_ana_index_get_max_index(gmx_ana_index_t *g)
 {
-    int  j;
-
-    for (j = 0; j < g->isize; ++j)
+    if (g->isize == 0)
     {
-        if (g->index[j] >= natoms)
-        {
-            gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) "
-                      "larger than number of atoms in trajectory (%d atoms)",
-                      g->index[j], g->name, g->isize, natoms);
-        }
-        else if (g->index[j] < 0)
-        {
-            gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) "
-                      "is less than zero",
-                      g->index[j], g->name, g->isize);
-        }
+        return 0;
+    }
+    else
+    {
+        return *std::max_element(g->index, g->index + g->isize);
     }
 }
 
@@ -487,6 +456,19 @@ gmx_ana_index_check_sorted(gmx_ana_index_t *g)
     return true;
 }
 
+bool
+gmx_ana_index_check_range(gmx_ana_index_t *g, int natoms)
+{
+    for (int i = 0; i < g->isize; ++i)
+    {
+        if (g->index[i] < 0 || g->index[i] >= natoms)
+        {
+            return false;
+        }
+    }
+    return true;
+}
+
 /********************************************************************
  * Set operations
  ********************************************************************/
@@ -512,7 +494,7 @@ cmp_atomid(const void *a, const void *b)
 void
 gmx_ana_index_sort(gmx_ana_index_t *g)
 {
-    qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
+    std::qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
 }
 
 /*!
@@ -823,7 +805,7 @@ gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
             srenew(t->a, g->isize);
             t->nalloc_a = g->isize;
         }
-        memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
+        std::memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
     }
 
     /* Allocate memory for the block index. We don't know in advance
@@ -892,7 +874,7 @@ gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
                         break;
 
                     default: /* Should not be reached */
-                        gmx_bug("internal error");
+                        GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
                         break;
                 }
             }
@@ -1005,6 +987,8 @@ gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b)
  * \returns   true if \p g consists of one or more complete elements of type
  *   \p type, false otherwise.
  *
+ * \p g is assumed to be sorted, otherwise may return false negatives.
+ *
  * If \p type is \ref INDEX_ATOM, the return value is always true.
  * If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
  * always false.
@@ -1013,6 +997,7 @@ bool
 gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
                                  t_topology *top)
 {
+    // TODO: Consider whether unsorted groups need to be supported better.
     switch (type)
     {
         case INDEX_UNKNOWN:
@@ -1069,12 +1054,14 @@ void
 gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
 {
     m->type              = INDEX_UNKNOWN;
-    m->nr                = 0;
     m->refid             = NULL;
     m->mapid             = NULL;
     m->mapb.nr           = 0;
     m->mapb.index        = NULL;
     m->mapb.nalloc_index = 0;
+    m->mapb.nra          = 0;
+    m->mapb.a            = NULL;
+    m->mapb.nalloc_a     = 0;
     m->orgid             = NULL;
     m->b.nr              = 0;
     m->b.index           = NULL;
@@ -1083,7 +1070,6 @@ gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
     m->b.nalloc_index    = 0;
     m->b.nalloc_a        = 0;
     m->bStatic           = true;
-    m->bMapStatic        = true;
 }
 
 /*!
@@ -1134,8 +1120,7 @@ gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     m->type   = type;
     gmx_ana_index_make_block(&m->b, top, g, type, false);
     gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
-    m->nr = m->b.nr;
-    for (i = mi = 0; i < m->nr; ++i)
+    for (i = mi = 0; i < m->b.nr; ++i)
     {
         ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
         switch (type)
@@ -1159,15 +1144,16 @@ gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                 break;
         }
     }
-    for (i = 0; i < m->nr; ++i)
+    for (i = 0; i < m->b.nr; ++i)
     {
         m->refid[i] = i;
         m->mapid[i] = m->orgid[i];
     }
-    m->mapb.nr = m->nr;
-    memcpy(m->mapb.index, m->b.index, (m->nr+1)*sizeof(*(m->mapb.index)));
-    m->bStatic    = true;
-    m->bMapStatic = true;
+    m->mapb.nr  = m->b.nr;
+    m->mapb.nra = m->b.nra;
+    m->mapb.a   = m->b.a;
+    std::memcpy(m->mapb.index, m->b.index, (m->b.nr+1)*sizeof(*(m->mapb.index)));
+    m->bStatic  = true;
 }
 
 /*!
@@ -1189,16 +1175,18 @@ void
 gmx_ana_indexmap_set_static(gmx_ana_indexmap_t *m, t_blocka *b)
 {
     sfree(m->mapid);
-    m->mapid = m->orgid;
-    sfree(m->b.index);
-    m->b.nalloc_index = 0;
-    m->b.index        = b->index;
     sfree(m->mapb.index);
-    m->mapb.nalloc_index = 0;
-    m->mapb.index        = m->b.index;
+    sfree(m->b.index);
     sfree(m->b.a);
-    m->b.nalloc_a = 0;
-    m->b.a        = b->a;
+    m->mapb.nalloc_index = 0;
+    m->mapb.nalloc_a     = 0;
+    m->b.nalloc_index    = 0;
+    m->b.nalloc_a        = 0;
+    m->mapid             = m->orgid;
+    m->mapb.index        = b->index;
+    m->mapb.a            = b->a;
+    m->b.index           = b->index;
+    m->b.a               = b->a;
 }
 
 /*!
@@ -1219,17 +1207,53 @@ gmx_ana_indexmap_copy(gmx_ana_indexmap_t *dest, gmx_ana_indexmap_t *src, bool bF
         dest->type       = src->type;
         dest->b.nr       = src->b.nr;
         dest->b.nra      = src->b.nra;
-        memcpy(dest->orgid,      src->orgid,      dest->b.nr*sizeof(*dest->orgid));
-        memcpy(dest->b.index,    src->b.index,   (dest->b.nr+1)*sizeof(*dest->b.index));
-        memcpy(dest->b.a,        src->b.a,        dest->b.nra*sizeof(*dest->b.a));
+        std::memcpy(dest->orgid,      src->orgid,      dest->b.nr*sizeof(*dest->orgid));
+        std::memcpy(dest->b.index,    src->b.index,   (dest->b.nr+1)*sizeof(*dest->b.index));
+        std::memcpy(dest->b.a,        src->b.a,        dest->b.nra*sizeof(*dest->b.a));
     }
-    dest->nr         = src->nr;
     dest->mapb.nr    = src->mapb.nr;
-    memcpy(dest->refid,      src->refid,      dest->nr*sizeof(*dest->refid));
-    memcpy(dest->mapid,      src->mapid,      dest->nr*sizeof(*dest->mapid));
-    memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
-    dest->bStatic    = src->bStatic;
-    dest->bMapStatic = src->bMapStatic;
+    dest->mapb.nra   = src->mapb.nra;
+    if (src->mapb.nalloc_a > 0)
+    {
+        if (bFirst)
+        {
+            snew(dest->mapb.a, src->mapb.nalloc_a);
+            dest->mapb.nalloc_a = src->mapb.nalloc_a;
+        }
+        std::memcpy(dest->mapb.a, src->mapb.a, dest->mapb.nra*sizeof(*dest->mapb.a));
+    }
+    else
+    {
+        dest->mapb.a = src->mapb.a;
+    }
+    std::memcpy(dest->refid,      src->refid,      dest->mapb.nr*sizeof(*dest->refid));
+    std::memcpy(dest->mapid,      src->mapid,      dest->mapb.nr*sizeof(*dest->mapid));
+    std::memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
+    dest->bStatic = src->bStatic;
+}
+
+/*! \brief
+ * Helper function to set the source atoms in an index map.
+ *
+ * \param[in,out] m     Mapping structure.
+ * \param[in]     isize Number of atoms in the \p index array.
+ * \param[in]     index List of atoms.
+ */
+static void
+set_atoms(gmx_ana_indexmap_t *m, int isize, int *index)
+{
+    m->mapb.nra = isize;
+    if (m->mapb.nalloc_a == 0)
+    {
+        m->mapb.a = index;
+    }
+    else
+    {
+        for (int i = 0; i < isize; ++i)
+        {
+            m->mapb.a[i] = index[i];
+        }
+    }
 }
 
 /*!
@@ -1247,7 +1271,6 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                         bool bMaskOnly)
 {
     int  i, j, bi, bj;
-    bool bStatic;
 
     /* Process the simple cases first */
     if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
@@ -1256,6 +1279,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     }
     if (m->type == INDEX_ALL)
     {
+        set_atoms(m, g->isize, g->index);
         if (m->b.nr > 0)
         {
             m->mapb.index[1] = g->isize;
@@ -1263,8 +1287,9 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
         return;
     }
     /* Reset the reference IDs and mapping if necessary */
-    bStatic = (g->isize == m->b.nra && m->nr == m->b.nr);
-    if (bStatic || bMaskOnly)
+    const bool bToFull  = (g->isize == m->b.nra);
+    const bool bWasFull = (m->mapb.nra == m->b.nra);
+    if (bToFull || bMaskOnly)
     {
         if (!m->bStatic)
         {
@@ -1273,7 +1298,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                 m->refid[bj] = bj;
             }
         }
-        if (!m->bMapStatic)
+        if (!bWasFull)
         {
             for (bj = 0; bj < m->b.nr; ++bj)
             {
@@ -1283,11 +1308,12 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
             {
                 m->mapb.index[bj] = m->b.index[bj];
             }
-            m->bMapStatic = true;
         }
+        set_atoms(m, m->b.nra, m->b.a);
+        m->mapb.nr = m->b.nr;
     }
     /* Exit immediately if the group is static */
-    if (bStatic)
+    if (bToFull)
     {
         m->bStatic = true;
         return;
@@ -1295,7 +1321,6 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
 
     if (bMaskOnly)
     {
-        m->nr = m->b.nr;
         for (i = j = bj = 0; i < g->isize; ++i, ++j)
         {
             /* Find the next atom in the block */
@@ -1322,6 +1347,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     }
     else
     {
+        set_atoms(m, g->isize, g->index);
         for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
         {
             /* Find the next atom in the block */
@@ -1345,10 +1371,8 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
         }
         /* Update the number of blocks */
         m->mapb.index[bi] = g->isize;
-        m->nr             = bi;
-        m->bMapStatic     = false;
+        m->mapb.nr        = bi;
     }
-    m->mapb.nr = m->nr;
     m->bStatic = false;
 }
 
@@ -1371,6 +1395,10 @@ gmx_ana_indexmap_deinit(gmx_ana_indexmap_t *m)
     {
         sfree(m->mapb.index);
     }
+    if (m->mapb.nalloc_a > 0)
+    {
+        sfree(m->mapb.a);
+    }
     sfree(m->orgid);
     if (m->b.nalloc_index > 0)
     {