1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2009,2010,2011,2012,2013 by the GROMACS development team.
5 * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
6 * Copyright (c) 2019,2020, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37 /*! \internal \file
38 * \brief
39 * Implements classes in selection.h.
40 *
41 * \author Teemu Murtola <teemu.murtola@gmail.com>
42 * \ingroup module_selection
43 */
44 #include "gmxpre.h"
45
46 #include "selection.h"
47
48 #include <string>
49
50 #include "gromacs/selection/nbsearch.h"
51 #include "gromacs/topology/mtop_lookup.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/exceptions.h"
54 #include "gromacs/utility/gmxassert.h"
55 #include "gromacs/utility/stringutil.h"
56 #include "gromacs/utility/textwriter.h"
57
58 #include "position.h"
59 #include "selelem.h"
60 #include "selvalue.h"
61
62 namespace gmx
63 {
64
65 namespace internal
66 {
67
68 /********************************************************************
69 * SelectionData
70 */
71
SelectionData(SelectionTreeElement * elem,const char * selstr)72 SelectionData::SelectionData(SelectionTreeElement* elem, const char* selstr) :
73 name_(elem->name()),
74 selectionText_(selstr),
75 rootElement_(*elem),
76 coveredFractionType_(CFRAC_NONE),
77 coveredFraction_(1.0),
78 averageCoveredFraction_(1.0),
79 bDynamic_(false),
80 bDynamicCoveredFraction_(false)
81 {
82 if (elem->child->type == SEL_CONST)
83 {
84 // TODO: This is not exception-safe if any called function throws.
85 gmx_ana_pos_copy(&rawPositions_, elem->child->v.u.p, true);
86 }
87 else
88 {
89 SelectionTreeElementPointer child = elem->child;
90 child->flags &= ~SEL_ALLOCVAL;
91 _gmx_selvalue_setstore(&child->v, &rawPositions_);
92 /* We should also skip any modifiers to determine the dynamic
93 * status. */
94 while (child->type == SEL_MODIFIER)
95 {
96 child = child->child;
97 if (!child)
98 {
99 break;
100 }
101 if (child->type == SEL_SUBEXPRREF)
102 {
103 child = child->child;
104 /* Because most subexpression elements are created
105 * during compilation, we need to check for them
106 * explicitly here.
107 */
108 if (child->type == SEL_SUBEXPR)
109 {
110 child = child->child;
111 }
112 }
113 }
114 if (child)
115 {
116 /* For variable references, we should skip the
117 * SEL_SUBEXPRREF and SEL_SUBEXPR elements. */
118 if (child->type == SEL_SUBEXPRREF)
119 {
120 child = child->child->child;
121 }
122 bDynamic_ = ((child->child->flags & SEL_DYNAMIC) != 0);
123 }
124 }
125 initCoveredFraction(CFRAC_NONE);
126 }
127
128
~SelectionData()129 SelectionData::~SelectionData() {}
130
131
initCoveredFraction(e_coverfrac_t type)132 bool SelectionData::initCoveredFraction(e_coverfrac_t type)
133 {
134 coveredFractionType_ = type;
135 if (type == CFRAC_NONE)
136 {
137 bDynamicCoveredFraction_ = false;
138 }
139 else if (!_gmx_selelem_can_estimate_cover(rootElement()))
140 {
141 coveredFractionType_ = CFRAC_NONE;
142 bDynamicCoveredFraction_ = false;
143 }
144 else
145 {
146 bDynamicCoveredFraction_ = true;
147 }
148 coveredFraction_ = bDynamicCoveredFraction_ ? 0.0 : 1.0;
149 averageCoveredFraction_ = coveredFraction_;
150 return type == CFRAC_NONE || coveredFractionType_ != CFRAC_NONE;
151 }
152
153 namespace
154 {
155
156 /*! \brief
157 * Helper function to compute total masses and charges for positions.
158 *
159 * \param[in] top Topology to take atom masses from.
160 * \param[in] pos Positions to compute masses and charges for.
161 * \param[out] masses Output masses.
162 * \param[out] charges Output charges.
163 *
164 * Does not throw if enough space has been reserved for the output vectors.
165 */
computeMassesAndCharges(const gmx_mtop_t * top,const gmx_ana_pos_t & pos,std::vector<real> * masses,std::vector<real> * charges)166 void computeMassesAndCharges(const gmx_mtop_t* top,
167 const gmx_ana_pos_t& pos,
168 std::vector<real>* masses,
169 std::vector<real>* charges)
170 {
171 GMX_ASSERT(top != nullptr, "Should not have been called with NULL topology");
172 masses->clear();
173 charges->clear();
174 int molb = 0;
175 for (int b = 0; b < pos.count(); ++b)
176 {
177 real mass = 0.0;
178 real charge = 0.0;
179 for (int i = pos.m.mapb.index[b]; i < pos.m.mapb.index[b + 1]; ++i)
180 {
181 const int index = pos.m.mapb.a[i];
182 const t_atom& atom = mtopGetAtomParameters(top, index, &molb);
183 mass += atom.m;
184 charge += atom.q;
185 }
186 masses->push_back(mass);
187 charges->push_back(charge);
188 }
189 }
190
191 } // namespace
192
hasSortedAtomIndices() const193 bool SelectionData::hasSortedAtomIndices() const
194 {
195 gmx_ana_index_t g;
196 gmx_ana_index_set(&g, rawPositions_.m.mapb.nra, rawPositions_.m.mapb.a, -1);
197 return gmx_ana_index_check_sorted(&g);
198 }
199
refreshName()200 void SelectionData::refreshName()
201 {
202 rootElement_.fillNameIfMissing(selectionText_.c_str());
203 name_ = rootElement_.name();
204 }
205
initializeMassesAndCharges(const gmx_mtop_t * top)206 void SelectionData::initializeMassesAndCharges(const gmx_mtop_t* top)
207 {
208 GMX_ASSERT(posMass_.empty() && posCharge_.empty(), "Should not be called more than once");
209 posMass_.reserve(posCount());
210 posCharge_.reserve(posCount());
211 if (top == nullptr)
212 {
213 posMass_.resize(posCount(), 1.0);
214 posCharge_.resize(posCount(), 0.0);
215 }
216 else
217 {
218 computeMassesAndCharges(top, rawPositions_, &posMass_, &posCharge_);
219 }
220 }
221
222
refreshMassesAndCharges(const gmx_mtop_t * top)223 void SelectionData::refreshMassesAndCharges(const gmx_mtop_t* top)
224 {
225 if (top != nullptr && isDynamic() && !hasFlag(efSelection_DynamicMask))
226 {
227 computeMassesAndCharges(top, rawPositions_, &posMass_, &posCharge_);
228 }
229 }
230
231
updateCoveredFractionForFrame()232 void SelectionData::updateCoveredFractionForFrame()
233 {
234 if (isCoveredFractionDynamic())
235 {
236 real cfrac = _gmx_selelem_estimate_coverfrac(rootElement());
237 coveredFraction_ = cfrac;
238 averageCoveredFraction_ += cfrac;
239 }
240 }
241
242
computeAverageCoveredFraction(int nframes)243 void SelectionData::computeAverageCoveredFraction(int nframes)
244 {
245 if (isCoveredFractionDynamic() && nframes > 0)
246 {
247 averageCoveredFraction_ /= nframes;
248 }
249 }
250
251
restoreOriginalPositions(const gmx_mtop_t * top)252 void SelectionData::restoreOriginalPositions(const gmx_mtop_t* top)
253 {
254 if (isDynamic())
255 {
256 gmx_ana_pos_t& p = rawPositions_;
257 gmx_ana_indexmap_update(&p.m, rootElement().v.u.g, hasFlag(gmx::efSelection_DynamicMask));
258 refreshMassesAndCharges(top);
259 }
260 }
261
262 } // namespace internal
263
264 /********************************************************************
265 * Selection
266 */
267
operator AnalysisNeighborhoodPositions() const268 Selection::operator AnalysisNeighborhoodPositions() const
269 {
270 AnalysisNeighborhoodPositions pos(data().rawPositions_.x, data().rawPositions_.count());
271 if (hasOnlyAtoms())
272 {
273 pos.exclusionIds(atomIndices());
274 }
275 return pos;
276 }
277
278
setOriginalId(int i,int id)279 void Selection::setOriginalId(int i, int id)
280 {
281 data().rawPositions_.m.mapid[i] = id;
282 data().rawPositions_.m.orgid[i] = id;
283 }
284
285
initOriginalIdsToGroup(const gmx_mtop_t * top,e_index_t type)286 int Selection::initOriginalIdsToGroup(const gmx_mtop_t* top, e_index_t type)
287 {
288 try
289 {
290 return gmx_ana_indexmap_init_orgid_group(&data().rawPositions_.m, top, type);
291 }
292 catch (const InconsistentInputError&)
293 {
294 GMX_ASSERT(type == INDEX_RES || type == INDEX_MOL,
295 "Expected that only grouping by residue/molecule would fail");
296 std::string message = formatString(
297 "Cannot group selection '%s' into %s, because some "
298 "positions have atoms from more than one such group.",
299 name(), type == INDEX_MOL ? "molecules" : "residues");
300 GMX_THROW(InconsistentInputError(message));
301 }
302 }
303
304
printInfo(FILE * fp) const305 void Selection::printInfo(FILE* fp) const
306 {
307 fprintf(fp, "\"%s\" (%d position%s, %d atom%s%s)", name(), posCount(), posCount() == 1 ? "" : "s",
308 atomCount(), atomCount() == 1 ? "" : "s", isDynamic() ? ", dynamic" : "");
309 fprintf(fp, "\n");
310 }
311
312
printDebugInfo(FILE * fp,int nmaxind) const313 void Selection::printDebugInfo(FILE* fp, int nmaxind) const
314 {
315 const gmx_ana_pos_t& p = data().rawPositions_;
316
317 fprintf(fp, " ");
318 printInfo(fp);
319 fprintf(fp, " Group ");
320 gmx_ana_index_t g;
321 gmx_ana_index_set(&g, p.m.mapb.nra, p.m.mapb.a, 0);
322 TextWriter writer(fp);
323 gmx_ana_index_dump(&writer, &g, nmaxind);
324
325 fprintf(fp, " Block (size=%d):", p.m.mapb.nr);
326 if (!p.m.mapb.index)
327 {
328 fprintf(fp, " (null)");
329 }
330 else
331 {
332 int n = p.m.mapb.nr;
333 if (nmaxind >= 0 && n > nmaxind)
334 {
335 n = nmaxind;
336 }
337 for (int i = 0; i <= n; ++i)
338 {
339 fprintf(fp, " %d", p.m.mapb.index[i]);
340 }
341 if (n < p.m.mapb.nr)
342 {
343 fprintf(fp, " ...");
344 }
345 }
346 fprintf(fp, "\n");
347
348 int n = posCount();
349 if (nmaxind >= 0 && n > nmaxind)
350 {
351 n = nmaxind;
352 }
353 fprintf(fp, " RefId:");
354 if (!p.m.refid)
355 {
356 fprintf(fp, " (null)");
357 }
358 else
359 {
360 for (int i = 0; i < n; ++i)
361 {
362 fprintf(fp, " %d", p.m.refid[i]);
363 }
364 if (n < posCount())
365 {
366 fprintf(fp, " ...");
367 }
368 }
369 fprintf(fp, "\n");
370
371 fprintf(fp, " MapId:");
372 if (!p.m.mapid)
373 {
374 fprintf(fp, " (null)");
375 }
376 else
377 {
378 for (int i = 0; i < n; ++i)
379 {
380 fprintf(fp, " %d", p.m.mapid[i]);
381 }
382 if (n < posCount())
383 {
384 fprintf(fp, " ...");
385 }
386 }
387 fprintf(fp, "\n");
388 }
389
390
391 /********************************************************************
392 * SelectionPosition
393 */
394
operator AnalysisNeighborhoodPositions() const395 SelectionPosition::operator AnalysisNeighborhoodPositions() const
396 {
397 AnalysisNeighborhoodPositions pos(sel_->rawPositions_.x, sel_->rawPositions_.count());
398 if (sel_->hasOnlyAtoms())
399 {
400 // TODO: Move atomIndices() such that it can be reused here as well.
401 pos.exclusionIds(constArrayRefFromArray<int>(sel_->rawPositions_.m.mapb.a,
402 sel_->rawPositions_.m.mapb.nra));
403 }
404 return pos.selectSingleFromArray(i_);
405 }
406
407 } // namespace gmx
408