src/gmxlib/trajana/displacement.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 /*! \page displacements Displacement calculation routines
  32  *
  33  * Functions to calculate particle displacements as a function of time are
  34  * defined in displacement.h.
  35  * These functions allow one to easily calculate the displacement during a
  36  * single pass over the trajectory.
  37  * A new calculation is initialized with gmx_ana_displ_create().
  38  * For each frame, gmx_ana_displ_start_frame() needs to be called, followed by
  39  * a series of gmx_ana_displ_store() or a single call to
  40  * gmx_ana_displ_store_array(), gmx_ana_displ_store_all() or
  41  * gmx_ana_displ_store_pos().
  42  * The displacements for the current frame can then be calculated with
  43  * gmx_ana_displ_vector(), gmx_ana_displ_vectors() and/or
  44  * gmx_ana_displ_vectors_all().
  45  * These calculate the displacements from old positions to the ones on the
  46  * current frame.
  47  * The input frames should be evenly spaced, otherwise a fatal error is
  48  * reported. A time interval can be converted to a number of steps for
  49  * the calculation functions using gmx_ana_displ_time_to_steps().
  50  * When the structure is no longer required, gmx_ana_displ_free() can be used
  51  * to free the memory.
  52  *
  53  * The functions support calculation for dynamic selections:
  54  * when storing the positions, it is possible to specify whether a particle
  55  * is present or not.
  56  * This data is then accessible in the \p pout array returned by the
  57  * calculation functions.
  58  * Note that if you only call gmx_ana_displ_store() for the particles that are
  59  * present, use gmx_ana_displ_store_array(), or use positions calculated without
  60  * \ref POS_MASKONLY for gmx_ana_displ_store_pos(), you should ensure that you
  61  * do not use the displacements for which \p pout is FALSE (the values cannot
  62  * be calculated based on the provided data, and are undefined).
  63  */
  64 /*! \internal \file
  65  * \brief Implementation of functions in displacement.h.
  66  */
  67 #ifdef HAVE_CONFIG_H
  68 #include <config.h>
  69 #endif
  70
  71 #include <math.h>
  72
  73 #include <pbc.h>
  74 #include <smalloc.h>
  75 #include <vec.h>
  76
  77 #include <displacement.h>
  78 #include <position.h>
  79
  80 /*! \internal \brief
  81  * Internal structure used by the displacement calculation routines in
  82  * displacement.h.
  83  */
  84 typedef struct gmx_ana_displpos_t
  85 {
  86     /** Stored position vector. */
  87     rvec                 x;
  88     /** TRUE if there is something stored. */
  89     gmx_bool                 bPres;
  90 } gmx_ana_displpos_t;
  91
  92 /*! \internal \brief
  93  * Stores displacement calculation data and parameters.
  94  *
  95  * There data can be accessed only through the functions defined in
  96  * displacement.h.
  97  */
  98 struct gmx_ana_displ_t
  99 {
 100     /** Maximum number of particles for which the displacements are calculated. */
 101     int                  nmax;
 102     /** Maximum time for which the displacements are needed. */
 103     real                 tmax;
 104
 105     /** TRUE if no frames have been read. */
 106     gmx_bool                 bFirst;
 107     /** Stores the time of the first frame. */
 108     real                 t0;
 109     /** Stores the time interval between frames. */
 110     real                 dt;
 111     /** Stores the time of the current frame. */
 112     real                 t;
 113     /** Stores the index in the store for the current positions. */
 114     int                  ci;
 115
 116     /** Maximum number of positions to store for a particle. */
 117     int                  max_store;
 118     /** The total number of positions ever stored (can be larger than \p max_store). */
 119     int                  nstored;
 120     /** Two-dimensional array of stored positions, second index is the particle position. */
 121     gmx_ana_displpos_t **p;
 122     /** Pointer to the memory allocated for \p p. */
 123     gmx_ana_displpos_t  *palloc;
 124 };
 125
 126 /*!
 127  * \param[out] data Displacement calculation data structure pointer to
 128  *   initialize.
 129  * \param[in]  nmax Number of particles
 130  *   for which displacements should be calculated.
 131  * \param[in]  tmax Maximum interval (in same units as in the trajectory)
 132  *   for which displacements are required.
 133  * \returns    0 on success.
 134  *
 135  * This function only allocates memory for the first frame;
 136  * gmx_ana_displ_start_frame() allocates a sufficient amount of memory for
 137  * storing all the required positions once it known the interval between
 138  * the trajectory frames.
 139  */
 140 int
 141 gmx_ana_displ_create(gmx_ana_displ_t **data, int nmax, real tmax)
 142 {
 143     gmx_ana_displ_t *d;
 144
 145     snew(d, 1);
 146     d->nmax        = nmax;
 147     d->tmax        = tmax;
 148
 149     d->bFirst      = TRUE;
 150     d->t0          = 0;
 151     d->dt          = 0;
 152     d->t           = 0;
 153     d->ci          = -1;
 154     d->nstored     = 0;
 155
 156     d->max_store   = -1;
 157     snew(d->palloc, nmax);
 158     snew(d->p,      1);
 159     d->p[0]        = d->palloc;
 160
 161     *data = d;
 162     return 0;
 163 }
 164
 165 /*!
 166  * \param[in,out] d   Displacement calculation data.
 167  * \param[in]     t   Time for the new frame.
 168  * \returns       0 on success, -1 on failure.
 169  *
 170  * This function should be called at the beginning of each frame for which
 171  * atomic positions should be stored.
 172  * It is required that the times for which this function is called are
 173  * evenly spaced; otherwise, it reports a fatal error.
 174  */
 175 int
 176 gmx_ana_displ_start_frame(gmx_ana_displ_t *d, real t)
 177 {
 178     int   i;
 179
 180     /* Initialize times */
 181     if (d->bFirst)
 182     {
 183         d->t0 = t;
 184     }
 185     else if (d->dt <= 0)
 186     {
 187         d->dt = t - d->t0;
 188     }
 189     else
 190     {
 191         if (!gmx_within_tol(t - d->t, d->dt, GMX_REAL_EPS))
 192         {
 193             gmx_input("Trajectory not evenly spaced");
 194             return -1;
 195         }
 196     }
 197     d->t = t;
 198
 199     /* Allocate memory for all the positions once it is possible */
 200     if (d->max_store == -1 && !d->bFirst)
 201     {
 202         d->max_store = (int)(d->tmax/d->dt + 1);
 203         srenew(d->palloc, d->nmax * d->max_store);
 204         sfree(d->p);
 205         snew(d->p, d->max_store);
 206         for (i = 0; i < d->max_store; ++i)
 207         {
 208             d->p[i] = &d->palloc[d->nmax*i];
 209         }
 210     }
 211
 212     /* Increment the index where current positions are stored */
 213     d->ci++;
 214     if (d->ci >= d->max_store)
 215     {
 216         d->ci = 0;
 217     }
 218
 219     for (i = 0; i < d->nmax; ++i)
 220     {
 221         d->p[d->ci][i].bPres = FALSE;
 222     }
 223     d->nstored++;
 224     d->bFirst = FALSE;
 225     return 0;
 226 }
 227
 228 /*!
 229  * \param[in]     d     Displacement calculation data.
 230  * \param[in]     time  Time to convert to steps.
 231  * \param[out]    steps The number of steps between the frames separated by
 232  *   \p time, 0 if the number of steps cannot be determined.
 233  * \returns       0 on success (also if \p *steps is 0), -1 on error.
 234  *
 235  * \p *steps is set to zero if gmx_ana_displ_start_frame() has not yet been
 236  * called twice.
 237  * -1 is returned if \p time is not an integer multiple of the interval
 238  * between frames.
 239  */
 240 int
 241 gmx_ana_displ_time_to_steps(gmx_ana_displ_t *d, real time, int *steps)
 242 {
 243     if (d->dt <= 0)
 244     {
 245         *steps = 0;
 246         return 0;
 247     }
 248     if (!gmx_within_tol(fmod(time, d->dt), 0, GMX_REAL_EPS))
 249     {
 250         gmx_input("Interval not multiple of trajectory time step");
 251         return -1;
 252     }
 253     *steps = (int)(time / d->dt + 0.5);
 254     return 0;
 255 }
 256
 257 /*!
 258  * \param[in,out] d     Displacement calculation data.
 259  * \param[in]     id    Particle ID.
 260  * \param[in]     x     Particle position.
 261  * \param[in]     bPres TRUE if the particle should be marked as present.
 262  * \returns       0 on success.
 263  */
 264 int
 265 gmx_ana_displ_store(gmx_ana_displ_t *d, atom_id id, rvec x, gmx_bool bPres)
 266 {
 267     copy_rvec(x, d->p[d->ci][id].x);
 268     d->p[d->ci][id].bPres = bPres;
 269     return 0;
 270 }
 271
 272 /*!
 273  * \param[in,out] d     Displacement calculation data.
 274  * \param[in]     n     Number of elements in the \p id and \p x arrays.
 275  * \param[in]     id    Particle IDs.
 276  * \param[in]     x     Particle positions.
 277  * \returns       0 on success.
 278  *
 279  * The positions of the given \p n particles are stored, and all of them
 280  * are marked as present.
 281  */
 282 int
 283 gmx_ana_displ_store_array(gmx_ana_displ_t *d, int n, atom_id id[], rvec x[])
 284 {
 285     int i;
 286
 287     for (i = 0; i < n; ++i)
 288     {
 289         gmx_ana_displ_store(d, id[i], x[i], TRUE);
 290     }
 291     return 0;
 292 }
 293
 294 /*!
 295  * \param[in,out] d     Displacement calculation data.
 296  * \param[in]     id    Particle IDs.
 297  * \param[in]     x     Particle positions.
 298  * \returns       0 on success.
 299  *
 300  * Both the \p id and \p x arrays should have
 301  * \p nmax items, where \p nmax was provided for gmx_ana_displ_create(), and
 302  * the \p id[i] should be either \p i or -1.
 303  * If \p id[i]==-1, then that particle is marked as not present,
 304  * but the position is still stored.
 305  */
 306 int
 307 gmx_ana_displ_store_all(gmx_ana_displ_t *d, atom_id id[], rvec x[])
 308 {
 309     int i;
 310
 311     for (i = 0; i < d->nmax; ++i)
 312     {
 313         gmx_ana_displ_store(d, i, x[i], id[i] >= 0);
 314     }
 315     return 0;
 316 }
 317
 318 /*!
 319  * \param[in,out] d     Displacement calculation data.
 320  * \param[in]     p     Particle positions.
 321  * \returns       0 on success.
 322  *
 323  * \p p should have at most \p nmax positions, where \p nmax was provided for
 324  * gmx_ana_displ_create().
 325  * If \p has exactly \p nmax positions, gmx_ana_displ_store_all() is called,
 326  * otherwise gmx_ana_displ_store_array() is used.
 327  */
 328 int
 329 gmx_ana_displ_store_pos(gmx_ana_displ_t *d, gmx_ana_pos_t *p)
 330 {
 331     if (p->nr == d->nmax)
 332     {
 333         gmx_ana_displ_store_all(d, p->m.refid, p->x);
 334     }
 335     else
 336     {
 337         gmx_ana_displ_store_array(d, p->nr, p->m.refid, p->x);
 338     }
 339     return 0;
 340 }
 341
 342 /*! \brief
 343  * Finds the index in the displacement position storage array.
 344  *
 345  * \param[in] d    Displacement calculation data.
 346  * \param[in] step Calculate the index of positions \p steps back.
 347  * \returns   Index into the \ref gmx_ana_displ_t::p "d->p" array
 348  *   that contains the positions \p steps back, -1 on error.
 349  *
 350  * If \p step is too large, -1 is returned.
 351  */
 352 static int
 353 find_store_index(gmx_ana_displ_t *d, int step)
 354 {
 355     int si;
 356
 357     si = d->ci - step;
 358     if (si < 0)
 359     {
 360         si += d->max_store;
 361     }
 362     if (si < 0)
 363     {
 364         gmx_incons("Displacement requested for an interval too long");
 365         return -1;
 366     }
 367
 368     return si;
 369 }
 370
 371 /*!
 372  * \param[in]     d     Displacement calculation data.
 373  * \param[in]     step  Displacement is calculated from the location
 374  *   \p step steps back.
 375  * \param[in]     pbc   Periodic boundary structure
 376  *   (if NULL, PBC are not used in distance calculation).
 377  * \param[in]     id    Particle ID for which the displacement is calculated.
 378  * \param[in]     x     Current particle position.
 379  * \param[out]    xout  Displacement of the particle.
 380  * \param[out]    pout  TRUE if the old particle position was marked present.
 381  *   If \p *pout is FALSE and the old particle position was not stored,
 382  *   the value of \p xout is undefined.
 383  *   Can be NULL, in which case the value is not stored.
 384  * \return        0 on success, -1 if less than \p step frames have been
 385  *   stored or \p step <= 0, or EINVAL if \p step is too large.
 386  */
 387 int
 388 gmx_ana_displ_vector(gmx_ana_displ_t *d, int step, t_pbc *pbc,
 389                      atom_id id, rvec x, rvec xout, gmx_bool *pout)
 390 {
 391     int si;
 392
 393     if (step >= d->nstored || step <= 0)
 394     {
 395         return -1;
 396     }
 397     si = find_store_index(d, step);
 398     if (si == -1)
 399     {
 400         return EINVAL;
 401     }
 402     if (pout)
 403     {
 404         *pout = d->p[si][id].bPres;
 405     }
 406     if (pbc)
 407     {
 408         pbc_dx(pbc, x, d->p[si][id].x, xout);
 409     }
 410     else
 411     {
 412         rvec_sub(x, d->p[si][id].x, xout);
 413     }
 414     return 0;
 415 }
 416
 417 /*!
 418  * \param[in]     d     Displacement calculation data.
 419  * \param[in]     step  Displacement is calculated from the location
 420  *   \p step steps back.
 421  * \param[in]     pbc   Periodic boundary structure
 422  *   (if NULL, PBC are not used in distance calculation).
 423  * \param[in]     n     Number of particles for which to calculate.
 424  * \param[in]     id    Particle IDs for which the displacement is calculated.
 425  * \param[in]     x     Current particle positions.
 426  * \param[out]    xout  Displacement of the particles.
 427  * \param[out]    pout  TRUE if the old particle position was marked present.
 428  *   If \p pout[i] is FALSE and the old particle position was not stored,
 429  *   the value of \p xout[i] is undefined.
 430  *   Can be NULL, in which case the value is not stored.
 431  * \return        0 on success, -1 if less than \p step frames have been
 432  *   stored or \p step <= 0, or EINVAL if \p step is too large.
 433  */
 434 int
 435 gmx_ana_displ_vectors(gmx_ana_displ_t *d, int step, t_pbc *pbc,
 436                       int n, atom_id id[], rvec x[], rvec xout[], gmx_bool *pout)
 437 {
 438     int si, i;
 439
 440     if (step >= d->nstored || step <= 0)
 441     {
 442         return -1;
 443     }
 444     si = find_store_index(d, step);
 445     if (si == -1)
 446     {
 447         return EINVAL;
 448     }
 449     for (i = 0; i < n; ++i)
 450     {
 451         if (pout)
 452         {
 453             pout[i] =  d->p[si][id[i]].bPres;
 454         }
 455         if (pbc)
 456         {
 457             pbc_dx(pbc, x[i], d->p[si][id[i]].x, xout[i]);
 458         }
 459         else
 460         {
 461             rvec_sub(x[i], d->p[si][id[i]].x, xout[i]);
 462         }
 463     }
 464     return 0;
 465 }
 466
 467 /*!
 468  * \param[in]     d     Displacement calculation data.
 469  * \param[in]     step  Displacement is calculated from the location
 470  *   \p step steps back.
 471  * \param[in]     pbc   Periodic boundary structure
 472  *   (if NULL, PBC are not used in distance calculation).
 473  * \param[in]     x     Current particle positions.
 474  * \param[out]    xout  Displacement of the particles.
 475  * \param[out]    pout  TRUE if the old particle position was marked present.
 476  *   If \p pout[i] is FALSE and the old particle position was not stored,
 477  *   the value of \p xout[i] is undefined.
 478  *   Can be NULL, in which case the value is not stored.
 479  * \return        0 on success, -1 if less than \p step frames have been
 480  *   stored or \p step <= 0, or EINVAL if \p step is too large.
 481  *
 482  * The \p x array should have \p nmax items, where \p nmax is the one provided
 483  * to init_displ_calc().
 484  */
 485 int
 486 gmx_ana_displ_vectors_all(gmx_ana_displ_t *d, int step, t_pbc *pbc,
 487                           rvec x[], rvec xout[], gmx_bool *pout)
 488 {
 489     int si, i;
 490
 491     if (step >= d->nstored || step <= 0)
 492     {
 493         return -1;
 494     }
 495     si = find_store_index(d, step);
 496     if (si == -1)
 497     {
 498         return EINVAL;
 499     }
 500     for (i = 0; i < d->nmax; ++i)
 501     {
 502         if (pout)
 503         {
 504             pout[i] =  d->p[si][i].bPres;
 505         }
 506         if (pbc)
 507         {
 508             pbc_dx(pbc, x[i], d->p[si][i].x, xout[i]);
 509         }
 510         else
 511         {
 512             rvec_sub(x[i], d->p[si][i].x, xout[i]);
 513         }
 514     }
 515     return 0;
 516 }
 517
 518 /*!
 519  * \param[in,out] d  Displacement calculation data.
 520  *
 521  * After the call, the pointer \p d is invalid.
 522  */
 523 void
 524 gmx_ana_displ_free(gmx_ana_displ_t *d)
 525 {
 526     sfree(d->p);
 527     sfree(d->palloc);
 528     sfree(d);
 529 }