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