src/gromacs/selection/position.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements functions in position.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "position.h"
  45
  46 #include <cstring>
  47
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/selection/indexutil.h"
  50 #include "gromacs/utility/gmxassert.h"
  51 #include "gromacs/utility/smalloc.h"
  52
  53 gmx_ana_pos_t::gmx_ana_pos_t()
  54 {
  55     x = nullptr;
  56     v = nullptr;
  57     f = nullptr;
  58     gmx_ana_indexmap_clear(&m);
  59     nalloc_x = 0;
  60 }
  61
  62 gmx_ana_pos_t::~gmx_ana_pos_t()
  63 {
  64     sfree(x);
  65     sfree(v);
  66     sfree(f);
  67     gmx_ana_indexmap_deinit(&m);
  68 }
  69
  70 /*!
  71  * \param[in,out] pos   Position data structure.
  72  * \param[in]     n     Maximum number of positions.
  73  * \param[in]     isize Maximum number of atoms.
  74  *
  75  * Ensures that enough memory is allocated in \p pos to calculate \p n
  76  * positions from \p isize atoms.
  77  */
  78 void
  79 gmx_ana_pos_reserve(gmx_ana_pos_t *pos, int n, int isize)
  80 {
  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.
  85     if (n == 0)
  86     {
  87         n = 1;
  88     }
  89     if (pos->nalloc_x < n)
  90     {
  91         pos->nalloc_x = n;
  92         srenew(pos->x, n);
  93         if (pos->v)
  94         {
  95             srenew(pos->v, n);
  96         }
  97         if (pos->f)
  98         {
  99             srenew(pos->f, n);
 100         }
 101     }
 102     if (isize >= 0)
 103     {
 104         gmx_ana_indexmap_reserve(&pos->m, n, isize);
 105     }
 106 }
 107
 108 /*!
 109  * \param[in,out] pos   Position data structure.
 110  *
 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.
 113  */
 114 void
 115 gmx_ana_pos_reserve_velocities(gmx_ana_pos_t *pos)
 116 {
 117     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 118                        "No memory reserved yet for positions");
 119     if (!pos->v)
 120     {
 121         snew(pos->v, pos->nalloc_x);
 122     }
 123 }
 124
 125 /*!
 126  * \param[in,out] pos   Position data structure.
 127  *
 128  * Currently, this function can only be called after gmx_ana_pos_reserve()
 129  * has been called at least once with a \p n >= 0.
 130  */
 131 void
 132 gmx_ana_pos_reserve_forces(gmx_ana_pos_t *pos)
 133 {
 134     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 135                        "No memory reserved yet for positions");
 136     if (!pos->f)
 137     {
 138         snew(pos->f, pos->nalloc_x);
 139     }
 140 }
 141
 142 /*!
 143  * \param[in,out] pos   Position data structure.
 144  * \param[in]     n     Maximum number of positions.
 145  * \param[in]     isize Maximum number of atoms.
 146  * \param[in]     bVelocities Whether to reserve space for velocities.
 147  * \param[in]     bForces     Whether to reserve space for forces.
 148  *
 149  * Ensures that enough memory is allocated in \p pos to calculate \p n
 150  * positions from \p isize atoms.
 151  *
 152  * This method needs to be called instead of gmx_ana_pos_reserve() if the
 153  * intent is to use gmx_ana_pos_append_init()/gmx_ana_pos_append().
 154  */
 155 void
 156 gmx_ana_pos_reserve_for_append(gmx_ana_pos_t *pos, int n, int isize,
 157                                bool bVelocities, bool bForces)
 158 {
 159     gmx_ana_pos_reserve(pos, n, isize);
 160     snew(pos->m.mapb.a, isize);
 161     pos->m.mapb.nalloc_a = isize;
 162     if (bVelocities)
 163     {
 164         gmx_ana_pos_reserve_velocities(pos);
 165     }
 166     if (bForces)
 167     {
 168         gmx_ana_pos_reserve_forces(pos);
 169     }
 170 }
 171
 172 /*!
 173  * \param[out]    pos  Position data structure to initialize.
 174  * \param[in]     x    Position vector to use.
 175  */
 176 void
 177 gmx_ana_pos_init_const(gmx_ana_pos_t *pos, const rvec x)
 178 {
 179     snew(pos->x, 1);
 180     snew(pos->v, 1);
 181     snew(pos->f, 1);
 182     pos->nalloc_x = 1;
 183     copy_rvec(x, pos->x[0]);
 184     clear_rvec(pos->v[0]);
 185     clear_rvec(pos->f[0]);
 186     gmx_ana_indexmap_init(&pos->m, nullptr, nullptr, INDEX_UNKNOWN);
 187 }
 188
 189 /*!
 190  * \param[in,out] dest   Destination positions.
 191  * \param[in]     src    Source positions.
 192  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
 193  *   copy is made; otherwise, only variable parts are copied, and no memory
 194  *   is allocated.
 195  *
 196  * \p dest should have been initialized somehow (calloc() is enough).
 197  */
 198 void
 199 gmx_ana_pos_copy(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, bool bFirst)
 200 {
 201     if (bFirst)
 202     {
 203         gmx_ana_pos_reserve(dest, src->count(), -1);
 204         if (src->v)
 205         {
 206             gmx_ana_pos_reserve_velocities(dest);
 207         }
 208         if (src->f)
 209         {
 210             gmx_ana_pos_reserve_forces(dest);
 211         }
 212     }
 213     memcpy(dest->x, src->x, src->count()*sizeof(*dest->x));
 214     if (dest->v)
 215     {
 216         GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
 217         memcpy(dest->v, src->v, src->count()*sizeof(*dest->v));
 218     }
 219     if (dest->f)
 220     {
 221         GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
 222         memcpy(dest->f, src->f, src->count()*sizeof(*dest->f));
 223     }
 224     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 225 }
 226
 227 /*!
 228  * \param[in,out] pos  Position data structure.
 229  * \param[in]     nr   Number of positions.
 230  */
 231 void
 232 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 233 {
 234     // TODO: This puts the mapping in a somewhat inconsistent state.
 235     pos->m.mapb.nr = nr;
 236 }
 237
 238 /*!
 239  * \param[in,out] pos   Position data structure.
 240  *
 241  * Sets the number of positions to 0.
 242  */
 243 void
 244 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 245 {
 246     pos->m.mapb.nr  = 0;
 247     pos->m.mapb.nra = 0;
 248     pos->m.b.nr     = 0;
 249     pos->m.b.nra    = 0;
 250     /* Initializing these should not really be necessary, but do it for
 251      * safety... */
 252     pos->m.mapb.index[0] = 0;
 253     pos->m.b.index[0]    = 0;
 254     /* This function should only be used to construct all the possible
 255      * positions, so the result should always be static. */
 256     pos->m.bStatic       = true;
 257 }
 258
 259 /*!
 260  * \param[in,out] pos   Position data structure.
 261  *
 262  * Sets the number of positions to 0.
 263  */
 264 void
 265 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 266 {
 267     pos->m.mapb.nr  = 0;
 268     pos->m.mapb.nra = 0;
 269     /* This should not really be necessary, but do it for safety... */
 270     pos->m.mapb.index[0] = 0;
 271     /* We set the flag to true, although really in the empty state it
 272      * should be false. This makes it possible to update the flag in
 273      * gmx_ana_pos_append(), and just make a simple check in
 274      * gmx_ana_pos_append_finish(). */
 275     pos->m.bStatic       = true;
 276 }
 277
 278 /*!
 279  * \param[in,out] dest  Data structure to which the new position is appended.
 280  * \param[in]     src   Data structure from which the position is copied.
 281  * \param[in]     i     Index in \p from to copy.
 282  */
 283 void
 284 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i)
 285 {
 286     int  j, k;
 287
 288     j = dest->count();
 289     copy_rvec(src->x[i], dest->x[j]);
 290     if (dest->v)
 291     {
 292         if (src->v)
 293         {
 294             copy_rvec(src->v[i], dest->v[j]);
 295         }
 296         else
 297         {
 298             clear_rvec(dest->v[j]);
 299         }
 300     }
 301     if (dest->f)
 302     {
 303         if (src->f)
 304         {
 305             copy_rvec(src->f[i], dest->f[j]);
 306         }
 307         else
 308         {
 309             clear_rvec(dest->f[j]);
 310         }
 311     }
 312     dest->m.refid[j] = j;
 313     dest->m.mapid[j] = src->m.mapid[i];
 314     dest->m.orgid[j] = src->m.orgid[i];
 315     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 316     {
 317         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 318         dest->m.b.a[dest->m.b.nra++]       = src->m.b.a[k];
 319     }
 320     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 321     dest->m.b.index[j+1]    = dest->m.mapb.nra;
 322     dest->m.mapb.nr++;
 323     dest->m.b.nr++;
 324 }
 325
 326 /*!
 327  * \param[in,out] dest  Data structure to which the new position is appended.
 328  * \param[in]     src   Data structure from which the position is copied.
 329  * \param[in]     i     Index in \p src to copy.
 330  * \param[in]     refid Reference ID in \p out
 331  *   (all negative values are treated as -1).
 332  */
 333 void
 334 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i, int refid)
 335 {
 336     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 337     {
 338         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 339     }
 340     const int j = dest->count();
 341     if (dest->v)
 342     {
 343         if (src->v)
 344         {
 345             copy_rvec(src->v[i], dest->v[j]);
 346         }
 347         else
 348         {
 349             clear_rvec(dest->v[j]);
 350         }
 351     }
 352     if (dest->f)
 353     {
 354         if (src->f)
 355         {
 356             copy_rvec(src->f[i], dest->f[j]);
 357         }
 358         else
 359         {
 360             clear_rvec(dest->f[j]);
 361         }
 362     }
 363     copy_rvec(src->x[i], dest->x[j]);
 364     if (refid < 0)
 365     {
 366         dest->m.refid[j] = -1;
 367         dest->m.bStatic  = false;
 368         /* If we are using masks, there is no need to alter the
 369          * mapid field. */
 370     }
 371     else
 372     {
 373         if (refid != j)
 374         {
 375             dest->m.bStatic = false;
 376         }
 377         dest->m.refid[j] = refid;
 378         /* Use the original IDs from the output structure to correctly
 379          * handle user customization. */
 380         dest->m.mapid[j] = dest->m.orgid[refid];
 381     }
 382     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 383     dest->m.mapb.nr++;
 384 }
 385
 386 /*!
 387  * \param[in,out] pos   Position data structure.
 388  *
 389  * After gmx_ana_pos_empty(), internal state of the position data structure
 390  * is not consistent before this function is called. This function should be
 391  * called after any gmx_ana_pos_append() calls have been made.
 392  */
 393 void
 394 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 395 {
 396     if (pos->m.mapb.nr != pos->m.b.nr)
 397     {
 398         pos->m.bStatic = false;
 399     }
 400 }
 401
 402 /*!
 403  * \param[in,out] g     Data structure to which the new atoms are appended.
 404  * \param[in]     src   Data structure from which the position is copied.
 405  * \param[in]     i     Index in \p src to copy.
 406  */
 407 void
 408 gmx_ana_pos_add_to_group(gmx_ana_index_t *g, gmx_ana_pos_t *src, int i)
 409 {
 410     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 411     {
 412         g->index[g->isize++] = src->m.mapb.a[k];
 413     }
 414 }