2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009-2018, The GROMACS development team.
5 * Copyright (c) 2019, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 * Implements functions in position.h.
40 * \author Teemu Murtola <teemu.murtola@gmail.com>
41 * \ingroup module_selection
49 #include "gromacs/math/vec.h"
50 #include "gromacs/selection/indexutil.h"
51 #include "gromacs/utility/gmxassert.h"
52 #include "gromacs/utility/smalloc.h"
54 gmx_ana_pos_t::gmx_ana_pos_t()
59 gmx_ana_indexmap_clear(&m);
63 gmx_ana_pos_t::~gmx_ana_pos_t()
68 gmx_ana_indexmap_deinit(&m);
72 * \param[in,out] pos Position data structure.
73 * \param[in] n Maximum number of positions.
74 * \param[in] isize Maximum number of atoms.
76 * Ensures that enough memory is allocated in \p pos to calculate \p n
77 * positions from \p isize atoms.
79 void gmx_ana_pos_reserve(gmx_ana_pos_t* pos, int n, int isize)
81 GMX_RELEASE_ASSERT(n >= 0, "Invalid position allocation count");
82 // Always reserve at least one entry to make NULL checks against pos->x
83 // and gmx_ana_pos_reserve_velocities/forces() work as expected in the case
84 // that there are actually no positions.
89 if (pos->nalloc_x < n)
104 gmx_ana_indexmap_reserve(&pos->m, n, isize);
109 * \param[in,out] pos Position data structure.
111 * Currently, this function can only be called after gmx_ana_pos_reserve()
112 * has been called at least once with a \p n >= 0.
114 void gmx_ana_pos_reserve_velocities(gmx_ana_pos_t* pos)
116 GMX_RELEASE_ASSERT(pos->nalloc_x > 0, "No memory reserved yet for positions");
119 snew(pos->v, pos->nalloc_x);
124 * \param[in,out] pos Position data structure.
126 * Currently, this function can only be called after gmx_ana_pos_reserve()
127 * has been called at least once with a \p n >= 0.
129 void gmx_ana_pos_reserve_forces(gmx_ana_pos_t* pos)
131 GMX_RELEASE_ASSERT(pos->nalloc_x > 0, "No memory reserved yet for positions");
134 snew(pos->f, pos->nalloc_x);
139 * \param[in,out] pos Position data structure.
140 * \param[in] n Maximum number of positions.
141 * \param[in] isize Maximum number of atoms.
142 * \param[in] bVelocities Whether to reserve space for velocities.
143 * \param[in] bForces Whether to reserve space for forces.
145 * Ensures that enough memory is allocated in \p pos to calculate \p n
146 * positions from \p isize atoms.
148 * This method needs to be called instead of gmx_ana_pos_reserve() if the
149 * intent is to use gmx_ana_pos_append_init()/gmx_ana_pos_append().
151 void gmx_ana_pos_reserve_for_append(gmx_ana_pos_t* pos, int n, int isize, bool bVelocities, bool bForces)
153 gmx_ana_pos_reserve(pos, n, isize);
154 snew(pos->m.mapb.a, isize);
155 pos->m.mapb.nalloc_a = isize;
158 gmx_ana_pos_reserve_velocities(pos);
162 gmx_ana_pos_reserve_forces(pos);
167 * \param[out] pos Position data structure to initialize.
168 * \param[in] x Position vector to use.
170 void gmx_ana_pos_init_const(gmx_ana_pos_t* pos, const rvec x)
176 copy_rvec(x, pos->x[0]);
177 clear_rvec(pos->v[0]);
178 clear_rvec(pos->f[0]);
179 gmx_ana_indexmap_init(&pos->m, nullptr, nullptr, INDEX_UNKNOWN);
183 * \param[in,out] dest Destination positions.
184 * \param[in] src Source positions.
185 * \param[in] bFirst If true, memory is allocated for \p dest and a full
186 * copy is made; otherwise, only variable parts are copied, and no memory
189 * \p dest should have been initialized somehow (calloc() is enough).
191 void gmx_ana_pos_copy(gmx_ana_pos_t* dest, gmx_ana_pos_t* src, bool bFirst)
195 gmx_ana_pos_reserve(dest, src->count(), -1);
198 gmx_ana_pos_reserve_velocities(dest);
202 gmx_ana_pos_reserve_forces(dest);
205 memcpy(dest->x, src->x, src->count() * sizeof(*dest->x));
208 GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
209 memcpy(dest->v, src->v, src->count() * sizeof(*dest->v));
213 GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
214 memcpy(dest->f, src->f, src->count() * sizeof(*dest->f));
216 gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
220 * \param[in,out] pos Position data structure.
221 * \param[in] nr Number of positions.
223 void gmx_ana_pos_set_nr(gmx_ana_pos_t* pos, int nr)
225 // TODO: This puts the mapping in a somewhat inconsistent state.
230 * \param[in,out] pos Position data structure.
232 * Sets the number of positions to 0.
234 void gmx_ana_pos_empty_init(gmx_ana_pos_t* pos)
240 /* Initializing these should not really be necessary, but do it for
242 pos->m.mapb.index[0] = 0;
243 pos->m.b.index[0] = 0;
244 /* This function should only be used to construct all the possible
245 * positions, so the result should always be static. */
246 pos->m.bStatic = true;
250 * \param[in,out] pos Position data structure.
252 * Sets the number of positions to 0.
254 void gmx_ana_pos_empty(gmx_ana_pos_t* pos)
258 /* This should not really be necessary, but do it for safety... */
259 pos->m.mapb.index[0] = 0;
260 /* We set the flag to true, although really in the empty state it
261 * should be false. This makes it possible to update the flag in
262 * gmx_ana_pos_append(), and just make a simple check in
263 * gmx_ana_pos_append_finish(). */
264 pos->m.bStatic = true;
268 * \param[in,out] dest Data structure to which the new position is appended.
269 * \param[in] src Data structure from which the position is copied.
270 * \param[in] i Index in \p from to copy.
272 void gmx_ana_pos_append_init(gmx_ana_pos_t* dest, gmx_ana_pos_t* src, int i)
277 copy_rvec(src->x[i], dest->x[j]);
282 copy_rvec(src->v[i], dest->v[j]);
286 clear_rvec(dest->v[j]);
293 copy_rvec(src->f[i], dest->f[j]);
297 clear_rvec(dest->f[j]);
300 dest->m.refid[j] = j;
301 dest->m.mapid[j] = src->m.mapid[i];
302 dest->m.orgid[j] = src->m.orgid[i];
303 for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i + 1]; ++k)
305 dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
306 dest->m.b.a[dest->m.b.nra++] = src->m.b.a[k];
308 dest->m.mapb.index[j + 1] = dest->m.mapb.nra;
309 dest->m.b.index[j + 1] = dest->m.mapb.nra;
315 * \param[in,out] dest Data structure to which the new position is appended.
316 * \param[in] src Data structure from which the position is copied.
317 * \param[in] i Index in \p src to copy.
318 * \param[in] refid Reference ID in \p out
319 * (all negative values are treated as -1).
321 void gmx_ana_pos_append(gmx_ana_pos_t* dest, gmx_ana_pos_t* src, int i, int refid)
323 for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i + 1]; ++k)
325 dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
327 const int j = dest->count();
332 copy_rvec(src->v[i], dest->v[j]);
336 clear_rvec(dest->v[j]);
343 copy_rvec(src->f[i], dest->f[j]);
347 clear_rvec(dest->f[j]);
350 copy_rvec(src->x[i], dest->x[j]);
353 dest->m.refid[j] = -1;
354 dest->m.bStatic = false;
355 /* If we are using masks, there is no need to alter the
362 dest->m.bStatic = false;
364 dest->m.refid[j] = refid;
365 /* Use the original IDs from the output structure to correctly
366 * handle user customization. */
367 dest->m.mapid[j] = dest->m.orgid[refid];
369 dest->m.mapb.index[j + 1] = dest->m.mapb.nra;
374 * \param[in,out] pos Position data structure.
376 * After gmx_ana_pos_empty(), internal state of the position data structure
377 * is not consistent before this function is called. This function should be
378 * called after any gmx_ana_pos_append() calls have been made.
380 void gmx_ana_pos_append_finish(gmx_ana_pos_t* pos)
382 if (pos->m.mapb.nr != pos->m.b.nr)
384 pos->m.bStatic = false;
389 * \param[in,out] g Data structure to which the new atoms are appended.
390 * \param[in] src Data structure from which the position is copied.
391 * \param[in] i Index in \p src to copy.
393 void gmx_ana_pos_add_to_group(gmx_ana_index_t* g, gmx_ana_pos_t* src, int i)
395 for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i + 1]; ++k)
397 g->index[g->isize++] = src->m.mapb.a[k];