src/gromacs/selection/position.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements functions in position.h.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #include <string.h>
  39
  40 #include "gromacs/legacyheaders/smalloc.h"
  41 #include "gromacs/legacyheaders/typedefs.h"
  42 #include "gromacs/legacyheaders/vec.h"
  43
  44 #include "gromacs/selection/indexutil.h"
  45 #include "gromacs/selection/position.h"
  46 #include "gromacs/utility/gmxassert.h"
  47
  48 /*!
  49  * \param[out] pos      Output structure.
  50  *
  51  * Any contents of \p pos are discarded without freeing.
  52  */
  53 void
  54 gmx_ana_pos_clear(gmx_ana_pos_t *pos)
  55 {
  56     pos->nr = 0;
  57     pos->x  = NULL;
  58     pos->v  = NULL;
  59     pos->f  = NULL;
  60     gmx_ana_indexmap_clear(&pos->m);
  61     pos->g        = NULL;
  62     pos->nalloc_x = 0;
  63 }
  64
  65 /*!
  66  * \param[in,out] pos   Position data structure.
  67  * \param[in]     n     Maximum number of positions.
  68  * \param[in]     isize Maximum number of atoms.
  69  *
  70  * Ensures that enough memory is allocated in \p pos to calculate \p n
  71  * positions from \p isize atoms.
  72  */
  73 void
  74 gmx_ana_pos_reserve(gmx_ana_pos_t *pos, int n, int isize)
  75 {
  76     GMX_RELEASE_ASSERT(n >= 0, "Invalid position allocation count");
  77     // Always reserve at least one entry to make NULL checks against pos->x
  78     // and gmx_ana_pos_reserve_velocities/forces() work as expected in the case
  79     // that there are actually no positions.
  80     if (n == 0)
  81     {
  82         n = 1;
  83     }
  84     if (pos->nalloc_x < n)
  85     {
  86         pos->nalloc_x = n;
  87         srenew(pos->x, n);
  88         if (pos->v)
  89         {
  90             srenew(pos->v, n);
  91         }
  92         if (pos->f)
  93         {
  94             srenew(pos->f, n);
  95         }
  96     }
  97     if (isize > 0)
  98     {
  99         gmx_ana_indexmap_reserve(&pos->m, n, isize);
 100     }
 101 }
 102
 103 /*!
 104  * \param[in,out] pos   Position data structure.
 105  *
 106  * Currently, this function can only be called after gmx_ana_pos_reserve()
 107  * has been called at least once with a \p n >= 0.
 108  */
 109 void
 110 gmx_ana_pos_reserve_velocities(gmx_ana_pos_t *pos)
 111 {
 112     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 113                        "No memory reserved yet for positions");
 114     if (!pos->v)
 115     {
 116         snew(pos->v, pos->nalloc_x);
 117     }
 118 }
 119
 120 /*!
 121  * \param[in,out] pos   Position data structure.
 122  *
 123  * Currently, this function can only be called after gmx_ana_pos_reserve()
 124  * has been called at least once with a \p n >= 0.
 125  */
 126 void
 127 gmx_ana_pos_reserve_forces(gmx_ana_pos_t *pos)
 128 {
 129     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 130                        "No memory reserved yet for positions");
 131     if (!pos->f)
 132     {
 133         snew(pos->f, pos->nalloc_x);
 134     }
 135 }
 136
 137 /*!
 138  * \param[out]    pos  Position data structure to initialize.
 139  * \param[in]     x    Position vector to use.
 140  */
 141 void
 142 gmx_ana_pos_init_const(gmx_ana_pos_t *pos, const rvec x)
 143 {
 144     gmx_ana_pos_clear(pos);
 145     pos->nr = 1;
 146     snew(pos->x, 1);
 147     snew(pos->v, 1);
 148     snew(pos->f, 1);
 149     pos->nalloc_x = 1;
 150     copy_rvec(x, pos->x[0]);
 151     clear_rvec(pos->v[0]);
 152     clear_rvec(pos->f[0]);
 153     gmx_ana_indexmap_init(&pos->m, NULL, NULL, INDEX_UNKNOWN);
 154 }
 155
 156 /*!
 157  * \param[in,out] pos   Position data structure.
 158  *
 159  * Frees any memory allocated within \p pos.
 160  * The pointer \p pos itself is not freed.
 161  *
 162  * \see gmx_ana_pos_free()
 163  */
 164 void
 165 gmx_ana_pos_deinit(gmx_ana_pos_t *pos)
 166 {
 167     pos->nr               = 0;
 168     sfree(pos->x); pos->x = NULL;
 169     sfree(pos->v); pos->v = NULL;
 170     sfree(pos->f); pos->f = NULL;
 171     pos->nalloc_x         = 0;
 172     gmx_ana_indexmap_deinit(&pos->m);
 173 }
 174
 175 /*!
 176  * \param[in,out] pos   Position data structure.
 177  *
 178  * Frees any memory allocated for \p pos.
 179  * The pointer \p pos is also freed, and is invalid after the call.
 180  *
 181  * \see gmx_ana_pos_deinit()
 182  */
 183 void
 184 gmx_ana_pos_free(gmx_ana_pos_t *pos)
 185 {
 186     gmx_ana_pos_deinit(pos);
 187     sfree(pos);
 188 }
 189
 190 /*!
 191  * \param[in,out] dest   Destination positions.
 192  * \param[in]     src    Source positions.
 193  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
 194  *   copy is made; otherwise, only variable parts are copied, and no memory
 195  *   is allocated.
 196  *
 197  * \p dest should have been initialized somehow (calloc() is enough).
 198  */
 199 void
 200 gmx_ana_pos_copy(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, bool bFirst)
 201 {
 202     if (bFirst)
 203     {
 204         gmx_ana_pos_reserve(dest, src->nr, 0);
 205         if (src->v)
 206         {
 207             gmx_ana_pos_reserve_velocities(dest);
 208         }
 209         if (src->f)
 210         {
 211             gmx_ana_pos_reserve_forces(dest);
 212         }
 213     }
 214     dest->nr = src->nr;
 215     memcpy(dest->x, src->x, dest->nr*sizeof(*dest->x));
 216     if (dest->v)
 217     {
 218         GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
 219         memcpy(dest->v, src->v, dest->nr*sizeof(*dest->v));
 220     }
 221     if (dest->f)
 222     {
 223         GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
 224         memcpy(dest->f, src->f, dest->nr*sizeof(*dest->f));
 225     }
 226     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 227     dest->g = src->g;
 228 }
 229
 230 /*!
 231  * \param[in,out] pos  Position data structure.
 232  * \param[in]     nr   Number of positions.
 233  */
 234 void
 235 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 236 {
 237     pos->nr = nr;
 238 }
 239
 240 /*!
 241  * \param[in,out] pos  Position data structure.
 242  * \param         g    Evaluation group.
 243  *
 244  * The old group, if any, is discarded.
 245  * Note that only a pointer to \p g is stored; it is the responsibility of
 246  * the caller to ensure that \p g is not freed while it can be accessed
 247  * through \p pos.
 248  */
 249 void
 250 gmx_ana_pos_set_evalgrp(gmx_ana_pos_t *pos, gmx_ana_index_t *g)
 251 {
 252     pos->g = g;
 253 }
 254
 255 /*!
 256  * \param[in,out] pos   Position data structure.
 257  *
 258  * Sets the number of positions to 0.
 259  */
 260 void
 261 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 262 {
 263     pos->nr        = 0;
 264     pos->m.nr      = 0;
 265     pos->m.mapb.nr = 0;
 266     pos->m.b.nr    = 0;
 267     pos->m.b.nra   = 0;
 268     /* This should not really be necessary, but do it for safety... */
 269     pos->m.mapb.index[0] = 0;
 270     pos->m.b.index[0]    = 0;
 271     /* This function should only be used to construct all the possible
 272      * positions, so the result should always be static. */
 273     pos->m.bStatic    = true;
 274     pos->m.bMapStatic = true;
 275 }
 276
 277 /*!
 278  * \param[in,out] pos   Position data structure.
 279  *
 280  * Sets the number of positions to 0.
 281  */
 282 void
 283 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 284 {
 285     pos->nr        = 0;
 286     pos->m.nr      = 0;
 287     pos->m.mapb.nr = 0;
 288     /* This should not really be necessary, but do it for safety... */
 289     pos->m.mapb.index[0] = 0;
 290     /* We set the flags to true, although really in the empty state they
 291      * should be false. This makes it possible to update the flags in
 292      * gmx_ana_pos_append(), and just make a simple check in
 293      * gmx_ana_pos_append_finish(). */
 294     pos->m.bStatic    = true;
 295     pos->m.bMapStatic = true;
 296 }
 297
 298 /*!
 299  * \param[in,out] dest  Data structure to which the new position is appended.
 300  * \param[in,out] g     Data structure to which the new atoms are appended.
 301  * \param[in]     src   Data structure from which the position is copied.
 302  * \param[in]     i     Index in \p from to copy.
 303  */
 304 void
 305 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 306                         gmx_ana_pos_t *src, int i)
 307 {
 308     int  j, k;
 309
 310     j = dest->nr;
 311     copy_rvec(src->x[i], dest->x[j]);
 312     if (dest->v)
 313     {
 314         if (src->v)
 315         {
 316             copy_rvec(src->v[i], dest->v[j]);
 317         }
 318         else
 319         {
 320             clear_rvec(dest->v[j]);
 321         }
 322     }
 323     if (dest->f)
 324     {
 325         if (src->f)
 326         {
 327             copy_rvec(src->f[i], dest->f[j]);
 328         }
 329         else
 330         {
 331             clear_rvec(dest->f[j]);
 332         }
 333     }
 334     dest->m.refid[j] = j;
 335     dest->m.mapid[j] = src->m.mapid[i];
 336     dest->m.orgid[j] = src->m.orgid[i];
 337     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 338     {
 339         g->index[g->isize++]         = src->g->index[k];
 340         dest->m.b.a[dest->m.b.nra++] = src->m.b.a[k];
 341     }
 342     dest->m.mapb.index[j+1] = g->isize;
 343     dest->m.b.index[j+1]    = g->isize;
 344     dest->nr++;
 345     dest->m.nr      = dest->nr;
 346     dest->m.mapb.nr = dest->nr;
 347     dest->m.b.nr    = dest->nr;
 348 }
 349
 350 /*!
 351  * \param[in,out] dest  Data structure to which the new position is appended
 352  *      (can be NULL, in which case only \p g is updated).
 353  * \param[in,out] g     Data structure to which the new atoms are appended.
 354  * \param[in]     src   Data structure from which the position is copied.
 355  * \param[in]     i     Index in \p src to copy.
 356  * \param[in]     refid Reference ID in \p out
 357  *   (all negative values are treated as -1).
 358  *
 359  * If \p dest is NULL, the value of \p refid is not used.
 360  */
 361 void
 362 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 363                    gmx_ana_pos_t *src, int i, int refid)
 364 {
 365     int  j, k;
 366
 367     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 368     {
 369         g->index[g->isize++] = src->g->index[k];
 370     }
 371     if (dest)
 372     {
 373         j = dest->nr;
 374         if (dest->v)
 375         {
 376             if (src->v)
 377             {
 378                 copy_rvec(src->v[i], dest->v[j]);
 379             }
 380             else
 381             {
 382                 clear_rvec(dest->v[j]);
 383             }
 384         }
 385         if (dest->f)
 386         {
 387             if (src->f)
 388             {
 389                 copy_rvec(src->f[i], dest->f[j]);
 390             }
 391             else
 392             {
 393                 clear_rvec(dest->f[j]);
 394             }
 395         }
 396         copy_rvec(src->x[i], dest->x[j]);
 397         if (refid < 0)
 398         {
 399             dest->m.refid[j] = -1;
 400             dest->m.bStatic  = false;
 401             /* If we are using masks, there is no need to alter the
 402              * mapid field. */
 403         }
 404         else
 405         {
 406             if (refid != j)
 407             {
 408                 dest->m.bStatic    = false;
 409                 dest->m.bMapStatic = false;
 410             }
 411             dest->m.refid[j] = refid;
 412             /* Use the original IDs from the output structure to correctly
 413              * handle user customization. */
 414             dest->m.mapid[j] = dest->m.orgid[refid];
 415         }
 416         dest->m.mapb.index[j+1] = g->isize;
 417         dest->nr++;
 418         dest->m.nr      = dest->nr;
 419         dest->m.mapb.nr = dest->nr;
 420     }
 421 }
 422
 423 /*!
 424  * \param[in,out] pos   Position data structure.
 425  *
 426  * After gmx_ana_pos_empty(), internal state of the position data structure
 427  * is not consistent before this function is called. This function should be
 428  * called after any gmx_ana_pos_append() calls have been made.
 429  */
 430 void
 431 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 432 {
 433     if (pos->m.nr != pos->m.b.nr)
 434     {
 435         pos->m.bStatic    = false;
 436         pos->m.bMapStatic = false;
 437     }
 438 }