2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
37 * Implements functions in indexutil.h.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_selection
44 #include "indexutil.h"
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"
64 /********************************************************************
65 * gmx_ana_indexgrps_t functions
66 ********************************************************************/
69 * Stores a set of index groups.
71 struct gmx_ana_indexgrps_t
73 //! Initializes an empty set of groups.
74 explicit gmx_ana_indexgrps_t(int nr) : nr(nr), g(nullptr)
79 ~gmx_ana_indexgrps_t()
81 for (int i = 0; i < nr; ++i)
83 gmx_ana_index_deinit(&g[i]);
88 /** Number of index groups. */
90 /** Array of index groups. */
93 std::vector<std::string> names;
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.
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.
110 gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, gmx_mtop_t *top,
113 t_blocka *block = nullptr;
114 char **names = nullptr;
118 block = init_index(fnm, &names);
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);
130 *g = new gmx_ana_indexgrps_t(0);
136 *g = new gmx_ana_indexgrps_t(block->nr);
137 for (int i = 0; i < block->nr; ++i)
139 gmx_ana_index_t *grp = &(*g)->g[i];
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)
145 grp->index[j] = block->a[block->index[i]+j];
147 grp->nalloc_index = grp->isize;
148 (*g)->names.emplace_back(names[i]);
153 for (int i = 0; i < block->nr; ++i)
162 for (int i = 0; i < block->nr; ++i)
172 * \param[in] g Index groups structure.
174 * The pointer \p g is invalid after the call.
177 gmx_ana_indexgrps_free(gmx_ana_indexgrps_t *g)
183 * \param[out] g Index group structure.
184 * \returns true if \p g is empty, i.e., has 0 index groups.
187 gmx_ana_indexgrps_is_empty(gmx_ana_indexgrps_t *g)
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.
197 * The returned pointer should not be freed.
200 gmx_ana_indexgrps_get_grp(gmx_ana_indexgrps_t *g, int n)
202 if (n < 0 || n >= g->nr)
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.
217 gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, std::string *destName,
218 gmx_ana_indexgrps_t *src, int n)
221 if (n < 0 || n >= src->nr)
227 if (destName != nullptr)
229 *destName = src->names[n];
231 gmx_ana_index_copy(dest, &src->g[n], true);
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.
242 * Uses the Gromacs routine find_group() to find the actual group;
243 * the comparison is case-insensitive.
246 gmx_ana_indexgrps_find(gmx_ana_index_t *dest, std::string *destName,
247 gmx_ana_indexgrps_t *src,
253 snew(names, src->nr);
254 for (int i = 0; i < src->nr; ++i)
256 names[i] = src->names[i].c_str();
258 int n = find_group(const_cast<char *>(name), src->nr,
259 const_cast<char **>(names));
267 return gmx_ana_indexgrps_extract(dest, destName, src, n);
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).
277 gmx_ana_indexgrps_print(gmx::TextWriter *writer, gmx_ana_indexgrps_t *g, int maxn)
279 for (int i = 0; i < g->nr; ++i)
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);
287 /********************************************************************
288 * gmx_ana_index_t functions
289 ********************************************************************/
292 * \param[in,out] g Index group structure.
293 * \param[in] isize Maximum number of atoms to reserve space for.
296 gmx_ana_index_reserve(gmx_ana_index_t *g, int isize)
298 if (g->nalloc_index < isize)
300 srenew(g->index, isize);
301 g->nalloc_index = isize;
306 * \param[in,out] g Index group structure.
308 * Resizes the memory allocated for holding the indices such that the
309 * current contents fit.
312 gmx_ana_index_squeeze(gmx_ana_index_t *g)
314 srenew(g->index, g->isize);
315 g->nalloc_index = g->isize;
319 * \param[out] g Output structure.
321 * Any contents of \p g are discarded without freeing.
324 gmx_ana_index_clear(gmx_ana_index_t *g)
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())
338 * No copy if \p index is made.
341 gmx_ana_index_set(gmx_ana_index_t *g, int isize, int *index, int nalloc)
345 g->nalloc_index = nalloc;
349 * \param[out] g Output structure.
350 * \param[in] natoms Number of atoms.
353 gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms)
358 snew(g->index, natoms);
359 for (i = 0; i < natoms; ++i)
363 g->nalloc_index = natoms;
367 * \param[in] g Index group structure.
369 * The pointer \p g is not freed.
372 gmx_ana_index_deinit(gmx_ana_index_t *g)
374 if (g->nalloc_index > 0)
378 gmx_ana_index_clear(g);
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.
388 gmx_ana_index_copy(gmx_ana_index_t *dest, gmx_ana_index_t *src, bool bAlloc)
390 dest->isize = src->isize;
393 snew(dest->index, dest->isize);
394 dest->nalloc_index = dest->isize;
398 std::memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
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).
408 gmx_ana_index_dump(gmx::TextWriter *writer, gmx_ana_index_t *g, int maxn)
410 writer->writeString(gmx::formatString("(%d atoms)", g->isize));
413 writer->writeString(":");
415 if (maxn >= 0 && n > maxn)
419 for (int j = 0; j < n; ++j)
421 writer->writeString(gmx::formatString(" %d", g->index[j]+1));
425 writer->writeString(" ...");
428 writer->ensureLineBreak();
432 gmx_ana_index_get_max_index(gmx_ana_index_t *g)
440 return *std::max_element(g->index, g->index + g->isize);
445 * \param[in] g Index group to check.
446 * \returns true if the index group is sorted and has no duplicates,
450 gmx_ana_index_check_sorted(gmx_ana_index_t *g)
454 for (i = 0; i < g->isize-1; ++i)
456 if (g->index[i+1] <= g->index[i])
465 gmx_ana_index_check_range(gmx_ana_index_t *g, int natoms)
467 for (int i = 0; i < g->isize; ++i)
469 if (g->index[i] < 0 || g->index[i] >= natoms)
477 /********************************************************************
479 ********************************************************************/
481 /** Helper function for gmx_ana_index_sort(). */
483 cmp_atomid(const void *a, const void *b)
485 if (*(int *)a < *(int *)b)
489 if (*(int *)a > *(int *)b)
497 * \param[in,out] g Index group to be sorted.
500 gmx_ana_index_sort(gmx_ana_index_t *g)
502 std::qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
506 gmx_ana_index_remove_duplicates(gmx_ana_index_t *g)
509 for (int i = 0; i < g->isize; ++i)
511 if (i == 0 || g->index[i-1] != g->index[i])
513 g->index[j] = g->index[i];
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.
526 gmx_ana_index_equals(gmx_ana_index_t *a, gmx_ana_index_t *b)
530 if (a->isize != b->isize)
534 for (i = 0; i < a->isize; ++i)
536 if (a->index[i] != b->index[i])
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,
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.
554 gmx_ana_index_contains(gmx_ana_index_t *a, gmx_ana_index_t *b)
558 for (i = j = 0; j < b->isize; ++i, ++j)
560 while (i < a->isize && a->index[i] != b->index[j])
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.
577 * \p dest can be the same as \p a or \p b.
580 gmx_ana_index_intersection(gmx_ana_index_t *dest,
581 gmx_ana_index_t *a, gmx_ana_index_t *b)
585 for (i = j = k = 0; i < a->isize && j < b->isize; ++i)
587 while (j < b->isize && b->index[j] < a->index[i])
591 if (j < b->isize && b->index[j] == a->index[i])
593 dest->index[k++] = b->index[j++];
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.
604 * \p dest can equal \p a, but not \p b.
607 gmx_ana_index_difference(gmx_ana_index_t *dest,
608 gmx_ana_index_t *a, gmx_ana_index_t *b)
612 for (i = j = k = 0; i < a->isize; ++i)
614 while (j < b->isize && b->index[j] < a->index[i])
618 if (j == b->isize || b->index[j] != a->index[i])
620 dest->index[k++] = a->index[i];
627 * \param[in] a First index group.
628 * \param[in] b Second index group.
629 * \returns Size of the difference \p a - \p b.
632 gmx_ana_index_difference_size(gmx_ana_index_t *a, gmx_ana_index_t *b)
636 for (i = j = k = 0; i < a->isize; ++i)
638 while (j < b->isize && b->index[j] < a->index[i])
642 if (j == b->isize || b->index[j] != a->index[i])
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.
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
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).
665 * The calculation can be performed in-place by setting \p dest1 equal to
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)
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)
678 while (j >= 0 && src->index[j] != g->index[i])
680 dest2->index[k--] = src->index[j--];
685 dest2->index[k--] = src->index[j--];
687 gmx_ana_index_copy(dest1, g, false);
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.
695 * \p a and \p b can have common items.
696 * \p dest can equal \p a or \p b.
698 * \see gmx_ana_index_merge()
701 gmx_ana_index_union(gmx_ana_index_t *dest,
702 gmx_ana_index_t *a, gmx_ana_index_t *b)
707 dsize = gmx_ana_index_difference_size(b, a);
710 dest->isize = a->isize + dsize;
711 for (k = dest->isize - 1; k >= 0; k--)
713 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
715 dest->index[k] = b->index[j--];
719 if (j >= 0 && a->index[i] == b->index[j])
723 dest->index[k] = a->index[i--];
729 gmx_ana_index_union_unsorted(gmx_ana_index_t *dest,
730 gmx_ana_index_t *a, gmx_ana_index_t *b)
732 if (gmx_ana_index_check_sorted(b))
734 gmx_ana_index_union(dest, a, b);
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);
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.
752 * \p a and \p b should not have common items.
753 * \p dest can equal \p a or \p b.
755 * \see gmx_ana_index_union()
758 gmx_ana_index_merge(gmx_ana_index_t *dest,
759 gmx_ana_index_t *a, gmx_ana_index_t *b)
765 dest->isize = a->isize + b->isize;
766 for (k = dest->isize - 1; k >= 0; k--)
768 if (i < 0 || (j >= 0 && a->index[i] < b->index[j]))
770 dest->index[k] = b->index[j--];
774 dest->index[k] = a->index[i--];
779 /********************************************************************
780 * gmx_ana_indexmap_t and related things
781 ********************************************************************/
784 * Helper for splitting a sequence of atom indices into groups.
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.
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.
796 * \ingroup module_selection
799 next_group_index(int atomIndex, const gmx_mtop_t *top,
800 e_index_t type, int *id)
810 int resind, molb = 0;
811 mtopGetAtomAndResidueName(top, atomIndex, &molb, nullptr, nullptr, nullptr, &resind);
818 *id = mtopGetMoleculeIndex(top, atomIndex, &molb);
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
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.
842 * \p m should have been initialized somehow (calloc() is enough).
843 * \p g should be sorted.
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)
849 if (type == INDEX_UNKNOWN)
863 // TODO: Check callers and either check these there as well, or turn these
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");
870 /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it
872 if (type != INDEX_RES && type != INDEX_MOL)
876 /* Allocate memory for the atom array and fill it unless we are using
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)
885 srenew(t->a, top->natoms);
886 t->nalloc_a = top->natoms;
892 if (t->nalloc_a < g->isize)
894 srenew(t->a, g->isize);
895 t->nalloc_a = g->isize;
897 std::memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
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
903 if (t->nalloc_index < g->isize + 1)
905 srenew(t->index, g->isize + 1);
906 t->nalloc_index = g->isize + 1;
912 for (int i = 0; i < g->isize; ++i)
914 const int ai = g->index[i];
915 /* Find the ID number of the atom/residue/molecule corresponding to
917 if (next_group_index(ai, top, type, &id))
919 /* If this is the first atom in a new block, initialize the block. */
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. */
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)
938 int first_atom = atnr_mol - 1;
939 while (first_atom >= 0
940 && mol_atoms.atom[first_atom].resind == id)
944 int first_mol_atom = top->moleculeBlockIndices[molb].globalAtomStart;
945 first_mol_atom += molnr*top->moleculeBlockIndices[molb].numAtomsPerMolecule;
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)
957 while (molb + 1 < top->molblock.size() && id >= top->moleculeBlockIndices[molb].moleculeIndexStart)
961 const MoleculeBlockIndices &blockIndices = top->moleculeBlockIndices[molb];
962 const int atomStart = blockIndices.globalAtomStart + (id - blockIndices.moleculeIndexStart)*blockIndices.numAtomsPerMolecule;
963 for (int j = 0; j < blockIndices.numAtomsPerMolecule; ++j)
965 t->a[t->nra++] = atomStart + j;
969 default: /* Should not be reached */
970 GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
976 /* If not using completion, simply store the start of the block. */
977 t->index[t->nr++] = i;
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;
988 srenew(t->a, t->nra);
989 t->nalloc_a = t->nra;
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,
999 * The atoms in \p g are assumed to be sorted.
1002 gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g,
1003 const gmx::RangePartitioning *b)
1008 /* Each round in the loop matches one block */
1009 while (i < g->isize)
1011 /* Find the block that begins with the first unmatched atom */
1012 while (bi < b->numBlocks() && b->block(bi).begin() != g->index[i])
1016 /* If not found, or if too large, return */
1017 if (bi == b->numBlocks() || i + b->block(bi).size() > g->isize)
1021 /* Check that the block matches the index */
1022 for (j = b->block(bi).begin(); j < b->block(bi).end(); ++j, ++i)
1024 if (g->index[i] != j)
1029 /* Move the search to the next block */
1036 * \param[in] g Index group to check.
1037 * \param[in] b Block data to check against.
1038 * \returns true if \p g consists of one or more complete blocks from \p b,
1041 * The atoms in \p g and \p b->a are assumed to be in the same order.
1044 gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b)
1049 /* Each round in the loop matches one block */
1050 while (i < g->isize)
1052 /* Find the block that begins with the first unmatched atom */
1053 while (bi < b->nr && b->a[b->index[bi]] != g->index[i])
1057 /* If not found, or if too large, return */
1058 if (bi == b->nr || i + b->index[bi+1] - b->index[bi] > g->isize)
1062 /* Check that the block matches the index */
1063 for (j = b->index[bi]; j < b->index[bi+1]; ++j, ++i)
1065 if (b->a[j] != g->index[i])
1070 /* Move the search to the next block */
1077 * \brief Returns if an atom is at a residue boundary.
1079 * \param[in] top Topology data.
1080 * \param[in] a Atom index to check, should be -1 <= \p a < top->natoms.
1081 * \param[in,out] molb The molecule block of atom a
1082 * \returns true if atoms \p a and \p a + 1 are in different residues, false otherwise.
1084 static bool is_at_residue_boundary(const gmx_mtop_t *top, int a, int *molb)
1086 if (a == -1 || a + 1 == top->natoms)
1091 mtopGetAtomAndResidueName(top, a, molb,
1092 nullptr, nullptr, nullptr, &resindA);
1094 mtopGetAtomAndResidueName(top, a + 1, molb,
1095 nullptr, nullptr, nullptr, &resindAPlusOne);
1096 return resindAPlusOne != resindA;
1100 * \param[in] g Index group to check.
1101 * \param[in] type Block data to check against.
1102 * \param[in] top Topology data.
1103 * \returns true if \p g consists of one or more complete elements of type
1104 * \p type, false otherwise.
1106 * \p g is assumed to be sorted, otherwise may return false negatives.
1108 * If \p type is \ref INDEX_ATOM, the return value is always true.
1109 * If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
1113 gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
1114 const gmx_mtop_t *top)
1121 // TODO: Consider whether unsorted groups need to be supported better.
1135 for (int i = 0; i < g->isize; ++i)
1137 const int a = g->index[i];
1138 // Check if a is consecutive or on a residue boundary
1141 if (!is_at_residue_boundary(top, aPrev, &molb))
1145 if (!is_at_residue_boundary(top, a - 1, &molb))
1152 GMX_ASSERT(g->isize > 0, "We return above when isize=0");
1153 const int a = g->index[g->isize - 1];
1154 if (!is_at_residue_boundary(top, a, &molb))
1163 auto molecules = gmx_mtop_molecules(*top);
1164 return gmx_ana_index_has_full_blocks(g, &molecules);
1171 * \param[out] m Output structure.
1173 * Any contents of \p m are discarded without freeing.
1176 gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
1178 m->type = INDEX_UNKNOWN;
1182 m->mapb.index = nullptr;
1183 m->mapb.nalloc_index = 0;
1185 m->mapb.a = nullptr;
1186 m->mapb.nalloc_a = 0;
1189 m->b.index = nullptr;
1192 m->b.nalloc_index = 0;
1198 * \param[in,out] m Mapping structure.
1199 * \param[in] nr Maximum number of blocks to reserve space for.
1200 * \param[in] isize Maximum number of atoms to reserve space for.
1203 gmx_ana_indexmap_reserve(gmx_ana_indexmap_t *m, int nr, int isize)
1205 if (m->mapb.nalloc_index < nr + 1)
1207 srenew(m->refid, nr);
1208 srenew(m->mapid, nr);
1209 srenew(m->orgid, nr);
1210 srenew(m->mapb.index, nr + 1);
1211 srenew(m->b.index, nr + 1);
1212 m->mapb.nalloc_index = nr + 1;
1213 m->b.nalloc_index = nr + 1;
1215 if (m->b.nalloc_a < isize)
1217 srenew(m->b.a, isize);
1218 m->b.nalloc_a = isize;
1223 * \param[in,out] m Mapping structure to initialize.
1224 * \param[in] g Index group to map
1225 * (can be NULL if \p type is \ref INDEX_UNKNOWN).
1226 * \param[in] top Topology structure
1227 * (can be NULL if \p type is not \ref INDEX_RES or \ref INDEX_MOL).
1228 * \param[in] type Type of mapping to construct.
1230 * Initializes a new index group mapping.
1231 * The index group provided to gmx_ana_indexmap_update() should always be a
1232 * subset of the \p g given here.
1234 * \p m should have been initialized somehow (calloc() is enough).
1237 gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1238 const gmx_mtop_t *top, e_index_t type)
1241 gmx_ana_index_make_block(&m->b, top, g, type, false);
1242 gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
1244 for (int i = 0; i < m->b.nr; ++i)
1246 const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1247 next_group_index(ii, top, type, &id);
1252 m->mapb.nr = m->b.nr;
1253 m->mapb.nra = m->b.nra;
1255 std::memcpy(m->mapb.index, m->b.index, (m->b.nr+1)*sizeof(*(m->mapb.index)));
1260 gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, const gmx_mtop_t *top,
1263 GMX_RELEASE_ASSERT(m->bStatic,
1264 "Changing original IDs is not supported after starting "
1265 "to use the mapping");
1266 GMX_RELEASE_ASSERT(top != nullptr || (type != INDEX_RES && type != INDEX_MOL),
1267 "Topology must be provided for residue or molecule blocks");
1268 // Check that all atoms in each block belong to the same group.
1269 // This is a separate loop for better error handling (no state is modified
1270 // if there is an error.
1271 if (type == INDEX_RES || type == INDEX_MOL)
1274 for (int i = 0; i < m->b.nr; ++i)
1276 const int ii = m->b.a[m->b.index[i]];
1277 if (next_group_index(ii, top, type, &id))
1279 for (int j = m->b.index[i] + 1; j < m->b.index[i+1]; ++j)
1281 if (next_group_index(m->b.a[j], top, type, &id))
1283 std::string message("Grouping into residues/molecules is ambiguous");
1284 GMX_THROW(gmx::InconsistentInputError(message));
1290 // Do a second loop, where things are actually set.
1293 for (int i = 0; i < m->b.nr; ++i)
1295 const int ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
1296 if (next_group_index(ii, top, type, &id))
1300 m->mapid[i] = group;
1301 m->orgid[i] = group;
1303 // Count also the last group.
1309 * \param[in,out] m Mapping structure to initialize.
1310 * \param[in] b Block information to use for data.
1312 * Frees some memory that is not necessary for static index group mappings.
1313 * Internal pointers are set to point to data in \p b; it is the responsibility
1314 * of the caller to ensure that the block information matches the contents of
1316 * After this function has been called, the index group provided to
1317 * gmx_ana_indexmap_update() should always be the same as \p g given here.
1319 * This function breaks modularity of the index group mapping interface in an
1320 * ugly way, but allows reducing memory usage of static selections by a
1321 * significant amount.
1324 gmx_ana_indexmap_set_static(gmx_ana_indexmap_t *m, t_blocka *b)
1327 sfree(m->mapb.index);
1330 m->mapb.nalloc_index = 0;
1331 m->mapb.nalloc_a = 0;
1332 m->b.nalloc_index = 0;
1334 m->mapid = m->orgid;
1335 m->mapb.index = b->index;
1337 m->b.index = b->index;
1342 * \param[in,out] dest Destination data structure.
1343 * \param[in] src Source mapping.
1344 * \param[in] bFirst If true, memory is allocated for \p dest and a full
1345 * copy is made; otherwise, only variable parts are copied, and no memory
1348 * \p dest should have been initialized somehow (calloc() is enough).
1351 gmx_ana_indexmap_copy(gmx_ana_indexmap_t *dest, gmx_ana_indexmap_t *src, bool bFirst)
1355 gmx_ana_indexmap_reserve(dest, src->b.nr, src->b.nra);
1356 dest->type = src->type;
1357 dest->b.nr = src->b.nr;
1358 dest->b.nra = src->b.nra;
1359 std::memcpy(dest->orgid, src->orgid, dest->b.nr*sizeof(*dest->orgid));
1360 std::memcpy(dest->b.index, src->b.index, (dest->b.nr+1)*sizeof(*dest->b.index));
1361 std::memcpy(dest->b.a, src->b.a, dest->b.nra*sizeof(*dest->b.a));
1363 dest->mapb.nr = src->mapb.nr;
1364 dest->mapb.nra = src->mapb.nra;
1365 if (src->mapb.nalloc_a > 0)
1369 snew(dest->mapb.a, src->mapb.nalloc_a);
1370 dest->mapb.nalloc_a = src->mapb.nalloc_a;
1372 std::memcpy(dest->mapb.a, src->mapb.a, dest->mapb.nra*sizeof(*dest->mapb.a));
1376 dest->mapb.a = src->mapb.a;
1378 std::memcpy(dest->refid, src->refid, dest->mapb.nr*sizeof(*dest->refid));
1379 std::memcpy(dest->mapid, src->mapid, dest->mapb.nr*sizeof(*dest->mapid));
1380 std::memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
1381 dest->bStatic = src->bStatic;
1385 * Helper function to set the source atoms in an index map.
1387 * \param[in,out] m Mapping structure.
1388 * \param[in] isize Number of atoms in the \p index array.
1389 * \param[in] index List of atoms.
1392 set_atoms(gmx_ana_indexmap_t *m, int isize, int *index)
1394 m->mapb.nra = isize;
1395 if (m->mapb.nalloc_a == 0)
1401 for (int i = 0; i < isize; ++i)
1403 m->mapb.a[i] = index[i];
1409 * \param[in,out] m Mapping structure.
1410 * \param[in] g Current index group.
1411 * \param[in] bMaskOnly true if the unused blocks should be masked with
1412 * -1 instead of removing them.
1414 * Updates the index group mapping with the new index group \p g.
1416 * \see gmx_ana_indexmap_t
1419 gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
1424 /* Process the simple cases first */
1425 if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
1429 if (m->type == INDEX_ALL)
1431 set_atoms(m, g->isize, g->index);
1434 m->mapb.index[1] = g->isize;
1438 /* Reset the reference IDs and mapping if necessary */
1439 const bool bToFull = (g->isize == m->b.nra);
1440 const bool bWasFull = (m->mapb.nra == m->b.nra);
1441 if (bToFull || bMaskOnly)
1445 for (bj = 0; bj < m->b.nr; ++bj)
1452 for (bj = 0; bj < m->b.nr; ++bj)
1454 m->mapid[bj] = m->orgid[bj];
1456 for (bj = 0; bj <= m->b.nr; ++bj)
1458 m->mapb.index[bj] = m->b.index[bj];
1461 set_atoms(m, m->b.nra, m->b.a);
1462 m->mapb.nr = m->b.nr;
1464 /* Exit immediately if the group is static */
1473 for (i = j = bj = 0; i < g->isize; ++i, ++j)
1475 /* Find the next atom in the block */
1476 while (m->b.a[j] != g->index[i])
1480 /* Mark blocks that did not contain any atoms */
1481 while (bj < m->b.nr && m->b.index[bj+1] <= j)
1483 m->refid[bj++] = -1;
1485 /* Advance the block index if we have reached the next block */
1486 if (m->b.index[bj] <= j)
1491 /* Mark the last blocks as not accessible */
1492 while (bj < m->b.nr)
1494 m->refid[bj++] = -1;
1499 set_atoms(m, g->isize, g->index);
1500 for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
1502 /* Find the next atom in the block */
1503 while (m->b.a[j] != g->index[i])
1507 /* If we have reached a new block, add it */
1508 if (m->b.index[bj+1] <= j)
1510 /* Skip any blocks in between */
1511 while (bj < m->b.nr && m->b.index[bj+1] <= j)
1516 m->mapid[bi] = m->orgid[bj];
1517 m->mapb.index[bi] = i;
1521 /* Update the number of blocks */
1522 m->mapb.index[bi] = g->isize;
1529 * \param[in,out] m Mapping structure to free.
1531 * All the memory allocated for the mapping structure is freed, and
1532 * the pointers set to NULL.
1533 * The pointer \p m is not freed.
1536 gmx_ana_indexmap_deinit(gmx_ana_indexmap_t *m)
1539 if (m->mapid != m->orgid)
1543 if (m->mapb.nalloc_index > 0)
1545 sfree(m->mapb.index);
1547 if (m->mapb.nalloc_a > 0)
1552 if (m->b.nalloc_index > 0)
1556 if (m->b.nalloc_a > 0)
1560 gmx_ana_indexmap_clear(m);