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