src/gromacs/selection/sm_merge.cpp

   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
  33  * Implements the \p merge and \p plus selection modifiers.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #include "gromacs/legacyheaders/macros.h"
  39 #include "gromacs/legacyheaders/smalloc.h"
  40 #include "gromacs/legacyheaders/vec.h"
  41
  42 #include "gromacs/selection/position.h"
  43 #include "gromacs/selection/selmethod.h"
  44 #include "gromacs/utility/exceptions.h"
  45
  46 /*! \internal \brief
  47  * Data structure for the merging selection modifiers.
  48  */
  49 typedef struct
  50 {
  51     /** Input positions. */
  52     gmx_ana_pos_t    p1;
  53     /** Other input positions. */
  54     gmx_ana_pos_t    p2;
  55     /** Group to store the output atom indices. */
  56     gmx_ana_index_t  g;
  57     /** Stride for merging (\c stride values from \c p1 for each in \c p2). */
  58     int              stride;
  59 } t_methoddata_merge;
  60
  61 /** Allocates data for the merging selection modifiers. */
  62 static void *
  63 init_data_merge(int npar, gmx_ana_selparam_t *param);
  64 /** Initializes data for the merging selection modifiers. */
  65 static void
  66 init_merge(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  67 /** Initializes output for the \p merge selection modifier. */
  68 static void
  69 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  70 /** Initializes output for the \p plus selection modifier. */
  71 static void
  72 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  73 /** Frees the memory allocated for the merging selection modifiers. */
  74 static void
  75 free_data_merge(void *data);
  76 /** Evaluates the \p merge selection modifier. */
  77 static void
  78 evaluate_merge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  79                gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
  80 /** Evaluates the \p plus selection modifier. */
  81 static void
  82 evaluate_plus(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  83               gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
  84
  85 /** Parameters for the merging selection modifiers. */
  86 static gmx_ana_selparam_t smparams_merge[] = {
  87     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
  88     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
  89     {"stride",   {INT_VALUE,  1, {NULL}}, NULL, SPAR_OPTIONAL},
  90 };
  91
  92 /** Help text for the merging selection modifiers. */
  93 static const char *help_merge[] = {
  94     "MERGING SELECTIONS[PAR]",
  95
  96     "[TT]POSEXPR merge POSEXPR [stride INT][tt][BR]",
  97     "[TT]POSEXPR merge POSEXPR [merge POSEXPR ...][tt][BR]",
  98     "[TT]POSEXPR plus POSEXPR [plus POSEXPR ...][tt][PAR]",
  99
 100     "Basic selection keywords can only create selections where each atom",
 101     "occurs at most once. The [TT]merge[tt] and [TT]plus[tt] selection",
 102     "keywords can be used to work around this limitation. Both create",
 103     "a selection that contains the positions from all the given position",
 104     "expressions, even if they contain duplicates.",
 105     "The difference between the two is that [TT]merge[tt] expects two or more",
 106     "selections with the same number of positions, and the output contains",
 107     "the input positions selected from each expression in turn, i.e.,",
 108     "the output is like A1 B1 A2 B2 and so on. It is also possible to merge",
 109     "selections of unequal size as long as the size of the first is a",
 110     "multiple of the second one. The [TT]stride[tt] parameter can be used",
 111     "to explicitly provide this multiplicity.",
 112     "[TT]plus[tt] simply concatenates the positions after each other, and",
 113     "can work also with selections of different sizes.",
 114     "These keywords are valid only at the selection level, not in any",
 115     "subexpressions.",
 116 };
 117
 118 /** \internal Selection method data for the \p plus modifier. */
 119 gmx_ana_selmethod_t sm_merge = {
 120     "merge", POS_VALUE, SMETH_MODIFIER,
 121     asize(smparams_merge), smparams_merge,
 122     &init_data_merge,
 123     NULL,
 124     &init_merge,
 125     &init_output_merge,
 126     &free_data_merge,
 127     NULL,
 128     NULL,
 129     &evaluate_merge,
 130     {"merge POSEXPR", asize(help_merge), help_merge},
 131 };
 132
 133 /** \internal Selection method data for the \p plus modifier. */
 134 gmx_ana_selmethod_t sm_plus = {
 135     "plus", POS_VALUE, SMETH_MODIFIER,
 136     asize(smparams_merge)-1, smparams_merge,
 137     &init_data_merge,
 138     NULL,
 139     &init_merge,
 140     &init_output_plus,
 141     &free_data_merge,
 142     NULL,
 143     NULL,
 144     &evaluate_plus,
 145     {"plus POSEXPR", asize(help_merge), help_merge},
 146 };
 147
 148 /*!
 149  * \param[in]     npar  Should be 2 for \c plus and 3 for \c merge.
 150  * \param[in,out] param Method parameters (should point to a copy of
 151  *   \ref smparams_merge).
 152  * \returns Pointer to the allocated data (\p t_methoddata_merge).
 153  *
 154  * Allocates memory for a \p t_methoddata_merge structure.
 155  */
 156 static void *
 157 init_data_merge(int npar, gmx_ana_selparam_t *param)
 158 {
 159     t_methoddata_merge *data;
 160
 161     snew(data, 1);
 162     data->stride     = 0;
 163     param[0].val.u.p = &data->p1;
 164     param[1].val.u.p = &data->p2;
 165     if (npar > 2)
 166     {
 167         param[2].val.u.i = &data->stride;
 168     }
 169     return data;
 170 }
 171
 172 /*!
 173  * \param[in] top   Not used.
 174  * \param[in] npar  Not used (should be 2 or 3).
 175  * \param[in] param Method parameters (should point to \ref smparams_merge).
 176  * \param[in] data  Should point to a \p t_methoddata_merge.
 177  * \returns   0 if everything is successful, -1 on error.
 178  */
 179 static void
 180 init_merge(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 181 {
 182     t_methoddata_merge *d = (t_methoddata_merge *)data;
 183
 184     if (d->stride < 0)
 185     {
 186         GMX_THROW(gmx::InvalidInputError("Stride for merging should be positive"));
 187     }
 188     /* If no stride given, deduce it from the input sizes */
 189     if (d->stride == 0)
 190     {
 191         d->stride = d->p1.nr / d->p2.nr;
 192     }
 193     if (d->p1.nr != d->stride*d->p2.nr)
 194     {
 195         GMX_THROW(gmx::InconsistentInputError("The number of positions to be merged are not compatible"));
 196     }
 197     /* We access the m.b.nra field instead of g->isize in the position
 198      * data structures to handle cases where g is NULL
 199      * (this occurs with constant positions. */
 200     gmx_ana_index_reserve(&d->g, d->p1.m.b.nra + d->p2.m.b.nra);
 201     d->g.isize = d->p1.m.b.nra + d->p2.m.b.nra;
 202 }
 203
 204 /*! \brief
 205  * Does common initialization to all merging modifiers.
 206  *
 207  * \param[in]     top   Topology data structure.
 208  * \param[in,out] out   Pointer to output data structure.
 209  * \param[in,out] data  Should point to \c t_methoddata_merge.
 210  */
 211 static void
 212 init_output_common(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 213 {
 214     t_methoddata_merge *d = (t_methoddata_merge *)data;
 215
 216     if (d->p1.m.type != d->p2.m.type)
 217     {
 218         /* TODO: Maybe we could pick something else here? */
 219         out->u.p->m.type = INDEX_UNKNOWN;
 220     }
 221     else
 222     {
 223         out->u.p->m.type = d->p1.m.type;
 224     }
 225     gmx_ana_pos_reserve(out->u.p, d->p1.nr + d->p2.nr, d->g.isize);
 226     if (d->p1.v)
 227     {
 228         gmx_ana_pos_reserve_velocities(out->u.p);
 229     }
 230     if (d->p1.f)
 231     {
 232         gmx_ana_pos_reserve_forces(out->u.p);
 233     }
 234     gmx_ana_pos_set_evalgrp(out->u.p, &d->g);
 235     gmx_ana_pos_empty_init(out->u.p);
 236     d->g.isize = 0;
 237 }
 238
 239 /*!
 240  * \param[in]     top   Topology data structure.
 241  * \param[in,out] out   Pointer to output data structure.
 242  * \param[in,out] data  Should point to \c t_methoddata_merge.
 243  */
 244 static void
 245 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 246 {
 247     t_methoddata_merge *d = (t_methoddata_merge *)data;
 248     int                 i, j;
 249
 250     init_output_common(top, out, data);
 251     for (i = 0; i < d->p2.nr; ++i)
 252     {
 253         for (j = 0; j < d->stride; ++j)
 254         {
 255             gmx_ana_pos_append_init(out->u.p, &d->g, &d->p1, d->stride*i+j);
 256         }
 257         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p2, i);
 258     }
 259 }
 260
 261 /*!
 262  * \param[in]     top   Topology data structure.
 263  * \param[in,out] out   Pointer to output data structure.
 264  * \param[in,out] data  Should point to \c t_methoddata_merge.
 265  */
 266 static void
 267 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 268 {
 269     t_methoddata_merge *d = (t_methoddata_merge *)data;
 270     int                 i;
 271
 272     init_output_common(top, out, data);
 273     for (i = 0; i < d->p1.nr; ++i)
 274     {
 275         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p1, i);
 276     }
 277     for (i = 0; i < d->p2.nr; ++i)
 278     {
 279         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p2, i);
 280     }
 281 }
 282
 283 /*!
 284  * \param data Data to free (should point to a \p t_methoddata_merge).
 285  *
 286  * Frees the memory allocated for \c t_methoddata_merge.
 287  */
 288 static void
 289 free_data_merge(void *data)
 290 {
 291     t_methoddata_merge *d = (t_methoddata_merge *)data;
 292
 293     gmx_ana_index_deinit(&d->g);
 294     sfree(d);
 295 }
 296
 297 /*!
 298  * \param[in]  top   Not used.
 299  * \param[in]  fr    Not used.
 300  * \param[in]  pbc   Not used.
 301  * \param[in]  p     Positions to merge (should point to \p data->p1).
 302  * \param[out] out   Output data structure (\p out->u.p is used).
 303  * \param[in]  data  Should point to a \p t_methoddata_merge.
 304  */
 305 static void
 306 evaluate_merge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 307                gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data)
 308 {
 309     t_methoddata_merge *d = (t_methoddata_merge *)data;
 310     int                 i, j;
 311     int                 refid;
 312
 313     if (d->p1.nr != d->stride*d->p2.nr)
 314     {
 315         GMX_THROW(gmx::InconsistentInputError("The number of positions to be merged are not compatible"));
 316     }
 317     d->g.isize = 0;
 318     gmx_ana_pos_empty(out->u.p);
 319     for (i = 0; i < d->p2.nr; ++i)
 320     {
 321         for (j = 0; j < d->stride; ++j)
 322         {
 323             refid = d->p1.m.refid[d->stride*i+j];
 324             if (refid != -1)
 325             {
 326                 refid = (d->stride+1) * (refid / d->stride) + (refid % d->stride);
 327             }
 328             gmx_ana_pos_append(out->u.p, &d->g, &d->p1, d->stride*i+j, refid);
 329         }
 330         refid = (d->stride+1)*d->p2.m.refid[i]+d->stride;
 331         gmx_ana_pos_append(out->u.p, &d->g, &d->p2, i, refid);
 332     }
 333     gmx_ana_pos_append_finish(out->u.p);
 334 }
 335
 336 /*!
 337  * \param[in]  top   Not used.
 338  * \param[in]  fr    Not used.
 339  * \param[in]  pbc   Not used.
 340  * \param[in]  p     Positions to merge (should point to \p data->p1).
 341  * \param[out] out   Output data structure (\p out->u.p is used).
 342  * \param[in]  data  Should point to a \p t_methoddata_merge.
 343  */
 344 static void
 345 evaluate_plus(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 346               gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data)
 347 {
 348     t_methoddata_merge *d = (t_methoddata_merge *)data;
 349     int                 i;
 350     int                 refid;
 351
 352     d->g.isize = 0;
 353     gmx_ana_pos_empty(out->u.p);
 354     for (i = 0; i < d->p1.nr; ++i)
 355     {
 356         refid = d->p1.m.refid[i];
 357         gmx_ana_pos_append(out->u.p, &d->g, &d->p1, i, refid);
 358     }
 359     for (i = 0; i < d->p2.nr; ++i)
 360     {
 361         refid = d->p2.m.refid[i];
 362         if (refid != -1)
 363         {
 364             refid += d->p1.m.b.nr;
 365         }
 366         gmx_ana_pos_append(out->u.p, &d->g, &d->p2, i, refid);
 367     }
 368     gmx_ana_pos_append_finish(out->u.p);
 369 }