src/gmxlib/trajana/position.c

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