2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, 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.
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.
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.
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.
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.
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.
37 * Declares gmx::Selection and supporting classes.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
41 * \ingroup module_selection
43 #ifndef GMX_SELECTION_SELECTION_H
44 #define GMX_SELECTION_SELECTION_H
49 #include "gromacs/selection/position.h"
50 #include "gromacs/selection/selectionenums.h"
51 #include "gromacs/utility/arrayref.h"
52 #include "gromacs/utility/classhelpers.h"
53 #include "gromacs/utility/gmxassert.h"
60 class SelectionOptionStorage
;
61 class SelectionTreeElement
;
63 class AnalysisNeighborhoodPositions
;
65 class SelectionPosition
;
67 //! Container of selections used in public selection interfaces.
68 typedef std::vector
<Selection
> SelectionList
;
75 * Internal data for a single selection.
77 * This class is internal to the selection module, but resides in a public
78 * header because of efficiency reasons: it allows frequently used access
79 * methods in \ref Selection to be inlined.
81 * Methods in this class do not throw unless otherwise specified.
83 * \ingroup module_selection
89 * Creates a new selection object.
91 * \param[in] elem Root of the evaluation tree for this selection.
92 * \param[in] selstr String that was parsed to produce this selection.
93 * \throws std::bad_alloc if out of memory.
95 SelectionData(SelectionTreeElement
*elem
, const char *selstr
);
98 //! Returns the name for this selection.
99 const char *name() const { return name_
.c_str(); }
100 //! Returns the string that was parsed to produce this selection.
101 const char *selectionText() const { return selectionText_
.c_str(); }
102 //! Returns true if the size of the selection (posCount()) is dynamic.
103 bool isDynamic() const { return bDynamic_
; }
104 //! Returns the type of positions in the selection.
105 e_index_t
type() const { return rawPositions_
.m
.type
; }
106 //! Returns true if the selection only contains positions with a single atom each.
107 bool hasOnlyAtoms() const { return type() == INDEX_ATOM
; }
108 //! Returns `true` if the atom indices in the selection are in ascending order.
109 bool hasSortedAtomIndices() const;
111 //! Number of positions in the selection.
112 int posCount() const { return rawPositions_
.count(); }
113 //! Returns the root of the evaluation tree for this selection.
114 SelectionTreeElement
&rootElement() { return rootElement_
; }
116 //! Returns whether the covered fraction can change between frames.
117 bool isCoveredFractionDynamic() const { return bDynamicCoveredFraction_
; }
119 //! Returns true if the given flag is set.
120 bool hasFlag(SelectionFlag flag
) const { return flags_
.test(flag
); }
121 //! Sets the flags for this selection.
122 void setFlags(SelectionFlags flags
) { flags_
= flags
; }
124 //! \copydoc Selection::initCoveredFraction()
125 bool initCoveredFraction(e_coverfrac_t type
);
128 * Updates the name of the selection if missing.
130 * \throws std::bad_alloc if out of memory.
132 * If selections get their value from a group reference that cannot be
133 * resolved during parsing, the name is final only after group
134 * references have been resolved.
136 * This function is called by SelectionCollection::setIndexGroups().
140 * Computes total masses and charges for all selection positions.
142 * \param[in] top Topology information.
143 * \throws std::bad_alloc if out of memory.
145 * For dynamic selections, the values need to be updated after each
146 * evaluation with refreshMassesAndCharges().
147 * This is done by SelectionEvaluator.
149 * This function is called by SelectionCompiler.
151 * Strong exception safety guarantee.
153 void initializeMassesAndCharges(const gmx_mtop_t
*top
);
155 * Updates masses and charges after dynamic selection has been
158 * \param[in] top Topology information.
160 * Called by SelectionEvaluator.
162 void refreshMassesAndCharges(const gmx_mtop_t
*top
);
164 * Updates the covered fraction after a selection has been evaluated.
166 * Called by SelectionEvaluator.
168 void updateCoveredFractionForFrame();
170 * Computes average covered fraction after all frames have been evaluated.
172 * \param[in] nframes Number of frames that have been evaluated.
174 * \p nframes should be equal to the number of calls to
175 * updateCoveredFractionForFrame().
176 * Called by SelectionEvaluator::evaluateFinal().
178 void computeAverageCoveredFraction(int nframes
);
180 * Restores position information to state it was in after compilation.
182 * \param[in] top Topology information.
184 * Depends on SelectionCompiler storing the original atoms in the
185 * \a rootElement_ object.
186 * Called by SelectionEvaluator::evaluateFinal().
188 void restoreOriginalPositions(const gmx_mtop_t
*top
);
191 //! Name of the selection.
193 //! The actual selection string.
194 std::string selectionText_
;
195 //! Low-level representation of selected positions.
196 gmx_ana_pos_t rawPositions_
;
197 //! Total masses for the current positions.
198 std::vector
<real
> posMass_
;
199 //! Total charges for the current positions.
200 std::vector
<real
> posCharge_
;
201 SelectionFlags flags_
;
202 //! Root of the selection evaluation tree.
203 SelectionTreeElement
&rootElement_
;
204 //! Type of the covered fraction.
205 e_coverfrac_t coveredFractionType_
;
206 //! Covered fraction of the selection for the current frame.
207 real coveredFraction_
;
208 //! The average covered fraction (over the trajectory).
209 real averageCoveredFraction_
;
210 //! true if the value can change as a function of time.
212 //! true if the covered fraction depends on the frame.
213 bool bDynamicCoveredFraction_
;
216 * Needed to wrap access to information.
218 friend class gmx::Selection
;
220 * Needed for proper access to position information.
222 friend class gmx::SelectionPosition
;
224 GMX_DISALLOW_COPY_AND_ASSIGN(SelectionData
);
227 } // namespace internal
230 * Provides access to a single selection.
232 * This class provides a public interface for accessing selection information.
233 * General information about the selection can be accessed with methods name(),
234 * selectionText(), isDynamic(), and type(). The first three can be accessed
235 * any time after the selection has been parsed, and type() can be accessed
236 * after the selection has been compiled.
238 * There are a few methods that can be used to change the behavior of the
239 * selection. setEvaluateVelocities() and setEvaluateForces() can be called
240 * before the selection is compiled to request evaluation of velocities and/or
241 * forces in addition to coordinates.
243 * Each selection is made of a set of positions. Each position has associated
244 * coordinates, and possibly velocities and forces if they have been requested
245 * and are available. It also has a set of atoms associated with it; typically
246 * the coordinates are the center-of-mass or center-of-geometry coordinates for
247 * that set of atoms. To access the number of positions in the selection, use
248 * posCount(). To access individual positions, use position().
249 * See SelectionPosition for details of how to use individual positions.
250 * setOriginalId() can be used to adjust the return value of
251 * SelectionPosition::mappedId(); see that method for details.
253 * It is also possible to access the list of atoms that make up all the
254 * positions directly: atomCount() returns the total number of atoms in the
255 * selection and atomIndices() an array of their indices.
256 * Similarly, it is possible to access the coordinates and other properties
257 * of the positions as continuous arrays through coordinates(), velocities(),
258 * forces(), masses(), charges(), refIds(), and mappedIds().
260 * Both positions and atoms can be accessed after the selection has been
261 * compiled. For dynamic selections, the return values of these methods change
262 * after each evaluation to reflect the situation for the current frame.
263 * Before any frame has been evaluated, these methods return the maximal set
264 * to which the selection can evaluate.
266 * There are two possible modes for how positions for dynamic selections are
267 * handled. In the default mode, posCount() can change, and for each frame,
268 * only the positions that are selected in that frame can be accessed. In a
269 * masked mode, posCount() remains constant, i.e., the positions are always
270 * evaluated for the maximal set, and SelectionPosition::selected() is used to
271 * determine whether a position is selected for a frame. The masked mode can
272 * be requested with SelectionOption::dynamicMask().
274 * The class also provides methods for printing out information: printInfo()
275 * and printDebugInfo(). These are mainly for internal use by Gromacs.
277 * This class works like a pointer type: copying and assignment is lightweight,
278 * and all copies work interchangeably, accessing the same internal data.
280 * Methods in this class do not throw.
282 * \see SelectionPosition
285 * \ingroup module_selection
291 * Creates a selection wrapper that has no associated selection.
293 * Any attempt to call methods in the object before a selection is
294 * assigned results in undefined behavior.
295 * isValid() returns `false` for the selection until it is initialized.
297 Selection() : sel_(nullptr) {}
299 * Creates a new selection object.
301 * \param sel Selection data to wrap.
303 * Only for internal use by the selection module.
305 explicit Selection(internal::SelectionData
*sel
) : sel_(sel
) {}
307 //! Returns whether the selection object is initialized.
308 bool isValid() const { return sel_
!= nullptr; }
310 //! Returns whether two selection objects wrap the same selection.
311 bool operator==(const Selection
&other
) const
313 return sel_
== other
.sel_
;
315 //! Returns whether two selection objects wrap different selections.
316 bool operator!=(const Selection
&other
) const
318 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
337 return data().rawPositions_
.m
.mapb
.nra
;
339 //! Returns atom indices of all atoms in the selection.
340 ArrayRef
<const int> atomIndices() const
342 return constArrayRefFromArray(sel_
->rawPositions_
.m
.mapb
.a
,
343 sel_
->rawPositions_
.m
.mapb
.nra
);
345 //! Number of positions in the selection.
346 int posCount() const { return data().posCount(); }
347 //! Access a single position.
348 SelectionPosition
position(int i
) const;
349 //! Returns coordinates for this selection as a continuous array.
350 ArrayRef
<const rvec
> coordinates() const
352 return constArrayRefFromArray(data().rawPositions_
.x
, posCount());
354 //! Returns whether velocities are available for this selection.
355 bool hasVelocities() const { return data().rawPositions_
.v
!= nullptr; }
357 * Returns velocities for this selection as a continuous array.
359 * Must not be called if hasVelocities() returns false.
361 ArrayRef
<const rvec
> velocities() const
363 GMX_ASSERT(hasVelocities(), "Velocities accessed, but unavailable");
364 return constArrayRefFromArray(data().rawPositions_
.v
, posCount());
366 //! Returns whether forces are available for this selection.
367 bool hasForces() const { return sel_
->rawPositions_
.f
!= nullptr; }
369 * Returns forces for this selection as a continuous array.
371 * Must not be called if hasForces() returns false.
373 ArrayRef
<const rvec
> forces() const
375 GMX_ASSERT(hasForces(), "Forces accessed, but unavailable");
376 return constArrayRefFromArray(data().rawPositions_
.f
, posCount());
378 //! Returns masses for this selection as a continuous array.
379 ArrayRef
<const real
> masses() const
381 // posMass_ may have more entries than posCount() in the case of
382 // dynamic selections that don't have a topology
383 // (and thus the masses and charges are fixed).
384 GMX_ASSERT(data().posMass_
.size() >= static_cast<size_t>(posCount()),
385 "Internal inconsistency");
386 return makeArrayRef(data().posMass_
).subArray(0, posCount());
388 //! Returns charges for this selection as a continuous array.
389 ArrayRef
<const real
> charges() const
391 // posCharge_ may have more entries than posCount() in the case of
392 // dynamic selections that don't have a topology
393 // (and thus the masses and charges are fixed).
394 GMX_ASSERT(data().posCharge_
.size() >= static_cast<size_t>(posCount()),
395 "Internal inconsistency");
396 return makeArrayRef(data().posCharge_
).subArray(0, posCount());
399 * Returns reference IDs for this selection as a continuous array.
401 * \see SelectionPosition::refId()
403 ArrayRef
<const int> refIds() const
405 return constArrayRefFromArray(data().rawPositions_
.m
.refid
, posCount());
408 * Returns mapped IDs for this selection as a continuous array.
410 * \see SelectionPosition::mappedId()
412 ArrayRef
<const int> mappedIds() const
414 return constArrayRefFromArray(data().rawPositions_
.m
.mapid
, posCount());
417 //! Returns whether the covered fraction can change between frames.
418 bool isCoveredFractionDynamic() const { return data().isCoveredFractionDynamic(); }
419 //! Returns the covered fraction for the current frame.
420 real
coveredFraction() const { return data().coveredFraction_
; }
423 * Allows passing a selection directly to neighborhood searching.
425 * When initialized this way, AnalysisNeighborhoodPair objects return
426 * indices that can be used to index the selection positions with
429 * Works exactly like if AnalysisNeighborhoodPositions had a
430 * constructor taking a Selection object as a parameter.
431 * See AnalysisNeighborhoodPositions for rationale and additional
434 operator AnalysisNeighborhoodPositions() const;
437 * Initializes information about covered fractions.
439 * \param[in] type Type of covered fraction required.
440 * \returns true if the covered fraction can be calculated for the
443 bool initCoveredFraction(e_coverfrac_t type
)
445 return data().initCoveredFraction(type
);
448 * Sets whether this selection evaluates velocities for positions.
450 * \param[in] bEnabled If true, velocities are evaluated.
452 * If you request the evaluation, but then evaluate the selection for
453 * a frame that does not contain velocity information, results are
457 * Implement it such that in the above case, hasVelocities() will
458 * return false for such frames.
462 void setEvaluateVelocities(bool bEnabled
)
464 data().flags_
.set(efSelection_EvaluateVelocities
, bEnabled
);
467 * Sets whether this selection evaluates forces for positions.
469 * \param[in] bEnabled If true, forces are evaluated.
471 * If you request the evaluation, but then evaluate the selection for
472 * a frame that does not contain force information, results are
477 void setEvaluateForces(bool bEnabled
)
479 data().flags_
.set(efSelection_EvaluateForces
, bEnabled
);
483 * Sets the ID for the \p i'th position for use with
484 * SelectionPosition::mappedId().
486 * \param[in] i Zero-based index
487 * \param[in] id Identifier to set.
489 * This method is not part of SelectionPosition because that interface
490 * only provides access to const data by design.
492 * This method can only be called after compilation, before the
493 * selection has been evaluated for any frame.
495 * \see SelectionPosition::mappedId()
497 void setOriginalId(int i
, int id
);
499 * Inits the IDs for use with SelectionPosition::mappedId() for
502 * \param[in] top Topology information
503 * (can be NULL if not required for \p type).
504 * \param[in] type Type of groups to generate.
505 * \returns Number of groups that were present in the selection.
506 * \throws InconsistentInputError if the selection positions cannot
507 * be assigned to groups of the given type.
509 * If `type == INDEX_ATOM`, the IDs are initialized to 0, 1, 2, ...,
510 * and the return value is the number of positions.
511 * If `type == INDEX_ALL`, all the IDs are initialized to 0, and the
512 * return value is one.
513 * If `type == INDEX_RES` or `type == INDEX_MOL`, the first position
514 * will get ID 0, and all following positions that belong to the same
515 * residue/molecule will get the same ID. The first position that
516 * belongs to a different residue/molecule will get ID 1, and so on.
517 * If some position contains atoms from multiple residues/molecules,
518 * i.e., the mapping is ambiguous, an exception is thrown.
519 * The return value is the number of residues/molecules that are
520 * present in the selection positions.
522 * This method is useful if the calling code needs to group the
523 * selection, e.g., for computing aggregate properties for each residue
524 * or molecule. It can then use this method to initialize the
525 * appropriate grouping, use the return value to allocate a
526 * sufficiently sized buffer to store the aggregated values, and then
527 * use SelectionPosition::mappedId() to identify the location where to
530 * \see setOriginalId()
531 * \see SelectionPosition::mappedId()
533 int initOriginalIdsToGroup(const gmx_mtop_t
*top
, e_index_t type
);
536 * Prints out one-line description of the selection.
538 * \param[in] fp Where to print the information.
540 * The output contains the name of the selection, the number of atoms
541 * and the number of positions, and indication of whether the selection
544 void printInfo(FILE *fp
) const;
546 * Prints out extended information about the selection for debugging.
548 * \param[in] fp Where to print the information.
549 * \param[in] nmaxind Maximum number of values to print in lists
552 void printDebugInfo(FILE *fp
, int nmaxind
) const;
555 internal::SelectionData
&data()
557 GMX_ASSERT(sel_
!= nullptr,
558 "Attempted to access uninitialized selection");
561 const internal::SelectionData
&data() const
563 GMX_ASSERT(sel_
!= nullptr,
564 "Attempted to access uninitialized selection");
569 * Pointer to internal data for the selection.
571 * The memory for this object is managed by a SelectionCollection
572 * object, and the \ref Selection class simply provides a public
573 * interface for accessing the data.
575 internal::SelectionData
*sel_
;
578 * Needed to access the data to adjust flags.
580 friend class SelectionOptionStorage
;
584 * Provides access to information about a single selected position.
586 * Each position has associated coordinates, and possibly velocities and forces
587 * if they have been requested and are available. It also has a set of atoms
588 * associated with it; typically the coordinates are the center-of-mass or
589 * center-of-geometry coordinates for that set of atoms. It is possible that
590 * there are not atoms associated if the selection has been provided as a fixed
593 * After the selection has been compiled, but not yet evaluated, the contents
594 * of the coordinate, velocity and force vectors are undefined.
596 * Default copy constructor and assignment operators are used, and work as
597 * intended: the copy references the same position and works identically.
599 * Methods in this class do not throw.
604 * \ingroup module_selection
606 class SelectionPosition
610 * Constructs a wrapper object for given selection position.
612 * \param[in] sel Selection from which the position is wrapped.
613 * \param[in] index Zero-based index of the position to wrap.
615 * Asserts if \p index is out of range.
617 * Only for internal use of the library. To obtain a SelectionPosition
618 * object in other code, use Selection::position().
620 SelectionPosition(const internal::SelectionData
&sel
, int index
)
621 : sel_(&sel
), i_(index
)
623 GMX_ASSERT(index
>= 0 && index
< sel
.posCount(),
624 "Invalid selection position index");
628 * Returns type of this position.
630 * Currently always returns the same as Selection::type().
632 e_index_t
type() const { return sel_
->type(); }
633 //! Returns coordinates for this position.
634 const rvec
&x() const
636 return sel_
->rawPositions_
.x
[i_
];
639 * Returns velocity for this position.
641 * Must not be called if Selection::hasVelocities() returns false.
643 const rvec
&v() const
645 GMX_ASSERT(sel_
->rawPositions_
.v
!= nullptr,
646 "Velocities accessed, but unavailable");
647 return sel_
->rawPositions_
.v
[i_
];
650 * Returns force for this position.
652 * Must not be called if Selection::hasForces() returns false.
654 const rvec
&f() const
656 GMX_ASSERT(sel_
->rawPositions_
.f
!= nullptr,
657 "Velocities accessed, but unavailable");
658 return sel_
->rawPositions_
.f
[i_
];
661 * Returns total mass for this position.
663 * Returns the total mass of atoms that make up this position.
664 * If there are no atoms associated or masses are not available,
669 return sel_
->posMass_
[i_
];
672 * Returns total charge for this position.
674 * Returns the sum of charges of atoms that make up this position.
675 * If there are no atoms associated or charges are not available,
680 return sel_
->posCharge_
[i_
];
682 //! Returns the number of atoms that make up this position.
683 int atomCount() const
685 return sel_
->rawPositions_
.m
.mapb
.index
[i_
+ 1]
686 - sel_
->rawPositions_
.m
.mapb
.index
[i_
];
688 //! Return atom indices that make up this position.
689 ArrayRef
<const int> atomIndices() const
691 const int *atoms
= sel_
->rawPositions_
.m
.mapb
.a
;
692 if (atoms
== nullptr)
694 return ArrayRef
<const int>();
696 const int first
= sel_
->rawPositions_
.m
.mapb
.index
[i_
];
697 return constArrayRefFromArray(&atoms
[first
], atomCount());
700 * Returns whether this position is selected in the current frame.
702 * The return value is equivalent to \c refid() == -1. Returns always
703 * true if SelectionOption::dynamicMask() has not been set.
707 bool selected() const
712 * Returns reference ID for this position.
714 * For dynamic selections, this provides means to associate positions
715 * across frames. After compilation, these IDs are consequently
716 * numbered starting from zero. For each frame, the ID then reflects
717 * the location of the position in the original array of positions.
718 * If SelectionOption::dynamicMask() has been set for the parent
719 * selection, the IDs for positions not present in the current
720 * selection are set to -1, otherwise they are removed completely.
723 * If a dynamic selection consists of at most three positions, after
724 * compilation refId() will return 0, 1, 2 for them, respectively.
725 * If for a particular frame, only the first and the third are present,
726 * refId() will return 0, 2.
727 * If SelectionOption::dynamicMask() has been set, all three positions
728 * can be accessed also for that frame and refId() will return 0, -1,
733 return sel_
->rawPositions_
.m
.refid
[i_
];
736 * Returns mapped ID for this position.
738 * Returns ID of the position that corresponds to that set with
739 * Selection::setOriginalId().
741 * If for an array \c id, \c setOriginalId(i, id[i]) has been called
742 * for each \c i, then it always holds that
743 * \c mappedId()==id[refId()].
745 * Selection::setOriginalId() has not been called, the default values
746 * are dependent on type():
747 * - ::INDEX_ATOM: atom indices
748 * - ::INDEX_RES: residue indices
749 * - ::INDEX_MOL: molecule indices
751 * All the default values are zero-based.
755 return sel_
->rawPositions_
.m
.mapid
[i_
];
759 * Allows passing a selection position directly to neighborhood searching.
761 * When initialized this way, AnalysisNeighborhoodPair objects return
762 * the index that can be used to access this position using
763 * Selection::position().
765 * Works exactly like if AnalysisNeighborhoodPositions had a
766 * constructor taking a SelectionPosition object as a parameter.
767 * See AnalysisNeighborhoodPositions for rationale and additional
770 operator AnalysisNeighborhoodPositions() const;
773 const internal::SelectionData
*sel_
;
778 inline SelectionPosition
779 Selection::position(int i
) const
781 return SelectionPosition(data(), i
);