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