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->nr = 0;
  61     pos->x  = NULL;
  62     pos->v  = NULL;
  63     pos->f  = NULL;
  64     gmx_ana_indexmap_clear(&pos->m);
  65     pos->nalloc_x = 0;
  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     gmx_ana_pos_clear(pos);
 178     pos->nr = 1;
 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] pos   Position data structure.
 191  *
 192  * Frees any memory allocated within \p pos.
 193  * The pointer \p pos itself is not freed.
 194  *
 195  * \see gmx_ana_pos_free()
 196  */
 197 void
 198 gmx_ana_pos_deinit(gmx_ana_pos_t *pos)
 199 {
 200     pos->nr               = 0;
 201     sfree(pos->x); pos->x = NULL;
 202     sfree(pos->v); pos->v = NULL;
 203     sfree(pos->f); pos->f = NULL;
 204     pos->nalloc_x         = 0;
 205     gmx_ana_indexmap_deinit(&pos->m);
 206 }
 207
 208 /*!
 209  * \param[in,out] pos   Position data structure.
 210  *
 211  * Frees any memory allocated for \p pos.
 212  * The pointer \p pos is also freed, and is invalid after the call.
 213  *
 214  * \see gmx_ana_pos_deinit()
 215  */
 216 void
 217 gmx_ana_pos_free(gmx_ana_pos_t *pos)
 218 {
 219     gmx_ana_pos_deinit(pos);
 220     sfree(pos);
 221 }
 222
 223 /*!
 224  * \param[in,out] dest   Destination positions.
 225  * \param[in]     src    Source positions.
 226  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
 227  *   copy is made; otherwise, only variable parts are copied, and no memory
 228  *   is allocated.
 229  *
 230  * \p dest should have been initialized somehow (calloc() is enough).
 231  */
 232 void
 233 gmx_ana_pos_copy(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, bool bFirst)
 234 {
 235     if (bFirst)
 236     {
 237         gmx_ana_pos_reserve(dest, src->nr, 0);
 238         if (src->v)
 239         {
 240             gmx_ana_pos_reserve_velocities(dest);
 241         }
 242         if (src->f)
 243         {
 244             gmx_ana_pos_reserve_forces(dest);
 245         }
 246     }
 247     dest->nr = src->nr;
 248     memcpy(dest->x, src->x, dest->nr*sizeof(*dest->x));
 249     if (dest->v)
 250     {
 251         GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
 252         memcpy(dest->v, src->v, dest->nr*sizeof(*dest->v));
 253     }
 254     if (dest->f)
 255     {
 256         GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
 257         memcpy(dest->f, src->f, dest->nr*sizeof(*dest->f));
 258     }
 259     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 260 }
 261
 262 /*!
 263  * \param[in,out] pos  Position data structure.
 264  * \param[in]     nr   Number of positions.
 265  */
 266 void
 267 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 268 {
 269     pos->nr = nr;
 270 }
 271
 272 /*!
 273  * \param[in,out] pos   Position data structure.
 274  *
 275  * Sets the number of positions to 0.
 276  */
 277 void
 278 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 279 {
 280     pos->nr         = 0;
 281     pos->m.nr       = 0;
 282     pos->m.mapb.nr  = 0;
 283     pos->m.mapb.nra = 0;
 284     pos->m.b.nr     = 0;
 285     pos->m.b.nra    = 0;
 286     /* This should not really be necessary, but do it for safety... */
 287     pos->m.mapb.index[0] = 0;
 288     pos->m.b.index[0]    = 0;
 289     /* This function should only be used to construct all the possible
 290      * positions, so the result should always be static. */
 291     pos->m.bStatic    = true;
 292     pos->m.bMapStatic = true;
 293 }
 294
 295 /*!
 296  * \param[in,out] pos   Position data structure.
 297  *
 298  * Sets the number of positions to 0.
 299  */
 300 void
 301 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 302 {
 303     pos->nr         = 0;
 304     pos->m.nr       = 0;
 305     pos->m.mapb.nr  = 0;
 306     pos->m.mapb.nra = 0;
 307     /* This should not really be necessary, but do it for safety... */
 308     pos->m.mapb.index[0] = 0;
 309     /* We set the flags to true, although really in the empty state they
 310      * should be false. This makes it possible to update the flags in
 311      * gmx_ana_pos_append(), and just make a simple check in
 312      * gmx_ana_pos_append_finish(). */
 313     pos->m.bStatic    = true;
 314     pos->m.bMapStatic = true;
 315 }
 316
 317 /*!
 318  * \param[in,out] dest  Data structure to which the new position is appended.
 319  * \param[in]     src   Data structure from which the position is copied.
 320  * \param[in]     i     Index in \p from to copy.
 321  */
 322 void
 323 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i)
 324 {
 325     int  j, k;
 326
 327     j = dest->nr;
 328     copy_rvec(src->x[i], dest->x[j]);
 329     if (dest->v)
 330     {
 331         if (src->v)
 332         {
 333             copy_rvec(src->v[i], dest->v[j]);
 334         }
 335         else
 336         {
 337             clear_rvec(dest->v[j]);
 338         }
 339     }
 340     if (dest->f)
 341     {
 342         if (src->f)
 343         {
 344             copy_rvec(src->f[i], dest->f[j]);
 345         }
 346         else
 347         {
 348             clear_rvec(dest->f[j]);
 349         }
 350     }
 351     dest->m.refid[j] = j;
 352     dest->m.mapid[j] = src->m.mapid[i];
 353     dest->m.orgid[j] = src->m.orgid[i];
 354     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 355     {
 356         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 357         dest->m.b.a[dest->m.b.nra++]       = src->m.b.a[k];
 358     }
 359     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 360     dest->m.b.index[j+1]    = dest->m.mapb.nra;
 361     dest->nr++;
 362     dest->m.nr      = dest->nr;
 363     dest->m.mapb.nr = dest->nr;
 364     dest->m.b.nr    = dest->nr;
 365 }
 366
 367 /*!
 368  * \param[in,out] dest  Data structure to which the new position is appended.
 369  * \param[in]     src   Data structure from which the position is copied.
 370  * \param[in]     i     Index in \p src to copy.
 371  * \param[in]     refid Reference ID in \p out
 372  *   (all negative values are treated as -1).
 373  */
 374 void
 375 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i, int refid)
 376 {
 377     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 378     {
 379         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 380     }
 381     const int j = dest->nr;
 382     if (dest->v)
 383     {
 384         if (src->v)
 385         {
 386             copy_rvec(src->v[i], dest->v[j]);
 387         }
 388         else
 389         {
 390             clear_rvec(dest->v[j]);
 391         }
 392     }
 393     if (dest->f)
 394     {
 395         if (src->f)
 396         {
 397             copy_rvec(src->f[i], dest->f[j]);
 398         }
 399         else
 400         {
 401             clear_rvec(dest->f[j]);
 402         }
 403     }
 404     copy_rvec(src->x[i], dest->x[j]);
 405     if (refid < 0)
 406     {
 407         dest->m.refid[j] = -1;
 408         dest->m.bStatic  = false;
 409         /* If we are using masks, there is no need to alter the
 410          * mapid field. */
 411     }
 412     else
 413     {
 414         if (refid != j)
 415         {
 416             dest->m.bStatic    = false;
 417             dest->m.bMapStatic = false;
 418         }
 419         dest->m.refid[j] = refid;
 420         /* Use the original IDs from the output structure to correctly
 421          * handle user customization. */
 422         dest->m.mapid[j] = dest->m.orgid[refid];
 423     }
 424     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 425     dest->nr++;
 426     dest->m.nr      = dest->nr;
 427     dest->m.mapb.nr = dest->nr;
 428 }
 429
 430 /*!
 431  * \param[in,out] pos   Position data structure.
 432  *
 433  * After gmx_ana_pos_empty(), internal state of the position data structure
 434  * is not consistent before this function is called. This function should be
 435  * called after any gmx_ana_pos_append() calls have been made.
 436  */
 437 void
 438 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 439 {
 440     if (pos->m.nr != pos->m.b.nr)
 441     {
 442         pos->m.bStatic    = false;
 443         pos->m.bMapStatic = false;
 444     }
 445 }
 446
 447 /*!
 448  * \param[in,out] g     Data structure to which the new atoms are appended.
 449  * \param[in]     src   Data structure from which the position is copied.
 450  * \param[in]     i     Index in \p src to copy.
 451  */
 452 void
 453 gmx_ana_pos_add_to_group(gmx_ana_index_t *g, gmx_ana_pos_t *src, int i)
 454 {
 455     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 456     {
 457         g->index[g->isize++] = src->m.mapb.a[k];
 458     }
 459 }