src/gmxlib/trajana/displacement.c

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