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     if (pos->nalloc_x < n)
  77     {
  78         pos->nalloc_x = n;
  79         srenew(pos->x, n);
  80         if (pos->v)
  81         {
  82             srenew(pos->v, n);
  83         }
  84         if (pos->f)
  85         {
  86             srenew(pos->f, n);
  87         }
  88     }
  89     if (isize > 0)
  90     {
  91         gmx_ana_indexmap_reserve(&pos->m, n, isize);
  92     }
  93 }
  94
  95 /*!
  96  * \param[in,out] pos   Position data structure.
  97  *
  98  * Currently, this function can only be called after gmx_ana_pos_reserve()
  99  * has been called at least once with a \p n > 0.
 100  */
 101 void
 102 gmx_ana_pos_reserve_velocities(gmx_ana_pos_t *pos)
 103 {
 104     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 105                        "No memory reserved yet for positions");
 106     if (!pos->v)
 107     {
 108         snew(pos->v, pos->nalloc_x);
 109     }
 110 }
 111
 112 /*!
 113  * \param[in,out] pos   Position data structure.
 114  *
 115  * Currently, this function can only be called after gmx_ana_pos_reserve()
 116  * has been called at least once with a \p n > 0.
 117  */
 118 void
 119 gmx_ana_pos_reserve_forces(gmx_ana_pos_t *pos)
 120 {
 121     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 122                        "No memory reserved yet for positions");
 123     if (!pos->f)
 124     {
 125         snew(pos->f, pos->nalloc_x);
 126     }
 127 }
 128
 129 /*!
 130  * \param[out]    pos  Position data structure to initialize.
 131  * \param[in]     x    Position vector to use.
 132  */
 133 void
 134 gmx_ana_pos_init_const(gmx_ana_pos_t *pos, const rvec x)
 135 {
 136     gmx_ana_pos_clear(pos);
 137     pos->nr = 1;
 138     snew(pos->x, 1);
 139     snew(pos->v, 1);
 140     snew(pos->f, 1);
 141     pos->nalloc_x = 1;
 142     copy_rvec(x, pos->x[0]);
 143     clear_rvec(pos->v[0]);
 144     clear_rvec(pos->f[0]);
 145     gmx_ana_indexmap_init(&pos->m, NULL, NULL, INDEX_UNKNOWN);
 146 }
 147
 148 /*!
 149  * \param[in,out] pos   Position data structure.
 150  *
 151  * Frees any memory allocated within \p pos.
 152  * The pointer \p pos itself is not freed.
 153  *
 154  * \see gmx_ana_pos_free()
 155  */
 156 void
 157 gmx_ana_pos_deinit(gmx_ana_pos_t *pos)
 158 {
 159     pos->nr = 0;
 160     sfree(pos->x); pos->x = NULL;
 161     sfree(pos->v); pos->v = NULL;
 162     sfree(pos->f); pos->f = NULL;
 163     pos->nalloc_x = 0;
 164     gmx_ana_indexmap_deinit(&pos->m);
 165 }
 166
 167 /*!
 168  * \param[in,out] pos   Position data structure.
 169  *
 170  * Frees any memory allocated for \p pos.
 171  * The pointer \p pos is also freed, and is invalid after the call.
 172  *
 173  * \see gmx_ana_pos_deinit()
 174  */
 175 void
 176 gmx_ana_pos_free(gmx_ana_pos_t *pos)
 177 {
 178     gmx_ana_pos_deinit(pos);
 179     sfree(pos);
 180 }
 181
 182 /*!
 183  * \param[in,out] dest   Destination positions.
 184  * \param[in]     src    Source positions.
 185  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
 186  *   copy is made; otherwise, only variable parts are copied, and no memory
 187  *   is allocated.
 188  *
 189  * \p dest should have been initialized somehow (calloc() is enough).
 190  */
 191 void
 192 gmx_ana_pos_copy(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, bool bFirst)
 193 {
 194     if (bFirst)
 195     {
 196         gmx_ana_pos_reserve(dest, src->nr, 0);
 197         if (src->v)
 198         {
 199             gmx_ana_pos_reserve_velocities(dest);
 200         }
 201         if (src->f)
 202         {
 203             gmx_ana_pos_reserve_forces(dest);
 204         }
 205     }
 206     dest->nr = src->nr;
 207     memcpy(dest->x, src->x, dest->nr*sizeof(*dest->x));
 208     if (dest->v)
 209     {
 210         GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
 211         memcpy(dest->v, src->v, dest->nr*sizeof(*dest->v));
 212     }
 213     if (dest->f)
 214     {
 215         GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
 216         memcpy(dest->f, src->f, dest->nr*sizeof(*dest->f));
 217     }
 218     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 219     dest->g = src->g;
 220 }
 221
 222 /*!
 223  * \param[in,out] pos  Position data structure.
 224  * \param[in]     nr   Number of positions.
 225  */
 226 void
 227 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 228 {
 229     pos->nr = nr;
 230 }
 231
 232 /*!
 233  * \param[in,out] pos  Position data structure.
 234  * \param         g    Evaluation group.
 235  *
 236  * The old group, if any, is discarded.
 237  * Note that only a pointer to \p g is stored; it is the responsibility of
 238  * the caller to ensure that \p g is not freed while it can be accessed
 239  * through \p pos.
 240  */
 241 void
 242 gmx_ana_pos_set_evalgrp(gmx_ana_pos_t *pos, gmx_ana_index_t *g)
 243 {
 244     pos->g = g;
 245 }
 246
 247 /*!
 248  * \param[in,out] pos   Position data structure.
 249  *
 250  * Sets the number of positions to 0.
 251  */
 252 void
 253 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 254 {
 255     pos->nr = 0;
 256     pos->m.nr = 0;
 257     pos->m.mapb.nr = 0;
 258     pos->m.b.nr = 0;
 259     pos->m.b.nra = 0;
 260     /* This should not really be necessary, but do it for safety... */
 261     pos->m.mapb.index[0] = 0;
 262     pos->m.b.index[0] = 0;
 263     /* This function should only be used to construct all the possible
 264      * positions, so the result should always be static. */
 265     pos->m.bStatic = true;
 266     pos->m.bMapStatic = true;
 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(gmx_ana_pos_t *pos)
 276 {
 277     pos->nr = 0;
 278     pos->m.nr = 0;
 279     pos->m.mapb.nr = 0;
 280     /* This should not really be necessary, but do it for safety... */
 281     pos->m.mapb.index[0] = 0;
 282     /* We set the flags to true, although really in the empty state they
 283      * should be false. This makes it possible to update the flags in
 284      * gmx_ana_pos_append(), and just make a simple check in
 285      * gmx_ana_pos_append_finish(). */
 286     pos->m.bStatic = true;
 287     pos->m.bMapStatic = true;
 288 }
 289
 290 /*!
 291  * \param[in,out] dest  Data structure to which the new position is appended.
 292  * \param[in,out] g     Data structure to which the new atoms are appended.
 293  * \param[in]     src   Data structure from which the position is copied.
 294  * \param[in]     i     Index in \p from to copy.
 295  */
 296 void
 297 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 298                         gmx_ana_pos_t *src, int i)
 299 {
 300     int  j, k;
 301
 302     j = dest->nr;
 303     copy_rvec(src->x[i], dest->x[j]);
 304     if (dest->v)
 305     {
 306         if (src->v)
 307         {
 308             copy_rvec(src->v[i], dest->v[j]);
 309         }
 310         else
 311         {
 312             clear_rvec(dest->v[j]);
 313         }
 314     }
 315     if (dest->f)
 316     {
 317         if (src->f)
 318         {
 319             copy_rvec(src->f[i], dest->f[j]);
 320         }
 321         else
 322         {
 323             clear_rvec(dest->f[j]);
 324         }
 325     }
 326     dest->m.refid[j] = j;
 327     dest->m.mapid[j] = src->m.mapid[i];
 328     dest->m.orgid[j] = src->m.orgid[i];
 329     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 330     {
 331         g->index[g->isize++]         = src->g->index[k];
 332         dest->m.b.a[dest->m.b.nra++] = src->m.b.a[k];
 333     }
 334     dest->m.mapb.index[j+1] = g->isize;
 335     dest->m.b.index[j+1]    = g->isize;
 336     dest->nr++;
 337     dest->m.nr = dest->nr;
 338     dest->m.mapb.nr = dest->nr;
 339     dest->m.b.nr = dest->nr;
 340 }
 341
 342 /*!
 343  * \param[in,out] dest  Data structure to which the new position is appended
 344  *      (can be NULL, in which case only \p g is updated).
 345  * \param[in,out] g     Data structure to which the new atoms are appended.
 346  * \param[in]     src   Data structure from which the position is copied.
 347  * \param[in]     i     Index in \p src to copy.
 348  * \param[in]     refid Reference ID in \p out
 349  *   (all negative values are treated as -1).
 350  *
 351  * If \p dest is NULL, the value of \p refid is not used.
 352  */
 353 void
 354 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 355                    gmx_ana_pos_t *src, int i, int refid)
 356 {
 357     int  j, k;
 358
 359     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 360     {
 361         g->index[g->isize++] = src->g->index[k];
 362     }
 363     if (dest)
 364     {
 365         j = dest->nr;
 366         if (dest->v)
 367         {
 368             if (src->v)
 369             {
 370                 copy_rvec(src->v[i], dest->v[j]);
 371             }
 372             else
 373             {
 374                 clear_rvec(dest->v[j]);
 375             }
 376         }
 377         if (dest->f)
 378         {
 379             if (src->f)
 380             {
 381                 copy_rvec(src->f[i], dest->f[j]);
 382             }
 383             else
 384             {
 385                 clear_rvec(dest->f[j]);
 386             }
 387         }
 388         copy_rvec(src->x[i], dest->x[j]);
 389         if (refid < 0)
 390         {
 391             dest->m.refid[j] = -1;
 392             dest->m.bStatic = false;
 393             /* If we are using masks, there is no need to alter the
 394              * mapid field. */
 395         }
 396         else
 397         {
 398             if (refid != j)
 399             {
 400                 dest->m.bStatic = false;
 401                 dest->m.bMapStatic = false;
 402             }
 403             dest->m.refid[j] = refid;
 404             /* Use the original IDs from the output structure to correctly
 405              * handle user customization. */
 406             dest->m.mapid[j] = dest->m.orgid[refid];
 407         }
 408         dest->m.mapb.index[j+1] = g->isize;
 409         dest->nr++;
 410         dest->m.nr = dest->nr;
 411         dest->m.mapb.nr = dest->nr;
 412     }
 413 }
 414
 415 /*!
 416  * \param[in,out] pos   Position data structure.
 417  *
 418  * After gmx_ana_pos_empty(), internal state of the position data structure
 419  * is not consistent before this function is called. This function should be
 420  * called after any gmx_ana_pos_append() calls have been made.
 421  */
 422 void
 423 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 424 {
 425     if (pos->m.nr != pos->m.b.nr)
 426     {
 427         pos->m.bStatic = false;
 428         pos->m.bMapStatic = false;
 429     }
 430 }