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