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