src/gromacs/selection/sm_permute.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 permute selection modifier.
  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/exceptions.h"
  49 #include "gromacs/utility/smalloc.h"
  50 #include "gromacs/utility/stringutil.h"
  51
  52 #include "position.h"
  53 #include "selmethod.h"
  54 #include "selmethod_impl.h"
  55
  56 /*! \internal \brief
  57  * Data structure for the \p permute selection modifier.
  58  */
  59 typedef struct methoddata_permute
  60 {
  61     /** Positions to permute. */
  62     gmx_ana_pos_t p;
  63     /** Number of elements in the permutation. */
  64     int n;
  65     /** Array describing the permutation. */
  66     int* perm;
  67     /** Array that has the permutation reversed. */
  68     int* rperm;
  69 } t_methoddata_permute;
  70
  71 /*! \brief
  72  * Allocates data for the \p permute selection modifier.
  73  *
  74  * \param[in]     npar  Not used (should be 2).
  75  * \param[in,out] param Method parameters (should point to a copy of
  76  *   \ref smparams_permute).
  77  * \returns Pointer to the allocated data (\p t_methoddata_permute).
  78  *
  79  * Allocates memory for a \p t_methoddata_permute structure.
  80  */
  81 static void* init_data_permute(int npar, gmx_ana_selparam_t* param);
  82 /*! \brief
  83  * Initializes data for the \p permute selection modifier.
  84  *
  85  * \param[in] top   Not used.
  86  * \param[in] npar  Not used (should be 2).
  87  * \param[in] param Method parameters (should point to \ref smparams_permute).
  88  * \param[in] data  Should point to a \p t_methoddata_permute.
  89  * \returns   0 if the input permutation is valid, -1 on error.
  90  */
  91 static void init_permute(const gmx_mtop_t* top, int npar, gmx_ana_selparam_t* param, void* data);
  92 /*! \brief
  93  * Initializes output for the \p permute selection modifier.
  94  *
  95  * \param[in]     top   Topology data structure.
  96  * \param[in,out] out   Pointer to output data structure.
  97  * \param[in,out] data  Should point to \c t_methoddata_permute.
  98  */
  99 static void init_output_permute(const gmx_mtop_t* top, gmx_ana_selvalue_t* out, void* data);
 100 /** Frees the memory allocated for the \p permute selection modifier. */
 101 static void free_data_permute(void* data);
 102 /*! \brief
 103  * Evaluates the \p permute selection modifier.
 104  *
 105  * \param[in]  context Not used.
 106  * \param[in]  p     Positions to permute (should point to \p data->p).
 107  * \param[out] out   Output data structure (\p out->u.p is used).
 108  * \param[in]  data  Should point to a \p t_methoddata_permute.
 109  *
 110  * Throws if the size of \p p is not divisible by the number of
 111  * elements in the permutation.
 112  */
 113 static void evaluate_permute(const gmx::SelMethodEvalContext& context,
 114                              gmx_ana_pos_t*                   p,
 115                              gmx_ana_selvalue_t*              out,
 116                              void*                            data);
 117
 118 /** Parameters for the \p permute selection modifier. */
 119 static gmx_ana_selparam_t smparams_permute[] = {
 120     { nullptr, { POS_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
 121     { nullptr, { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_VARNUM },
 122 };
 123
 124 /** Help text for the \p permute selection modifier. */
 125 static const char* const help_permute[] = {
 126     "::",
 127     "",
 128     "  permute P1 ... PN",
 129     "",
 130     "By default, all selections are evaluated such that the atom indices are",
 131     "returned in ascending order. This can be changed by appending",
 132     "[TT]permute P1 P2 ... PN[tt] to an expression.",
 133     "The [TT]Pi[tt] should form a permutation of the numbers 1 to N.",
 134     "This keyword permutes each N-position block in the selection such that",
 135     "the i'th position in the block becomes Pi'th.",
 136     "Note that it is the positions that are permuted, not individual atoms.",
 137     "A fatal error occurs if the size of the selection is not a multiple of n.",
 138     "It is only possible to permute the whole selection expression, not any",
 139     "subexpressions, i.e., the [TT]permute[tt] keyword should appear last in",
 140     "a selection.",
 141 };
 142
 143 /** Selection method data for the \p permute modifier. */
 144 gmx_ana_selmethod_t sm_permute = {
 145     "permute",
 146     POS_VALUE,
 147     SMETH_MODIFIER,
 148     asize(smparams_permute),
 149     smparams_permute,
 150     &init_data_permute,
 151     nullptr,
 152     &init_permute,
 153     &init_output_permute,
 154     &free_data_permute,
 155     nullptr,
 156     nullptr,
 157     &evaluate_permute,
 158     { "POSEXPR permute P1 ... PN", "Permuting selections", asize(help_permute), help_permute },
 159 };
 160
 161 static void* init_data_permute(int /* npar */, gmx_ana_selparam_t* param)
 162 {
 163     t_methoddata_permute* data = new t_methoddata_permute();
 164     data->n                    = 0;
 165     data->perm                 = nullptr;
 166     data->rperm                = nullptr;
 167     param[0].val.u.p           = &data->p;
 168     return data;
 169 }
 170
 171 static void init_permute(const gmx_mtop_t* /* top */, int /* npar */, gmx_ana_selparam_t* param, void* data)
 172 {
 173     t_methoddata_permute* d = static_cast<t_methoddata_permute*>(data);
 174     int                   i;
 175
 176     d->n    = param[1].val.nr;
 177     d->perm = param[1].val.u.i;
 178     if (d->p.count() % d->n != 0)
 179     {
 180         GMX_THROW(gmx::InconsistentInputError(gmx::formatString(
 181                 "The number of positions to be permuted is not divisible by %d", d->n)));
 182     }
 183     snew(d->rperm, d->n);
 184     for (i = 0; i < d->n; ++i)
 185     {
 186         d->rperm[i] = -1;
 187     }
 188     for (i = 0; i < d->n; ++i)
 189     {
 190         d->perm[i]--;
 191         if (d->perm[i] < 0 || d->perm[i] >= d->n)
 192         {
 193             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 194         }
 195         if (d->rperm[d->perm[i]] >= 0)
 196         {
 197             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 198         }
 199         d->rperm[d->perm[i]] = i;
 200     }
 201 }
 202
 203 static void init_output_permute(const gmx_mtop_t* /* top */, gmx_ana_selvalue_t* out, void* data)
 204 {
 205     t_methoddata_permute* d = static_cast<t_methoddata_permute*>(data);
 206     int                   i, j, b;
 207
 208     out->u.p->m.type = d->p.m.type;
 209     gmx_ana_pos_reserve_for_append(out->u.p, d->p.count(), d->p.m.b.nra, d->p.v != nullptr,
 210                                    d->p.f != nullptr);
 211     gmx_ana_pos_empty_init(out->u.p);
 212     for (i = 0; i < d->p.count(); i += d->n)
 213     {
 214         for (j = 0; j < d->n; ++j)
 215         {
 216             b = i + d->rperm[j];
 217             gmx_ana_pos_append_init(out->u.p, &d->p, b);
 218         }
 219     }
 220 }
 221
 222 /*!
 223  * \param data Data to free (should point to a \p t_methoddata_permute).
 224  *
 225  * Frees the memory allocated for \c t_methoddata_permute.
 226  */
 227 static void free_data_permute(void* data)
 228 {
 229     t_methoddata_permute* d = static_cast<t_methoddata_permute*>(data);
 230
 231     sfree(d->rperm);
 232     delete d;
 233 }
 234
 235 static void evaluate_permute(const gmx::SelMethodEvalContext& /*context*/,
 236                              gmx_ana_pos_t* /*p*/,
 237                              gmx_ana_selvalue_t* out,
 238                              void*               data)
 239 {
 240     t_methoddata_permute* d = static_cast<t_methoddata_permute*>(data);
 241     int                   i, j, b;
 242     int                   refid;
 243
 244     if (d->p.count() % d->n != 0)
 245     {
 246         GMX_THROW(gmx::InconsistentInputError(gmx::formatString(
 247                 "The number of positions to be permuted is not divisible by %d", d->n)));
 248     }
 249     gmx_ana_pos_empty(out->u.p);
 250     for (i = 0; i < d->p.count(); i += d->n)
 251     {
 252         for (j = 0; j < d->n; ++j)
 253         {
 254             b     = i + d->rperm[j];
 255             refid = d->p.m.refid[b];
 256             if (refid != -1)
 257             {
 258                 /* De-permute the reference ID */
 259                 refid = refid - (refid % d->n) + d->perm[refid % d->n];
 260             }
 261             gmx_ana_pos_append(out->u.p, &d->p, b, refid);
 262         }
 263     }
 264     gmx_ana_pos_append_finish(out->u.p);
 265 }