src/gromacs/selection/poscalc.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 gmx::PositionCalculationCollection and functions in poscalc.h.
  40  *
  41  * \todo
  42  * There is probably some room for optimization in the calculation of
  43  * positions with bases.
  44  * In particular, the current implementation may do a lot of unnecessary
  45  * copying.
  46  * The interface would need to be changed to make it possible to use the
  47  * same output positions for several calculations.
  48  *
  49  * \todo
  50  * The current algorithm for setting up base calculations could be improved
  51  * in cases when there are calculations that cannot use a common base but
  52  * still overlap partially (e.g., with three calculations A, B, and C
  53  * such that A could use both B and C as a base, but B and C cannot use the
  54  * same base).
  55  * Setting up the bases in an optimal manner in every possible situation can
  56  * be quite difficult unless several bases are allowed for one calculation,
  57  * but better heuristics could probably be implemented.
  58  * For best results, the setup should probably be postponed (at least
  59  * partially) to gmx_ana_poscalc_init_eval().
  60  *
  61  * \author Teemu Murtola <teemu.murtola@gmail.com>
  62  * \ingroup module_selection
  63  */
  64 #include "gmxpre.h"
  65
  66 #include "poscalc.h"
  67
  68 #include <cstring>
  69
  70 #include <algorithm>
  71 #include <vector>
  72
  73 #include "gromacs/math/vec.h"
  74 #include "gromacs/selection/indexutil.h"
  75 #include "gromacs/trajectory/trajectoryframe.h"
  76 #include "gromacs/utility/arrayref.h"
  77 #include "gromacs/utility/exceptions.h"
  78 #include "gromacs/utility/gmxassert.h"
  79 #include "gromacs/utility/smalloc.h"
  80
  81 #include "centerofmass.h"
  82 #include "position.h"
  83
  84 namespace gmx
  85 {
  86
  87 /*! \internal \brief
  88  * Private implementation class for PositionCalculationCollection.
  89  *
  90  * \ingroup module_selection
  91  */
  92 class PositionCalculationCollection::Impl
  93 {
  94 public:
  95     Impl();
  96     ~Impl();
  97
  98     /*! \brief
  99      * Inserts a position calculation structure into its collection.
 100      *
 101      * \param pc     Data structure to insert.
 102      * \param before Data structure before which to insert
 103      *   (NULL = insert at end).
 104      *
 105      * Inserts \p pc to its collection before \p before.
 106      * If \p before is NULL, \p pc is appended to the list.
 107      */
 108     void insertCalculation(gmx_ana_poscalc_t* pc, gmx_ana_poscalc_t* before);
 109     /*! \brief
 110      * Removes a position calculation structure from its collection.
 111      *
 112      * \param pc    Data structure to remove.
 113      *
 114      * Removes \p pc from its collection.
 115      */
 116     void removeCalculation(gmx_ana_poscalc_t* pc);
 117
 118     /*! \copydoc PositionCalculationCollection::createCalculation()
 119      *
 120      * This method contains the actual implementation of the similarly
 121      * named method in the parent class.
 122      */
 123     gmx_ana_poscalc_t* createCalculation(e_poscalc_t type, int flags);
 124
 125     /*! \brief
 126      * Maps given topology indices into frame indices.
 127      *
 128      * Only one position calculation at a time needs to access this (and
 129      * there are also other thread-unsafe constructs here), so a temporary
 130      * array is used to avoid repeated memory allocation.
 131      */
 132     ArrayRef<const int> getFrameIndices(int size, int index[])
 133     {
 134         if (mapToFrameAtoms_.empty())
 135         {
 136             return constArrayRefFromArray(index, size);
 137         }
 138         tmpFrameAtoms_.resize(size);
 139         for (int i = 0; i < size; ++i)
 140         {
 141             const int ii = index[i];
 142             GMX_ASSERT(ii >= 0 && ii <= ssize(mapToFrameAtoms_) && mapToFrameAtoms_[ii] != -1,
 143                        "Invalid input atom index");
 144             tmpFrameAtoms_[i] = mapToFrameAtoms_[ii];
 145         }
 146         return tmpFrameAtoms_;
 147     }
 148
 149     /*! \brief
 150      * Topology data.
 151      *
 152      * Can be NULL if none of the calculations require topology data or if
 153      * setTopology() has not been called.
 154      */
 155     const gmx_mtop_t* top_;
 156     //! Pointer to the first data structure.
 157     gmx_ana_poscalc_t* first_;
 158     //! Pointer to the last data structure.
 159     gmx_ana_poscalc_t* last_;
 160     //! Whether the collection has been initialized for evaluation.
 161     bool bInit_;
 162     //! Mapping from topology atoms to frame atoms (one index for each topology atom).
 163     std::vector<int> mapToFrameAtoms_;
 164     //! Working array for updating positions.
 165     std::vector<int> tmpFrameAtoms_;
 166 };
 167
 168 } // namespace gmx
 169
 170 /*! \internal \brief
 171  * Data structure for position calculation.
 172  */
 173 struct gmx_ana_poscalc_t
 174 {
 175     /*! \brief
 176      * Type of calculation.
 177      *
 178      * This field may differ from the type requested by the user, because
 179      * it is changed internally to the most effective calculation.
 180      * For example, if the user requests a COM calculation for residues
 181      * consisting of single atoms, it is simply set to POS_ATOM.
 182      * To provide a consistent interface to the user, the field \p itype
 183      * should be used when information should be given out.
 184      */
 185     e_poscalc_t type;
 186     /*! \brief
 187      * Flags for calculation options.
 188      *
 189      * See \ref poscalc_flags "documentation of the flags".
 190      */
 191     int flags;
 192
 193     /*! \brief
 194      * Type for the created indices.
 195      *
 196      * This field always agrees with the type that the user requested, but
 197      * may differ from \p type.
 198      */
 199     e_index_t itype;
 200     /*! \brief
 201      * Block data for the calculation.
 202      */
 203     t_blocka b;
 204     /*! \brief
 205      * Mapping from the blocks to the blocks of \p sbase.
 206      *
 207      * If \p sbase is NULL, this field is also.
 208      */
 209     int* baseid;
 210     /*! \brief
 211      * Maximum evaluation group.
 212      */
 213     gmx_ana_index_t gmax;
 214
 215     /** Position storage for calculations that are used as a base. */
 216     gmx_ana_pos_t* p;
 217
 218     /** true if the positions have been evaluated for the current frame. */
 219     bool bEval;
 220     /*! \brief
 221      * Base position data for this calculation.
 222      *
 223      * If not NULL, the centers required by this calculation have
 224      * already been calculated in \p sbase.
 225      * The structure pointed by \p sbase is always a static calculation.
 226      */
 227     gmx_ana_poscalc_t* sbase;
 228     /** Next structure in the linked list of calculations. */
 229     gmx_ana_poscalc_t* next;
 230     /** Previous structure in the linked list of calculations. */
 231     gmx_ana_poscalc_t* prev;
 232     /** Number of references to this structure. */
 233     int refcount;
 234     /** Collection this calculation belongs to. */
 235     gmx::PositionCalculationCollection::Impl* coll;
 236 };
 237
 238 const char* const gmx::PositionCalculationCollection::typeEnumValues[] = {
 239     "atom",          "res_com",       "res_cog",       "mol_com",       "mol_cog",
 240     "whole_res_com", "whole_res_cog", "whole_mol_com", "whole_mol_cog", "part_res_com",
 241     "part_res_cog",  "part_mol_com",  "part_mol_cog",  "dyn_res_com",   "dyn_res_cog",
 242     "dyn_mol_com",   "dyn_mol_cog",   nullptr,
 243 };
 244
 245 /*! \brief
 246  * Returns the partition type for a given position type.
 247  *
 248  * \param [in] type  \c e_poscalc_t value to convert.
 249  * \returns    Corresponding \c e_indet_t.
 250  */
 251 static e_index_t index_type_for_poscalc(e_poscalc_t type)
 252 {
 253     switch (type)
 254     {
 255         case POS_ATOM: return INDEX_ATOM;
 256         case POS_RES: return INDEX_RES;
 257         case POS_MOL: return INDEX_MOL;
 258         case POS_ALL: // Intended fall through
 259         case POS_ALL_PBC: return INDEX_ALL;
 260     }
 261     return INDEX_UNKNOWN;
 262 }
 263
 264 namespace gmx
 265 {
 266
 267 namespace
 268 {
 269
 270 //! Helper function for determining required topology information.
 271 PositionCalculationCollection::RequiredTopologyInfo requiredTopologyInfo(e_poscalc_t type, int flags)
 272 {
 273     if (type != POS_ATOM)
 274     {
 275         if ((flags & POS_MASS) || (flags & POS_FORCES))
 276         {
 277             return PositionCalculationCollection::RequiredTopologyInfo::TopologyAndMasses;
 278         }
 279         if (type == POS_RES || type == POS_MOL)
 280         {
 281             return PositionCalculationCollection::RequiredTopologyInfo::Topology;
 282         }
 283     }
 284     return PositionCalculationCollection::RequiredTopologyInfo::None;
 285 }
 286
 287 } // namespace
 288
 289 // static
 290 void PositionCalculationCollection::typeFromEnum(const char* post, e_poscalc_t* type, int* flags)
 291 {
 292     if (post[0] == 'a')
 293     {
 294         *type = POS_ATOM;
 295         *flags &= ~(POS_MASS | POS_COMPLMAX | POS_COMPLWHOLE);
 296         return;
 297     }
 298
 299     /* Process the prefix */
 300     const char* ptr = post;
 301     if (post[0] == 'w')
 302     {
 303         *flags &= ~POS_COMPLMAX;
 304         *flags |= POS_COMPLWHOLE;
 305         ptr = post + 6;
 306     }
 307     else if (post[0] == 'p')
 308     {
 309         *flags &= ~POS_COMPLWHOLE;
 310         *flags |= POS_COMPLMAX;
 311         ptr = post + 5;
 312     }
 313     else if (post[0] == 'd')
 314     {
 315         *flags &= ~(POS_COMPLMAX | POS_COMPLWHOLE);
 316         ptr = post + 4;
 317     }
 318
 319     if (ptr[0] == 'r')
 320     {
 321         *type = POS_RES;
 322     }
 323     else if (ptr[0] == 'm')
 324     {
 325         *type = POS_MOL;
 326     }
 327     else
 328     {
 329         GMX_THROW(InternalError("Unknown position calculation type"));
 330     }
 331     if (strlen(ptr) < 7)
 332     {
 333         GMX_THROW(InternalError("Unknown position calculation type"));
 334     }
 335     if (ptr[6] == 'm')
 336     {
 337         *flags |= POS_MASS;
 338     }
 339     else if (ptr[6] == 'g')
 340     {
 341         *flags &= ~POS_MASS;
 342     }
 343     else
 344     {
 345         GMX_THROW(InternalError("Unknown position calculation type"));
 346     }
 347 }
 348
 349 // static
 350 PositionCalculationCollection::RequiredTopologyInfo
 351 PositionCalculationCollection::requiredTopologyInfoForType(const char* post, bool forces)
 352 {
 353     e_poscalc_t type;
 354     int         flags = (forces ? POS_FORCES : 0);
 355     PositionCalculationCollection::typeFromEnum(post, &type, &flags);
 356     return requiredTopologyInfo(type, flags);
 357 }
 358
 359 /********************************************************************
 360  * PositionCalculationCollection::Impl
 361  */
 362
 363 PositionCalculationCollection::Impl::Impl() :
 364     top_(nullptr), first_(nullptr), last_(nullptr), bInit_(false)
 365 {
 366 }
 367
 368 PositionCalculationCollection::Impl::~Impl()
 369 {
 370     // Loop backwards, because there can be internal references in that are
 371     // correctly handled by this direction.
 372     while (last_ != nullptr)
 373     {
 374         GMX_ASSERT(last_->refcount == 1, "Dangling references to position calculations");
 375         gmx_ana_poscalc_free(last_);
 376     }
 377 }
 378
 379 void PositionCalculationCollection::Impl::insertCalculation(gmx_ana_poscalc_t* pc, gmx_ana_poscalc_t* before)
 380 {
 381     GMX_RELEASE_ASSERT(pc->coll == this, "Inconsistent collections");
 382     if (before == nullptr)
 383     {
 384         pc->next = nullptr;
 385         pc->prev = last_;
 386         if (last_ != nullptr)
 387         {
 388             last_->next = pc;
 389         }
 390         last_ = pc;
 391     }
 392     else
 393     {
 394         pc->prev = before->prev;
 395         pc->next = before;
 396         if (before->prev)
 397         {
 398             before->prev->next = pc;
 399         }
 400         before->prev = pc;
 401     }
 402     if (pc->prev == nullptr)
 403     {
 404         first_ = pc;
 405     }
 406 }
 407
 408 void PositionCalculationCollection::Impl::removeCalculation(gmx_ana_poscalc_t* pc)
 409 {
 410     GMX_RELEASE_ASSERT(pc->coll == this, "Inconsistent collections");
 411     if (pc->prev != nullptr)
 412     {
 413         pc->prev->next = pc->next;
 414     }
 415     else if (pc == first_)
 416     {
 417         first_ = pc->next;
 418     }
 419     if (pc->next != nullptr)
 420     {
 421         pc->next->prev = pc->prev;
 422     }
 423     else if (pc == last_)
 424     {
 425         last_ = pc->prev;
 426     }
 427     pc->prev = pc->next = nullptr;
 428 }
 429
 430 gmx_ana_poscalc_t* PositionCalculationCollection::Impl::createCalculation(e_poscalc_t type, int flags)
 431 {
 432     gmx_ana_poscalc_t* pc;
 433
 434     snew(pc, 1);
 435     pc->type  = type;
 436     pc->itype = index_type_for_poscalc(type);
 437     gmx_ana_poscalc_set_flags(pc, flags);
 438     pc->refcount = 1;
 439     pc->coll     = this;
 440     insertCalculation(pc, nullptr);
 441     return pc;
 442 }
 443
 444
 445 /********************************************************************
 446  * PositionCalculationCollection
 447  */
 448
 449 PositionCalculationCollection::PositionCalculationCollection() : impl_(new Impl) {}
 450
 451 PositionCalculationCollection::~PositionCalculationCollection() {}
 452
 453 void PositionCalculationCollection::setTopology(const gmx_mtop_t* top)
 454 {
 455     impl_->top_ = top;
 456 }
 457
 458 void PositionCalculationCollection::printTree(FILE* fp) const
 459 {
 460     gmx_ana_poscalc_t* pc;
 461     int                i, j;
 462
 463     fprintf(fp, "Position calculations:\n");
 464     i  = 1;
 465     pc = impl_->first_;
 466     while (pc)
 467     {
 468         fprintf(fp, "%2d ", i);
 469         switch (pc->type)
 470         {
 471             case POS_ATOM: fprintf(fp, "ATOM"); break;
 472             case POS_RES: fprintf(fp, "RES"); break;
 473             case POS_MOL: fprintf(fp, "MOL"); break;
 474             case POS_ALL: fprintf(fp, "ALL"); break;
 475             case POS_ALL_PBC: fprintf(fp, "ALL_PBC"); break;
 476         }
 477         if (pc->itype != index_type_for_poscalc(pc->type))
 478         {
 479             fprintf(fp, " (");
 480             switch (pc->itype)
 481             {
 482                 case INDEX_UNKNOWN: fprintf(fp, "???"); break;
 483                 case INDEX_ATOM: fprintf(fp, "ATOM"); break;
 484                 case INDEX_RES: fprintf(fp, "RES"); break;
 485                 case INDEX_MOL: fprintf(fp, "MOL"); break;
 486                 case INDEX_ALL: fprintf(fp, "ALL"); break;
 487             }
 488             fprintf(fp, ")");
 489         }
 490         fprintf(fp, " flg=");
 491         if (pc->flags & POS_MASS)
 492         {
 493             fprintf(fp, "M");
 494         }
 495         if (pc->flags & POS_DYNAMIC)
 496         {
 497             fprintf(fp, "D");
 498         }
 499         if (pc->flags & POS_MASKONLY)
 500         {
 501             fprintf(fp, "A");
 502         }
 503         if (pc->flags & POS_COMPLMAX)
 504         {
 505             fprintf(fp, "Cm");
 506         }
 507         if (pc->flags & POS_COMPLWHOLE)
 508         {
 509             fprintf(fp, "Cw");
 510         }
 511         if (!pc->flags)
 512         {
 513             fprintf(fp, "0");
 514         }
 515         fprintf(fp, " nr=%d nra=%d", pc->b.nr, pc->b.nra);
 516         fprintf(fp, " refc=%d", pc->refcount);
 517         fprintf(fp, "\n");
 518         if (pc->gmax.nalloc_index > 0)
 519         {
 520             fprintf(fp, "   Group: ");
 521             if (pc->gmax.isize > 20)
 522             {
 523                 fprintf(fp, " %d atoms", pc->gmax.isize);
 524             }
 525             else
 526             {
 527                 for (j = 0; j < pc->gmax.isize; ++j)
 528                 {
 529                     fprintf(fp, " %d", pc->gmax.index[j] + 1);
 530                 }
 531             }
 532             fprintf(fp, "\n");
 533         }
 534         if (pc->b.nalloc_a > 0)
 535         {
 536             fprintf(fp, "   Atoms: ");
 537             if (pc->b.nra > 20)
 538             {
 539                 fprintf(fp, " %d atoms", pc->b.nra);
 540             }
 541             else
 542             {
 543                 for (j = 0; j < pc->b.nra; ++j)
 544                 {
 545                     fprintf(fp, " %d", pc->b.a[j] + 1);
 546                 }
 547             }
 548             fprintf(fp, "\n");
 549         }
 550         if (pc->b.nalloc_index > 0)
 551         {
 552             fprintf(fp, "   Blocks:");
 553             if (pc->b.nr > 20)
 554             {
 555                 fprintf(fp, " %d pcs", pc->b.nr);
 556             }
 557             else
 558             {
 559                 for (j = 0; j <= pc->b.nr; ++j)
 560                 {
 561                     fprintf(fp, " %d", pc->b.index[j]);
 562                 }
 563             }
 564             fprintf(fp, "\n");
 565         }
 566         if (pc->sbase)
 567         {
 568             gmx_ana_poscalc_t* base;
 569
 570             fprintf(fp, "   Base: ");
 571             j    = 1;
 572             base = impl_->first_;
 573             while (base && base != pc->sbase)
 574             {
 575                 ++j;
 576                 base = base->next;
 577             }
 578             fprintf(fp, "%d", j);
 579             if (pc->baseid && pc->b.nr <= 20)
 580             {
 581                 fprintf(fp, " id:");
 582                 for (j = 0; j < pc->b.nr; ++j)
 583                 {
 584                     fprintf(fp, " %d", pc->baseid[j] + 1);
 585                 }
 586             }
 587             fprintf(fp, "\n");
 588         }
 589         ++i;
 590         pc = pc->next;
 591     }
 592 }
 593
 594 gmx_ana_poscalc_t* PositionCalculationCollection::createCalculation(e_poscalc_t type, int flags)
 595 {
 596     return impl_->createCalculation(type, flags);
 597 }
 598
 599 gmx_ana_poscalc_t* PositionCalculationCollection::createCalculationFromEnum(const char* post, int flags)
 600 {
 601     e_poscalc_t type;
 602     int         cflags = flags;
 603     typeFromEnum(post, &type, &cflags);
 604     return impl_->createCalculation(type, cflags);
 605 }
 606
 607 void PositionCalculationCollection::getRequiredAtoms(gmx_ana_index_t* out) const
 608 {
 609     gmx_ana_poscalc_t* pc = impl_->first_;
 610     while (pc)
 611     {
 612         // Calculations with a base just copy positions from the base, so
 613         // those do not need to be considered in the check.
 614         if (!pc->sbase)
 615         {
 616             gmx_ana_index_t g;
 617             gmx_ana_index_set(&g, pc->b.nra, pc->b.a, 0);
 618             gmx_ana_index_union_unsorted(out, out, &g);
 619         }
 620         pc = pc->next;
 621     }
 622 }
 623
 624 void PositionCalculationCollection::initEvaluation()
 625 {
 626     if (impl_->bInit_)
 627     {
 628         return;
 629     }
 630     gmx_ana_poscalc_t* pc = impl_->first_;
 631     while (pc)
 632     {
 633         /* Initialize position storage for base calculations */
 634         if (pc->p)
 635         {
 636             gmx_ana_poscalc_init_pos(pc, pc->p);
 637         }
 638         /* Construct the mapping of the base positions */
 639         if (pc->sbase)
 640         {
 641             int bi, bj;
 642
 643             snew(pc->baseid, pc->b.nr);
 644             for (bi = bj = 0; bi < pc->b.nr; ++bi, ++bj)
 645             {
 646                 while (pc->sbase->b.a[pc->sbase->b.index[bj]] != pc->b.a[pc->b.index[bi]])
 647                 {
 648                     ++bj;
 649                 }
 650                 pc->baseid[bi] = bj;
 651             }
 652         }
 653         /* Free the block data for dynamic calculations */
 654         if (pc->flags & POS_DYNAMIC)
 655         {
 656             if (pc->b.nalloc_index > 0)
 657             {
 658                 sfree(pc->b.index);
 659                 pc->b.nalloc_index = 0;
 660             }
 661             if (pc->b.nalloc_a > 0)
 662             {
 663                 sfree(pc->b.a);
 664                 pc->b.nalloc_a = 0;
 665             }
 666         }
 667         pc = pc->next;
 668     }
 669     impl_->bInit_ = true;
 670 }
 671
 672 void PositionCalculationCollection::initFrame(const t_trxframe* fr)
 673 {
 674     if (!impl_->bInit_)
 675     {
 676         initEvaluation();
 677     }
 678     /* Clear the evaluation flags */
 679     gmx_ana_poscalc_t* pc = impl_->first_;
 680     while (pc)
 681     {
 682         pc->bEval = false;
 683         pc        = pc->next;
 684     }
 685     if (fr->bIndex && fr->natoms > 0)
 686     {
 687         const int highestAtom = *std::max_element(fr->index, fr->index + fr->natoms);
 688         impl_->mapToFrameAtoms_.resize(highestAtom + 1);
 689         std::fill(impl_->mapToFrameAtoms_.begin(), impl_->mapToFrameAtoms_.end(), -1);
 690         for (int i = 0; i < fr->natoms; ++i)
 691         {
 692             impl_->mapToFrameAtoms_[fr->index[i]] = i;
 693         }
 694     }
 695     else
 696     {
 697         impl_->mapToFrameAtoms_.clear();
 698     }
 699 }
 700
 701 } // namespace gmx
 702
 703 /*! \brief
 704  * Initializes position calculation using the maximum possible input index.
 705  *
 706  * \param[in,out] pc  Position calculation data structure.
 707  * \param[in]     g   Maximum index group for the calculation.
 708  * \param[in]     bBase Whether \p pc will be used as a base or not.
 709  *
 710  * \p bBase affects on how the \p pc->gmax field is initialized.
 711  */
 712 static void set_poscalc_maxindex(gmx_ana_poscalc_t* pc, gmx_ana_index_t* g, bool bBase)
 713 {
 714     const gmx_mtop_t* top = pc->coll->top_;
 715     gmx_ana_index_make_block(&pc->b, top, g, pc->itype, (pc->flags & POS_COMPLWHOLE) != 0);
 716     /* Set the type to POS_ATOM if the calculation in fact is such. */
 717     if (pc->b.nr == pc->b.nra)
 718     {
 719         pc->type = POS_ATOM;
 720         pc->flags &= ~(POS_MASS | POS_COMPLMAX | POS_COMPLWHOLE);
 721     }
 722     /* Set the POS_COMPLWHOLE flag if the calculation in fact always uses
 723      * complete residues and molecules. */
 724     if (!(pc->flags & POS_COMPLWHOLE) && (!(pc->flags & POS_DYNAMIC) || (pc->flags & POS_COMPLMAX))
 725         && (pc->type == POS_RES || pc->type == POS_MOL)
 726         && gmx_ana_index_has_complete_elems(g, pc->itype, top))
 727     {
 728         pc->flags &= ~POS_COMPLMAX;
 729         pc->flags |= POS_COMPLWHOLE;
 730     }
 731     /* Setup the gmax field */
 732     if ((pc->flags & POS_COMPLWHOLE) && !bBase && pc->b.nra > g->isize)
 733     {
 734         gmx_ana_index_copy(&pc->gmax, g, true);
 735     }
 736     else
 737     {
 738         gmx_ana_index_set(&pc->gmax, pc->b.nra, pc->b.a, 0);
 739     }
 740 }
 741
 742 /*! \brief
 743  * Checks whether a position calculation should use a base at all.
 744  *
 745  * \param[in] pc   Position calculation data to check.
 746  * \returns   true if \p pc can use a base and gets some benefit out of it,
 747  *   false otherwise.
 748  */
 749 static bool can_use_base(gmx_ana_poscalc_t* pc)
 750 {
 751     /* For atoms, it should be faster to do a simple copy, so don't use a
 752      * base. */
 753     if (pc->type == POS_ATOM)
 754     {
 755         return false;
 756     }
 757     /* For dynamic selections that do not use completion, it is not possible
 758      * to use a base. */
 759     if ((pc->type == POS_RES || pc->type == POS_MOL) && (pc->flags & POS_DYNAMIC)
 760         && !(pc->flags & (POS_COMPLMAX | POS_COMPLWHOLE)))
 761     {
 762         return false;
 763     }
 764     /* Dynamic calculations for a single position cannot use a base. */
 765     if ((pc->type == POS_ALL || pc->type == POS_ALL_PBC) && (pc->flags & POS_DYNAMIC))
 766     {
 767         return false;
 768     }
 769     return true;
 770 }
 771
 772 /*! \brief
 773  * Checks whether two position calculations should use a common base.
 774  *
 775  * \param[in]     pc1 Calculation 1 to check for.
 776  * \param[in]     pc2 Calculation 2 to check for.
 777  * \param[in]     g1  Index group structure that contains the atoms from
 778  *   \p pc1.
 779  * \param[in,out] g   Working space, should have enough allocated memory to
 780  *   contain the intersection of the atoms in \p pc1 and \p pc2.
 781  * \returns   true if the two calculations should be merged to use a common
 782  *   base, false otherwise.
 783  */
 784 static bool should_merge(gmx_ana_poscalc_t* pc1, gmx_ana_poscalc_t* pc2, gmx_ana_index_t* g1, gmx_ana_index_t* g)
 785 {
 786     gmx_ana_index_t g2;
 787
 788     /* Do not merge calculations with different mass weighting. */
 789     if ((pc1->flags & POS_MASS) != (pc2->flags & POS_MASS))
 790     {
 791         return false;
 792     }
 793     /* Avoid messing up complete calculations. */
 794     if ((pc1->flags & POS_COMPLWHOLE) != (pc2->flags & POS_COMPLWHOLE))
 795     {
 796         return false;
 797     }
 798     /* Find the overlap between the calculations. */
 799     gmx_ana_index_set(&g2, pc2->b.nra, pc2->b.a, 0);
 800     gmx_ana_index_intersection(g, g1, &g2);
 801     /* Do not merge if there is no overlap. */
 802     if (g->isize == 0)
 803     {
 804         return false;
 805     }
 806     /* Full completion calculations always match if the type is correct. */
 807     if ((pc1->flags & POS_COMPLWHOLE) && (pc2->flags & POS_COMPLWHOLE) && pc1->type == pc2->type)
 808     {
 809         return true;
 810     }
 811     /* The calculations also match if the intersection consists of full
 812      * blocks. */
 813     if (gmx_ana_index_has_full_ablocks(g, &pc1->b) && gmx_ana_index_has_full_ablocks(g, &pc2->b))
 814     {
 815         return true;
 816     }
 817     return false;
 818 }
 819
 820 /*! \brief
 821  * Creates a static base for position calculation.
 822  *
 823  * \param     pc  Data structure to copy.
 824  * \returns   Pointer to a newly allocated base for \p pc.
 825  *
 826  * Creates and returns a deep copy of \p pc, but clears the
 827  * \ref POS_DYNAMIC and \ref POS_MASKONLY flags.
 828  * The newly created structure is set as the base (\c gmx_ana_poscalc_t::sbase)
 829  * of \p pc and inserted in the collection before \p pc.
 830  */
 831 static gmx_ana_poscalc_t* create_simple_base(gmx_ana_poscalc_t* pc)
 832 {
 833     gmx_ana_poscalc_t* base;
 834     int                flags;
 835
 836     flags = pc->flags & ~(POS_DYNAMIC | POS_MASKONLY);
 837     base  = pc->coll->createCalculation(pc->type, flags);
 838     set_poscalc_maxindex(base, &pc->gmax, true);
 839
 840     base->p = new gmx_ana_pos_t();
 841
 842     pc->sbase = base;
 843     pc->coll->removeCalculation(base);
 844     pc->coll->insertCalculation(base, pc);
 845
 846     return base;
 847 }
 848
 849 /*! \brief
 850  * Merges a calculation into another calculation such that the new calculation
 851  * can be used as a base.
 852  *
 853  * \param[in,out] base Base calculation to merge to.
 854  * \param[in,out] pc   Position calculation to merge to \p base.
 855  *
 856  * After the call, \p base can be used as a base for \p pc (or any calculation
 857  * that used it as a base).
 858  * It is assumed that any overlap between \p base and \p pc is in complete
 859  * blocks, i.e., that the merge is possible.
 860  */
 861 static void merge_to_base(gmx_ana_poscalc_t* base, gmx_ana_poscalc_t* pc)
 862 {
 863     gmx_ana_index_t gp, gb, g;
 864     int             isize, bnr;
 865     int             i, bi, bj, bo;
 866
 867     base->flags |= pc->flags & (POS_VELOCITIES | POS_FORCES);
 868     gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, 0);
 869     gmx_ana_index_set(&gb, base->b.nra, base->b.a, 0);
 870     isize = gmx_ana_index_difference_size(&gp, &gb);
 871     if (isize > 0)
 872     {
 873         gmx_ana_index_clear(&g);
 874         gmx_ana_index_reserve(&g, base->b.nra + isize);
 875         /* Find the new blocks */
 876         gmx_ana_index_difference(&g, &gp, &gb);
 877         /* Count the blocks in g */
 878         i = bi = bnr = 0;
 879         while (i < g.isize)
 880         {
 881             while (pc->b.a[pc->b.index[bi]] != g.index[i])
 882             {
 883                 ++bi;
 884             }
 885             i += pc->b.index[bi + 1] - pc->b.index[bi];
 886             ++bnr;
 887             ++bi;
 888         }
 889         /* Merge the atoms into a temporary structure */
 890         gmx_ana_index_merge(&g, &gb, &g);
 891         /* Merge the blocks */
 892         srenew(base->b.index, base->b.nr + bnr + 1);
 893         i                     = g.isize - 1;
 894         bi                    = base->b.nr - 1;
 895         bj                    = pc->b.nr - 1;
 896         bo                    = base->b.nr + bnr - 1;
 897         base->b.index[bo + 1] = i + 1;
 898         while (bo >= 0)
 899         {
 900             if (bi < 0 || base->b.a[base->b.index[bi + 1] - 1] != g.index[i])
 901             {
 902                 i -= pc->b.index[bj + 1] - pc->b.index[bj];
 903                 --bj;
 904             }
 905             else
 906             {
 907                 if (bj >= 0 && pc->b.a[pc->b.index[bj + 1] - 1] == g.index[i])
 908                 {
 909                     --bj;
 910                 }
 911                 i -= base->b.index[bi + 1] - base->b.index[bi];
 912                 --bi;
 913             }
 914             base->b.index[bo] = i + 1;
 915             --bo;
 916         }
 917         base->b.nr += bnr;
 918         base->b.nalloc_index += bnr;
 919         sfree(base->b.a);
 920         base->b.nra      = g.isize;
 921         base->b.a        = g.index;
 922         base->b.nalloc_a = g.isize;
 923         /* Refresh the gmax field */
 924         gmx_ana_index_set(&base->gmax, base->b.nra, base->b.a, 0);
 925     }
 926 }
 927
 928 /*! \brief
 929  * Merges two bases into one.
 930  *
 931  * \param[in,out] tbase Base calculation to merge to.
 932  * \param[in]     mbase Base calculation to merge to \p tbase.
 933  *
 934  * After the call, \p mbase has been freed and \p tbase is used as the base
 935  * for all calculations that previously had \p mbase as their base.
 936  * It is assumed that any overlap between \p tbase and \p mbase is in complete
 937  * blocks, i.e., that the merge is possible.
 938  */
 939 static void merge_bases(gmx_ana_poscalc_t* tbase, gmx_ana_poscalc_t* mbase)
 940 {
 941     gmx_ana_poscalc_t* pc;
 942
 943     merge_to_base(tbase, mbase);
 944     mbase->coll->removeCalculation(mbase);
 945     /* Set tbase as the base for all calculations that had mbase */
 946     pc = tbase->coll->first_;
 947     while (pc)
 948     {
 949         if (pc->sbase == mbase)
 950         {
 951             pc->sbase = tbase;
 952             tbase->refcount++;
 953         }
 954         pc = pc->next;
 955     }
 956     /* Free mbase */
 957     mbase->refcount = 0;
 958     gmx_ana_poscalc_free(mbase);
 959 }
 960
 961 /*! \brief
 962  * Setups the static base calculation for a position calculation.
 963  *
 964  * \param[in,out] pc   Position calculation to setup the base for.
 965  */
 966 static void setup_base(gmx_ana_poscalc_t* pc)
 967 {
 968     gmx_ana_poscalc_t *base, *pbase, *next;
 969     gmx_ana_index_t    gp, g;
 970
 971     /* Exit immediately if pc should not have a base. */
 972     if (!can_use_base(pc))
 973     {
 974         return;
 975     }
 976
 977     gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, 0);
 978     gmx_ana_index_clear(&g);
 979     gmx_ana_index_reserve(&g, pc->b.nra);
 980     pbase = pc;
 981     base  = pc->coll->first_;
 982     while (base)
 983     {
 984         /* Save the next calculation so that we can safely delete base */
 985         next = base->next;
 986         /* Skip pc, calculations that already have a base (we should match the
 987          * base instead), as well as calculations that should not have a base.
 988          * If the above conditions are met, check whether we should do a
 989          * merge.
 990          */
 991         if (base != pc && !base->sbase && can_use_base(base) && should_merge(pbase, base, &gp, &g))
 992         {
 993             /* Check whether this is the first base found */
 994             if (pbase == pc)
 995             {
 996                 /* Create a real base if one is not present */
 997                 if (!base->p)
 998                 {
 999                     pbase = create_simple_base(base);
1000                 }
1001                 else
1002                 {
1003                     pbase = base;
1004                 }
1005                 /* Make it a base for pc as well */
1006                 merge_to_base(pbase, pc);
1007                 pc->sbase = pbase;
1008                 pbase->refcount++;
1009             }
1010             else /* This was not the first base */
1011             {
1012                 if (!base->p)
1013                 {
1014                     /* If it is not a real base, just make the new base as
1015                      * the base for it as well. */
1016                     merge_to_base(pbase, base);
1017                     base->sbase = pbase;
1018                     pbase->refcount++;
1019                 }
1020                 else
1021                 {
1022                     /* If base is a real base, merge it with the new base
1023                      * and delete it. */
1024                     merge_bases(pbase, base);
1025                 }
1026             }
1027             gmx_ana_index_set(&gp, pbase->b.nra, pbase->b.a, 0);
1028             gmx_ana_index_reserve(&g, pc->b.nra);
1029         }
1030         /* Proceed to the next unchecked calculation */
1031         base = next;
1032     }
1033
1034     gmx_ana_index_deinit(&g);
1035
1036     /* If no base was found, create one if one is required */
1037     if (!pc->sbase && (pc->flags & POS_DYNAMIC) && (pc->flags & (POS_COMPLMAX | POS_COMPLWHOLE)))
1038     {
1039         create_simple_base(pc);
1040     }
1041 }
1042
1043 /*!
1044  * \param[in,out] pc    Position calculation data structure.
1045  * \param[in]     flags New flags.
1046  *
1047  * \p flags are added to the old flags.
1048  * If calculation type is \ref POS_ATOM, \ref POS_MASS is automatically
1049  * cleared.
1050  * If both \ref POS_DYNAMIC and \ref POS_MASKONLY are provided,
1051  * \ref POS_DYNAMIC is cleared.
1052  * If calculation type is not \ref POS_RES or \ref POS_MOL,
1053  * \ref POS_COMPLMAX and \ref POS_COMPLWHOLE are automatically cleared.
1054  */
1055 void gmx_ana_poscalc_set_flags(gmx_ana_poscalc_t* pc, int flags)
1056 {
1057     if (pc->type == POS_ATOM)
1058     {
1059         flags &= ~POS_MASS;
1060     }
1061     if (flags & POS_MASKONLY)
1062     {
1063         flags &= ~POS_DYNAMIC;
1064     }
1065     if (pc->type != POS_RES && pc->type != POS_MOL)
1066     {
1067         flags &= ~(POS_COMPLMAX | POS_COMPLWHOLE);
1068     }
1069     pc->flags |= flags;
1070 }
1071
1072 /*!
1073  * \param[in,out] pc  Position calculation data structure.
1074  * \param[in]     g   Maximum index group for the calculation.
1075  *
1076  * Subsequent calls to gmx_ana_poscalc_update() should use only subsets of
1077  * \p g for evaluation.
1078  *
1079  * The topology should have been set for the collection of which \p pc is
1080  * a member.
1081  */
1082 void gmx_ana_poscalc_set_maxindex(gmx_ana_poscalc_t* pc, gmx_ana_index_t* g)
1083 {
1084     set_poscalc_maxindex(pc, g, false);
1085     setup_base(pc);
1086 }
1087
1088 /*!
1089  * \param[in]  pc  Position calculation data structure.
1090  * \param[out] p   Output positions.
1091  *
1092  * Calls to gmx_ana_poscalc_update() using \p pc should use only positions
1093  * initialized with this function.
1094  * The \c p->g pointer is initialized to point to an internal group that
1095  * contains the maximum index group set with gmx_ana_poscalc_set_maxindex().
1096  */
1097 void gmx_ana_poscalc_init_pos(gmx_ana_poscalc_t* pc, gmx_ana_pos_t* p)
1098 {
1099     gmx_ana_indexmap_init(&p->m, &pc->gmax, pc->coll->top_, pc->itype);
1100     /* Only do the static optimization when there is no completion */
1101     if (!(pc->flags & POS_DYNAMIC) && pc->b.nra == pc->gmax.isize)
1102     {
1103         gmx_ana_indexmap_set_static(&p->m, &pc->b);
1104     }
1105     gmx_ana_pos_reserve(p, p->m.mapb.nr, -1);
1106     if (pc->flags & POS_VELOCITIES)
1107     {
1108         gmx_ana_pos_reserve_velocities(p);
1109     }
1110     if (pc->flags & POS_FORCES)
1111     {
1112         gmx_ana_pos_reserve_forces(p);
1113     }
1114 }
1115
1116 /*!
1117  * \param  pc  Position calculation data to be freed.
1118  *
1119  * The \p pc pointer is invalid after the call.
1120  */
1121 void gmx_ana_poscalc_free(gmx_ana_poscalc_t* pc)
1122 {
1123     if (!pc)
1124     {
1125         return;
1126     }
1127
1128     pc->refcount--;
1129     if (pc->refcount > 0)
1130     {
1131         return;
1132     }
1133
1134     pc->coll->removeCalculation(pc);
1135     if (pc->b.nalloc_index > 0)
1136     {
1137         sfree(pc->b.index);
1138     }
1139     if (pc->b.nalloc_a > 0)
1140     {
1141         sfree(pc->b.a);
1142     }
1143     if (pc->flags & POS_COMPLWHOLE)
1144     {
1145         gmx_ana_index_deinit(&pc->gmax);
1146     }
1147     delete pc->p;
1148     if (pc->sbase)
1149     {
1150         gmx_ana_poscalc_free(pc->sbase);
1151         sfree(pc->baseid);
1152     }
1153     sfree(pc);
1154 }
1155
1156 gmx::PositionCalculationCollection::RequiredTopologyInfo gmx_ana_poscalc_required_topology_info(gmx_ana_poscalc_t* pc)
1157 {
1158     return gmx::requiredTopologyInfo(pc->type, pc->flags);
1159 }
1160
1161 /*!
1162  * \param[in]     pc   Position calculation data.
1163  * \param[in,out] p    Output positions, initialized previously with
1164  *   gmx_ana_poscalc_init_pos() using \p pc.
1165  * \param[in]     g    Index group to use for the update.
1166  * \param[in]     fr   Current frame.
1167  * \param[in]     pbc  PBC data, or NULL if no PBC should be used.
1168  *
1169  * gmx_ana_poscalc_init_frame() should be called for each frame before calling
1170  * this function.
1171  */
1172 void gmx_ana_poscalc_update(gmx_ana_poscalc_t* pc,
1173                             gmx_ana_pos_t*     p,
1174                             gmx_ana_index_t*   g,
1175                             t_trxframe*        fr,
1176                             const t_pbc*       pbc)
1177 {
1178     int i, bi, bj;
1179
1180     if (pc->bEval && !(pc->flags & POS_MASKONLY))
1181     {
1182         return;
1183     }
1184     if (pc->sbase)
1185     {
1186         gmx_ana_poscalc_update(pc->sbase, nullptr, nullptr, fr, pbc);
1187     }
1188     if (!p)
1189     {
1190         p = pc->p;
1191     }
1192     if (!g)
1193     {
1194         g = &pc->gmax;
1195     }
1196
1197     /* Update the index map */
1198     if (pc->flags & POS_DYNAMIC)
1199     {
1200         gmx_ana_indexmap_update(&p->m, g, false);
1201     }
1202     else if (pc->flags & POS_MASKONLY)
1203     {
1204         gmx_ana_indexmap_update(&p->m, g, true);
1205         if (pc->bEval)
1206         {
1207             return;
1208         }
1209     }
1210     if (!(pc->flags & POS_DYNAMIC))
1211     {
1212         pc->bEval = true;
1213     }
1214
1215     /* Evaluate the positions */
1216     if (pc->sbase)
1217     {
1218         /* TODO: It might be faster to evaluate the positions within this
1219          * loop instead of in the beginning. */
1220         if (pc->flags & POS_DYNAMIC)
1221         {
1222             for (bi = 0; bi < p->count(); ++bi)
1223             {
1224                 bj = pc->baseid[p->m.refid[bi]];
1225                 copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
1226             }
1227             if (p->v)
1228             {
1229                 for (bi = 0; bi < p->count(); ++bi)
1230                 {
1231                     bj = pc->baseid[p->m.refid[bi]];
1232                     copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
1233                 }
1234             }
1235             if (p->f)
1236             {
1237                 for (bi = 0; bi < p->count(); ++bi)
1238                 {
1239                     bj = pc->baseid[p->m.refid[bi]];
1240                     copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
1241                 }
1242             }
1243         }
1244         else
1245         {
1246             for (bi = 0; bi < p->count(); ++bi)
1247             {
1248                 bj = pc->baseid[bi];
1249                 copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
1250             }
1251             if (p->v)
1252             {
1253                 for (bi = 0; bi < p->count(); ++bi)
1254                 {
1255                     bj = pc->baseid[bi];
1256                     copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
1257                 }
1258             }
1259             if (p->f)
1260             {
1261                 for (bi = 0; bi < p->count(); ++bi)
1262                 {
1263                     bj = pc->baseid[bi];
1264                     copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
1265                 }
1266             }
1267         }
1268     }
1269     else /* pc->sbase is NULL */
1270     {
1271         if (pc->flags & POS_DYNAMIC)
1272         {
1273             pc->b.nr    = p->m.mapb.nr;
1274             pc->b.index = p->m.mapb.index;
1275             pc->b.nra   = g->isize;
1276             pc->b.a     = g->index;
1277         }
1278         if (p->v && !fr->bV)
1279         {
1280             for (i = 0; i < pc->b.nra; ++i)
1281             {
1282                 clear_rvec(p->v[i]);
1283             }
1284         }
1285         if (p->f && !fr->bF)
1286         {
1287             for (i = 0; i < pc->b.nra; ++i)
1288             {
1289                 clear_rvec(p->f[i]);
1290             }
1291         }
1292         gmx::ArrayRef<const int> index = pc->coll->getFrameIndices(pc->b.nra, pc->b.a);
1293         const gmx_mtop_t*        top   = pc->coll->top_;
1294         const bool               bMass = (pc->flags & POS_MASS) != 0;
1295         switch (pc->type)
1296         {
1297             case POS_ATOM:
1298                 for (i = 0; i < pc->b.nra; ++i)
1299                 {
1300                     copy_rvec(fr->x[index[i]], p->x[i]);
1301                 }
1302                 if (p->v && fr->bV)
1303                 {
1304                     for (i = 0; i < pc->b.nra; ++i)
1305                     {
1306                         copy_rvec(fr->v[index[i]], p->v[i]);
1307                     }
1308                 }
1309                 if (p->f && fr->bF)
1310                 {
1311                     for (i = 0; i < pc->b.nra; ++i)
1312                     {
1313                         copy_rvec(fr->f[index[i]], p->f[i]);
1314                     }
1315                 }
1316                 break;
1317             case POS_ALL:
1318                 gmx_calc_comg(top, fr->x, index.size(), index.data(), bMass, p->x[0]);
1319                 if (p->v && fr->bV)
1320                 {
1321                     gmx_calc_comg(top, fr->v, index.size(), index.data(), bMass, p->v[0]);
1322                 }
1323                 if (p->f && fr->bF)
1324                 {
1325                     gmx_calc_comg_f(top, fr->f, index.size(), index.data(), bMass, p->f[0]);
1326                 }
1327                 break;
1328             case POS_ALL_PBC:
1329                 gmx_calc_comg_pbc(top, fr->x, pbc, index.size(), index.data(), bMass, p->x[0]);
1330                 if (p->v && fr->bV)
1331                 {
1332                     gmx_calc_comg(top, fr->v, index.size(), index.data(), bMass, p->v[0]);
1333                 }
1334                 if (p->f && fr->bF)
1335                 {
1336                     gmx_calc_comg_f(top, fr->f, index.size(), index.data(), bMass, p->f[0]);
1337                 }
1338                 break;
1339             default:
1340                 // TODO: It would probably be better to do this without the type casts.
1341                 gmx_calc_comg_block(
1342                         top, fr->x, reinterpret_cast<t_block*>(&pc->b), index.data(), bMass, p->x);
1343                 if (p->v && fr->bV)
1344                 {
1345                     gmx_calc_comg_block(
1346                             top, fr->v, reinterpret_cast<t_block*>(&pc->b), index.data(), bMass, p->v);
1347                 }
1348                 if (p->f && fr->bF)
1349                 {
1350                     gmx_calc_comg_f_block(
1351                             top, fr->f, reinterpret_cast<t_block*>(&pc->b), index.data(), bMass, p->f);
1352                 }
1353                 break;
1354         }
1355     }
1356 }