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, 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         memcpy(dest->v, src->v, dest->nr*sizeof(*dest->v));
 209     }
 210     if (dest->f)
 211     {
 212         memcpy(dest->f, src->f, dest->nr*sizeof(*dest->f));
 213     }
 214     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 215     dest->g = src->g;
 216 }
 217
 218 /*!
 219  * \param[in,out] pos  Position data structure.
 220  * \param[in]     nr   Number of positions.
 221  */
 222 void
 223 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 224 {
 225     pos->nr = nr;
 226 }
 227
 228 /*!
 229  * \param[in,out] pos  Position data structure.
 230  * \param         g    Evaluation group.
 231  *
 232  * The old group, if any, is discarded.
 233  * Note that only a pointer to \p g is stored; it is the responsibility of
 234  * the caller to ensure that \p g is not freed while it can be accessed
 235  * through \p pos.
 236  */
 237 void
 238 gmx_ana_pos_set_evalgrp(gmx_ana_pos_t *pos, gmx_ana_index_t *g)
 239 {
 240     pos->g = g;
 241 }
 242
 243 /*!
 244  * \param[in,out] pos   Position data structure.
 245  *
 246  * Sets the number of positions to 0.
 247  */
 248 void
 249 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 250 {
 251     pos->nr = 0;
 252     pos->m.nr = 0;
 253     pos->m.mapb.nr = 0;
 254     pos->m.b.nr = 0;
 255     pos->m.b.nra = 0;
 256     /* This should not really be necessary, but do it for safety... */
 257     pos->m.mapb.index[0] = 0;
 258     pos->m.b.index[0] = 0;
 259     /* This function should only be used to construct all the possible
 260      * positions, so the result should always be static. */
 261     pos->m.bStatic = TRUE;
 262     pos->m.bMapStatic = TRUE;
 263 }
 264
 265 /*!
 266  * \param[in,out] pos   Position data structure.
 267  *
 268  * Sets the number of positions to 0.
 269  */
 270 void
 271 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 272 {
 273     pos->nr = 0;
 274     pos->m.nr = 0;
 275     pos->m.mapb.nr = 0;
 276     /* This should not really be necessary, but do it for safety... */
 277     pos->m.mapb.index[0] = 0;
 278     /* We set the flags to TRUE, although really in the empty state they
 279      * should be FALSE. This makes it possible to update the flags in
 280      * gmx_ana_pos_append(), and just make a simple check in
 281      * gmx_ana_pos_append_finish(). */
 282     pos->m.bStatic = TRUE;
 283     pos->m.bMapStatic = TRUE;
 284 }
 285
 286 /*!
 287  * \param[in,out] dest  Data structure to which the new position is appended.
 288  * \param[in,out] g     Data structure to which the new atoms are appended.
 289  * \param[in]     src   Data structure from which the position is copied.
 290  * \param[in]     i     Index in \p from to copy.
 291  */
 292 void
 293 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 294                         gmx_ana_pos_t *src, int i)
 295 {
 296     int  j, k;
 297
 298     j = dest->nr;
 299     copy_rvec(src->x[i], dest->x[j]);
 300     if (dest->v)
 301     {
 302         if (src->v)
 303         {
 304             copy_rvec(src->v[i], dest->v[j]);
 305         }
 306         else
 307         {
 308             clear_rvec(dest->v[j]);
 309         }
 310     }
 311     if (dest->f)
 312     {
 313         if (src->f)
 314         {
 315             copy_rvec(src->f[i], dest->f[j]);
 316         }
 317         else
 318         {
 319             clear_rvec(dest->f[j]);
 320         }
 321     }
 322     dest->m.refid[j] = j;
 323     dest->m.mapid[j] = src->m.mapid[i];
 324     dest->m.orgid[j] = src->m.orgid[i];
 325     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 326     {
 327         g->index[g->isize++]         = src->g->index[k];
 328         dest->m.b.a[dest->m.b.nra++] = src->m.b.a[k];
 329     }
 330     dest->m.mapb.index[j+1] = g->isize;
 331     dest->m.b.index[j+1]    = g->isize;
 332     dest->nr++;
 333     dest->m.nr = dest->nr;
 334     dest->m.mapb.nr = dest->nr;
 335     dest->m.b.nr = dest->nr;
 336 }
 337
 338 /*!
 339  * \param[in,out] dest  Data structure to which the new position is appended
 340  *      (can be NULL, in which case only \p g is updated).
 341  * \param[in,out] g     Data structure to which the new atoms are appended.
 342  * \param[in]     src   Data structure from which the position is copied.
 343  * \param[in]     i     Index in \p src to copy.
 344  * \param[in]     refid Reference ID in \p out
 345  *   (all negative values are treated as -1).
 346  *
 347  * If \p dest is NULL, the value of \p refid is not used.
 348  */
 349 void
 350 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 351                    gmx_ana_pos_t *src, int i, int refid)
 352 {
 353     int  j, k;
 354
 355     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 356     {
 357         g->index[g->isize++] = src->g->index[k];
 358     }
 359     if (dest)
 360     {
 361         j = dest->nr;
 362         if (dest->v)
 363         {
 364             if (src->v)
 365             {
 366                 copy_rvec(src->v[i], dest->v[j]);
 367             }
 368             else
 369             {
 370                 clear_rvec(dest->v[j]);
 371             }
 372         }
 373         if (dest->f)
 374         {
 375             if (src->f)
 376             {
 377                 copy_rvec(src->f[i], dest->f[j]);
 378             }
 379             else
 380             {
 381                 clear_rvec(dest->f[j]);
 382             }
 383         }
 384         copy_rvec(src->x[i], dest->x[j]);
 385         if (refid < 0)
 386         {
 387             dest->m.refid[j] = -1;
 388             dest->m.bStatic = FALSE;
 389             /* If we are using masks, there is no need to alter the
 390              * mapid field. */
 391         }
 392         else
 393         {
 394             if (refid != j)
 395             {
 396                 dest->m.bStatic = FALSE;
 397                 dest->m.bMapStatic = FALSE;
 398             }
 399             dest->m.refid[j] = refid;
 400             /* Use the original IDs from the output structure to correctly
 401              * handle user customization. */
 402             dest->m.mapid[j] = dest->m.orgid[refid];
 403         }
 404         dest->m.mapb.index[j+1] = g->isize;
 405         dest->nr++;
 406         dest->m.nr = dest->nr;
 407         dest->m.mapb.nr = dest->nr;
 408     }
 409 }
 410
 411 /*!
 412  * \param[in,out] pos   Position data structure.
 413  *
 414  * After gmx_ana_pos_empty(), internal state of the position data structure
 415  * is not consistent before this function is called. This function should be
 416  * called after any gmx_ana_pos_append() calls have been made.
 417  */
 418 void
 419 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 420 {
 421     if (pos->m.nr != pos->m.b.nr)
 422     {
 423         pos->m.bStatic = FALSE;
 424         pos->m.bMapStatic = FALSE;
 425     }
 426 }