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