src/gromacs/selection/selection.h

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \file
  32  * \brief
  33  * Declares gmx::Selection and supporting classes.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \inpublicapi
  37  * \ingroup module_selection
  38  */
  39 #ifndef GMX_SELECTION_SELECTION_H
  40 #define GMX_SELECTION_SELECTION_H
  41
  42 #include <string>
  43 #include <vector>
  44
  45 #include "../legacyheaders/typedefs.h"
  46
  47 #include "../utility/arrayref.h"
  48 #include "../utility/common.h"
  49 #include "../utility/gmxassert.h"
  50
  51 #include "position.h"
  52 #include "indexutil.h"
  53 #include "selectionenums.h"
  54
  55 namespace gmx
  56 {
  57
  58 class SelectionOptionStorage;
  59 class SelectionTreeElement;
  60
  61 class Selection;
  62 class SelectionPosition;
  63
  64 //! Container of selections used in public selection interfaces.
  65 typedef std::vector<Selection> SelectionList;
  66
  67 namespace internal
  68 {
  69
  70 /*! \internal \brief
  71  * Internal data for a single selection.
  72  *
  73  * This class is internal to the selection module, but resides in a public
  74  * header because of efficiency reasons: it allows frequently used access
  75  * methods in \ref Selection to be inlined.
  76  *
  77  * Methods in this class do not throw unless otherwise specified.
  78  *
  79  * \ingroup module_selection
  80  */
  81 class SelectionData
  82 {
  83     public:
  84         /*! \brief
  85          * Creates a new selection object.
  86          *
  87          * \param[in] elem   Root of the evaluation tree for this selection.
  88          * \param[in] selstr String that was parsed to produce this selection.
  89          * \throws    std::bad_alloc if out of memory.
  90          */
  91         SelectionData(SelectionTreeElement *elem, const char *selstr);
  92         ~SelectionData();
  93
  94         //! Returns the string that was parsed to produce this selection.
  95         const char *selectionText() const { return selectionText_.c_str(); }
  96         //! Returns true if the size of the selection (posCount()) is dynamic.
  97         bool isDynamic() const { return bDynamic_; }
  98         //! Number of positions in the selection.
  99         int posCount() const { return rawPositions_.nr; }
 100         //! Returns the root of the evaluation tree for this selection.
 101         SelectionTreeElement &rootElement() { return rootElement_; }
 102
 103         //! Returns whether the covered fraction can change between frames.
 104         bool isCoveredFractionDynamic() const { return bDynamicCoveredFraction_; }
 105
 106         //! Returns true if the given flag is set.
 107         bool hasFlag(SelectionFlag flag) const { return flags_.test(flag); }
 108         //! Sets the flags for this selection.
 109         void setFlags(SelectionFlags flags) { flags_ = flags; }
 110
 111         //! \copydoc Selection::initCoveredFraction()
 112         bool initCoveredFraction(e_coverfrac_t type);
 113
 114         /*! \brief
 115          * Computes total masses and charges for all selection positions.
 116          *
 117          * \param[in] top   Topology information.
 118          * \throws    std::bad_alloc if out of memory.
 119          *
 120          * For dynamic selections, the values need to be updated after each
 121          * evaluation with refreshMassesAndCharges().
 122          * This is done by SelectionEvaluator.
 123          *
 124          * This function is called by SelectionCompiler.
 125          *
 126          * Strong exception safety guarantee.
 127          */
 128         void initializeMassesAndCharges(const t_topology *top);
 129         /*! \brief
 130          * Updates masses and charges after dynamic selection has been
 131          * evaluated.
 132          *
 133          * \param[in] top   Topology information.
 134          *
 135          * Called by SelectionEvaluator.
 136          */
 137         void refreshMassesAndCharges(const t_topology *top);
 138         /*! \brief
 139          * Updates the covered fraction after a selection has been evaluated.
 140          *
 141          * Called by SelectionEvaluator.
 142          */
 143         void updateCoveredFractionForFrame();
 144         /*! \brief
 145          * Computes average covered fraction after all frames have been evaluated.
 146          *
 147          * \param[in] nframes  Number of frames that have been evaluated.
 148          *
 149          * \p nframes should be equal to the number of calls to
 150          * updateCoveredFractionForFrame().
 151          * Called by SelectionEvaluator::evaluateFinal().
 152          */
 153         void computeAverageCoveredFraction(int nframes);
 154         /*! \brief
 155          * Restores position information to state it was in after compilation.
 156          *
 157          * \param[in] top   Topology information.
 158          *
 159          * Depends on SelectionCompiler storing the original atoms in the
 160          * \a rootElement_ object.
 161          * Called by SelectionEvaluator::evaluateFinal().
 162          */
 163         void restoreOriginalPositions(const t_topology *top);
 164
 165     private:
 166         //! Name of the selection.
 167         std::string               name_;
 168         //! The actual selection string.
 169         std::string               selectionText_;
 170         //! Low-level representation of selected positions.
 171         gmx_ana_pos_t             rawPositions_;
 172         //! Total masses for the current positions.
 173         std::vector<real>         posMass_;
 174         //! Total charges for the current positions.
 175         std::vector<real>         posCharge_;
 176         SelectionFlags            flags_;
 177         //! Root of the selection evaluation tree.
 178         SelectionTreeElement     &rootElement_;
 179         //! Type of the covered fraction.
 180         e_coverfrac_t             coveredFractionType_;
 181         //! Covered fraction of the selection for the current frame.
 182         real                      coveredFraction_;
 183         //! The average covered fraction (over the trajectory).
 184         real                      averageCoveredFraction_;
 185         //! true if the value can change as a function of time.
 186         bool                      bDynamic_;
 187         //! true if the covered fraction depends on the frame.
 188         bool                      bDynamicCoveredFraction_;
 189
 190         /*! \brief
 191          * Needed to wrap access to information.
 192          */
 193         friend class gmx::Selection;
 194         /*! \brief
 195          * Needed for proper access to position information.
 196          */
 197         friend class gmx::SelectionPosition;
 198
 199         GMX_DISALLOW_COPY_AND_ASSIGN(SelectionData);
 200 };
 201
 202 }   // namespace internal
 203
 204 /*! \brief
 205  * Provides access to a single selection.
 206  *
 207  * This class provides a public interface for accessing selection information.
 208  * General information about the selection can be accessed with methods name(),
 209  * selectionText(), isDynamic(), and type().  The first three can be accessed
 210  * any time after the selection has been parsed, and type() can be accessed
 211  * after the selection has been compiled.
 212  *
 213  * Each selection is made of a set of positions.  Each position has associated
 214  * coordinates, and possibly velocities and forces if they have been requested
 215  * and are available.  It also has a set of atoms associated with it; typically
 216  * the coordinates are the center-of-mass or center-of-geometry coordinates for
 217  * that set of atoms.  To access the number of positions in the selection, use
 218  * posCount().  To access individual positions, use position().
 219  * See SelectionPosition for details of how to use individual positions.
 220  * setOriginalId() can be used to adjust the return value of
 221  * SelectionPosition::mappedId(); see that method for details.
 222  *
 223  * It is also possible to access the list of atoms that make up all the
 224  * positions directly: atomCount() returns the total number of atoms in the
 225  * selection and atomIndices() an array of their indices.
 226  * Similarly, it is possible to access the coordinates and other properties
 227  * of the positions as continuous arrays through coordinates(), velocities(),
 228  * forces(), masses(), charges(), refIds(), and mappedIds().
 229  *
 230  * Both positions and atoms can be accessed after the selection has been
 231  * compiled.  For dynamic selections, the return values of these methods change
 232  * after each evaluation to reflect the situation for the current frame.
 233  * Before any frame has been evaluated, these methods return the maximal set
 234  * to which the selection can evaluate.
 235  *
 236  * There are two possible modes for how positions for dynamic selections are
 237  * handled.  In the default mode, posCount() can change, and for each frame,
 238  * only the positions that are selected in that frame can be accessed.  In a
 239  * masked mode, posCount() remains constant, i.e., the positions are always
 240  * evaluated for the maximal set, and SelectionPosition::selected() is used to
 241  * determine whether a position is selected for a frame.  The masked mode can
 242  * be requested with SelectionOption::dynamicMask().
 243  *
 244  * The class also provides methods for printing out information: printInfo()
 245  * and printDebugInfo().  These are mainly for internal use by Gromacs.
 246  *
 247  * This class works like a pointer type: copying and assignment is lightweight,
 248  * and all copies work interchangeably, accessing the same internal data.
 249  *
 250  * Methods in this class do not throw.
 251  *
 252  * \see SelectionPosition
 253  *
 254  * \inpublicapi
 255  * \ingroup module_selection
 256  */
 257 class Selection
 258 {
 259     public:
 260         /*! \brief
 261          * Creates a selection wrapper that has no associated selection.
 262          *
 263          * Any attempt to call methods in the object before a selection is
 264          * assigned results in undefined behavior.
 265          */
 266         Selection() : sel_(NULL) {}
 267         /*! \brief
 268          * Creates a new selection object.
 269          *
 270          * \param  sel  Selection data to wrap.
 271          *
 272          * Only for internal use by the selection module.
 273          */
 274         explicit Selection(internal::SelectionData *sel) : sel_(sel) {}
 275
 276         //! Returns the name of the selection.
 277         const char *name() const  { return data().name_.c_str(); }
 278         //! Returns the string that was parsed to produce this selection.
 279         const char *selectionText() const { return data().selectionText(); }
 280         //! Returns true if the size of the selection (posCount()) is dynamic.
 281         bool isDynamic() const { return data().isDynamic(); }
 282         //! Returns the type of positions in the selection.
 283         e_index_t type() const { return data().rawPositions_.m.type; }
 284
 285         //! Total number of atoms in the selection.
 286         int atomCount() const
 287         {
 288             return data().rawPositions_.g != NULL ? data().rawPositions_.g->isize : 0;
 289         }
 290         //! Returns atom indices of all atoms in the selection.
 291         ConstArrayRef<int> atomIndices() const
 292         {
 293             if (data().rawPositions_.g == NULL)
 294             {
 295                 return ConstArrayRef<int>();
 296             }
 297             return ConstArrayRef<int>(data().rawPositions_.g->index,
 298                                       data().rawPositions_.g->isize);
 299         }
 300         //! Number of positions in the selection.
 301         int posCount() const { return data().posCount(); }
 302         //! Access a single position.
 303         SelectionPosition position(int i) const;
 304         //! Returns coordinates for this selection as a continuous array.
 305         ConstArrayRef<rvec> coordinates() const
 306         {
 307             return ConstArrayRef<rvec>(data().rawPositions_.x, posCount());
 308         }
 309         //! Returns whether velocities are available for this selection.
 310         bool hasVelocities() const { return data().rawPositions_.v != NULL; }
 311         /*! \brief
 312          * Returns velocities for this selection as a continuous array.
 313          *
 314          * Must not be called if hasVelocities() returns false.
 315          */
 316         ConstArrayRef<rvec> velocities() const
 317         {
 318             GMX_ASSERT(hasVelocities(), "Velocities accessed, but unavailable");
 319             return ConstArrayRef<rvec>(data().rawPositions_.v, posCount());
 320         }
 321         //! Returns whether forces are available for this selection.
 322         bool hasForces() const { return sel_->rawPositions_.f != NULL; }
 323         /*! \brief
 324          * Returns forces for this selection as a continuous array.
 325          *
 326          * Must not be called if hasForces() returns false.
 327          */
 328         ConstArrayRef<rvec> forces() const
 329         {
 330             GMX_ASSERT(hasForces(), "Forces accessed, but unavailable");
 331             return ConstArrayRef<rvec>(data().rawPositions_.f, posCount());
 332         }
 333         //! Returns masses for this selection as a continuous array.
 334         ConstArrayRef<real> masses() const
 335         {
 336             // posMass_ may have more entries than posCount() in the case of
 337             // dynamic selections that don't have a topology
 338             // (and thus the masses and charges are fixed).
 339             GMX_ASSERT(data().posMass_.size() >= static_cast<size_t>(posCount()),
 340                        "Internal inconsistency");
 341             return ConstArrayRef<real>(data().posMass_.begin(),
 342                                        data().posMass_.begin() + posCount());
 343         }
 344         //! Returns charges for this selection as a continuous array.
 345         ConstArrayRef<real> charges() const
 346         {
 347             // posCharge_ may have more entries than posCount() in the case of
 348             // dynamic selections that don't have a topology
 349             // (and thus the masses and charges are fixed).
 350             GMX_ASSERT(data().posCharge_.size() >= static_cast<size_t>(posCount()),
 351                        "Internal inconsistency");
 352             return ConstArrayRef<real>(data().posCharge_.begin(),
 353                                        data().posCharge_.begin() + posCount());
 354         }
 355         /*! \brief
 356          * Returns reference IDs for this selection as a continuous array.
 357          *
 358          * \see SelectionPosition::refId()
 359          */
 360         ConstArrayRef<int> refIds() const
 361         {
 362             return ConstArrayRef<int>(data().rawPositions_.m.refid, posCount());
 363         }
 364         /*! \brief
 365          * Returns mapped IDs for this selection as a continuous array.
 366          *
 367          * \see SelectionPosition::mappedId()
 368          */
 369         ConstArrayRef<int> mappedIds() const
 370         {
 371             return ConstArrayRef<int>(data().rawPositions_.m.mapid, posCount());
 372         }
 373         /*! \brief
 374          * Sets the ID for the \p i'th position for use with
 375          * SelectionPosition::mappedId().
 376          *
 377          * \param[in] i  Zero-based index
 378          * \param[in] id Identifier to set.
 379          *
 380          * This method is not part of SelectionPosition because that interface
 381          * only provides access to const data by design.
 382          *
 383          * This method can only be called after compilation, before the
 384          * selection has been evaluated for any frame.
 385          *
 386          * \see SelectionPosition::mappedId()
 387          */
 388         void setOriginalId(int i, int id) { data().rawPositions_.m.orgid[i] = id; }
 389
 390         //! Deprecated method for direct access to position data.
 391         const gmx_ana_pos_t *positions() const { return &data().rawPositions_; }
 392
 393         //! Returns whether the covered fraction can change between frames.
 394         bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
 395         //! Returns the covered fraction for the current frame.
 396         real coveredFraction() const { return data().coveredFraction_; }
 397         /*! \brief
 398          * Initializes information about covered fractions.
 399          *
 400          * \param[in] type Type of covered fraction required.
 401          * \returns   true if the covered fraction can be calculated for the
 402          *      selection.
 403          */
 404         bool initCoveredFraction(e_coverfrac_t type)
 405         {
 406             return data().initCoveredFraction(type);
 407         }
 408
 409         /*! \brief
 410          * Prints out one-line description of the selection.
 411          *
 412          * \param[in] fp      Where to print the information.
 413          *
 414          * The output contains the name of the selection, the number of atoms
 415          * and the number of positions, and indication of whether the selection
 416          * is dynamic.
 417          */
 418         void printInfo(FILE *fp) const;
 419         /*! \brief
 420          * Prints out extended information about the selection for debugging.
 421          *
 422          * \param[in] fp      Where to print the information.
 423          * \param[in] nmaxind Maximum number of values to print in lists
 424          *      (-1 = print all).
 425          */
 426         void printDebugInfo(FILE *fp, int nmaxind) const;
 427
 428     private:
 429         internal::SelectionData &data()
 430         {
 431             GMX_ASSERT(sel_ != NULL,
 432                        "Attempted to access uninitialized selection");
 433             return *sel_;
 434         }
 435         const internal::SelectionData &data() const
 436         {
 437             GMX_ASSERT(sel_ != NULL,
 438                        "Attempted to access uninitialized selection");
 439             return *sel_;
 440         }
 441
 442         /*! \brief
 443          * Pointer to internal data for the selection.
 444          *
 445          * The memory for this object is managed by a SelectionCollection
 446          * object, and the \ref Selection class simply provides a public
 447          * interface for accessing the data.
 448          */
 449         internal::SelectionData *sel_;
 450
 451         /*! \brief
 452          * Needed to access the data to adjust flags.
 453          */
 454         friend class SelectionOptionStorage;
 455 };
 456
 457 /*! \brief
 458  * Provides access to information about a single selected position.
 459  *
 460  * Each position has associated coordinates, and possibly velocities and forces
 461  * if they have been requested and are available.  It also has a set of atoms
 462  * associated with it; typically the coordinates are the center-of-mass or
 463  * center-of-geometry coordinates for that set of atoms.  It is possible that
 464  * there are not atoms associated if the selection has been provided as a fixed
 465  * position.
 466  *
 467  * After the selection has been compiled, but not yet evaluated, the contents
 468  * of the coordinate, velocity and force vectors are undefined.
 469  *
 470  * Default copy constructor and assignment operators are used, and work as
 471  * intended: the copy references the same position and works identically.
 472  *
 473  * Methods in this class do not throw.
 474  *
 475  * \see Selection
 476  *
 477  * \inpublicapi
 478  * \ingroup module_selection
 479  */
 480 class SelectionPosition
 481 {
 482     public:
 483         /*! \brief
 484          * Constructs a wrapper object for given selection position.
 485          *
 486          * \param[in] sel    Selection from which the position is wrapped.
 487          * \param[in] index  Zero-based index of the position to wrap.
 488          *
 489          * Asserts if \p index is out of range.
 490          *
 491          * Only for internal use of the library.  To obtain a SelectionPosition
 492          * object in other code, use Selection::position().
 493          */
 494         SelectionPosition(const internal::SelectionData &sel, int index)
 495             : sel_(&sel), i_(index)
 496         {
 497             GMX_ASSERT(index >= 0 && index < sel.posCount(),
 498                        "Invalid selection position index");
 499         }
 500
 501         /*! \brief
 502          * Returns type of this position.
 503          *
 504          * Currently always returns the same as Selection::type().
 505          */
 506         e_index_t type() const { return sel_->rawPositions_.m.type; }
 507         //! Returns coordinates for this position.
 508         const rvec &x() const
 509         {
 510             return sel_->rawPositions_.x[i_];
 511         }
 512         /*! \brief
 513          * Returns velocity for this position.
 514          *
 515          * Must not be called if Selection::hasVelocities() returns false.
 516          */
 517         const rvec &v() const
 518         {
 519             GMX_ASSERT(sel_->rawPositions_.v != NULL,
 520                        "Velocities accessed, but unavailable");
 521             return sel_->rawPositions_.v[i_];
 522         }
 523         /*! \brief
 524          * Returns force for this position.
 525          *
 526          * Must not be called if Selection::hasForces() returns false.
 527          */
 528         const rvec &f() const
 529         {
 530             GMX_ASSERT(sel_->rawPositions_.f != NULL,
 531                        "Velocities accessed, but unavailable");
 532             return sel_->rawPositions_.f[i_];
 533         }
 534         /*! \brief
 535          * Returns total mass for this position.
 536          *
 537          * Returns the total mass of atoms that make up this position.
 538          * If there are no atoms associated or masses are not available,
 539          * returns unity.
 540          */
 541         real mass() const
 542         {
 543             return sel_->posMass_[i_];
 544         }
 545         /*! \brief
 546          * Returns total charge for this position.
 547          *
 548          * Returns the sum of charges of atoms that make up this position.
 549          * If there are no atoms associated or charges are not available,
 550          * returns zero.
 551          */
 552         real charge() const
 553         {
 554             return sel_->posCharge_[i_];
 555         }
 556         //! Returns the number of atoms that make up this position.
 557         int atomCount() const
 558         {
 559             return sel_->rawPositions_.m.mapb.index[i_ + 1]
 560                    - sel_->rawPositions_.m.mapb.index[i_];
 561         }
 562         //! Return atom indices that make up this position.
 563         ConstArrayRef<int> atomIndices() const
 564         {
 565             if (sel_->rawPositions_.g == NULL)
 566             {
 567                 return ConstArrayRef<int>();
 568             }
 569             int first = sel_->rawPositions_.m.mapb.index[i_];
 570             return ConstArrayRef<int>(&sel_->rawPositions_.g->index[first],
 571                                       atomCount());
 572         }
 573         /*! \brief
 574          * Returns whether this position is selected in the current frame.
 575          *
 576          * The return value is equivalent to \c refid() == -1.  Returns always
 577          * true if SelectionOption::dynamicMask() has not been set.
 578          *
 579          * \see refId()
 580          */
 581         bool selected() const
 582         {
 583             return refId() >= 0;
 584         }
 585         /*! \brief
 586          * Returns reference ID for this position.
 587          *
 588          * For dynamic selections, this provides means to associate positions
 589          * across frames.  After compilation, these IDs are consequently
 590          * numbered starting from zero.  For each frame, the ID then reflects
 591          * the location of the position in the original array of positions.
 592          * If SelectionOption::dynamicMask() has been set for the parent
 593          * selection, the IDs for positions not present in the current
 594          * selection are set to -1, otherwise they are removed completely.
 595          *
 596          * Example:
 597          * If a dynamic selection consists of at most three positions, after
 598          * compilation refId() will return 0, 1, 2 for them, respectively.
 599          * If for a particular frame, only the first and the third are present,
 600          * refId() will return 0, 2.
 601          * If SelectionOption::dynamicMask() has been set, all three positions
 602          * can be accessed also for that frame and refId() will return 0, -1,
 603          * 2.
 604          */
 605         int refId() const
 606         {
 607             return sel_->rawPositions_.m.refid[i_];
 608         }
 609         /*! \brief
 610          * Returns mapped ID for this position.
 611          *
 612          * Returns ID of the position that corresponds to that set with
 613          * Selection::setOriginalId().
 614          *
 615          * If for an array \c id, \c setOriginalId(i, id[i]) has been called
 616          * for each \c i, then it always holds that
 617          * \c mappedId()==id[refId()].
 618          *
 619          * Selection::setOriginalId() has not been called, the default values
 620          * are dependent on type():
 621          *  - ::INDEX_ATOM: atom indices
 622          *  - ::INDEX_RES:  residue numbers
 623          *  - ::INDEX_MOL:  molecule numbers
 624          *  .
 625          * All the default values are zero-based
 626          */
 627         int mappedId() const
 628         {
 629             return sel_->rawPositions_.m.mapid[i_];
 630         }
 631
 632     private:
 633         const internal::SelectionData  *sel_;
 634         int                             i_;
 635 };
 636
 637
 638 inline SelectionPosition
 639 Selection::position(int i) const
 640 {
 641     return SelectionPosition(data(), i);
 642 }
 643
 644 } // namespace gmx
 645
 646 #endif