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