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