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