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