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, 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 /*! \brief Determines up to what local atom count a pull group gets processed single-threaded.
  62  *
  63  * We set this limit to 1 with debug to catch bugs.
  64  * On Haswell with GCC 5 the cross-over point is around 400 atoms,
  65  * independent of thread count and hyper-threading.
  66  */
  67 #ifdef NDEBUG
  68 static const int c_pullMaxNumLocalAtomsSingleThreaded = 100;
  69 #else
  70 static const int c_pullMaxNumLocalAtomsSingleThreaded = 1;
  71 #endif
  72
  73 class PullHistory;
  74 enum class PbcType : int;
  75
  76 enum
  77 {
  78     epgrppbcNONE,
  79     epgrppbcREFAT,
  80     epgrppbcCOS,
  81     epgrppbcPREVSTEPCOM
  82 };
  83
  84 /*! \internal
  85  * \brief Pull group data used during pulling
  86  */
  87 struct pull_group_work_t
  88 {
  89     /*! \brief Constructor
  90      *
  91      * \param[in] params                  The group parameters set by the user
  92      * \param[in] atomSet                 The global to local atom set manager
  93      * \param[in] setPbcRefToPrevStepCOM Does this pull group use the COM from the previous step as reference position?
  94      */
  95     pull_group_work_t(const t_pull_group& params, gmx::LocalAtomSet atomSet, bool setPbcRefToPrevStepCOM);
  96
  97     /* Data only modified at initialization */
  98     const t_pull_group params;   /**< The pull group parameters */
  99     const int          epgrppbc; /**< The type of pbc for this pull group, see enum above */
 100     bool               needToCalcCom; /**< Do we need to calculate the COM? (Not for group 0 or if only used as cylinder group) */
 101     std::vector<real>  globalWeights; /**< Weights per atom set by the user and/or mass/friction coefficients, if empty all weights are equal */
 102
 103     /* Data modified only at init or at domain decomposition */
 104     gmx::LocalAtomSet                  atomSet;      /**< Global to local atom set mapper */
 105     std::vector<real>                  localWeights; /**< Weights for the local atoms */
 106     std::unique_ptr<gmx::LocalAtomSet> pbcAtomSet;   /**< Keeps index of the pbc reference atom.
 107                                                           The stored LocalAtomSet consists of exactly   one atom when pbc reference atom is required.
 108                                                           When no pbc refence atom is used, this   pointer   shall be null. */
 109
 110     /* Data, potentially, changed at every pull call */
 111     real mwscale; /**< mass*weight scaling factor 1/sum w m */
 112     real wscale;  /**< scaling factor for the weights: sum w m/sum w w m */
 113     real invtm;   /**< inverse total mass of the group: 1/wscale sum w m */
 114     std::vector<gmx::BasicVector<double>> mdw; /**< mass*gradient(weight) for atoms */
 115     std::vector<double>                   dv;  /**< distance to the other group(s) along vec */
 116     dvec                                  x;   /**< COM before update */
 117     dvec                                  xp;  /**< COM after update before constraining */
 118     dvec                                  x_prev_step; /**< center of mass of the previous step */
 119 };
 120
 121 /* Struct describing the instantaneous spatial layout of a pull coordinate */
 122 struct PullCoordSpatialData
 123 {
 124     dvec   dr01;       /* The direction vector of group 1 relative to group 0 */
 125     dvec   dr23;       /* The direction vector of group 3 relative to group 2 */
 126     dvec   dr45;       /* The direction vector of group 5 relative to group 4 */
 127     dvec   vec;        /* The pull direction */
 128     double vec_len;    /* Length of vec for direction-relative */
 129     dvec   ffrad;      /* conversion factor from vec to radial force */
 130     double cyl_dev;    /* The deviation from the reference position */
 131     dvec   planevec_m; /* Normal of plane for groups 0, 1, 2, 3 for geometry dihedral */
 132     dvec   planevec_n; /* Normal of plane for groups 2, 3, 4, 5 for geometry dihedral */
 133
 134     double value; /* The current value of the coordinate, units of nm or rad */
 135 };
 136
 137 /* Struct with parameters and force evaluation local data for a pull coordinate */
 138 struct pull_coord_work_t
 139 {
 140     /* Constructor */
 141     pull_coord_work_t(const t_pull_coord& params) :
 142         params(params),
 143         value_ref(0),
 144         spatialData(),
 145         scalarForce(0),
 146         bExternalPotentialProviderHasBeenRegistered(false)
 147     {
 148     }
 149
 150     const t_pull_coord params; /* Pull coordinate parameters */
 151
 152     double value_ref; /* The reference value, usually init+rate*t, units of nm or rad */
 153
 154     PullCoordSpatialData spatialData; /* Data defining the current geometry */
 155
 156     double scalarForce; /* Scalar force for this cooordinate */
 157
 158     /* For external-potential coordinates only, for checking if a provider has been registered */
 159     bool bExternalPotentialProviderHasBeenRegistered;
 160 };
 161
 162 /* Struct for storing vectorial forces for a pull coordinate */
 163 struct PullCoordVectorForces
 164 {
 165     dvec force01; /* Force due to the pulling/constraining for groups 0, 1 */
 166     dvec force23; /* Force for groups 2 and 3 */
 167     dvec force45; /* Force for groups 4 and 5 */
 168 };
 169
 170 /* Struct for sums over (local) atoms in a pull group */
 171 struct ComSums
 172 {
 173     /* For normal weighting */
 174     double sum_wm;   /* Sum of weight*mass        */
 175     double sum_wwm;  /* Sum of weight*weight*mass */
 176     dvec   sum_wmx;  /* Sum of weight*mass*x      */
 177     dvec   sum_wmxp; /* Sum of weight*mass*xp     */
 178
 179     /* For cosine weighting */
 180     double sum_cm;  /* Sum of cos(x)*mass          */
 181     double sum_sm;  /* Sum of sin(x)*mass          */
 182     double sum_ccm; /* Sum of cos(x)*cos(x)*mass   */
 183     double sum_csm; /* Sum of cos(x)*sin(x)*mass   */
 184     double sum_ssm; /* Sum of sin(x)*sin(x)*mass   */
 185     double sum_cmp; /* Sum of cos(xp)*sin(xp)*mass */
 186     double sum_smp; /* Sum of sin(xp)*sin(xp)*mass */
 187
 188     /* Dummy data to ensure adjacent elements in an array are separated
 189      * by a cache line size, max 128 bytes.
 190      * TODO: Replace this by some automated mechanism.
 191      */
 192     int dummy[32];
 193 };
 194
 195 /*! \brief The normal COM buffer needs 3 elements per group */
 196 static constexpr int c_comBufferStride = 3;
 197
 198 /*! \brief The cylinder buffer needs 9 elements per group */
 199 static constexpr int c_cylinderBufferStride = 9;
 200
 201 struct pull_comm_t
 202 {
 203     gmx_bool bParticipateAll; /* Do all ranks always participate in pulling? */
 204     gmx_bool bParticipate;    /* Does our rank participate in pulling? */
 205 #if GMX_MPI
 206     MPI_Comm mpi_comm_com; /* Communicator for pulling */
 207 #endif
 208     int  nparticipate; /* The number of ranks participating */
 209     bool isMasterRank; /* Tells whether our rank is the master rank and thus should add the pull virial */
 210
 211     int64_t setup_count; /* The number of decomposition calls */
 212     int64_t must_count;  /* The last count our rank needed to be part */
 213
 214     /* Buffers for parallel reductions */
 215     std::vector<gmx::RVec>                pbcAtomBuffer; /* COM calculation buffer */
 216     std::vector<gmx::BasicVector<double>> comBuffer;     /* COM calculation buffer */
 217     std::vector<double> cylinderBuffer; /* cylinder ref. groups calculation buffer */
 218 };
 219
 220 // The COM pull force calculation data structure
 221 // TODO Convert this into a ForceProvider
 222 struct pull_t
 223 {
 224     /* Global parameters */
 225     pull_params_t params; /* The pull parameters, from inputrec */
 226
 227     gmx_bool bPotential;  /* Are there coordinates with potential? */
 228     gmx_bool bConstraint; /* Are there constrained coordinates? */
 229     gmx_bool bAngle;      /* Are there angle geometry coordinates? */
 230
 231     PbcType  pbcType;   /* the boundary conditions */
 232     int      npbcdim;   /* do pbc in dims 0 <= dim < npbcdim */
 233     gmx_bool bRefAt;    /* do we need reference atoms for a group COM ? */
 234     int      cosdim;    /* dimension for cosine weighting, -1 if none */
 235     gmx_bool bCylinder; /* Is group 0 a cylinder group? */
 236
 237     /* Parameters + dynamic data for groups */
 238     std::vector<pull_group_work_t> group; /* The pull group param and work data */
 239     std::vector<pull_group_work_t> dyna;  /* Dynamic groups for geom=cylinder */
 240
 241     /* Parameters + dynamic data for coordinates */
 242     std::vector<pull_coord_work_t> coord; /* The pull group param and work data */
 243
 244     /* Global dynamic data */
 245     gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
 246
 247     int                  nthreads; /* Number of threads used by the pull code */
 248     std::vector<ComSums> comSums;  /* Work array for summing for COM, 1 entry per thread */
 249
 250     pull_comm_t comm; /* Communication parameters, communicator and buffers */
 251
 252     FILE* out_x; /* Output file for pull data */
 253     FILE* out_f; /* Output file for pull data */
 254
 255     bool bXOutAverage; /* Output average pull coordinates */
 256     bool bFOutAverage; /* Output average pull forces */
 257
 258     PullHistory* coordForceHistory; /* Pull coordinate and force history */
 259
 260     /* The number of coordinates using an external potential */
 261     int numCoordinatesWithExternalPotential;
 262     /* Counter for checking external potential registration */
 263     int numUnregisteredExternalPotentials;
 264     /* */
 265     int numExternalPotentialsStillToBeAppliedThisStep;
 266 };
 267
 268 /*! \brief Copies the pull group COM of the previous step from the checkpoint state to the pull state
 269  *
 270  * \param[in]   pull  The COM pull force calculation data structure
 271  * \param[in]   state The global state container
 272  */
 273 void setPrevStepPullComFromState(struct pull_t* pull, const t_state* state);
 274
 275 /*! \brief Resizes the vector, in the state container, containing the COMs from the previous step
 276  *
 277  * \param[in]   state The global state container
 278  * \param[in]   pull  The COM pull force calculation data structure
 279  */
 280 void allocStatePrevStepPullCom(t_state* state, const pull_t* pull);
 281
 282
 283 #endif