2 * This file is part of the GROMACS molecular simulation package.
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,2021, 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.
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.
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.
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.
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.
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.
39 * Declares gmx::Selection and supporting classes.
41 * \author Teemu Murtola <teemu.murtola@gmail.com>
43 * \ingroup module_selection
45 #ifndef GMX_SELECTION_SELECTION_H
46 #define GMX_SELECTION_SELECTION_H
52 #include "gromacs/utility/arrayref.h"
53 #include "gromacs/utility/classhelpers.h"
54 #include "gromacs/utility/gmxassert.h"
57 #include "selectionenums.h"
64 class SelectionOptionStorage;
65 class SelectionTreeElement;
67 class AnalysisNeighborhoodPositions;
69 class SelectionPosition;
71 //! Container of selections used in public selection interfaces.
72 typedef std::vector<Selection> SelectionList;
79 * Internal data for a single selection.
81 * This class is internal to the selection module, but resides in a public
82 * header because of efficiency reasons: it allows frequently used access
83 * methods in \ref Selection to be inlined.
85 * Methods in this class do not throw unless otherwise specified.
87 * \ingroup module_selection
93 * Creates a new selection object.
95 * \param[in] elem Root of the evaluation tree for this selection.
96 * \param[in] selstr String that was parsed to produce this selection.
97 * \throws std::bad_alloc if out of memory.
99 SelectionData(SelectionTreeElement* elem, const char* selstr);
102 //! Returns the name for this selection.
103 const char* name() const { return name_.c_str(); }
104 //! Returns the string that was parsed to produce this selection.
105 const char* selectionText() const { return selectionText_.c_str(); }
106 //! Returns true if the size of the selection (posCount()) is dynamic.
107 bool isDynamic() const { return bDynamic_; }
108 //! Returns the type of positions in the selection.
109 e_index_t type() const { return rawPositions_.m.type; }
110 //! Returns true if the selection only contains positions with a single atom each.
111 bool hasOnlyAtoms() const { return type() == INDEX_ATOM; }
112 //! Returns `true` if the atom indices in the selection are in ascending order.
113 bool hasSortedAtomIndices() const;
115 //! Number of positions in the selection.
116 int posCount() const { return rawPositions_.count(); }
117 //! Returns the root of the evaluation tree for this selection.
118 SelectionTreeElement& rootElement() { return rootElement_; }
120 //! Returns whether the covered fraction can change between frames.
121 bool isCoveredFractionDynamic() const { return bDynamicCoveredFraction_; }
123 //! Returns true if the given flag is set.
124 bool hasFlag(SelectionFlag flag) const { return flags_.test(flag); }
125 //! Sets the flags for this selection.
126 void setFlags(SelectionFlags flags) { flags_ = flags; }
127 //! Returns the current flags.
128 SelectionFlags flags() const { return flags_; }
130 //! \copydoc Selection::initCoveredFraction()
131 bool initCoveredFraction(e_coverfrac_t type);
134 * Updates the name of the selection if missing.
136 * \throws std::bad_alloc if out of memory.
138 * If selections get their value from a group reference that cannot be
139 * resolved during parsing, the name is final only after group
140 * references have been resolved.
142 * This function is called by SelectionCollection::setIndexGroups().
146 * Computes total masses and charges for all selection positions.
148 * \param[in] top Topology information.
149 * \throws std::bad_alloc if out of memory.
151 * For dynamic selections, the values need to be updated after each
152 * evaluation with refreshMassesAndCharges().
153 * This is done by SelectionEvaluator.
155 * This function is called by SelectionCompiler.
157 * Strong exception safety guarantee.
159 void initializeMassesAndCharges(const gmx_mtop_t* top);
161 * Updates masses and charges after dynamic selection has been
164 * \param[in] top Topology information.
166 * Called by SelectionEvaluator.
168 void refreshMassesAndCharges(const gmx_mtop_t* top);
170 * Updates the covered fraction after a selection has been evaluated.
172 * Called by SelectionEvaluator.
174 void updateCoveredFractionForFrame();
176 * Computes average covered fraction after all frames have been evaluated.
178 * \param[in] nframes Number of frames that have been evaluated.
180 * \p nframes should be equal to the number of calls to
181 * updateCoveredFractionForFrame().
182 * Called by SelectionEvaluator::evaluateFinal().
184 void computeAverageCoveredFraction(int nframes);
186 * Restores position information to state it was in after compilation.
188 * \param[in] top Topology information.
190 * Depends on SelectionCompiler storing the original atoms in the
191 * \a rootElement_ object.
192 * Called by SelectionEvaluator::evaluateFinal().
194 void restoreOriginalPositions(const gmx_mtop_t* top);
197 //! Name of the selection.
199 //! The actual selection string.
200 std::string selectionText_;
201 //! Low-level representation of selected positions.
202 gmx_ana_pos_t rawPositions_;
203 //! Total masses for the current positions.
204 std::vector<real> posMass_;
205 //! Total charges for the current positions.
206 std::vector<real> posCharge_;
207 SelectionFlags flags_;
208 //! Root of the selection evaluation tree.
209 SelectionTreeElement& rootElement_;
210 //! Type of the covered fraction.
211 e_coverfrac_t coveredFractionType_;
212 //! Covered fraction of the selection for the current frame.
213 real coveredFraction_;
214 //! The average covered fraction (over the trajectory).
215 real averageCoveredFraction_;
216 //! true if the value can change as a function of time.
218 //! true if the covered fraction depends on the frame.
219 bool bDynamicCoveredFraction_;
222 * Needed to wrap access to information.
224 friend class gmx::Selection;
226 * Needed for proper access to position information.
228 friend class gmx::SelectionPosition;
230 GMX_DISALLOW_COPY_AND_ASSIGN(SelectionData);
233 } // namespace internal
236 * Provides access to a single selection.
238 * This class provides a public interface for accessing selection information.
239 * General information about the selection can be accessed with methods name(),
240 * selectionText(), isDynamic(), and type(). The first three can be accessed
241 * any time after the selection has been parsed, and type() can be accessed
242 * after the selection has been compiled.
244 * There are a few methods that can be used to change the behavior of the
245 * selection. setEvaluateVelocities() and setEvaluateForces() can be called
246 * before the selection is compiled to request evaluation of velocities and/or
247 * forces in addition to coordinates.
249 * Each selection is made of a set of positions. Each position has associated
250 * coordinates, and possibly velocities and forces if they have been requested
251 * and are available. It also has a set of atoms associated with it; typically
252 * the coordinates are the center-of-mass or center-of-geometry coordinates for
253 * that set of atoms. To access the number of positions in the selection, use
254 * posCount(). To access individual positions, use position().
255 * See SelectionPosition for details of how to use individual positions.
256 * setOriginalId() can be used to adjust the return value of
257 * SelectionPosition::mappedId(); see that method for details.
259 * It is also possible to access the list of atoms that make up all the
260 * positions directly: atomCount() returns the total number of atoms in the
261 * selection and atomIndices() an array of their indices.
262 * Similarly, it is possible to access the coordinates and other properties
263 * of the positions as continuous arrays through coordinates(), velocities(),
264 * forces(), masses(), charges(), refIds(), and mappedIds().
266 * Both positions and atoms can be accessed after the selection has been
267 * compiled. For dynamic selections, the return values of these methods change
268 * after each evaluation to reflect the situation for the current frame.
269 * Before any frame has been evaluated, these methods return the maximal set
270 * to which the selection can evaluate.
272 * There are two possible modes for how positions for dynamic selections are
273 * handled. In the default mode, posCount() can change, and for each frame,
274 * only the positions that are selected in that frame can be accessed. In a
275 * masked mode, posCount() remains constant, i.e., the positions are always
276 * evaluated for the maximal set, and SelectionPosition::selected() is used to
277 * determine whether a position is selected for a frame. The masked mode can
278 * be requested with SelectionOption::dynamicMask().
280 * The class also provides methods for printing out information: printInfo()
281 * and printDebugInfo(). These are mainly for internal use by Gromacs.
283 * This class works like a pointer type: copying and assignment is lightweight,
284 * and all copies work interchangeably, accessing the same internal data.
286 * Methods in this class do not throw.
288 * \see SelectionPosition
291 * \ingroup module_selection
297 * Creates a selection wrapper that has no associated selection.
299 * Any attempt to call methods in the object before a selection is
300 * assigned results in undefined behavior.
301 * isValid() returns `false` for the selection until it is initialized.
303 Selection() : sel_(nullptr) {}
305 * Creates a new selection object.
307 * \param sel Selection data to wrap.
309 * Only for internal use by the selection module.
311 explicit Selection(internal::SelectionData* sel) : sel_(sel) {}
313 //! Returns whether the selection object is initialized.
314 bool isValid() const { return sel_ != nullptr; }
316 //! Returns whether two selection objects wrap the same selection.
317 bool operator==(const Selection& other) const { return sel_ == other.sel_; }
318 //! Returns whether two selection objects wrap different selections.
319 bool operator!=(const Selection& other) const { return !operator==(other); }
321 //! Returns the name of the selection.
322 const char* name() const { return data().name(); }
323 //! Returns the string that was parsed to produce this selection.
324 const char* selectionText() const { return data().selectionText(); }
325 //! Returns true if the size of the selection (posCount()) is dynamic.
326 bool isDynamic() const { return data().isDynamic(); }
327 //! Returns the type of positions in the selection.
328 e_index_t type() const { return data().type(); }
329 //! Returns true if the selection only contains positions with a single atom each.
330 bool hasOnlyAtoms() const { return data().hasOnlyAtoms(); }
331 //! Returns `true` if the atom indices in the selection are in ascending order.
332 bool hasSortedAtomIndices() const { return data().hasSortedAtomIndices(); }
334 //! Total number of atoms in the selection.
335 int atomCount() const { return data().rawPositions_.m.mapb.nra; }
336 //! Returns atom indices of all atoms in the selection.
337 ArrayRef<const int> atomIndices() const
339 return constArrayRefFromArray(sel_->rawPositions_.m.mapb.a, sel_->rawPositions_.m.mapb.nra);
341 //! Number of positions in the selection.
342 int posCount() const { return data().posCount(); }
343 //! Access a single position.
344 SelectionPosition position(int i) const;
345 //! Returns coordinates for this selection as a continuous array.
346 ArrayRef<const rvec> coordinates() const
348 return constArrayRefFromArray(data().rawPositions_.x, posCount());
350 //! Returns whether velocities are available for this selection.
351 bool hasVelocities() const { return data().rawPositions_.v != nullptr; }
353 * Returns velocities for this selection as a continuous array.
355 * Must not be called if hasVelocities() returns false.
357 ArrayRef<const rvec> velocities() const
359 GMX_ASSERT(hasVelocities(), "Velocities accessed, but unavailable");
360 return constArrayRefFromArray(data().rawPositions_.v, posCount());
362 //! Returns whether forces are available for this selection.
363 bool hasForces() const { return sel_->rawPositions_.f != nullptr; }
365 * Returns forces for this selection as a continuous array.
367 * Must not be called if hasForces() returns false.
369 ArrayRef<const rvec> forces() const
371 GMX_ASSERT(hasForces(), "Forces accessed, but unavailable");
372 return constArrayRefFromArray(data().rawPositions_.f, posCount());
374 //! Returns masses for this selection as a continuous array.
375 ArrayRef<const real> masses() const
377 // posMass_ may have more entries than posCount() in the case of
378 // dynamic selections that don't have a topology
379 // (and thus the masses and charges are fixed).
380 GMX_ASSERT(data().posMass_.size() >= static_cast<size_t>(posCount()),
381 "Internal inconsistency");
382 return makeArrayRef(data().posMass_).subArray(0, posCount());
384 //! Returns charges for this selection as a continuous array.
385 ArrayRef<const real> charges() const
387 // posCharge_ may have more entries than posCount() in the case of
388 // dynamic selections that don't have a topology
389 // (and thus the masses and charges are fixed).
390 GMX_ASSERT(data().posCharge_.size() >= static_cast<size_t>(posCount()),
391 "Internal inconsistency");
392 return makeArrayRef(data().posCharge_).subArray(0, posCount());
395 * Returns reference IDs for this selection as a continuous array.
397 * \see SelectionPosition::refId()
399 ArrayRef<const int> refIds() const
401 return constArrayRefFromArray(data().rawPositions_.m.refid, posCount());
404 * Returns mapped IDs for this selection as a continuous array.
406 * \see SelectionPosition::mappedId()
408 ArrayRef<const int> mappedIds() const
410 return constArrayRefFromArray(data().rawPositions_.m.mapid, posCount());
413 //! Returns whether the covered fraction can change between frames.
414 bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
415 //! Returns the covered fraction for the current frame.
416 real coveredFraction() const { return data().coveredFraction_; }
419 * Allows passing a selection directly to neighborhood searching.
421 * When initialized this way, AnalysisNeighborhoodPair objects return
422 * indices that can be used to index the selection positions with
425 * Works exactly like if AnalysisNeighborhoodPositions had a
426 * constructor taking a Selection object as a parameter.
427 * See AnalysisNeighborhoodPositions for rationale and additional
430 operator AnalysisNeighborhoodPositions() const;
433 * Initializes information about covered fractions.
435 * \param[in] type Type of covered fraction required.
436 * \returns true if the covered fraction can be calculated for the
439 bool initCoveredFraction(e_coverfrac_t type) { return data().initCoveredFraction(type); }
441 * Sets whether this selection evaluates velocities for positions.
443 * \param[in] bEnabled If true, velocities are evaluated.
445 * If you request the evaluation, but then evaluate the selection for
446 * a frame that does not contain velocity information, results are
450 * Implement it such that in the above case, hasVelocities() will
451 * return false for such frames.
455 void setEvaluateVelocities(bool bEnabled)
457 data().flags_.set(efSelection_EvaluateVelocities, bEnabled);
460 * Sets whether this selection evaluates forces for positions.
462 * \param[in] bEnabled If true, forces are evaluated.
464 * If you request the evaluation, but then evaluate the selection for
465 * a frame that does not contain force information, results are
470 void setEvaluateForces(bool bEnabled)
472 data().flags_.set(efSelection_EvaluateForces, bEnabled);
476 * Sets the ID for the \p i'th position for use with
477 * SelectionPosition::mappedId().
479 * \param[in] i Zero-based index
480 * \param[in] id Identifier to set.
482 * This method is not part of SelectionPosition because that interface
483 * only provides access to const data by design.
485 * This method can only be called after compilation, before the
486 * selection has been evaluated for any frame.
488 * \see SelectionPosition::mappedId()
490 void setOriginalId(int i, int id);
492 * Inits the IDs for use with SelectionPosition::mappedId() for
495 * \param[in] top Topology information
496 * (can be NULL if not required for \p type).
497 * \param[in] type Type of groups to generate.
498 * \returns Number of groups that were present in the selection.
499 * \throws InconsistentInputError if the selection positions cannot
500 * be assigned to groups of the given type.
502 * If `type == INDEX_ATOM`, the IDs are initialized to 0, 1, 2, ...,
503 * and the return value is the number of positions.
504 * If `type == INDEX_ALL`, all the IDs are initialized to 0, and the
505 * return value is one.
506 * If `type == INDEX_RES` or `type == INDEX_MOL`, the first position
507 * will get ID 0, and all following positions that belong to the same
508 * residue/molecule will get the same ID. The first position that
509 * belongs to a different residue/molecule will get ID 1, and so on.
510 * If some position contains atoms from multiple residues/molecules,
511 * i.e., the mapping is ambiguous, an exception is thrown.
512 * The return value is the number of residues/molecules that are
513 * present in the selection positions.
515 * This method is useful if the calling code needs to group the
516 * selection, e.g., for computing aggregate properties for each residue
517 * or molecule. It can then use this method to initialize the
518 * appropriate grouping, use the return value to allocate a
519 * sufficiently sized buffer to store the aggregated values, and then
520 * use SelectionPosition::mappedId() to identify the location where to
523 * \see setOriginalId()
524 * \see SelectionPosition::mappedId()
526 int initOriginalIdsToGroup(const gmx_mtop_t* top, e_index_t type);
529 * Prints out one-line description of the selection.
531 * \param[in] fp Where to print the information.
533 * The output contains the name of the selection, the number of atoms
534 * and the number of positions, and indication of whether the selection
537 void printInfo(FILE* fp) const;
539 * Prints out extended information about the selection for debugging.
541 * \param[in] fp Where to print the information.
542 * \param[in] nmaxind Maximum number of values to print in lists
545 void printDebugInfo(FILE* fp, int nmaxind) const;
548 internal::SelectionData& data()
550 GMX_ASSERT(sel_ != nullptr, "Attempted to access uninitialized selection");
553 const internal::SelectionData& data() const
555 GMX_ASSERT(sel_ != nullptr, "Attempted to access uninitialized selection");
560 * Pointer to internal data for the selection.
562 * The memory for this object is managed by a SelectionCollection
563 * object, and the \ref Selection class simply provides a public
564 * interface for accessing the data.
566 internal::SelectionData* sel_;
569 * Needed to access the data to adjust flags.
571 friend class SelectionOptionStorage;
575 * Provides access to information about a single selected position.
577 * Each position has associated coordinates, and possibly velocities and forces
578 * if they have been requested and are available. It also has a set of atoms
579 * associated with it; typically the coordinates are the center-of-mass or
580 * center-of-geometry coordinates for that set of atoms. It is possible that
581 * there are not atoms associated if the selection has been provided as a fixed
584 * After the selection has been compiled, but not yet evaluated, the contents
585 * of the coordinate, velocity and force vectors are undefined.
587 * Default copy constructor and assignment operators are used, and work as
588 * intended: the copy references the same position and works identically.
590 * Methods in this class do not throw.
595 * \ingroup module_selection
597 class SelectionPosition
601 * Constructs a wrapper object for given selection position.
603 * \param[in] sel Selection from which the position is wrapped.
604 * \param[in] index Zero-based index of the position to wrap.
606 * Asserts if \p index is out of range.
608 * Only for internal use of the library. To obtain a SelectionPosition
609 * object in other code, use Selection::position().
611 SelectionPosition(const internal::SelectionData& sel, int index) : sel_(&sel), i_(index)
613 GMX_ASSERT(index >= 0 && index < sel.posCount(), "Invalid selection position index");
617 * Returns type of this position.
619 * Currently always returns the same as Selection::type().
621 e_index_t type() const { return sel_->type(); }
622 //! Returns coordinates for this position.
623 const rvec& x() const { return sel_->rawPositions_.x[i_]; }
625 * Returns velocity for this position.
627 * Must not be called if Selection::hasVelocities() returns false.
629 const rvec& v() const
631 GMX_ASSERT(sel_->rawPositions_.v != nullptr, "Velocities accessed, but unavailable");
632 return sel_->rawPositions_.v[i_];
635 * Returns force for this position.
637 * Must not be called if Selection::hasForces() returns false.
639 const rvec& f() const
641 GMX_ASSERT(sel_->rawPositions_.f != nullptr, "Forces accessed, but unavailable");
642 return sel_->rawPositions_.f[i_];
645 * Returns total mass for this position.
647 * Returns the total mass of atoms that make up this position.
648 * If there are no atoms associated or masses are not available,
651 real mass() const { return sel_->posMass_[i_]; }
653 * Returns total charge for this position.
655 * Returns the sum of charges of atoms that make up this position.
656 * If there are no atoms associated or charges are not available,
659 real charge() const { return sel_->posCharge_[i_]; }
660 //! Returns the number of atoms that make up this position.
661 int atomCount() const
663 return sel_->rawPositions_.m.mapb.index[i_ + 1] - sel_->rawPositions_.m.mapb.index[i_];
665 //! Return atom indices that make up this position.
666 ArrayRef<const int> atomIndices() const
668 const int* atoms = sel_->rawPositions_.m.mapb.a;
669 if (atoms == nullptr)
671 return ArrayRef<const int>();
673 const int first = sel_->rawPositions_.m.mapb.index[i_];
674 return constArrayRefFromArray(&atoms[first], atomCount());
677 * Returns whether this position is selected in the current frame.
679 * The return value is equivalent to \c refid() == -1. Returns always
680 * true if SelectionOption::dynamicMask() has not been set.
684 bool selected() const { return refId() >= 0; }
686 * Returns reference ID for this position.
688 * For dynamic selections, this provides means to associate positions
689 * across frames. After compilation, these IDs are consequently
690 * numbered starting from zero. For each frame, the ID then reflects
691 * the location of the position in the original array of positions.
692 * If SelectionOption::dynamicMask() has been set for the parent
693 * selection, the IDs for positions not present in the current
694 * selection are set to -1, otherwise they are removed completely.
697 * If a dynamic selection consists of at most three positions, after
698 * compilation refId() will return 0, 1, 2 for them, respectively.
699 * If for a particular frame, only the first and the third are present,
700 * refId() will return 0, 2.
701 * If SelectionOption::dynamicMask() has been set, all three positions
702 * can be accessed also for that frame and refId() will return 0, -1,
705 int refId() const { return sel_->rawPositions_.m.refid[i_]; }
707 * Returns mapped ID for this position.
709 * Returns ID of the position that corresponds to that set with
710 * Selection::setOriginalId().
712 * If for an array \c id, \c setOriginalId(i, id[i]) has been called
713 * for each \c i, then it always holds that
714 * \c mappedId()==id[refId()].
716 * Selection::setOriginalId() has not been called, the default values
717 * are dependent on type():
718 * - ::INDEX_ATOM: atom indices
719 * - ::INDEX_RES: residue indices
720 * - ::INDEX_MOL: molecule indices
722 * All the default values are zero-based.
724 int mappedId() const { return sel_->rawPositions_.m.mapid[i_]; }
727 * Allows passing a selection position directly to neighborhood searching.
729 * When initialized this way, AnalysisNeighborhoodPair objects return
730 * the index that can be used to access this position using
731 * Selection::position().
733 * Works exactly like if AnalysisNeighborhoodPositions had a
734 * constructor taking a SelectionPosition object as a parameter.
735 * See AnalysisNeighborhoodPositions for rationale and additional
738 operator AnalysisNeighborhoodPositions() const;
741 const internal::SelectionData* sel_;
746 inline SelectionPosition Selection::position(int i) const
748 return SelectionPosition(data(), i);