src/gromacs/selection/indexutil.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements functions in indexutil.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "indexutil.h"
  45
  46 #include <cstdlib>
  47 #include <cstring>
  48
  49 #include <algorithm>
  50 #include <string>
  51 #include <vector>
  52
  53 #include "gromacs/topology/block.h"
  54 #include "gromacs/topology/index.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/exceptions.h"
  57 #include "gromacs/utility/gmxassert.h"
  58 #include "gromacs/utility/smalloc.h"
  59 #include "gromacs/utility/stringutil.h"
  60 #include "gromacs/utility/textwriter.h"
  61
  62 /********************************************************************
  63  * gmx_ana_indexgrps_t functions
  64  ********************************************************************/
  65
  66 /*! \internal \brief
  67  * Stores a set of index groups.
  68  */
  69 struct gmx_ana_indexgrps_t
  70 {
  71     //! Initializes an empty set of groups.
  72     explicit gmx_ana_indexgrps_t(int nr) : nr(nr), g(NULL)
  73     {
  74         names.reserve(nr);
  75         snew(g, nr);
  76     }
  77     ~gmx_ana_indexgrps_t()
  78     {
  79         for (int i = 0; i < nr; ++i)
  80         {
  81             gmx_ana_index_deinit(&g[i]);
  82         }
  83         sfree(g);
  84     }
  85
  86     /** Number of index groups. */
  87     int                       nr;
  88     /** Array of index groups. */
  89     gmx_ana_index_t          *g;
  90     /** Group names. */
  91     std::vector<std::string>  names;
  92 };
  93
  94 /*!
  95  * \param[out] g     Index group structure.
  96  * \param[in]  top   Topology structure.
  97  * \param[in]  fnm   File name for the index file.
  98  *   Memory is automatically allocated.
  99  *
 100  * One or both of \p top or \p fnm can be NULL.
 101  * If \p top is NULL, an index file is required and the groups are read
 102  * from the file (uses Gromacs routine init_index()).
 103  * If \p fnm is NULL, default groups are constructed based on the
 104  * topology (uses Gromacs routine analyse()).
 105  * If both are null, the index group structure is initialized empty.
 106  */
 107 void
 108 gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
 109                        const char *fnm)
 110 {
 111     t_blocka *block = NULL;
 112     char    **names = NULL;
 113
 114     if (fnm)
 115     {
 116         block = init_index(fnm, &names);
 117     }
 118     else if (top)
 119     {
 120         block = new_blocka();
 121         analyse(&top->atoms, block, &names, FALSE, FALSE);
 122     }
 123     else
 124     {
 125         *g = new gmx_ana_indexgrps_t(0);
 126         return;
 127     }
 128
 129     try
 130     {
 131         *g = new gmx_ana_indexgrps_t(block->nr);
 132         for (int i = 0; i < block->nr; ++i)
 133         {
 134             gmx_ana_index_t *grp = &(*g)->g[i];
 135
 136             grp->isize = block->index[i+1] - block->index[i];
 137             snew(grp->index, grp->isize);
 138             for (int j = 0; j < grp->isize; ++j)
 139             {
 140                 grp->index[j] = block->a[block->index[i]+j];
 141             }
 142             grp->nalloc_index = grp->isize;
 143             (*g)->names.push_back(names[i]);
 144         }
 145     }
 146     catch (...)
 147     {
 148         for (int i = 0; i < block->nr; ++i)
 149         {
 150             sfree(names[i]);
 151         }
 152         sfree(names);
 153         done_blocka(block);
 154         sfree(block);
 155         throw;
 156     }
 157     for (int i = 0; i < block->nr; ++i)
 158     {
 159         sfree(names[i]);
 160     }
 161     sfree(names);
 162     done_blocka(block);
 163     sfree(block);
 164 }
 165
 166 /*!
 167  * \param[in] g  Index groups structure.
 168  *
 169  * The pointer \p g is invalid after the call.
 170  */
 171 void
 172 gmx_ana_indexgrps_free(gmx_ana_indexgrps_t *g)
 173 {
 174     delete g;
 175 }
 176
 177 /*!
 178  * \param[out] g     Index group structure.
 179  * \returns    true if \p g is empty, i.e., has 0 index groups.
 180  */
 181 bool
 182 gmx_ana_indexgrps_is_empty(gmx_ana_indexgrps_t *g)
 183 {
 184     return g->nr == 0;
 185 }
 186
 187 /*!
 188  * \param[in]  g     Index groups structure.
 189  * \param[in]  n     Index group number to get.
 190  * \returns    Pointer to the \p n'th index group in \p g.
 191  *
 192  * The returned pointer should not be freed.
 193  */
 194 gmx_ana_index_t *
 195 gmx_ana_indexgrps_get_grp(gmx_ana_indexgrps_t *g, int n)
 196 {
 197     if (n < 0 || n >= g->nr)
 198     {
 199         return NULL;
 200     }
 201     return &g->g[n];
 202 }
 203
 204 /*!
 205  * \param[out] dest     Output structure.
 206  * \param[out] destName Receives the name of the group if found.
 207  * \param[in]  src      Input index groups.
 208  * \param[in]  n        Number of the group to extract.
 209  * \returns true if \p n is a valid group in \p src, false otherwise.
 210  */
 211 bool
 212 gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, std::string *destName,
 213                           gmx_ana_indexgrps_t *src, int n)
 214 {
 215     destName->clear();
 216     if (n < 0 || n >= src->nr)
 217     {
 218         dest->isize = 0;
 219         return false;
 220     }
 221
 222     if (destName != NULL)
 223     {
 224         *destName = src->names[n];
 225     }
 226     gmx_ana_index_copy(dest, &src->g[n], true);
 227     return true;
 228 }
 229
 230 /*!
 231  * \param[out] dest     Output structure.
 232  * \param[out] destName Receives the name of the group if found.
 233  * \param[in]  src      Input index groups.
 234  * \param[in]  name     Name (or part of the name) of the group to extract.
 235  * \returns true if \p name is a valid group in \p src, false otherwise.
 236  *
 237  * Uses the Gromacs routine find_group() to find the actual group;
 238  * the comparison is case-insensitive.
 239  */
 240 bool
 241 gmx_ana_indexgrps_find(gmx_ana_index_t *dest, std::string *destName,
 242                        gmx_ana_indexgrps_t *src,
 243                        const char *name)
 244 {
 245     const char **names;
 246
 247     destName->clear();
 248     snew(names, src->nr);
 249     for (int i = 0; i < src->nr; ++i)
 250     {
 251         names[i] = src->names[i].c_str();
 252     }
 253     int n = find_group(const_cast<char *>(name), src->nr,
 254                        const_cast<char **>(names));
 255     sfree(names);
 256     if (n < 0)
 257     {
 258         dest->isize = 0;
 259         return false;
 260     }
 261
 262     return gmx_ana_indexgrps_extract(dest, destName, src, n);
 263 }
 264
 265 /*!
 266  * \param[in]  writer Writer to use for output.
 267  * \param[in]  g      Index groups to print.
 268  * \param[in]  maxn   Maximum number of indices to print
 269  *      (-1 = print all, 0 = print only names).
 270  */
 271 void
 272 gmx_ana_indexgrps_print(gmx::TextWriter *writer, gmx_ana_indexgrps_t *g, int maxn)
 273 {
 274     for (int i = 0; i < g->nr; ++i)
 275     {
 276         writer->writeString(gmx::formatString(" Group %2d \"%s\" ",
 277                                               i, g->names[i].c_str()));
 278         gmx_ana_index_dump(writer, &g->g[i], maxn);
 279     }
 280 }
 281
 282 /********************************************************************
 283  * gmx_ana_index_t functions
 284  ********************************************************************/
 285
 286 /*!
 287  * \param[in,out] g      Index group structure.
 288  * \param[in]     isize  Maximum number of atoms to reserve space for.
 289  */
 290 void
 291 gmx_ana_index_reserve(gmx_ana_index_t *g, int isize)
 292 {
 293     if (g->nalloc_index < isize)
 294     {
 295         srenew(g->index, isize);
 296         g->nalloc_index = isize;
 297     }
 298 }
 299
 300 /*!
 301  * \param[in,out] g      Index group structure.
 302  *
 303  * Resizes the memory allocated for holding the indices such that the
 304  * current contents fit.
 305  */
 306 void
 307 gmx_ana_index_squeeze(gmx_ana_index_t *g)
 308 {
 309     srenew(g->index, g->isize);
 310     g->nalloc_index = g->isize;
 311 }
 312
 313 /*!
 314  * \param[out] g      Output structure.
 315  *
 316  * Any contents of \p g are discarded without freeing.
 317  */
 318 void
 319 gmx_ana_index_clear(gmx_ana_index_t *g)
 320 {
 321     g->isize        = 0;
 322     g->index        = NULL;
 323     g->nalloc_index = 0;
 324 }
 325
 326 /*!
 327  * \param[out] g      Output structure.
 328  * \param[in]  isize  Number of atoms in the new group.
 329  * \param[in]  index  Array of \p isize atoms (can be NULL if \p isize is 0).
 330  * \param[in]  nalloc Number of elements allocated for \p index
 331  *   (if 0, \p index is not freed in gmx_ana_index_deinit())
 332  *
 333  * No copy if \p index is made.
 334  */
 335 void
 336 gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, int nalloc)
 337 {
 338     g->isize        = isize;
 339     g->index        = index;
 340     g->nalloc_index = nalloc;
 341 }
 342
 343 /*!
 344  * \param[out] g      Output structure.
 345  * \param[in]  natoms Number of atoms.
 346  */
 347 void
 348 gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms)
 349 {
 350     int  i;
 351
 352     g->isize = natoms;
 353     snew(g->index, natoms);
 354     for (i = 0; i < natoms; ++i)
 355     {
 356         g->index[i] = i;
 357     }
 358     g->nalloc_index = natoms;
 359 }
 360
 361 /*!
 362  * \param[in] g  Index group structure.
 363  *
 364  * The pointer \p g is not freed.
 365  */
 366 void
 367 gmx_ana_index_deinit(gmx_ana_index_t *g)
 368 {
 369     if (g->nalloc_index > 0)
 370     {
 371         sfree(g->index);
 372     }
 373     gmx_ana_index_clear(g);
 374 }
 375
 376 /*!
 377  * \param[out] dest   Destination index group.
 378  * \param[in]  src    Source index group.
 379  * \param[in]  bAlloc If true, memory is allocated at \p dest; otherwise,
 380  *   it is assumed that enough memory has been allocated for index.
 381  *
 382  * A deep copy of the name is only made if \p bAlloc is true.
 383  */
 384 void
 385 gmx_ana_index_copy(gmx_ana_index_t *dest, gmx_ana_index_t *src, bool bAlloc)
 386 {
 387     dest->isize = src->isize;
 388     if (dest->isize > 0)
 389     {
 390         if (bAlloc)
 391         {
 392             snew(dest->index, dest->isize);
 393             dest->nalloc_index = dest->isize;
 394         }
 395         std::memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
 396     }
 397 }
 398
 399 /*!
 400  * \param[in]  writer Writer to use for output.
 401  * \param[in]  g      Index group to print.
 402  * \param[in]  maxn   Maximum number of indices to print (-1 = print all).
 403  */
 404 void
 405 gmx_ana_index_dump(gmx::TextWriter *writer, gmx_ana_index_t *g, int maxn)
 406 {
 407     writer->writeString(gmx::formatString("(%d atoms)", g->isize));
 408     if (maxn != 0)
 409     {
 410         writer->writeString(":");
 411         int n = g->isize;
 412         if (maxn >= 0 && n > maxn)
 413         {
 414             n = maxn;
 415         }
 416         for (int j = 0; j < n; ++j)
 417         {
 418             writer->writeString(gmx::formatString(" %d", g->index[j]+1));
 419         }
 420         if (n < g->isize)
 421         {
 422             writer->writeString(" ...");
 423         }
 424     }
 425     writer->writeLine();
 426 }
 427
 428 int
 429 gmx_ana_index_get_max_index(gmx_ana_index_t *g)
 430 {
 431     if (g->isize == 0)
 432     {
 433         return 0;
 434     }
 435     else
 436     {
 437         return *std::max_element(g->index, g->index + g->isize);
 438     }
 439 }
 440
 441 /*!
 442  * \param[in]  g      Index group to check.
 443  * \returns    true if the index group is sorted and has no duplicates,
 444  *   false otherwise.
 445  */
 446 bool
 447 gmx_ana_index_check_sorted(gmx_ana_index_t *g)
 448 {
 449     int  i;
 450
 451     for (i = 0; i < g->isize-1; ++i)
 452     {
 453         if (g->index[i+1] <= g->index[i])
 454         {
 455             return false;
 456         }
 457     }
 458     return true;
 459 }
 460
 461 bool
 462 gmx_ana_index_check_range(gmx_ana_index_t *g, int natoms)
 463 {
 464     for (int i = 0; i < g->isize; ++i)
 465     {
 466         if (g->index[i] < 0 || g->index[i] >= natoms)
 467         {
 468             return false;
 469         }
 470     }
 471     return true;
 472 }
 473
 474 /********************************************************************
 475  * Set operations
 476  ********************************************************************/
 477
 478 /** Helper function for gmx_ana_index_sort(). */
 479 static int
 480 cmp_atomid(const void *a, const void *b)
 481 {
 482     if (*(atom_id *)a < *(atom_id *)b)
 483     {
 484         return -1;
 485     }
 486     if (*(atom_id *)a > *(atom_id *)b)
 487     {
 488         return 1;
 489     }
 490     return 0;
 491 }
 492
 493 /*!
 494  * \param[in,out] g  Index group to be sorted.
 495  */
 496 void
 497 gmx_ana_index_sort(gmx_ana_index_t *g)
 498 {
 499     std::qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
 500 }
 501
 502 /*!
 503  * \param[in]  a      Index group to check.
 504  * \param[in]  b      Index group to check.
 505  * \returns    true if \p a and \p b are equal, false otherwise.
 506  */
 507 bool
 508 gmx_ana_index_equals(gmx_ana_index_t *a, gmx_ana_index_t *b)
 509 {
 510     int  i;
 511
 512     if (a->isize != b->isize)
 513     {
 514         return false;
 515     }
 516     for (i = 0; i < a->isize; ++i)
 517     {
 518         if (a->index[i] != b->index[i])
 519         {
 520             return false;
 521         }
 522     }
 523     return true;
 524 }
 525
 526 /*!
 527  * \param[in]  a      Index group to check against.
 528  * \param[in]  b      Index group to check.
 529  * \returns    true if \p b is contained in \p a,
 530  *   false otherwise.
 531  *
 532  * If the elements are not in the same order in both groups, the function
 533  * fails. However, the groups do not need to be sorted.
 534  */
 535 bool
 536 gmx_ana_index_contains(gmx_ana_index_t *a, gmx_ana_index_t *b)
 537 {
 538     int  i, j;
 539
 540     for (i = j = 0; j < b->isize; ++i, ++j)
 541     {
 542         while (i < a->isize && a->index[i] != b->index[j])
 543         {
 544             ++i;
 545         }
 546         if (i == a->isize)
 547         {
 548             return false;
 549         }
 550     }
 551     return true;
 552 }
 553
 554 /*!
 555  * \param[out] dest Output index group (the intersection of \p a and \p b).
 556  * \param[in]  a    First index group.
 557  * \param[in]  b    Second index group.
 558  *
 559  * \p dest can be the same as \p a or \p b.
 560  */
 561 void
 562 gmx_ana_index_intersection(gmx_ana_index_t *dest,
 563                            gmx_ana_index_t *a, gmx_ana_index_t *b)
 564 {
 565     int i, j, k;
 566
 567     for (i = j = k = 0; i < a->isize && j < b->isize; ++i)
 568     {
 569         while (j < b->isize && b->index[j] < a->index[i])
 570         {
 571             ++j;
 572         }
 573         if (j < b->isize && b->index[j] == a->index[i])
 574         {
 575             dest->index[k++] = b->index[j++];
 576         }
 577     }
 578     dest->isize = k;
 579 }
 580
 581 /*!
 582  * \param[out] dest Output index group (the difference \p a - \p b).
 583  * \param[in]  a    First index group.
 584  * \param[in]  b    Second index group.
 585  *
 586  * \p dest can equal \p a, but not \p b.
 587  */
 588 void
 589 gmx_ana_index_difference(gmx_ana_index_t *dest,
 590                          gmx_ana_index_t *a, gmx_ana_index_t *b)
 591 {
 592     int i, j, k;
 593
 594     for (i = j = k = 0; i < a->isize; ++i)
 595     {
 596         while (j < b->isize && b->index[j] < a->index[i])
 597         {
 598             ++j;
 599         }
 600         if (j == b->isize || b->index[j] != a->index[i])
 601         {
 602             dest->index[k++] = a->index[i];
 603         }
 604     }
 605     dest->isize = k;
 606 }
 607
 608 /*!
 609  * \param[in]  a    First index group.
 610  * \param[in]  b    Second index group.
 611  * \returns    Size of the difference \p a - \p b.
 612  */
 613 int
 614 gmx_ana_index_difference_size(gmx_ana_index_t *a, gmx_ana_index_t *b)
 615 {
 616     int i, j, k;
 617
 618     for (i = j = k = 0; i < a->isize; ++i)
 619     {
 620         while (j < b->isize && b->index[j] < a->index[i])
 621         {
 622             ++j;
 623         }
 624         if (j == b->isize || b->index[j] != a->index[i])
 625         {
 626             ++k;
 627         }
 628     }
 629     return k;
 630 }
 631
 632 /*!
 633  * \param[out] dest1 Output group 1 (will equal \p g).
 634  * \param[out] dest2 Output group 2 (will equal \p src - \p g).
 635  * \param[in]  src   Group to be partitioned.
 636  * \param[in]  g     One partition.
 637  *
 638  * \pre \p g is a subset of \p src and both sets are sorted
 639  * \pre \p dest1 has allocated storage to store \p src
 640  * \post \p dest1 == \p g
 641  * \post \p dest2 == \p src - \p g
 642  *
 643  * No storage should be allocated for \p dest2; after the call,
 644  * \p dest2->index points to the memory allocated for \p dest1
 645  * (to a part that is not used by \p dest1).
 646  *
 647  * The calculation can be performed in-place by setting \p dest1 equal to
 648  * \p src.
 649  */
 650 void
 651 gmx_ana_index_partition(gmx_ana_index_t *dest1, gmx_ana_index_t *dest2,
 652                         gmx_ana_index_t *src, gmx_ana_index_t *g)
 653 {
 654     int i, j, k;
 655
 656     dest2->index = dest1->index + g->isize;
 657     dest2->isize = src->isize - g->isize;
 658     for (i = g->isize-1, j = src->isize-1, k = dest2->isize-1; i >= 0; --i, --j)
 659     {
 660         while (j >= 0 && src->index[j] != g->index[i])
 661         {
 662             dest2->index[k--] = src->index[j--];
 663         }
 664     }
 665     while (j >= 0)
 666     {
 667         dest2->index[k--] = src->index[j--];
 668     }
 669     gmx_ana_index_copy(dest1, g, false);
 670 }
 671
 672 /*!
 673  * \param[out] dest Output index group (the union of \p a and \p b).
 674  * \param[in]  a    First index group.
 675  * \param[in]  b    Second index group.
 676  *
 677  * \p a and \p b can have common items.
 678  * \p dest can equal \p a or \p b.
 679  *
 680  * \see gmx_ana_index_merge()
 681  */
 682 void
 683 gmx_ana_index_union(gmx_ana_index_t *dest,
 684                     gmx_ana_index_t *a, gmx_ana_index_t *b)
 685 {
 686     int dsize;
 687     int i, j, k;
 688
 689     dsize       = gmx_ana_index_difference_size(b, a);
 690     i           = a->isize - 1;
 691     j           = b->isize - 1;
 692     dest->isize = a->isize + dsize;
 693     for (k = dest->isize - 1; k >= 0; k--)
 694     {
 695         if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
 696         {
 697             dest->index[k] = b->index[j--];
 698         }
 699         else
 700         {
 701             if (j >= 0 && a->index[i] == b->index[j])
 702             {
 703                 --j;
 704             }
 705             dest->index[k] = a->index[i--];
 706         }
 707     }
 708 }
 709
 710 /*!
 711  * \param[out] dest Output index group (the union of \p a and \p b).
 712  * \param[in]  a    First index group.
 713  * \param[in]  b    Second index group.
 714  *
 715  * \p a and \p b should not have common items.
 716  * \p dest can equal \p a or \p b.
 717  *
 718  * \see gmx_ana_index_union()
 719  */
 720 void
 721 gmx_ana_index_merge(gmx_ana_index_t *dest,
 722                     gmx_ana_index_t *a, gmx_ana_index_t *b)
 723 {
 724     int i, j, k;
 725
 726     i           = a->isize - 1;
 727     j           = b->isize - 1;
 728     dest->isize = a->isize + b->isize;
 729     for (k = dest->isize - 1; k >= 0; k--)
 730     {
 731         if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
 732         {
 733             dest->index[k] = b->index[j--];
 734         }
 735         else
 736         {
 737             dest->index[k] = a->index[i--];
 738         }
 739     }
 740 }
 741
 742 /********************************************************************
 743  * gmx_ana_indexmap_t and related things
 744  ********************************************************************/
 745
 746 /*! \brief
 747  * Helper for splitting a sequence of atom indices into groups.
 748  *
 749  * \param[in]     atomIndex  Index of the next atom in the sequence.
 750  * \param[in]     top        Topology structure.
 751  * \param[in]     type       Type of group to split into.
 752  * \param[in,out] id         Variable to receive the next group id.
 753  * \returns  `true` if \p atomIndex starts a new group in the sequence, i.e.,
 754  *     if \p *id was changed.
 755  *
 756  * \p *id should be initialized to `-1` before first call of this function, and
 757  * then each atom index in the sequence passed to the function in turn.
 758  *
 759  * \ingroup module_selection
 760  */
 761 static bool
 762 next_group_index(int atomIndex, t_topology *top, e_index_t type, int *id)
 763 {
 764     int prev = *id;
 765     switch (type)
 766     {
 767         case INDEX_ATOM:
 768             *id = atomIndex;
 769             break;
 770         case INDEX_RES:
 771             *id = top->atoms.atom[atomIndex].resind;
 772             break;
 773         case INDEX_MOL:
 774             if (*id >= 0 && top->mols.index[*id] > atomIndex)
 775             {
 776                 *id = 0;
 777             }
 778             while (*id < top->mols.nr && atomIndex >= top->mols.index[*id+1])
 779             {
 780                 ++*id;
 781             }
 782             GMX_ASSERT(*id < top->mols.nr, "Molecules do not span all the atoms");
 783             break;
 784         case INDEX_UNKNOWN:
 785         case INDEX_ALL:
 786             *id = 0;
 787             break;
 788     }
 789     return prev != *id;
 790 }
 791
 792 /*!
 793  * \param[in,out] t    Output block.
 794  * \param[in]  top  Topology structure
 795  *   (only used if \p type is \ref INDEX_RES or \ref INDEX_MOL, can be NULL
 796  *   otherwise).
 797  * \param[in]  g    Index group
 798  *   (can be NULL if \p type is \ref INDEX_UNKNOWN).
 799  * \param[in]  type Type of partitioning to make.
 800  * \param[in]  bComplete
 801  *   If true, the index group is expanded to include any residue/molecule
 802  *   (depending on \p type) that is partially contained in the group.
 803  *   If \p type is not INDEX_RES or INDEX_MOL, this has no effect.
 804  *
 805  * \p m should have been initialized somehow (calloc() is enough).
 806  * \p g should be sorted.
 807  */
 808 void
 809 gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
 810                          e_index_t type, bool bComplete)
 811 {
 812     if (type == INDEX_UNKNOWN)
 813     {
 814         sfree(t->a);
 815         srenew(t->index, 2);
 816         t->nr           = 1;
 817         t->nalloc_index = 2;
 818         t->index[0]     = 0;
 819         t->index[1]     = 0;
 820         t->nra          = 0;
 821         t->a            = NULL;
 822         t->nalloc_a     = 0;
 823         return;
 824     }
 825
 826     // TODO: Check callers and either check these there as well, or turn these
 827     // into exceptions.
 828     GMX_RELEASE_ASSERT(top != NULL || (type != INDEX_RES && type != INDEX_MOL),
 829                        "Topology must be provided for residue or molecule blocks");
 830     GMX_RELEASE_ASSERT(!(type == INDEX_MOL && top->mols.nr == 0),
 831                        "Molecule information must be present for molecule blocks");
 832
 833     /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it
 834      * off otherwise. */
 835     if (type != INDEX_RES && type != INDEX_MOL)
 836     {
 837         bComplete = false;
 838     }
 839     /* Allocate memory for the atom array and fill it unless we are using
 840      * completion. */
 841     if (bComplete)
 842     {
 843         t->nra = 0;
 844         /* We may allocate some extra memory here because we don't know in
 845          * advance how much will be needed. */
 846         if (t->nalloc_a < top->atoms.nr)
 847         {
 848             srenew(t->a, top->atoms.nr);
 849             t->nalloc_a = top->atoms.nr;
 850         }
 851     }
 852     else
 853     {
 854         t->nra      = g->isize;
 855         if (t->nalloc_a < g->isize)
 856         {
 857             srenew(t->a, g->isize);
 858             t->nalloc_a = g->isize;
 859         }
 860         std::memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
 861     }
 862
 863     /* Allocate memory for the block index. We don't know in advance
 864      * how much will be needed, so we allocate some extra and free it in the
 865      * end. */
 866     if (t->nalloc_index < g->isize + 1)
 867     {
 868         srenew(t->index, g->isize + 1);
 869         t->nalloc_index = g->isize + 1;
 870     }
 871     /* Clear counters */
 872     t->nr  = 0;
 873     int j  = 0; /* j is used by residue completion for the first atom not stored */
 874     int id = -1;
 875     for (int i = 0; i < g->isize; ++i)
 876     {
 877         /* Find the ID number of the atom/residue/molecule corresponding to
 878          * the atom. */
 879         if (next_group_index(g->index[i], top, type, &id))
 880         {
 881             /* If this is the first atom in a new block, initialize the block. */
 882             if (bComplete)
 883             {
 884                 /* For completion, we first set the start of the block. */
 885                 t->index[t->nr++] = t->nra;
 886                 /* And then we find all the atoms that should be included. */
 887                 switch (type)
 888                 {
 889                     case INDEX_RES:
 890                         while (top->atoms.atom[j].resind != id)
 891                         {
 892                             ++j;
 893                         }
 894                         while (j < top->atoms.nr && top->atoms.atom[j].resind == id)
 895                         {
 896                             t->a[t->nra++] = j;
 897                             ++j;
 898                         }
 899                         break;
 900
 901                     case INDEX_MOL:
 902                         for (j = top->mols.index[id]; j < top->mols.index[id+1]; ++j)
 903                         {
 904                             t->a[t->nra++] = j;
 905                         }
 906                         break;
 907
 908                     default: /* Should not be reached */
 909                         GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
 910                         break;
 911                 }
 912             }
 913             else
 914             {
 915                 /* If not using completion, simply store the start of the block. */
 916                 t->index[t->nr++] = i;
 917             }
 918         }
 919     }
 920     /* Set the end of the last block */
 921     t->index[t->nr] = t->nra;
 922     /* Free any unnecessary memory */
 923     srenew(t->index, t->nr+1);
 924     t->nalloc_index = t->nr+1;
 925     if (bComplete)
 926     {
 927         srenew(t->a, t->nra);
 928         t->nalloc_a = t->nra;
 929     }
 930 }
 931
 932 /*!
 933  * \param[in] g   Index group to check.
 934  * \param[in] b   Block data to check against.
 935  * \returns   true if \p g consists of one or more complete blocks from \p b,
 936  *   false otherwise.
 937  *
 938  * The atoms in \p g are assumed to be sorted.
 939  */
 940 bool
 941 gmx_ana_index_has_full_blocks(gmx_ana_index_t *g, t_block *b)
 942 {
 943     int  i, j, bi;
 944
 945     i = bi = 0;
 946     /* Each round in the loop matches one block */
 947     while (i < g->isize)
 948     {
 949         /* Find the block that begins with the first unmatched atom */
 950         while (bi < b->nr && b->index[bi] != g->index[i])
 951         {
 952             ++bi;
 953         }
 954         /* If not found, or if too large, return */
 955         if (bi == b->nr || i + b->index[bi+1] -  b->index[bi] > g->isize)
 956         {
 957             return false;
 958         }
 959         /* Check that the block matches the index */
 960         for (j = b->index[bi]; j < b->index[bi+1]; ++j, ++i)
 961         {
 962             if (g->index[i] != j)
 963             {
 964                 return false;
 965             }
 966         }
 967         /* Move the search to the next block */
 968         ++bi;
 969     }
 970     return true;
 971 }
 972
 973 /*!
 974  * \param[in] g   Index group to check.
 975  * \param[in] b   Block data to check against.
 976  * \returns   true if \p g consists of one or more complete blocks from \p b,
 977  *   false otherwise.
 978  *
 979  * The atoms in \p g and \p b->a are assumed to be in the same order.
 980  */
 981 bool
 982 gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b)
 983 {
 984     int  i, j, bi;
 985
 986     i = bi = 0;
 987     /* Each round in the loop matches one block */
 988     while (i < g->isize)
 989     {
 990         /* Find the block that begins with the first unmatched atom */
 991         while (bi < b->nr && b->a[b->index[bi]] != g->index[i])
 992         {
 993             ++bi;
 994         }
 995         /* If not found, or if too large, return */
 996         if (bi == b->nr || i + b->index[bi+1] -  b->index[bi] > g->isize)
 997         {
 998             return false;
 999         }
1000         /* Check that the block matches the index */
1001         for (j = b->index[bi]; j < b->index[bi+1]; ++j, ++i)
1002         {
1003             if (b->a[j] != g->index[i])
1004             {
1005                 return false;
1006             }
1007         }
1008         /* Move the search to the next block */
1009         ++bi;
1010     }
1011     return true;
1012 }
1013
1014 /*!
1015  * \param[in] g     Index group to check.
1016  * \param[in] type  Block data to check against.
1017  * \param[in] top   Topology data.
1018  * \returns   true if \p g consists of one or more complete elements of type
1019  *   \p type, false otherwise.
1020  *
1021  * \p g is assumed to be sorted, otherwise may return false negatives.
1022  *
1023  * If \p type is \ref INDEX_ATOM, the return value is always true.
1024  * If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
1025  * always false.
1026  */
1027 bool
1028 gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
1029                                  t_topology *top)
1030 {
1031     // TODO: Consider whether unsorted groups need to be supported better.
1032     switch (type)
1033     {
1034         case INDEX_UNKNOWN:
1035         case INDEX_ALL:
1036             return false;
1037
1038         case INDEX_ATOM:
1039             return true;
1040
1041         case INDEX_RES:
1042         {
1043             int      i, ai;
1044             int      id, prev;
1045
1046             prev = -1;
1047             for (i = 0; i < g->isize; ++i)
1048             {
1049                 ai = g->index[i];
1050                 id = top->atoms.atom[ai].resind;
1051                 if (id != prev)
1052                 {
1053                     if (ai > 0 && top->atoms.atom[ai-1].resind == id)
1054                     {
1055                         return false;
1056                     }
1057                     if (i > 0 && g->index[i-1] < top->atoms.nr - 1
1058                         && top->atoms.atom[g->index[i-1]+1].resind == prev)
1059                     {
1060                         return false;
1061                     }
1062                 }
1063                 prev = id;
1064             }
1065             if (g->index[i-1] < top->atoms.nr - 1
1066                 && top->atoms.atom[g->index[i-1]+1].resind == prev)
1067             {
1068                 return false;
1069             }
1070             break;
1071         }
1072
1073         case INDEX_MOL:
1074             return gmx_ana_index_has_full_blocks(g, &top->mols);
1075     }
1076     return true;
1077 }
1078
1079 /*!
1080  * \param[out] m      Output structure.
1081  *
1082  * Any contents of \p m are discarded without freeing.
1083  */
1084 void
1085 gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
1086 {
1087     m->type              = INDEX_UNKNOWN;
1088     m->refid             = NULL;
1089     m->mapid             = NULL;
1090     m->mapb.nr           = 0;
1091     m->mapb.index        = NULL;
1092     m->mapb.nalloc_index = 0;
1093     m->mapb.nra          = 0;
1094     m->mapb.a            = NULL;
1095     m->mapb.nalloc_a     = 0;
1096     m->orgid             = NULL;
1097     m->b.nr              = 0;
1098     m->b.index           = NULL;
1099     m->b.nra             = 0;
1100     m->b.a               = NULL;
1101     m->b.nalloc_index    = 0;
1102     m->b.nalloc_a        = 0;
1103     m->bStatic           = true;
1104 }
1105
1106 /*!
1107  * \param[in,out] m      Mapping structure.
1108  * \param[in]     nr     Maximum number of blocks to reserve space for.
1109  * \param[in]     isize  Maximum number of atoms to reserve space for.
1110  */
1111 void
1112 gmx_ana_indexmap_reserve(gmx_ana_indexmap_t *m, int nr, int isize)
1113 {
1114     if (m->mapb.nalloc_index < nr + 1)
1115     {
1116         srenew(m->refid,      nr);
1117         srenew(m->mapid,      nr);
1118         srenew(m->orgid,      nr);
1119         srenew(m->mapb.index, nr + 1);
1120         srenew(m->b.index,    nr + 1);
1121         m->mapb.nalloc_index = nr + 1;
1122         m->b.nalloc_index    = nr + 1;
1123     }
1124     if (m->b.nalloc_a < isize)
1125     {
1126         srenew(m->b.a,        isize);
1127         m->b.nalloc_a = isize;
1128     }
1129 }
1130
1131 /*!
1132  * \param[in,out] m    Mapping structure to initialize.
1133  * \param[in]     g    Index group to map
1134  *   (can be NULL if \p type is \ref INDEX_UNKNOWN).
1135  * \param[in]     top  Topology structure
1136  *   (can be NULL if \p type is not \ref INDEX_RES or \ref INDEX_MOL).
1137  * \param[in]     type Type of mapping to construct.
1138  *
1139  * Initializes a new index group mapping.
1140  * The index group provided to gmx_ana_indexmap_update() should always be a
1141  * subset of the \p g given here.
1142  *
1143  * \p m should have been initialized somehow (calloc() is enough).
1144  */
1145 void
1146 gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1147                       t_topology *top, e_index_t type)
1148 {
1149     m->type   = type;
1150     gmx_ana_index_make_block(&m->b, top, g, type, false);
1151     gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
1152     int id = -1;
1153     for (int i = 0; i < m->b.nr; ++i)
1154     {
1155         const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1156         next_group_index(ii, top, type, &id);
1157         m->refid[i] = i;
1158         m->mapid[i] = id;
1159         m->orgid[i] = id;
1160     }
1161     m->mapb.nr  = m->b.nr;
1162     m->mapb.nra = m->b.nra;
1163     m->mapb.a   = m->b.a;
1164     std::memcpy(m->mapb.index, m->b.index, (m->b.nr+1)*sizeof(*(m->mapb.index)));
1165     m->bStatic  = true;
1166 }
1167
1168 int
1169 gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, t_topology *top,
1170                                   e_index_t type)
1171 {
1172     GMX_RELEASE_ASSERT(m->bStatic,
1173                        "Changing original IDs is not supported after starting "
1174                        "to use the mapping");
1175     GMX_RELEASE_ASSERT(top != NULL || (type != INDEX_RES && type != INDEX_MOL),
1176                        "Topology must be provided for residue or molecule blocks");
1177     GMX_RELEASE_ASSERT(!(type == INDEX_MOL && top->mols.nr == 0),
1178                        "Molecule information must be present for molecule blocks");
1179     // Check that all atoms in each block belong to the same group.
1180     // This is a separate loop for better error handling (no state is modified
1181     // if there is an error.
1182     if (type == INDEX_RES || type == INDEX_MOL)
1183     {
1184         int id = -1;
1185         for (int i = 0; i < m->b.nr; ++i)
1186         {
1187             const int ii = m->b.a[m->b.index[i]];
1188             if (next_group_index(ii, top, type, &id))
1189             {
1190                 for (int j = m->b.index[i] + 1; j < m->b.index[i+1]; ++j)
1191                 {
1192                     if (next_group_index(m->b.a[j], top, type, &id))
1193                     {
1194                         std::string message("Grouping into residues/molecules is ambiguous");
1195                         GMX_THROW(gmx::InconsistentInputError(message));
1196                     }
1197                 }
1198             }
1199         }
1200     }
1201     // Do a second loop, where things are actually set.
1202     int id    = -1;
1203     int group = -1;
1204     for (int i = 0; i < m->b.nr; ++i)
1205     {
1206         const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1207         if (next_group_index(ii, top, type, &id))
1208         {
1209             ++group;
1210         }
1211         m->mapid[i] = group;
1212         m->orgid[i] = group;
1213     }
1214     // Count also the last group.
1215     ++group;
1216     return group;
1217 }
1218
1219 /*!
1220  * \param[in,out] m    Mapping structure to initialize.
1221  * \param[in]     b    Block information to use for data.
1222  *
1223  * Frees some memory that is not necessary for static index group mappings.
1224  * Internal pointers are set to point to data in \p b; it is the responsibility
1225  * of the caller to ensure that the block information matches the contents of
1226  * the mapping.
1227  * After this function has been called, the index group provided to
1228  * gmx_ana_indexmap_update() should always be the same as \p g given here.
1229  *
1230  * This function breaks modularity of the index group mapping interface in an
1231  * ugly way, but allows reducing memory usage of static selections by a
1232  * significant amount.
1233  */
1234 void
1235 gmx_ana_indexmap_set_static(gmx_ana_indexmap_t *m, t_blocka *b)
1236 {
1237     sfree(m->mapid);
1238     sfree(m->mapb.index);
1239     sfree(m->b.index);
1240     sfree(m->b.a);
1241     m->mapb.nalloc_index = 0;
1242     m->mapb.nalloc_a     = 0;
1243     m->b.nalloc_index    = 0;
1244     m->b.nalloc_a        = 0;
1245     m->mapid             = m->orgid;
1246     m->mapb.index        = b->index;
1247     m->mapb.a            = b->a;
1248     m->b.index           = b->index;
1249     m->b.a               = b->a;
1250 }
1251
1252 /*!
1253  * \param[in,out] dest Destination data structure.
1254  * \param[in]     src  Source mapping.
1255  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
1256  *   copy is made; otherwise, only variable parts are copied, and no memory
1257  *   is allocated.
1258  *
1259  * \p dest should have been initialized somehow (calloc() is enough).
1260  */
1261 void
1262 gmx_ana_indexmap_copy(gmx_ana_indexmap_t *dest, gmx_ana_indexmap_t *src, bool bFirst)
1263 {
1264     if (bFirst)
1265     {
1266         gmx_ana_indexmap_reserve(dest, src->b.nr, src->b.nra);
1267         dest->type       = src->type;
1268         dest->b.nr       = src->b.nr;
1269         dest->b.nra      = src->b.nra;
1270         std::memcpy(dest->orgid,      src->orgid,      dest->b.nr*sizeof(*dest->orgid));
1271         std::memcpy(dest->b.index,    src->b.index,   (dest->b.nr+1)*sizeof(*dest->b.index));
1272         std::memcpy(dest->b.a,        src->b.a,        dest->b.nra*sizeof(*dest->b.a));
1273     }
1274     dest->mapb.nr    = src->mapb.nr;
1275     dest->mapb.nra   = src->mapb.nra;
1276     if (src->mapb.nalloc_a > 0)
1277     {
1278         if (bFirst)
1279         {
1280             snew(dest->mapb.a, src->mapb.nalloc_a);
1281             dest->mapb.nalloc_a = src->mapb.nalloc_a;
1282         }
1283         std::memcpy(dest->mapb.a, src->mapb.a, dest->mapb.nra*sizeof(*dest->mapb.a));
1284     }
1285     else
1286     {
1287         dest->mapb.a = src->mapb.a;
1288     }
1289     std::memcpy(dest->refid,      src->refid,      dest->mapb.nr*sizeof(*dest->refid));
1290     std::memcpy(dest->mapid,      src->mapid,      dest->mapb.nr*sizeof(*dest->mapid));
1291     std::memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
1292     dest->bStatic = src->bStatic;
1293 }
1294
1295 /*! \brief
1296  * Helper function to set the source atoms in an index map.
1297  *
1298  * \param[in,out] m     Mapping structure.
1299  * \param[in]     isize Number of atoms in the \p index array.
1300  * \param[in]     index List of atoms.
1301  */
1302 static void
1303 set_atoms(gmx_ana_indexmap_t *m, int isize, int *index)
1304 {
1305     m->mapb.nra = isize;
1306     if (m->mapb.nalloc_a == 0)
1307     {
1308         m->mapb.a = index;
1309     }
1310     else
1311     {
1312         for (int i = 0; i < isize; ++i)
1313         {
1314             m->mapb.a[i] = index[i];
1315         }
1316     }
1317 }
1318
1319 /*!
1320  * \param[in,out] m         Mapping structure.
1321  * \param[in]     g         Current index group.
1322  * \param[in]     bMaskOnly true if the unused blocks should be masked with
1323  *   -1 instead of removing them.
1324  *
1325  * Updates the index group mapping with the new index group \p g.
1326  *
1327  * \see gmx_ana_indexmap_t
1328  */
1329 void
1330 gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1331                         bool bMaskOnly)
1332 {
1333     int  i, j, bi, bj;
1334
1335     /* Process the simple cases first */
1336     if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
1337     {
1338         return;
1339     }
1340     if (m->type == INDEX_ALL)
1341     {
1342         set_atoms(m, g->isize, g->index);
1343         if (m->b.nr > 0)
1344         {
1345             m->mapb.index[1] = g->isize;
1346         }
1347         return;
1348     }
1349     /* Reset the reference IDs and mapping if necessary */
1350     const bool bToFull  = (g->isize == m->b.nra);
1351     const bool bWasFull = (m->mapb.nra == m->b.nra);
1352     if (bToFull || bMaskOnly)
1353     {
1354         if (!m->bStatic)
1355         {
1356             for (bj = 0; bj < m->b.nr; ++bj)
1357             {
1358                 m->refid[bj] = bj;
1359             }
1360         }
1361         if (!bWasFull)
1362         {
1363             for (bj = 0; bj < m->b.nr; ++bj)
1364             {
1365                 m->mapid[bj] = m->orgid[bj];
1366             }
1367             for (bj = 0; bj <= m->b.nr; ++bj)
1368             {
1369                 m->mapb.index[bj] = m->b.index[bj];
1370             }
1371         }
1372         set_atoms(m, m->b.nra, m->b.a);
1373         m->mapb.nr = m->b.nr;
1374     }
1375     /* Exit immediately if the group is static */
1376     if (bToFull)
1377     {
1378         m->bStatic = true;
1379         return;
1380     }
1381
1382     if (bMaskOnly)
1383     {
1384         for (i = j = bj = 0; i < g->isize; ++i, ++j)
1385         {
1386             /* Find the next atom in the block */
1387             while (m->b.a[j] != g->index[i])
1388             {
1389                 ++j;
1390             }
1391             /* Mark blocks that did not contain any atoms */
1392             while (bj < m->b.nr && m->b.index[bj+1] <= j)
1393             {
1394                 m->refid[bj++] = -1;
1395             }
1396             /* Advance the block index if we have reached the next block */
1397             if (m->b.index[bj] <= j)
1398             {
1399                 ++bj;
1400             }
1401         }
1402         /* Mark the last blocks as not accessible */
1403         while (bj < m->b.nr)
1404         {
1405             m->refid[bj++] = -1;
1406         }
1407     }
1408     else
1409     {
1410         set_atoms(m, g->isize, g->index);
1411         for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
1412         {
1413             /* Find the next atom in the block */
1414             while (m->b.a[j] != g->index[i])
1415             {
1416                 ++j;
1417             }
1418             /* If we have reached a new block, add it */
1419             if (m->b.index[bj+1] <= j)
1420             {
1421                 /* Skip any blocks in between */
1422                 while (bj < m->b.nr && m->b.index[bj+1] <= j)
1423                 {
1424                     ++bj;
1425                 }
1426                 m->refid[bi]      = bj;
1427                 m->mapid[bi]      = m->orgid[bj];
1428                 m->mapb.index[bi] = i;
1429                 bi++;
1430             }
1431         }
1432         /* Update the number of blocks */
1433         m->mapb.index[bi] = g->isize;
1434         m->mapb.nr        = bi;
1435     }
1436     m->bStatic = false;
1437 }
1438
1439 /*!
1440  * \param[in,out] m         Mapping structure to free.
1441  *
1442  * All the memory allocated for the mapping structure is freed, and
1443  * the pointers set to NULL.
1444  * The pointer \p m is not freed.
1445  */
1446 void
1447 gmx_ana_indexmap_deinit(gmx_ana_indexmap_t *m)
1448 {
1449     sfree(m->refid);
1450     if (m->mapid != m->orgid)
1451     {
1452         sfree(m->mapid);
1453     }
1454     if (m->mapb.nalloc_index > 0)
1455     {
1456         sfree(m->mapb.index);
1457     }
1458     if (m->mapb.nalloc_a > 0)
1459     {
1460         sfree(m->mapb.a);
1461     }
1462     sfree(m->orgid);
1463     if (m->b.nalloc_index > 0)
1464     {
1465         sfree(m->b.index);
1466     }
1467     if (m->b.nalloc_a > 0)
1468     {
1469         sfree(m->b.a);
1470     }
1471     gmx_ana_indexmap_clear(m);
1472 }