src/gromacs/selection/selection.h

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