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