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,2013 by the GROMACS development team.
   5  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief
  39  * Implements the \p merge and \p plus selection modifiers.
  40  *
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  * \ingroup module_selection
  43  */
  44 #include "gmxpre.h"
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/utility/arraysize.h"
  48 #include "gromacs/utility/basedefinitions.h"
  49 #include "gromacs/utility/exceptions.h"
  50 #include "gromacs/utility/smalloc.h"
  51
  52 #include "position.h"
  53 #include "selmethod.h"
  54 #include "selmethod_impl.h"
  55
  56 /*! \internal \brief
  57  * Data structure for the merging selection modifiers.
  58  */
  59 typedef struct methoddata_merge
  60 {
  61     /** Input positions. */
  62     gmx_ana_pos_t p1;
  63     /** Other input positions. */
  64     gmx_ana_pos_t p2;
  65     /** Stride for merging (\c stride values from \c p1 for each in \c p2). */
  66     int stride;
  67 } t_methoddata_merge;
  68
  69 /** Allocates data for the merging selection modifiers. */
  70 static void* init_data_merge(int npar, gmx_ana_selparam_t* param);
  71 /*! \brief
  72  * Initializes data for the merging selection modifiers.
  73  *
  74  * \param[in] top   Not used.
  75  * \param[in] npar  Not used (should be 2 or 3).
  76  * \param[in] param Method parameters (should point to \ref smparams_merge).
  77  * \param[in] data  Should point to a \p t_methoddata_merge.
  78  * \returns   0 if everything is successful, -1 on error.
  79  */
  80 static void init_merge(const gmx_mtop_t* top, int npar, gmx_ana_selparam_t* param, void* data);
  81 /** Initializes output for the \p merge selection modifier. */
  82 static void init_output_merge(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data);
  83 /** Initializes output for the \p plus selection modifier. */
  84 static void init_output_plus(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data);
  85 /** Frees the memory allocated for the merging selection modifiers. */
  86 static void free_data_merge(void* data);
  87 /*! \brief
  88  * Evaluates the \p merge selection modifier.
  89  *
  90  * \param[in]  context Not used.
  91  * \param[in]  p     Positions to merge (should point to \p data->p1).
  92  * \param[out] out   Output data structure (\p out->u.p is used).
  93  * \param[in]  data  Should point to a \p t_methoddata_merge.
  94  */
  95 static void evaluate_merge(const gmx::SelMethodEvalContext& context,
  96                            gmx_ana_pos_t*                   p,
  97                            gmx_ana_selvalue_t*              out,
  98                            void*                            data);
  99 /*! \brief
 100  * Evaluates the \p plus selection modifier.
 101  *
 102  * \param[in]  context Not used.
 103  * \param[in]  p     Positions to merge (should point to \p data->p1).
 104  * \param[out] out   Output data structure (\p out->u.p is used).
 105  * \param[in]  data  Should point to a \p t_methoddata_merge.
 106  */
 107 static void evaluate_plus(const gmx::SelMethodEvalContext& context,
 108                           gmx_ana_pos_t*                   p,
 109                           gmx_ana_selvalue_t*              out,
 110                           void*                            data);
 111
 112 /** Parameters for the merging selection modifiers. */
 113 static gmx_ana_selparam_t smparams_merge[] = {
 114     { nullptr, { POS_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
 115     { nullptr, { POS_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
 116     { "stride", { INT_VALUE, 1, { nullptr } }, nullptr, SPAR_OPTIONAL },
 117 };
 118
 119 //! Help title for the merging selection modifiers.
 120 static const char helptitle_merge[] = "Merging selections";
 121 //! Help text for the merging selection modifiers.
 122 static const char* const help_merge[] = {
 123     "::",
 124     "",
 125     "  POSEXPR merge POSEXPR [stride INT]",
 126     "  POSEXPR merge POSEXPR [merge POSEXPR ...]",
 127     "  POSEXPR plus POSEXPR [plus POSEXPR ...]",
 128     "",
 129     "Basic selection keywords can only create selections where each atom",
 130     "occurs at most once. The [TT]merge[tt] and [TT]plus[tt] selection",
 131     "keywords can be used to work around this limitation. Both create",
 132     "a selection that contains the positions from all the given position",
 133     "expressions, even if they contain duplicates.",
 134     "The difference between the two is that [TT]merge[tt] expects two or more",
 135     "selections with the same number of positions, and the output contains",
 136     "the input positions selected from each expression in turn, i.e.,",
 137     "the output is like A1 B1 A2 B2 and so on. It is also possible to merge",
 138     "selections of unequal size as long as the size of the first is a",
 139     "multiple of the second one. The [TT]stride[tt] parameter can be used",
 140     "to explicitly provide this multiplicity.",
 141     "[TT]plus[tt] simply concatenates the positions after each other, and",
 142     "can work also with selections of different sizes.",
 143     "These keywords are valid only at the selection level, not in any",
 144     "subexpressions.",
 145 };
 146
 147 /** Selection method data for the \p plus modifier. */
 148 gmx_ana_selmethod_t sm_merge = {
 149     "merge",
 150     POS_VALUE,
 151     SMETH_MODIFIER,
 152     asize(smparams_merge),
 153     smparams_merge,
 154     &init_data_merge,
 155     nullptr,
 156     &init_merge,
 157     &init_output_merge,
 158     &free_data_merge,
 159     nullptr,
 160     nullptr,
 161     &evaluate_merge,
 162     { "merge POSEXPR", helptitle_merge, asize(help_merge), help_merge },
 163 };
 164
 165 /** Selection method data for the \p plus modifier. */
 166 gmx_ana_selmethod_t sm_plus = {
 167     "plus",
 168     POS_VALUE,
 169     SMETH_MODIFIER,
 170     asize(smparams_merge) - 1,
 171     smparams_merge,
 172     &init_data_merge,
 173     nullptr,
 174     &init_merge,
 175     &init_output_plus,
 176     &free_data_merge,
 177     nullptr,
 178     nullptr,
 179     &evaluate_plus,
 180     { "plus POSEXPR", helptitle_merge, asize(help_merge), help_merge },
 181 };
 182
 183 /*!
 184  * \param[in]     npar  Should be 2 for \c plus and 3 for \c merge.
 185  * \param[in,out] param Method parameters (should point to a copy of
 186  *   \ref smparams_merge).
 187  * \returns Pointer to the allocated data (\p t_methoddata_merge).
 188  *
 189  * Allocates memory for a \p t_methoddata_merge structure.
 190  */
 191 static void* init_data_merge(int npar, gmx_ana_selparam_t* param)
 192 {
 193     t_methoddata_merge* data = new t_methoddata_merge();
 194     data->stride             = 0;
 195     param[0].val.u.p         = &data->p1;
 196     param[1].val.u.p         = &data->p2;
 197     if (npar > 2)
 198     {
 199         param[2].val.u.i = &data->stride;
 200     }
 201     return data;
 202 }
 203
 204 static void init_merge(const gmx_mtop_t* /* top */, int /* npar */, gmx_ana_selparam_t* /* param */, void* data)
 205 {
 206     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 207
 208     if (d->stride < 0)
 209     {
 210         GMX_THROW(gmx::InvalidInputError("Stride for merging should be positive"));
 211     }
 212     /* If no stride given, deduce it from the input sizes */
 213     if (d->stride == 0)
 214     {
 215         d->stride = d->p1.count() / d->p2.count();
 216     }
 217     if (d->p1.count() != d->stride * d->p2.count())
 218     {
 219         GMX_THROW(gmx::InconsistentInputError(
 220                 "The number of positions to be merged are not compatible"));
 221     }
 222 }
 223
 224 /*! \brief
 225  * Does common initialization to all merging modifiers.
 226  *
 227  * \param[in]     top   Topology data structure.
 228  * \param[in,out] out   Pointer to output data structure.
 229  * \param[in,out] data  Should point to \c t_methoddata_merge.
 230  */
 231 static void init_output_common(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data)
 232 {
 233     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 234
 235     GMX_UNUSED_VALUE(top);
 236     if (d->p1.m.type != d->p2.m.type)
 237     {
 238         /* TODO: Maybe we could pick something else here? */
 239         out->u.p->m.type = INDEX_UNKNOWN;
 240     }
 241     else
 242     {
 243         out->u.p->m.type = d->p1.m.type;
 244     }
 245     gmx_ana_pos_reserve_for_append(out->u.p,
 246                                    d->p1.count() + d->p2.count(),
 247                                    d->p1.m.b.nra + d->p2.m.b.nra,
 248                                    d->p1.v != nullptr,
 249                                    d->p1.f != nullptr);
 250     gmx_ana_pos_empty_init(out->u.p);
 251 }
 252
 253 /*!
 254  * \param[in]     top   Topology data structure.
 255  * \param[in,out] out   Pointer to output data structure.
 256  * \param[in,out] data  Should point to \c t_methoddata_merge.
 257  */
 258 static void init_output_merge(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data)
 259 {
 260     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 261     int                 i, j;
 262
 263     init_output_common(top, out, data);
 264     for (i = 0; i < d->p2.count(); ++i)
 265     {
 266         for (j = 0; j < d->stride; ++j)
 267         {
 268             gmx_ana_pos_append_init(out->u.p, &d->p1, d->stride * i + j);
 269         }
 270         gmx_ana_pos_append_init(out->u.p, &d->p2, i);
 271     }
 272 }
 273
 274 /*!
 275  * \param[in]     top   Topology data structure.
 276  * \param[in,out] out   Pointer to output data structure.
 277  * \param[in,out] data  Should point to \c t_methoddata_merge.
 278  */
 279 static void init_output_plus(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data)
 280 {
 281     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 282     int                 i;
 283
 284     init_output_common(top, out, data);
 285     for (i = 0; i < d->p1.count(); ++i)
 286     {
 287         gmx_ana_pos_append_init(out->u.p, &d->p1, i);
 288     }
 289     for (i = 0; i < d->p2.count(); ++i)
 290     {
 291         gmx_ana_pos_append_init(out->u.p, &d->p2, i);
 292     }
 293 }
 294
 295 /*!
 296  * \param data Data to free (should point to a \p t_methoddata_merge).
 297  *
 298  * Frees the memory allocated for \c t_methoddata_merge.
 299  */
 300 static void free_data_merge(void* data)
 301 {
 302     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 303     delete d;
 304 }
 305
 306 static void evaluate_merge(const gmx::SelMethodEvalContext& /*context*/,
 307                            gmx_ana_pos_t* /* p */,
 308                            gmx_ana_selvalue_t* out,
 309                            void*               data)
 310 {
 311     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 312     int                 i, j;
 313     int                 refid;
 314
 315     if (d->p1.count() != d->stride * d->p2.count())
 316     {
 317         GMX_THROW(gmx::InconsistentInputError(
 318                 "The number of positions to be merged are not compatible"));
 319     }
 320     gmx_ana_pos_empty(out->u.p);
 321     for (i = 0; i < d->p2.count(); ++i)
 322     {
 323         for (j = 0; j < d->stride; ++j)
 324         {
 325             refid = d->p1.m.refid[d->stride * i + j];
 326             if (refid != -1)
 327             {
 328                 refid = (d->stride + 1) * (refid / d->stride) + (refid % d->stride);
 329             }
 330             gmx_ana_pos_append(out->u.p, &d->p1, d->stride * i + j, refid);
 331         }
 332         refid = (d->stride + 1) * d->p2.m.refid[i] + d->stride;
 333         gmx_ana_pos_append(out->u.p, &d->p2, i, refid);
 334     }
 335     gmx_ana_pos_append_finish(out->u.p);
 336 }
 337
 338 static void evaluate_plus(const gmx::SelMethodEvalContext& /*context*/,
 339                           gmx_ana_pos_t* /* p */,
 340                           gmx_ana_selvalue_t* out,
 341                           void*               data)
 342 {
 343     t_methoddata_merge* d = static_cast<t_methoddata_merge*>(data);
 344     int                 i;
 345     int                 refid;
 346
 347     gmx_ana_pos_empty(out->u.p);
 348     for (i = 0; i < d->p1.count(); ++i)
 349     {
 350         refid = d->p1.m.refid[i];
 351         gmx_ana_pos_append(out->u.p, &d->p1, i, refid);
 352     }
 353     for (i = 0; i < d->p2.count(); ++i)
 354     {
 355         refid = d->p2.m.refid[i];
 356         if (refid != -1)
 357         {
 358             refid += d->p1.m.b.nr;
 359         }
 360         gmx_ana_pos_append(out->u.p, &d->p2, i, refid);
 361     }
 362     gmx_ana_pos_append_finish(out->u.p);
 363 }