2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, 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 * Implements classes in selection.h.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_selection
44 #include "selection.h"
48 #include "gromacs/selection/nbsearch.h"
49 #include "gromacs/selection/position.h"
50 #include "gromacs/topology/topology.h"
51 #include "gromacs/utility/exceptions.h"
52 #include "gromacs/utility/gmxassert.h"
53 #include "gromacs/utility/stringutil.h"
54 #include "gromacs/utility/textwriter.h"
65 /********************************************************************
69 SelectionData::SelectionData(SelectionTreeElement
*elem
,
71 : name_(elem
->name()), selectionText_(selstr
),
72 rootElement_(*elem
), coveredFractionType_(CFRAC_NONE
),
73 coveredFraction_(1.0), averageCoveredFraction_(1.0),
74 bDynamic_(false), bDynamicCoveredFraction_(false)
76 if (elem
->child
->type
== SEL_CONST
)
78 // TODO: This is not exception-safe if any called function throws.
79 gmx_ana_pos_copy(&rawPositions_
, elem
->child
->v
.u
.p
, true);
83 SelectionTreeElementPointer child
= elem
->child
;
84 child
->flags
&= ~SEL_ALLOCVAL
;
85 _gmx_selvalue_setstore(&child
->v
, &rawPositions_
);
86 /* We should also skip any modifiers to determine the dynamic
88 while (child
->type
== SEL_MODIFIER
)
95 if (child
->type
== SEL_SUBEXPRREF
)
98 /* Because most subexpression elements are created
99 * during compilation, we need to check for them
102 if (child
->type
== SEL_SUBEXPR
)
104 child
= child
->child
;
110 /* For variable references, we should skip the
111 * SEL_SUBEXPRREF and SEL_SUBEXPR elements. */
112 if (child
->type
== SEL_SUBEXPRREF
)
114 child
= child
->child
->child
;
116 bDynamic_
= (child
->child
->flags
& SEL_DYNAMIC
);
119 initCoveredFraction(CFRAC_NONE
);
123 SelectionData::~SelectionData()
129 SelectionData::initCoveredFraction(e_coverfrac_t type
)
131 coveredFractionType_
= type
;
132 if (type
== CFRAC_NONE
)
134 bDynamicCoveredFraction_
= false;
136 else if (!_gmx_selelem_can_estimate_cover(rootElement()))
138 coveredFractionType_
= CFRAC_NONE
;
139 bDynamicCoveredFraction_
= false;
143 bDynamicCoveredFraction_
= true;
145 coveredFraction_
= bDynamicCoveredFraction_
? 0.0 : 1.0;
146 averageCoveredFraction_
= coveredFraction_
;
147 return type
== CFRAC_NONE
|| coveredFractionType_
!= CFRAC_NONE
;
154 * Helper function to compute total masses and charges for positions.
156 * \param[in] top Topology to take atom masses from.
157 * \param[in] pos Positions to compute masses and charges for.
158 * \param[out] masses Output masses.
159 * \param[out] charges Output charges.
161 * Does not throw if enough space has been reserved for the output vectors.
163 void computeMassesAndCharges(const t_topology
*top
, const gmx_ana_pos_t
&pos
,
164 std::vector
<real
> *masses
,
165 std::vector
<real
> *charges
)
167 GMX_ASSERT(top
!= NULL
, "Should not have been called with NULL topology");
170 for (int b
= 0; b
< pos
.count(); ++b
)
174 for (int i
= pos
.m
.mapb
.index
[b
]; i
< pos
.m
.mapb
.index
[b
+1]; ++i
)
176 const int index
= pos
.m
.mapb
.a
[i
];
177 mass
+= top
->atoms
.atom
[index
].m
;
178 charge
+= top
->atoms
.atom
[index
].q
;
180 masses
->push_back(mass
);
181 charges
->push_back(charge
);
188 SelectionData::hasSortedAtomIndices() const
191 gmx_ana_index_set(&g
, rawPositions_
.m
.mapb
.nra
, rawPositions_
.m
.mapb
.a
, -1);
192 return gmx_ana_index_check_sorted(&g
);
196 SelectionData::refreshName()
198 rootElement_
.fillNameIfMissing(selectionText_
.c_str());
199 name_
= rootElement_
.name();
203 SelectionData::initializeMassesAndCharges(const t_topology
*top
)
205 GMX_ASSERT(posMass_
.empty() && posCharge_
.empty(),
206 "Should not be called more than once");
207 posMass_
.reserve(posCount());
208 posCharge_
.reserve(posCount());
211 posMass_
.resize(posCount(), 1.0);
212 posCharge_
.resize(posCount(), 0.0);
216 computeMassesAndCharges(top
, rawPositions_
, &posMass_
, &posCharge_
);
222 SelectionData::refreshMassesAndCharges(const t_topology
*top
)
224 if (top
!= NULL
&& isDynamic() && !hasFlag(efSelection_DynamicMask
))
226 computeMassesAndCharges(top
, rawPositions_
, &posMass_
, &posCharge_
);
232 SelectionData::updateCoveredFractionForFrame()
234 if (isCoveredFractionDynamic())
236 real cfrac
= _gmx_selelem_estimate_coverfrac(rootElement());
237 coveredFraction_
= cfrac
;
238 averageCoveredFraction_
+= cfrac
;
244 SelectionData::computeAverageCoveredFraction(int nframes
)
246 if (isCoveredFractionDynamic() && nframes
> 0)
248 averageCoveredFraction_
/= nframes
;
254 SelectionData::restoreOriginalPositions(const t_topology
*top
)
258 gmx_ana_pos_t
&p
= rawPositions_
;
259 gmx_ana_indexmap_update(&p
.m
, rootElement().v
.u
.g
,
260 hasFlag(gmx::efSelection_DynamicMask
));
261 refreshMassesAndCharges(top
);
265 } // namespace internal
267 /********************************************************************
271 Selection::operator AnalysisNeighborhoodPositions() const
273 AnalysisNeighborhoodPositions
pos(data().rawPositions_
.x
,
274 data().rawPositions_
.count());
277 pos
.exclusionIds(atomIndices());
284 Selection::setOriginalId(int i
, int id
)
286 data().rawPositions_
.m
.mapid
[i
] = id
;
287 data().rawPositions_
.m
.orgid
[i
] = id
;
292 Selection::initOriginalIdsToGroup(const gmx_mtop_t
*top
, e_index_t type
)
296 return gmx_ana_indexmap_init_orgid_group(&data().rawPositions_
.m
, top
, type
);
298 catch (const InconsistentInputError
&)
300 GMX_ASSERT(type
== INDEX_RES
|| type
== INDEX_MOL
,
301 "Expected that only grouping by residue/molecule would fail");
302 std::string message
=
303 formatString("Cannot group selection '%s' into %s, because some "
304 "positions have atoms from more than one such group.",
305 name(), type
== INDEX_MOL
? "molecules" : "residues");
306 GMX_THROW(InconsistentInputError(message
));
312 Selection::printInfo(FILE *fp
) const
314 fprintf(fp
, "\"%s\" (%d position%s, %d atom%s%s)", name(),
315 posCount(), posCount() == 1 ? "" : "s",
316 atomCount(), atomCount() == 1 ? "" : "s",
317 isDynamic() ? ", dynamic" : "");
323 Selection::printDebugInfo(FILE *fp
, int nmaxind
) const
325 const gmx_ana_pos_t
&p
= data().rawPositions_
;
329 fprintf(fp
, " Group ");
331 gmx_ana_index_set(&g
, p
.m
.mapb
.nra
, p
.m
.mapb
.a
, 0);
332 TextWriter
writer(fp
);
333 gmx_ana_index_dump(&writer
, &g
, nmaxind
);
335 fprintf(fp
, " Block (size=%d):", p
.m
.mapb
.nr
);
338 fprintf(fp
, " (null)");
343 if (nmaxind
>= 0 && n
> nmaxind
)
347 for (int i
= 0; i
<= n
; ++i
)
349 fprintf(fp
, " %d", p
.m
.mapb
.index
[i
]);
359 if (nmaxind
>= 0 && n
> nmaxind
)
363 fprintf(fp
, " RefId:");
366 fprintf(fp
, " (null)");
370 for (int i
= 0; i
< n
; ++i
)
372 fprintf(fp
, " %d", p
.m
.refid
[i
]);
381 fprintf(fp
, " MapId:");
384 fprintf(fp
, " (null)");
388 for (int i
= 0; i
< n
; ++i
)
390 fprintf(fp
, " %d", p
.m
.mapid
[i
]);
401 /********************************************************************
405 SelectionPosition::operator AnalysisNeighborhoodPositions() const
407 AnalysisNeighborhoodPositions
pos(sel_
->rawPositions_
.x
,
408 sel_
->rawPositions_
.count());
409 if (sel_
->hasOnlyAtoms())
411 // TODO: Move atomIndices() such that it can be reused here as well.
412 pos
.exclusionIds(constArrayRefFromArray
<int>(sel_
->rawPositions_
.m
.mapb
.a
,
413 sel_
->rawPositions_
.m
.mapb
.nra
));
415 return pos
.selectSingleFromArray(i_
);