src/gromacs/selection/sm_merge.cpp

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