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