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