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 }
 293
 294 /*!
 295  * \param[in,out] pos   Position data structure.
 296  *
 297  * Sets the number of positions to 0.
 298  */
 299 void
 300 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 301 {
 302     pos->nr         = 0;
 303     pos->m.nr       = 0;
 304     pos->m.mapb.nr  = 0;
 305     pos->m.mapb.nra = 0;
 306     /* This should not really be necessary, but do it for safety... */
 307     pos->m.mapb.index[0] = 0;
 308     /* We set the flag to true, although really in the empty state it
 309      * should be false. This makes it possible to update the flag in
 310      * gmx_ana_pos_append(), and just make a simple check in
 311      * gmx_ana_pos_append_finish(). */
 312     pos->m.bStatic       = true;
 313 }
 314
 315 /*!
 316  * \param[in,out] dest  Data structure to which the new position is appended.
 317  * \param[in]     src   Data structure from which the position is copied.
 318  * \param[in]     i     Index in \p from to copy.
 319  */
 320 void
 321 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i)
 322 {
 323     int  j, k;
 324
 325     j = dest->nr;
 326     copy_rvec(src->x[i], dest->x[j]);
 327     if (dest->v)
 328     {
 329         if (src->v)
 330         {
 331             copy_rvec(src->v[i], dest->v[j]);
 332         }
 333         else
 334         {
 335             clear_rvec(dest->v[j]);
 336         }
 337     }
 338     if (dest->f)
 339     {
 340         if (src->f)
 341         {
 342             copy_rvec(src->f[i], dest->f[j]);
 343         }
 344         else
 345         {
 346             clear_rvec(dest->f[j]);
 347         }
 348     }
 349     dest->m.refid[j] = j;
 350     dest->m.mapid[j] = src->m.mapid[i];
 351     dest->m.orgid[j] = src->m.orgid[i];
 352     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 353     {
 354         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 355         dest->m.b.a[dest->m.b.nra++]       = src->m.b.a[k];
 356     }
 357     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 358     dest->m.b.index[j+1]    = dest->m.mapb.nra;
 359     dest->nr++;
 360     dest->m.nr      = dest->nr;
 361     dest->m.mapb.nr = dest->nr;
 362     dest->m.b.nr    = dest->nr;
 363 }
 364
 365 /*!
 366  * \param[in,out] dest  Data structure to which the new position is appended.
 367  * \param[in]     src   Data structure from which the position is copied.
 368  * \param[in]     i     Index in \p src to copy.
 369  * \param[in]     refid Reference ID in \p out
 370  *   (all negative values are treated as -1).
 371  */
 372 void
 373 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, int i, int refid)
 374 {
 375     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 376     {
 377         dest->m.mapb.a[dest->m.mapb.nra++] = src->m.mapb.a[k];
 378     }
 379     const int j = dest->nr;
 380     if (dest->v)
 381     {
 382         if (src->v)
 383         {
 384             copy_rvec(src->v[i], dest->v[j]);
 385         }
 386         else
 387         {
 388             clear_rvec(dest->v[j]);
 389         }
 390     }
 391     if (dest->f)
 392     {
 393         if (src->f)
 394         {
 395             copy_rvec(src->f[i], dest->f[j]);
 396         }
 397         else
 398         {
 399             clear_rvec(dest->f[j]);
 400         }
 401     }
 402     copy_rvec(src->x[i], dest->x[j]);
 403     if (refid < 0)
 404     {
 405         dest->m.refid[j] = -1;
 406         dest->m.bStatic  = false;
 407         /* If we are using masks, there is no need to alter the
 408          * mapid field. */
 409     }
 410     else
 411     {
 412         if (refid != j)
 413         {
 414             dest->m.bStatic = false;
 415         }
 416         dest->m.refid[j] = refid;
 417         /* Use the original IDs from the output structure to correctly
 418          * handle user customization. */
 419         dest->m.mapid[j] = dest->m.orgid[refid];
 420     }
 421     dest->m.mapb.index[j+1] = dest->m.mapb.nra;
 422     dest->nr++;
 423     dest->m.nr      = dest->nr;
 424     dest->m.mapb.nr = dest->nr;
 425 }
 426
 427 /*!
 428  * \param[in,out] pos   Position data structure.
 429  *
 430  * After gmx_ana_pos_empty(), internal state of the position data structure
 431  * is not consistent before this function is called. This function should be
 432  * called after any gmx_ana_pos_append() calls have been made.
 433  */
 434 void
 435 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 436 {
 437     if (pos->m.nr != pos->m.b.nr)
 438     {
 439         pos->m.bStatic = false;
 440     }
 441 }
 442
 443 /*!
 444  * \param[in,out] g     Data structure to which the new atoms are appended.
 445  * \param[in]     src   Data structure from which the position is copied.
 446  * \param[in]     i     Index in \p src to copy.
 447  */
 448 void
 449 gmx_ana_pos_add_to_group(gmx_ana_index_t *g, gmx_ana_pos_t *src, int i)
 450 {
 451     for (int k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 452     {
 453         g->index[g->isize++] = src->m.mapb.a[k];
 454     }
 455 }