src/gromacs/selection/position.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
   5  * David van der Spoel, Berk Hess, Erik Lindahl, and including many
   6  * others, as listed in the AUTHORS file in the top-level source
   7  * directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements functions in position.h.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include <string.h>
  43
  44 #include "gromacs/legacyheaders/smalloc.h"
  45 #include "gromacs/legacyheaders/typedefs.h"
  46 #include "gromacs/legacyheaders/vec.h"
  47
  48 #include "gromacs/selection/indexutil.h"
  49 #include "gromacs/selection/position.h"
  50 #include "gromacs/utility/gmxassert.h"
  51
  52 /*!
  53  * \param[out] pos      Output structure.
  54  *
  55  * Any contents of \p pos are discarded without freeing.
  56  */
  57 void
  58 gmx_ana_pos_clear(gmx_ana_pos_t *pos)
  59 {
  60     pos->nr = 0;
  61     pos->x  = NULL;
  62     pos->v  = NULL;
  63     pos->f  = NULL;
  64     gmx_ana_indexmap_clear(&pos->m);
  65     pos->g        = NULL;
  66     pos->nalloc_x = 0;
  67 }
  68
  69 /*!
  70  * \param[in,out] pos   Position data structure.
  71  * \param[in]     n     Maximum number of positions.
  72  * \param[in]     isize Maximum number of atoms.
  73  *
  74  * Ensures that enough memory is allocated in \p pos to calculate \p n
  75  * positions from \p isize atoms.
  76  */
  77 void
  78 gmx_ana_pos_reserve(gmx_ana_pos_t *pos, int n, int isize)
  79 {
  80     GMX_RELEASE_ASSERT(n >= 0, "Invalid position allocation count");
  81     // Always reserve at least one entry to make NULL checks against pos->x
  82     // and gmx_ana_pos_reserve_velocities/forces() work as expected in the case
  83     // that there are actually no positions.
  84     if (n == 0)
  85     {
  86         n = 1;
  87     }
  88     if (pos->nalloc_x < n)
  89     {
  90         pos->nalloc_x = n;
  91         srenew(pos->x, n);
  92         if (pos->v)
  93         {
  94             srenew(pos->v, n);
  95         }
  96         if (pos->f)
  97         {
  98             srenew(pos->f, n);
  99         }
 100     }
 101     if (isize > 0)
 102     {
 103         gmx_ana_indexmap_reserve(&pos->m, n, isize);
 104     }
 105 }
 106
 107 /*!
 108  * \param[in,out] pos   Position data structure.
 109  *
 110  * Currently, this function can only be called after gmx_ana_pos_reserve()
 111  * has been called at least once with a \p n >= 0.
 112  */
 113 void
 114 gmx_ana_pos_reserve_velocities(gmx_ana_pos_t *pos)
 115 {
 116     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 117                        "No memory reserved yet for positions");
 118     if (!pos->v)
 119     {
 120         snew(pos->v, pos->nalloc_x);
 121     }
 122 }
 123
 124 /*!
 125  * \param[in,out] pos   Position data structure.
 126  *
 127  * Currently, this function can only be called after gmx_ana_pos_reserve()
 128  * has been called at least once with a \p n >= 0.
 129  */
 130 void
 131 gmx_ana_pos_reserve_forces(gmx_ana_pos_t *pos)
 132 {
 133     GMX_RELEASE_ASSERT(pos->nalloc_x > 0,
 134                        "No memory reserved yet for positions");
 135     if (!pos->f)
 136     {
 137         snew(pos->f, pos->nalloc_x);
 138     }
 139 }
 140
 141 /*!
 142  * \param[out]    pos  Position data structure to initialize.
 143  * \param[in]     x    Position vector to use.
 144  */
 145 void
 146 gmx_ana_pos_init_const(gmx_ana_pos_t *pos, const rvec x)
 147 {
 148     gmx_ana_pos_clear(pos);
 149     pos->nr = 1;
 150     snew(pos->x, 1);
 151     snew(pos->v, 1);
 152     snew(pos->f, 1);
 153     pos->nalloc_x = 1;
 154     copy_rvec(x, pos->x[0]);
 155     clear_rvec(pos->v[0]);
 156     clear_rvec(pos->f[0]);
 157     gmx_ana_indexmap_init(&pos->m, NULL, NULL, INDEX_UNKNOWN);
 158 }
 159
 160 /*!
 161  * \param[in,out] pos   Position data structure.
 162  *
 163  * Frees any memory allocated within \p pos.
 164  * The pointer \p pos itself is not freed.
 165  *
 166  * \see gmx_ana_pos_free()
 167  */
 168 void
 169 gmx_ana_pos_deinit(gmx_ana_pos_t *pos)
 170 {
 171     pos->nr               = 0;
 172     sfree(pos->x); pos->x = NULL;
 173     sfree(pos->v); pos->v = NULL;
 174     sfree(pos->f); pos->f = NULL;
 175     pos->nalloc_x         = 0;
 176     gmx_ana_indexmap_deinit(&pos->m);
 177 }
 178
 179 /*!
 180  * \param[in,out] pos   Position data structure.
 181  *
 182  * Frees any memory allocated for \p pos.
 183  * The pointer \p pos is also freed, and is invalid after the call.
 184  *
 185  * \see gmx_ana_pos_deinit()
 186  */
 187 void
 188 gmx_ana_pos_free(gmx_ana_pos_t *pos)
 189 {
 190     gmx_ana_pos_deinit(pos);
 191     sfree(pos);
 192 }
 193
 194 /*!
 195  * \param[in,out] dest   Destination positions.
 196  * \param[in]     src    Source positions.
 197  * \param[in]     bFirst If true, memory is allocated for \p dest and a full
 198  *   copy is made; otherwise, only variable parts are copied, and no memory
 199  *   is allocated.
 200  *
 201  * \p dest should have been initialized somehow (calloc() is enough).
 202  */
 203 void
 204 gmx_ana_pos_copy(gmx_ana_pos_t *dest, gmx_ana_pos_t *src, bool bFirst)
 205 {
 206     if (bFirst)
 207     {
 208         gmx_ana_pos_reserve(dest, src->nr, 0);
 209         if (src->v)
 210         {
 211             gmx_ana_pos_reserve_velocities(dest);
 212         }
 213         if (src->f)
 214         {
 215             gmx_ana_pos_reserve_forces(dest);
 216         }
 217     }
 218     dest->nr = src->nr;
 219     memcpy(dest->x, src->x, dest->nr*sizeof(*dest->x));
 220     if (dest->v)
 221     {
 222         GMX_ASSERT(src->v, "src velocities should be non-null if dest velocities are allocated");
 223         memcpy(dest->v, src->v, dest->nr*sizeof(*dest->v));
 224     }
 225     if (dest->f)
 226     {
 227         GMX_ASSERT(src->f, "src forces should be non-null if dest forces are allocated");
 228         memcpy(dest->f, src->f, dest->nr*sizeof(*dest->f));
 229     }
 230     gmx_ana_indexmap_copy(&dest->m, &src->m, bFirst);
 231     dest->g = src->g;
 232 }
 233
 234 /*!
 235  * \param[in,out] pos  Position data structure.
 236  * \param[in]     nr   Number of positions.
 237  */
 238 void
 239 gmx_ana_pos_set_nr(gmx_ana_pos_t *pos, int nr)
 240 {
 241     pos->nr = nr;
 242 }
 243
 244 /*!
 245  * \param[in,out] pos  Position data structure.
 246  * \param         g    Evaluation group.
 247  *
 248  * The old group, if any, is discarded.
 249  * Note that only a pointer to \p g is stored; it is the responsibility of
 250  * the caller to ensure that \p g is not freed while it can be accessed
 251  * through \p pos.
 252  */
 253 void
 254 gmx_ana_pos_set_evalgrp(gmx_ana_pos_t *pos, gmx_ana_index_t *g)
 255 {
 256     pos->g = g;
 257 }
 258
 259 /*!
 260  * \param[in,out] pos   Position data structure.
 261  *
 262  * Sets the number of positions to 0.
 263  */
 264 void
 265 gmx_ana_pos_empty_init(gmx_ana_pos_t *pos)
 266 {
 267     pos->nr        = 0;
 268     pos->m.nr      = 0;
 269     pos->m.mapb.nr = 0;
 270     pos->m.b.nr    = 0;
 271     pos->m.b.nra   = 0;
 272     /* This should not really be necessary, but do it for safety... */
 273     pos->m.mapb.index[0] = 0;
 274     pos->m.b.index[0]    = 0;
 275     /* This function should only be used to construct all the possible
 276      * positions, so the result should always be static. */
 277     pos->m.bStatic    = true;
 278     pos->m.bMapStatic = true;
 279 }
 280
 281 /*!
 282  * \param[in,out] pos   Position data structure.
 283  *
 284  * Sets the number of positions to 0.
 285  */
 286 void
 287 gmx_ana_pos_empty(gmx_ana_pos_t *pos)
 288 {
 289     pos->nr        = 0;
 290     pos->m.nr      = 0;
 291     pos->m.mapb.nr = 0;
 292     /* This should not really be necessary, but do it for safety... */
 293     pos->m.mapb.index[0] = 0;
 294     /* We set the flags to true, although really in the empty state they
 295      * should be false. This makes it possible to update the flags in
 296      * gmx_ana_pos_append(), and just make a simple check in
 297      * gmx_ana_pos_append_finish(). */
 298     pos->m.bStatic    = true;
 299     pos->m.bMapStatic = true;
 300 }
 301
 302 /*!
 303  * \param[in,out] dest  Data structure to which the new position is appended.
 304  * \param[in,out] g     Data structure to which the new atoms are appended.
 305  * \param[in]     src   Data structure from which the position is copied.
 306  * \param[in]     i     Index in \p from to copy.
 307  */
 308 void
 309 gmx_ana_pos_append_init(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 310                         gmx_ana_pos_t *src, int i)
 311 {
 312     int  j, k;
 313
 314     j = dest->nr;
 315     copy_rvec(src->x[i], dest->x[j]);
 316     if (dest->v)
 317     {
 318         if (src->v)
 319         {
 320             copy_rvec(src->v[i], dest->v[j]);
 321         }
 322         else
 323         {
 324             clear_rvec(dest->v[j]);
 325         }
 326     }
 327     if (dest->f)
 328     {
 329         if (src->f)
 330         {
 331             copy_rvec(src->f[i], dest->f[j]);
 332         }
 333         else
 334         {
 335             clear_rvec(dest->f[j]);
 336         }
 337     }
 338     dest->m.refid[j] = j;
 339     dest->m.mapid[j] = src->m.mapid[i];
 340     dest->m.orgid[j] = src->m.orgid[i];
 341     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 342     {
 343         g->index[g->isize++]         = src->g->index[k];
 344         dest->m.b.a[dest->m.b.nra++] = src->m.b.a[k];
 345     }
 346     dest->m.mapb.index[j+1] = g->isize;
 347     dest->m.b.index[j+1]    = g->isize;
 348     dest->nr++;
 349     dest->m.nr      = dest->nr;
 350     dest->m.mapb.nr = dest->nr;
 351     dest->m.b.nr    = dest->nr;
 352 }
 353
 354 /*!
 355  * \param[in,out] dest  Data structure to which the new position is appended
 356  *      (can be NULL, in which case only \p g is updated).
 357  * \param[in,out] g     Data structure to which the new atoms are appended.
 358  * \param[in]     src   Data structure from which the position is copied.
 359  * \param[in]     i     Index in \p src to copy.
 360  * \param[in]     refid Reference ID in \p out
 361  *   (all negative values are treated as -1).
 362  *
 363  * If \p dest is NULL, the value of \p refid is not used.
 364  */
 365 void
 366 gmx_ana_pos_append(gmx_ana_pos_t *dest, gmx_ana_index_t *g,
 367                    gmx_ana_pos_t *src, int i, int refid)
 368 {
 369     int  j, k;
 370
 371     for (k = src->m.mapb.index[i]; k < src->m.mapb.index[i+1]; ++k)
 372     {
 373         g->index[g->isize++] = src->g->index[k];
 374     }
 375     if (dest)
 376     {
 377         j = dest->nr;
 378         if (dest->v)
 379         {
 380             if (src->v)
 381             {
 382                 copy_rvec(src->v[i], dest->v[j]);
 383             }
 384             else
 385             {
 386                 clear_rvec(dest->v[j]);
 387             }
 388         }
 389         if (dest->f)
 390         {
 391             if (src->f)
 392             {
 393                 copy_rvec(src->f[i], dest->f[j]);
 394             }
 395             else
 396             {
 397                 clear_rvec(dest->f[j]);
 398             }
 399         }
 400         copy_rvec(src->x[i], dest->x[j]);
 401         if (refid < 0)
 402         {
 403             dest->m.refid[j] = -1;
 404             dest->m.bStatic  = false;
 405             /* If we are using masks, there is no need to alter the
 406              * mapid field. */
 407         }
 408         else
 409         {
 410             if (refid != j)
 411             {
 412                 dest->m.bStatic    = false;
 413                 dest->m.bMapStatic = false;
 414             }
 415             dest->m.refid[j] = refid;
 416             /* Use the original IDs from the output structure to correctly
 417              * handle user customization. */
 418             dest->m.mapid[j] = dest->m.orgid[refid];
 419         }
 420         dest->m.mapb.index[j+1] = g->isize;
 421         dest->nr++;
 422         dest->m.nr      = dest->nr;
 423         dest->m.mapb.nr = dest->nr;
 424     }
 425 }
 426
 427 /*!
 428  * \param[in,out] pos   Position data structure.
 429  *
 430  * After gmx_ana_pos_empty(), internal state of the position data structure
 431  * is not consistent before this function is called. This function should be
 432  * called after any gmx_ana_pos_append() calls have been made.
 433  */
 434 void
 435 gmx_ana_pos_append_finish(gmx_ana_pos_t *pos)
 436 {
 437     if (pos->m.nr != pos->m.b.nr)
 438     {
 439         pos->m.bStatic    = false;
 440         pos->m.bMapStatic = false;
 441     }
 442 }