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