src/gromacs/selection/poscalc.cpp

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