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, 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 <string.h>
  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 = NULL;
  56     v = NULL;
  57     f = NULL;
  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, NULL, NULL, 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(), 0);
 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     /* This should not really be necessary, but do it for safety... */
 251     pos->m.mapb.index[0] = 0;
 252     pos->m.b.index[0]    = 0;
 253     /* This function should only be used to construct all the possible
 254      * positions, so the result should always be static. */
 255     pos->m.bStatic       = true;
 256 }
 257
 258 /*!
 259  * \param[in,out] pos   Position data structure.
 260  *
 261  * Sets the number of positions to 0.
 262  */
 263 void
 264 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 265 {
 266     pos->m.mapb.nr  = 0;
 267     pos->m.mapb.nra = 0;
 268     /* This should not really be necessary, but do it for safety... */
 269     pos->m.mapb.index[0] = 0;
 270     /* We set the flag to true, although really in the empty state it
 271      * should be false. This makes it possible to update the flag in
 272      * gmx_ana_pos_append(), and just make a simple check in
 273      * gmx_ana_pos_append_finish(). */
 274     pos->m.bStatic       = true;
 275 }
 276
 277 /*!
 278  * \param[in,out] dest  Data structure to which the new position is appended.
 279  * \param[in]     src   Data structure from which the position is copied.
 280  * \param[in]     i     Index in \p from to copy.
 281  */
 282 void
 283 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i)
 284 {
 285     int  j, k;
 286
 287     j = dest->count();
 288     copy_rvec(src->x[i], dest->x[j]);
 289     if (dest->v)
 290     {
 291         if (src->v)
 292         {
 293             copy_rvec(src->v[i], dest->v[j]);
 294         }
 295         else
 296         {
 297             clear_rvec(dest->v[j]);
 298         }
 299     }
 300     if (dest->f)
 301     {
 302         if (src->f)
 303         {
 304             copy_rvec(src->f[i], dest->f[j]);
 305         }
 306         else
 307         {
 308             clear_rvec(dest->f[j]);
 309         }
 310     }
 311     dest->m.refid[j] = j;
 312     dest->m.mapid[j] = src->m.mapid[i];
 313     dest->m.orgid[j] = src->m.orgid[i];
 314     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 315     {
 316         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 317         dest->m.b.a[dest->m.b.nra++]       = src->m.b.a[k];
 318     }
 319     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 320     dest->m.b.index[j+1]    = dest->m.mapb.nra;
 321     dest->m.mapb.nr++;
 322     dest->m.b.nr++;
 323 }
 324
 325 /*!
 326  * \param[in,out] dest  Data structure to which the new position is appended.
 327  * \param[in]     src   Data structure from which the position is copied.
 328  * \param[in]     i     Index in \p src to copy.
 329  * \param[in]     refid Reference ID in \p out
 330  *   (all negative values are treated as -1).
 331  */
 332 void
 333 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i, int refid)
 334 {
 335     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 336     {
 337         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 338     }
 339     const int j = dest->count();
 340     if (dest->v)
 341     {
 342         if (src->v)
 343         {
 344             copy_rvec(src->v[i], dest->v[j]);
 345         }
 346         else
 347         {
 348             clear_rvec(dest->v[j]);
 349         }
 350     }
 351     if (dest->f)
 352     {
 353         if (src->f)
 354         {
 355             copy_rvec(src->f[i], dest->f[j]);
 356         }
 357         else
 358         {
 359             clear_rvec(dest->f[j]);
 360         }
 361     }
 362     copy_rvec(src->x[i], dest->x[j]);
 363     if (refid < 0)
 364     {
 365         dest->m.refid[j] = -1;
 366         dest->m.bStatic  = false;
 367         /* If we are using masks, there is no need to alter the
 368          * mapid field. */
 369     }
 370     else
 371     {
 372         if (refid != j)
 373         {
 374             dest->m.bStatic = false;
 375         }
 376         dest->m.refid[j] = refid;
 377         /* Use the original IDs from the output structure to correctly
 378          * handle user customization. */
 379         dest->m.mapid[j] = dest->m.orgid[refid];
 380     }
 381     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 382     dest->m.mapb.nr++;
 383 }
 384
 385 /*!
 386  * \param[in,out] pos   Position data structure.
 387  *
 388  * After gmx_ana_pos_empty(), internal state of the position data structure
 389  * is not consistent before this function is called. This function should be
 390  * called after any gmx_ana_pos_append() calls have been made.
 391  */
 392 void
 393 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 394 {
 395     if (pos->m.mapb.nr != pos->m.b.nr)
 396     {
 397         pos->m.bStatic = false;
 398     }
 399 }
 400
 401 /*!
 402  * \param[in,out] g     Data structure to which the new atoms are appended.
 403  * \param[in]     src   Data structure from which the position is copied.
 404  * \param[in]     i     Index in \p src to copy.
 405  */
 406 void
 407 gmx_ana_pos_add_to_group(gmx_ana_index_t *g, gmx_ana_pos_t *src, int i)
 408 {
 409     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 410     {
 411         g->index[g->isize++] = src->m.mapb.a[k];
 412     }
 413 }