src/gromacs/pulling/pull_internal.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2018 by the GROMACS development team.
   7  * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38
  39 /*! \internal \file
  40  *
  41  *
  42  * \brief
  43  * This file contains datatypes and function declarations for internal
  44    use in the pull code.
  45  *
  46  * \author Berk Hess
  47  */
  48
  49 #ifndef GMX_PULLING_PULL_INTERNAL_H
  50 #define GMX_PULLING_PULL_INTERNAL_H
  51
  52 #include "config.h"
  53
  54 #include <memory>
  55 #include <vector>
  56
  57 #include "gromacs/domdec/localatomset.h"
  58 #include "gromacs/mdtypes/pull_params.h"
  59 #include "gromacs/utility/gmxmpi.h"
  60
  61 #include "pullcoordexpressionparser.h"
  62
  63 /*! \brief Determines up to what local atom count a pull group gets processed single-threaded.
  64  *
  65  * We set this limit to 1 with debug to catch bugs.
  66  * On Haswell with GCC 5 the cross-over point is around 400 atoms,
  67  * independent of thread count and hyper-threading.
  68  */
  69 #ifdef NDEBUG
  70 static const int c_pullMaxNumLocalAtomsSingleThreaded = 100;
  71 #else
  72 static const int c_pullMaxNumLocalAtomsSingleThreaded = 1;
  73 #endif
  74
  75 class PullHistory;
  76 enum class PbcType : int;
  77
  78 class t_state;
  79
  80 enum
  81 {
  82     epgrppbcNONE,
  83     epgrppbcREFAT,
  84     epgrppbcCOS,
  85     epgrppbcPREVSTEPCOM
  86 };
  87
  88 /*! \internal
  89  * \brief Pull group data used during pulling
  90  */
  91 struct pull_group_work_t
  92 {
  93     /*! \brief Constructor
  94      *
  95      * \param[in] params                  The group parameters set by the user
  96      * \param[in] atomSet                 The global to local atom set manager
  97      * \param[in] setPbcRefToPrevStepCOM Does this pull group use the COM from the previous step as reference position?
  98      * \param[in] maxNumThreads           Use either this number of threads of 1 for operations on x and f
  99      */
 100     pull_group_work_t(const t_pull_group& params,
 101                       gmx::LocalAtomSet   atomSet,
 102                       bool                setPbcRefToPrevStepCOM,
 103                       int                 maxNumThreads);
 104
 105     //! Returns the number of threads to use for local atom operations based on the local atom count
 106     int numThreads() const
 107     {
 108         return atomSet.numAtomsLocal() <= c_pullMaxNumLocalAtomsSingleThreaded ? 1 : maxNumThreads;
 109     }
 110
 111     /* Data only modified at initialization */
 112     const t_pull_group params;   /**< The pull group parameters */
 113     const int          epgrppbc; /**< The type of pbc for this pull group, see enum above */
 114     const int maxNumThreads; /**< The maximum number of threads to use for operations on x and f */
 115     bool      needToCalcCom; /**< Do we need to calculate the COM? (Not for group 0 or if only used as cylinder group) */
 116     std::vector<real> globalWeights; /**< Weights per atom set by the user and/or mass/friction coefficients, if empty all weights are equal */
 117
 118     /* Data modified only at init or at domain decomposition */
 119     gmx::LocalAtomSet                  atomSet;      /**< Global to local atom set mapper */
 120     std::vector<real>                  localWeights; /**< Weights for the local atoms */
 121     std::unique_ptr<gmx::LocalAtomSet> pbcAtomSet;   /**< Keeps index of the pbc reference atom.
 122                                                           The stored LocalAtomSet consists of exactly   one atom when pbc reference atom is required.
 123                                                           When no pbc refence atom is used, this   pointer   shall be null. */
 124
 125     /* Data, potentially, changed at every pull call */
 126     real mwscale; /**< mass*weight scaling factor 1/sum w m */
 127     real wscale;  /**< scaling factor for the weights: sum w m/sum w w m */
 128     real invtm;   /**< inverse total mass of the group: 1/wscale sum w m */
 129     std::vector<gmx::BasicVector<double>> mdw; /**< mass*gradient(weight) for atoms */
 130     std::vector<double>                   dv;  /**< distance to the other group(s) along vec */
 131     dvec                                  x;   /**< COM before update */
 132     dvec                                  xp;  /**< COM after update before constraining */
 133     dvec                                  x_prev_step; /**< center of mass of the previous step */
 134 };
 135
 136 /* Struct describing the instantaneous spatial layout of a pull coordinate */
 137 struct PullCoordSpatialData
 138 {
 139     dvec   dr01;       /* The direction vector of group 1 relative to group 0 */
 140     dvec   dr23;       /* The direction vector of group 3 relative to group 2 */
 141     dvec   dr45;       /* The direction vector of group 5 relative to group 4 */
 142     dvec   vec;        /* The pull direction */
 143     double vec_len;    /* Length of vec for direction-relative */
 144     dvec   ffrad;      /* conversion factor from vec to radial force */
 145     double cyl_dev;    /* The deviation from the reference position */
 146     dvec   planevec_m; /* Normal of plane for groups 0, 1, 2, 3 for geometry dihedral */
 147     dvec   planevec_n; /* Normal of plane for groups 2, 3, 4, 5 for geometry dihedral */
 148
 149     double value; /* The current value of the coordinate, units of nm or rad */
 150 };
 151
 152 //! \brief Struct with parameters and force evaluation local data for a pull coordinate
 153 struct pull_coord_work_t
 154 {
 155     //! Constructor
 156     pull_coord_work_t(const t_pull_coord& params) :
 157         params(params),
 158         value_ref(0),
 159         spatialData(),
 160         scalarForce(0),
 161         bExternalPotentialProviderHasBeenRegistered(false),
 162         expressionParser(params.eGeom == PullGroupGeometry::Transformation ? params.expression : "",
 163                          params.coordIndex),
 164         transformationVariables(params.eGeom == PullGroupGeometry::Transformation ? params.coordIndex : 0)
 165     {
 166     }
 167
 168     //! Pull coordinate parameters
 169     const t_pull_coord params;
 170
 171     //! Dynamic pull group 0 for this coordinate with dynamic weights, only present when needed */
 172     std::unique_ptr<pull_group_work_t> dynamicGroup0;
 173     //! The reference value, usually init+rate*t, units of nm or rad.
 174     double value_ref;
 175
 176     //! Data defining the current geometry
 177     PullCoordSpatialData spatialData;
 178
 179     //! Scalar force for this cooordinate
 180     double scalarForce;
 181
 182     //! For external-potential coordinates only, for checking if a provider has been registered
 183     bool bExternalPotentialProviderHasBeenRegistered;
 184
 185     //! The expression parser for a transformation coordinate
 186     gmx::PullCoordExpressionParser expressionParser;
 187     //! Variables from other pull coordinates for a transformation coordinate
 188     std::vector<double> transformationVariables;
 189 };
 190
 191 /* Struct for storing vectorial forces for a pull coordinate */
 192 struct PullCoordVectorForces
 193 {
 194     dvec force01; /* Force due to the pulling/constraining for groups 0, 1 */
 195     dvec force23; /* Force for groups 2 and 3 */
 196     dvec force45; /* Force for groups 4 and 5 */
 197 };
 198
 199 /* Struct for sums over (local) atoms in a pull group */
 200 struct ComSums
 201 {
 202     /* For normal weighting */
 203     double sum_wm;   /* Sum of weight*mass        */
 204     double sum_wwm;  /* Sum of weight*weight*mass */
 205     dvec   sum_wmx;  /* Sum of weight*mass*x      */
 206     dvec   sum_wmxp; /* Sum of weight*mass*xp     */
 207
 208     /* For cosine weighting */
 209     double sum_cm;  /* Sum of cos(x)*mass          */
 210     double sum_sm;  /* Sum of sin(x)*mass          */
 211     double sum_ccm; /* Sum of cos(x)*cos(x)*mass   */
 212     double sum_csm; /* Sum of cos(x)*sin(x)*mass   */
 213     double sum_ssm; /* Sum of sin(x)*sin(x)*mass   */
 214     double sum_cmp; /* Sum of cos(xp)*sin(xp)*mass */
 215     double sum_smp; /* Sum of sin(xp)*sin(xp)*mass */
 216
 217     /* Dummy data to ensure adjacent elements in an array are separated
 218      * by a cache line size, max 128 bytes.
 219      * TODO: Replace this by some automated mechanism.
 220      */
 221     int dummy[32];
 222 };
 223
 224 /*! \brief The normal COM buffer needs 3 elements per group */
 225 static constexpr int c_comBufferStride = 3;
 226
 227 /*! \brief The cylinder buffer needs 9 elements per group */
 228 static constexpr int c_cylinderBufferStride = 9;
 229
 230 struct pull_comm_t
 231 {
 232     gmx_bool bParticipateAll; /* Do all ranks always participate in pulling? */
 233     gmx_bool bParticipate;    /* Does our rank participate in pulling? */
 234 #if GMX_MPI
 235     MPI_Comm mpi_comm_com; /* Communicator for pulling */
 236 #endif
 237     int nparticipate; /* The number of ranks participating */
 238     bool isMasterRank; /* Tells whether our rank is the master rank and thus should add the pull virial */
 239
 240     int64_t setup_count; /* The number of decomposition calls */
 241     int64_t must_count;  /* The last count our rank needed to be part */
 242
 243     /* Buffers for parallel reductions */
 244     std::vector<gmx::RVec>                pbcAtomBuffer; /* COM calculation buffer */
 245     std::vector<gmx::BasicVector<double>> comBuffer;     /* COM calculation buffer */
 246     std::vector<double> cylinderBuffer; /* cylinder ref. groups calculation buffer */
 247 };
 248
 249 // The COM pull force calculation data structure
 250 // TODO Convert this into a ForceProvider
 251 struct pull_t
 252 {
 253     /* Global parameters */
 254     pull_params_t params; /* The pull parameters, from inputrec */
 255
 256     gmx_bool bPotential;  /* Are there coordinates with potential? */
 257     gmx_bool bConstraint; /* Are there constrained coordinates? */
 258     gmx_bool bAngle;      /* Are there angle geometry coordinates? */
 259
 260     PbcType  pbcType;   /* the boundary conditions */
 261     int      npbcdim;   /* do pbc in dims 0 <= dim < npbcdim */
 262     gmx_bool bRefAt;    /* do we need reference atoms for a group COM ? */
 263     int      cosdim;    /* dimension for cosine weighting, -1 if none */
 264     gmx_bool bCylinder; /* Is group 0 a cylinder group? */
 265
 266     /* Parameters + dynamic data for groups */
 267     std::vector<pull_group_work_t> group; /* The pull group param and work data */
 268
 269     /* Parameters + dynamic data for coordinates */
 270     std::vector<pull_coord_work_t> coord; /* The pull group param and work data */
 271
 272     /* Global dynamic data */
 273     gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
 274
 275     std::vector<ComSums> comSums; /* Work array for summing for COM, 1 entry per thread */
 276
 277     pull_comm_t comm; /* Communication parameters, communicator and buffers */
 278
 279     FILE* out_x; /* Output file for pull data */
 280     FILE* out_f; /* Output file for pull data */
 281
 282     bool bXOutAverage; /* Output average pull coordinates */
 283     bool bFOutAverage; /* Output average pull forces */
 284
 285     PullHistory* coordForceHistory; /* Pull coordinate and force history */
 286
 287     /* The number of coordinates using an external potential */
 288     int numCoordinatesWithExternalPotential;
 289     /* Counter for checking external potential registration */
 290     int numUnregisteredExternalPotentials;
 291     /* */
 292     int numExternalPotentialsStillToBeAppliedThisStep;
 293 };
 294
 295 /*! \brief Copies the pull group COM of the previous step from the checkpoint state to the pull state
 296  *
 297  * \param[in]   pull  The COM pull force calculation data structure
 298  * \param[in]   state The global state container
 299  */
 300 void setPrevStepPullComFromState(struct pull_t* pull, const t_state* state);
 301
 302 /*! \brief Resizes the vector, in the state container, containing the COMs from the previous step
 303  *
 304  * \param[in]   state The global state container
 305  * \param[in]   pull  The COM pull force calculation data structure
 306  */
 307 void allocStatePrevStepPullCom(t_state* state, const pull_t* pull);
 308
 309
 310 #endif