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          * Computed values are cached, and need to be updated for dynamic
 121          * selections with refreshMassesAndCharges() after the selection has
 122          * been evaluated.  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          * Called by SelectionEvaluator.
 134          */
 135         void refreshMassesAndCharges();
 136         /*! \brief
 137          * Updates the covered fraction after a selection has been evaluated.
 138          *
 139          * Called by SelectionEvaluator.
 140          */
 141         void updateCoveredFractionForFrame();
 142         /*! \brief
 143          * Computes average covered fraction after all frames have been evaluated.
 144          *
 145          * \param[in] nframes  Number of frames that have been evaluated.
 146          *
 147          * \p nframes should be equal to the number of calls to
 148          * updateCoveredFractionForFrame().
 149          * Called by SelectionEvaluator::evaluateFinal().
 150          */
 151         void computeAverageCoveredFraction(int nframes);
 152         /*! \brief
 153          * Restores position information to state it was in after compilation.
 154          *
 155          * Depends on SelectionCompiler storing the original atoms in the
 156          * \a rootElement_ object.
 157          * Called by SelectionEvaluator::evaluateFinal().
 158          */
 159         void restoreOriginalPositions();
 160
 161     private:
 162         /*! \brief
 163          * Additional information about positions.
 164          *
 165          * This structure contains information about positions in the
 166          * selection that is not stored in ::gmx_ana_pos_t.
 167          *
 168          * Methods in this class do not throw.
 169          */
 170         struct PositionInfo
 171         {
 172             //! Construct position information with unit mass and no charge.
 173             PositionInfo() : mass(1.0), charge(0.0) {}
 174             //! Construct position information with the given information.
 175             PositionInfo(real mass, real charge) : mass(mass), charge(charge) {}
 176
 177             //! Total mass of atoms that make up the position.
 178             real                mass;
 179             //! Total charge of atoms that make up the position.
 180             real                charge;
 181         };
 182
 183         //! Name of the selection.
 184         std::string             name_;
 185         //! The actual selection string.
 186         std::string             selectionText_;
 187         //! Low-level representation of selected positions.
 188         gmx_ana_pos_t           rawPositions_;
 189         //! Information associated with the current positions.
 190         std::vector<PositionInfo> posInfo_;
 191         //! Information for all possible positions.
 192         std::vector<PositionInfo> originalPosInfo_;
 193         SelectionFlags          flags_;
 194         //! Root of the selection evaluation tree.
 195         SelectionTreeElement   &rootElement_;
 196         //! Type of the covered fraction.
 197         e_coverfrac_t           coveredFractionType_;
 198         //! Covered fraction of the selection for the current frame.
 199         real                    coveredFraction_;
 200         //! The average covered fraction (over the trajectory).
 201         real                    averageCoveredFraction_;
 202         //! true if the value can change as a function of time.
 203         bool                    bDynamic_;
 204         //! true if the covered fraction depends on the frame.
 205         bool                    bDynamicCoveredFraction_;
 206
 207         /*! \brief
 208          * Needed to wrap access to information.
 209          */
 210         friend class gmx::Selection;
 211         /*! \brief
 212          * Needed for proper access to position information.
 213          */
 214         friend class gmx::SelectionPosition;
 215
 216         GMX_DISALLOW_COPY_AND_ASSIGN(SelectionData);
 217 };
 218
 219 } // namespace internal
 220
 221 /*! \brief
 222  * Provides access to a single selection.
 223  *
 224  * This class provides a public interface for accessing selection information.
 225  * General information about the selection can be accessed with methods name(),
 226  * selectionText(), isDynamic(), and type().  The first three can be accessed
 227  * any time after the selection has been parsed, and type() can be accessed
 228  * after the selection has been compiled.
 229  *
 230  * Each selection is made of a set of positions.  Each position has associated
 231  * coordinates, and possibly velocities and forces if they have been requested
 232  * and are available.  It also has a set of atoms associated with it; typically
 233  * the coordinates are the center-of-mass or center-of-geometry coordinates for
 234  * that set of atoms.  To access the number of positions in the selection, use
 235  * posCount().  To access individual positions, use position().
 236  * See SelectionPosition for details of how to use individual positions.
 237  * setOriginalId() can be used to adjust the return value of
 238  * SelectionPosition::mappedId(); see that method for details.
 239  *
 240  * It is also possible to access the list of atoms that make up all the
 241  * positions directly: atomCount() returns the total number of atoms in the
 242  * selection and atomIndices() an array of their indices.
 243  *
 244  * Both positions and atoms can be accessed after the selection has been
 245  * compiled.  For dynamic selections, the return values of these methods change
 246  * after each evaluation to reflect the situation for the current frame.
 247  * Before any frame has been evaluated, these methods return the maximal set
 248  * to which the selection can evaluate.
 249  *
 250  * There are two possible modes for how positions for dynamic selections are
 251  * handled.  In the default mode, posCount() can change, and for each frame,
 252  * only the positions that are selected in that frame can be accessed.  In a
 253  * masked mode, posCount() remains constant, i.e., the positions are always
 254  * evaluated for the maximal set, and SelectionPosition::selected() is used to
 255  * determine whether a position is selected for a frame.  The masked mode can
 256  * be requested with SelectionOption::dynamicMask().
 257  *
 258  * The class also provides methods for printing out information: printInfo()
 259  * and printDebugInfo().  These are mainly for internal use by Gromacs.
 260  *
 261  * This class works like a pointer type: copying and assignment is lightweight,
 262  * and all copies work interchangeably, accessing the same internal data.
 263  *
 264  * Methods in this class do not throw.
 265  *
 266  * \see SelectionPosition
 267  *
 268  * \inpublicapi
 269  * \ingroup module_selection
 270  */
 271 class Selection
 272 {
 273     public:
 274         /*! \brief
 275          * Creates a selection wrapper that has no associated selection.
 276          *
 277          * Any attempt to call methods in the object before a selection is
 278          * assigned results in undefined behavior.
 279          */
 280         Selection() : sel_(NULL) {}
 281         /*! \brief
 282          * Creates a new selection object.
 283          *
 284          * \param  sel  Selection data to wrap.
 285          *
 286          * Only for internal use by the selection module.
 287          */
 288         explicit Selection(internal::SelectionData *sel) : sel_(sel) {}
 289
 290         //! Returns the name of the selection.
 291         const char *name() const  { return data().name_.c_str(); }
 292         //! Returns the string that was parsed to produce this selection.
 293         const char *selectionText() const { return data().selectionText(); }
 294         //! Returns true if the size of the selection (posCount()) is dynamic.
 295         bool isDynamic() const { return data().isDynamic(); }
 296         //! Returns the type of positions in the selection.
 297         e_index_t type() const { return data().rawPositions_.m.type; }
 298
 299         //! Total number of atoms in the selection.
 300         int atomCount() const
 301         {
 302             return data().rawPositions_.g != NULL ? data().rawPositions_.g->isize : 0;
 303         }
 304         //! Returns atom indices of all atoms in the selection.
 305         ConstArrayRef<int> atomIndices() const
 306         {
 307             if (data().rawPositions_.g == NULL)
 308             {
 309                 return ConstArrayRef<int>();
 310             }
 311             return ConstArrayRef<int>(data().rawPositions_.g->isize,
 312                                       data().rawPositions_.g->index);
 313         }
 314         //! Number of positions in the selection.
 315         int posCount() const { return data().posCount(); }
 316         //! Access a single position.
 317         SelectionPosition position(int i) const;
 318         /*! \brief
 319          * Sets the ID for the \p i'th position for use with
 320          * SelectionPosition::mappedId().
 321          *
 322          * \param[in] i  Zero-based index
 323          * \param[in] id Identifier to set.
 324          *
 325          * This method is not part of SelectionPosition because that interface
 326          * only provides access to const data by design.
 327          *
 328          * This method can only be called after compilation, before the
 329          * selection has been evaluated for any frame.
 330          *
 331          * \see SelectionPosition::mappedId()
 332          */
 333         void setOriginalId(int i, int id) { data().rawPositions_.m.orgid[i] = id; }
 334
 335         //! Deprecated method for direct access to position data.
 336         const gmx_ana_pos_t *positions() const { return &data().rawPositions_; }
 337
 338         //! Returns whether the covered fraction can change between frames.
 339         bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
 340         //! Returns the covered fraction for the current frame.
 341         real coveredFraction() const { return data().coveredFraction_; }
 342         /*! \brief
 343          * Initializes information about covered fractions.
 344          *
 345          * \param[in] type Type of covered fraction required.
 346          * \returns   true if the covered fraction can be calculated for the
 347          *      selection.
 348          */
 349         bool initCoveredFraction(e_coverfrac_t type)
 350         {
 351             return data().initCoveredFraction(type);
 352         }
 353
 354         /*! \brief
 355          * Prints out one-line description of the selection.
 356          *
 357          * \param[in] fp      Where to print the information.
 358          *
 359          * The output contains the name of the selection, the number of atoms
 360          * and the number of positions, and indication of whether the selection
 361          * is dynamic.
 362          */
 363         void printInfo(FILE *fp) const;
 364         /*! \brief
 365          * Prints out extended information about the selection for debugging.
 366          *
 367          * \param[in] fp      Where to print the information.
 368          * \param[in] nmaxind Maximum number of values to print in lists
 369          *      (-1 = print all).
 370          */
 371         void printDebugInfo(FILE *fp, int nmaxind) const;
 372
 373     private:
 374         internal::SelectionData &data()
 375         {
 376             GMX_ASSERT(sel_ != NULL,
 377                        "Attempted to access uninitialized selection");
 378             return *sel_;
 379         }
 380         const internal::SelectionData &data() const
 381         {
 382             GMX_ASSERT(sel_ != NULL,
 383                        "Attempted to access uninitialized selection");
 384             return *sel_;
 385         }
 386
 387         /*! \brief
 388          * Pointer to internal data for the selection.
 389          *
 390          * The memory for this object is managed by a SelectionCollection
 391          * object, and the \ref Selection class simply provides a public
 392          * interface for accessing the data.
 393          */
 394         internal::SelectionData *sel_;
 395
 396         /*! \brief
 397          * Needed to access the data to adjust flags.
 398          */
 399         friend class SelectionOptionStorage;
 400 };
 401
 402 /*! \brief
 403  * Provides access to information about a single selected position.
 404  *
 405  * Each position has associated coordinates, and possibly velocities and forces
 406  * if they have been requested and are available.  It also has a set of atoms
 407  * associated with it; typically the coordinates are the center-of-mass or
 408  * center-of-geometry coordinates for that set of atoms.  It is possible that
 409  * there are not atoms associated if the selection has been provided as a fixed
 410  * position.
 411  *
 412  * After the selection has been compiled, but not yet evaluated, the contents
 413  * of the coordinate, velocity and force vectors are undefined.
 414  *
 415  * Default copy constructor and assignment operators are used, and work as
 416  * intended: the copy references the same position and works identically.
 417  *
 418  * Methods in this class do not throw.
 419  *
 420  * \see Selection
 421  *
 422  * \inpublicapi
 423  * \ingroup module_selection
 424  */
 425 class SelectionPosition
 426 {
 427     public:
 428         /*! \brief
 429          * Constructs a wrapper object for given selection position.
 430          *
 431          * \param[in] sel    Selection from which the position is wrapped.
 432          * \param[in] index  Zero-based index of the position to wrap.
 433          *
 434          * Asserts if \p index is out of range.
 435          *
 436          * Only for internal use of the library.  To obtain a SelectionPosition
 437          * object in other code, use Selection::position().
 438          */
 439         SelectionPosition(const internal::SelectionData &sel, int index)
 440             : sel_(&sel), i_(index)
 441         {
 442             GMX_ASSERT(index >= 0 && index < sel.posCount(),
 443                        "Invalid selection position index");
 444         }
 445
 446         /*! \brief
 447          * Returns type of this position.
 448          *
 449          * Currently always returns the same as Selection::type().
 450          */
 451         e_index_t type() const { return sel_->rawPositions_.m.type; }
 452         //! Returns coordinates for this position.
 453         const rvec &x() const
 454         {
 455             return sel_->rawPositions_.x[i_];
 456         }
 457         //! Returns whether velocity is available for this position.
 458         bool hasVelocity() const { return sel_->rawPositions_.v != NULL; }
 459         /*! \brief
 460          * Returns velocity for this position.
 461          *
 462          * Must not be called if hasVelocity() returns false.
 463          */
 464         const rvec &v() const
 465         {
 466             GMX_ASSERT(hasVelocity(), "Velocities accessed, but unavailable");
 467             return sel_->rawPositions_.v[i_];
 468         }
 469         //! Returns whether force is available for this position.
 470         bool hasForce() const { return sel_->rawPositions_.f != NULL; }
 471         /*! \brief
 472          * Returns force for this position.
 473          *
 474          * Must not be called if hasForce() returns false.
 475          */
 476         const rvec &f() const
 477         {
 478             GMX_ASSERT(hasForce(), "Forces accessed, but unavailable");
 479             return sel_->rawPositions_.f[i_];
 480         }
 481         /*! \brief
 482          * Returns total mass for this position.
 483          *
 484          * Returns the total mass of atoms that make up this position.
 485          * If there are no atoms associated or masses are not available,
 486          * returns unity.
 487          */
 488         real mass() const
 489         {
 490             return sel_->posInfo_[i_].mass;
 491         }
 492         /*! \brief
 493          * Returns total charge for this position.
 494          *
 495          * Returns the sum of charges of atoms that make up this position.
 496          * If there are no atoms associated or charges are not available,
 497          * returns zero.
 498          */
 499         real charge() const
 500         {
 501             return sel_->posInfo_[i_].charge;
 502         }
 503         //! Returns the number of atoms that make up this position.
 504         int atomCount() const
 505         {
 506             return sel_->rawPositions_.m.mapb.index[i_ + 1]
 507                  - sel_->rawPositions_.m.mapb.index[i_];
 508         }
 509         //! Return atom indices that make up this position.
 510         ConstArrayRef<int> atomIndices() const
 511         {
 512             if (sel_->rawPositions_.g == NULL)
 513             {
 514                 return ConstArrayRef<int>();
 515             }
 516             int first = sel_->rawPositions_.m.mapb.index[i_];
 517             return ConstArrayRef<int>(atomCount(),
 518                                       &sel_->rawPositions_.g->index[first]);
 519         }
 520         /*! \brief
 521          * Returns whether this position is selected in the current frame.
 522          *
 523          * The return value is equivalent to \c refid() == -1.  Returns always
 524          * true if SelectionOption::dynamicMask() has not been set.
 525          *
 526          * \see refId()
 527          */
 528         bool selected() const
 529         {
 530             return refId() >= 0;
 531         }
 532         /*! \brief
 533          * Returns reference ID for this position.
 534          *
 535          * For dynamic selections, this provides means to associate positions
 536          * across frames.  After compilation, these IDs are consequently
 537          * numbered starting from zero.  For each frame, the ID then reflects
 538          * the location of the position in the original array of positions.
 539          * If SelectionOption::dynamicMask() has been set for the parent
 540          * selection, the IDs for positions not present in the current
 541          * selection are set to -1, otherwise they are removed completely.
 542          *
 543          * Example:
 544          * If a dynamic selection consists of at most three positions, after
 545          * compilation refId() will return 0, 1, 2 for them, respectively.
 546          * If for a particular frame, only the first and the third are present,
 547          * refId() will return 0, 2.
 548          * If SelectionOption::dynamicMask() has been set, all three positions
 549          * can be accessed also for that frame and refId() will return 0, -1,
 550          * 2.
 551          */
 552         int refId() const
 553         {
 554             return sel_->rawPositions_.m.refid[i_];
 555         }
 556         /*! \brief
 557          * Returns mapped ID for this position.
 558          *
 559          * Returns ID of the position that corresponds to that set with
 560          * Selection::setOriginalId().
 561          *
 562          * If for an array \c id, \c setOriginalId(i, id[i]) has been called
 563          * for each \c i, then it always holds that
 564          * \c mappedId()==id[refId()].
 565          *
 566          * Selection::setOriginalId() has not been called, the default values
 567          * are dependent on type():
 568          *  - ::INDEX_ATOM: atom indices
 569          *  - ::INDEX_RES:  residue numbers
 570          *  - ::INDEX_MOL:  molecule numbers
 571          *  .
 572          * All the default values are zero-based
 573          */
 574         int mappedId() const
 575         {
 576             return sel_->rawPositions_.m.mapid[i_];
 577         }
 578
 579     private:
 580         const internal::SelectionData  *sel_;
 581         int                             i_;
 582 };
 583
 584
 585 inline SelectionPosition
 586 Selection::position(int i) const
 587 {
 588     return SelectionPosition(data(), i);
 589 }
 590
 591 } // namespace gmx
 592
 593 #endif