src/gromacs/mdlib/tests/constr_impl.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2018,2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief SHAKE and LINCS tests header.
  37  *
  38  * Contains description and constructor for the test data accumulating object,
  39  * declares CPU- and GPU-based functions used to apply SHAKE or LINCS on the
  40  * test data.
  41  *
  42  * \author Artem Zhmurov <zhmurov@gmail.com>
  43  * \ingroup module_mdlib
  44  */
  45
  46 #ifndef GMX_MDLIB_TESTS_CONSTR_IMPL_H
  47 #define GMX_MDLIB_TESTS_CONSTR_IMPL_H
  48
  49 #include "gmxpre.h"
  50
  51 #include <assert.h>
  52
  53 #include <cmath>
  54
  55 #include <algorithm>
  56 #include <unordered_map>
  57 #include <vector>
  58
  59 #include "gromacs/fileio/gmxfio.h"
  60 #include "gromacs/gmxlib/nrnb.h"
  61 #include "gromacs/gmxlib/nonbonded/nonbonded.h"
  62 #include "gromacs/gpu_utils/gpu_utils.h"
  63 #include "gromacs/math/paddedvector.h"
  64 #include "gromacs/math/vec.h"
  65 #include "gromacs/math/vectypes.h"
  66 #include "gromacs/mdlib/constr.h"
  67 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  68 #include "gromacs/mdlib/lincs.h"
  69 #include "gromacs/mdlib/shake.h"
  70 #include "gromacs/mdrunutility/multisim.h"
  71 #include "gromacs/mdtypes/commrec.h"
  72 #include "gromacs/mdtypes/inputrec.h"
  73 #include "gromacs/mdtypes/mdatom.h"
  74 #include "gromacs/pbcutil/pbc.h"
  75 #include "gromacs/topology/block.h"
  76 #include "gromacs/topology/idef.h"
  77 #include "gromacs/topology/ifunc.h"
  78 #include "gromacs/topology/topology.h"
  79 #include "gromacs/utility/smalloc.h"
  80 #include "gromacs/utility/stringutil.h"
  81 #include "gromacs/utility/unique_cptr.h"
  82
  83 namespace gmx
  84 {
  85 namespace test
  86 {
  87
  88 /* \brief
  89  * Constraints test data structure.
  90  *
  91  * Structure to collect all the necessary data, including system coordinates and topology,
  92  * constraints information, etc. The structure can be reset and reused.
  93  */
  94 class ConstraintsTestData
  95 {
  96     public:
  97         //! Human-friendly name for a system
  98         std::string           title_;
  99         //! Number of atoms
 100         int                   numAtoms_;
 101         //! Topology
 102         gmx_mtop_t            mtop_;
 103         //! Masses
 104         std::vector<real>     masses_;
 105         //! Inverse masses
 106         std::vector<real>     invmass_;
 107         //! Communication record
 108         t_commrec             cr_;
 109         //! Input record (info that usually in .mdp file)
 110         t_inputrec            ir_;
 111         //! Local topology
 112         t_idef                idef_;
 113         //! MD atoms
 114         t_mdatoms             md_;
 115         //! Multisim data
 116         gmx_multisim_t        ms_;
 117         //! Computational time array (normally used to benchmark performance)
 118         t_nrnb                nrnb_;
 119
 120         //! Inverse timestep
 121         real                  invdt_;
 122         //! Number of flexible constraints
 123         int                   nflexcon_   = 0;
 124         //! Whether the virial should be computed
 125         bool                  computeVirial_;
 126         //! Scaled virial
 127         tensor                virialScaled_;
 128         //! Scaled virial (reference values)
 129         tensor                virialScaledRef_;
 130         //! If the free energy is computed
 131         bool                  compute_dHdLambda_;
 132         //! For free energy computation
 133         real                  dHdLambda_;
 134         //! For free energy computation (reference value)
 135         real                  dHdLambdaRef_;
 136
 137         //! Coordinates before the timestep
 138         PaddedVector<RVec>    x_;
 139         //! Coordinates after timestep, output for the constraints
 140         PaddedVector<RVec>    xPrime_;
 141         //! Backup for coordinates (for reset)
 142         PaddedVector<RVec>    xPrime0_;
 143         //! Intermediate set of coordinates (normally used for projection correction)
 144         PaddedVector<RVec>    xPrime2_;
 145         //! Velocities
 146         PaddedVector<RVec>    v_;
 147         //! Backup for velocities (for reset)
 148         PaddedVector<RVec>    v0_;
 149
 150         //! Constraints data (type1-i1-j1-type2-i2-j2-...)
 151         std::vector<int>      constraints_;
 152         //! Target lengths for all constraint types
 153         std::vector<real>     constraintsR0_;
 154
 155         /*! \brief
 156          * Constructor for the object with all parameters and variables needed by constraints algorithms.
 157          *
 158          * This constructor assembles stubs for all the data structures, required to initialize
 159          * and apply LINCS and SHAKE constraints. The coordinates and velocities before constraining
 160          * are saved to allow for reset. The constraints data are stored for testing after constraints
 161          * were applied.
 162          *
 163          * \param[in]  title                Human-friendly name of the system.
 164          * \param[in]  numAtoms             Number of atoms in the system.
 165          * \param[in]  masses               Atom masses. Size of this vector should be equal to numAtoms.
 166          * \param[in]  constraints          List of constraints, organized in triples of integers.
 167          *                                  First integer is the index of type for a constraint, second
 168          *                                  and third are the indices of constrained atoms. The types
 169          *                                  of constraints should be sequential but not necessarily
 170          *                                  start from zero (which is the way they normally are in
 171          *                                  GROMACS).
 172          * \param[in]  constraintsR0        Target values for bond lengths for bonds of each type. The
 173          *                                  size of this vector should be equal to the total number of
 174          *                                  unique types in constraints vector.
 175          * \param[in]  computeVirial        Whether the virial should be computed.
 176          * \param[in]  virialScaledRef      Reference values for scaled virial tensor.
 177          * \param[in]  compute_dHdLambda    Whether free energy should be computed.
 178          * \param[in]  dHdLambdaRef         Reference value for dHdLambda.
 179          * \param[in]  initialTime          Initial time.
 180          * \param[in]  timestep             Timestep.
 181          * \param[in]  x                    Coordinates before integration step.
 182          * \param[in]  xPrime               Coordinates after integration step, but before constraining.
 183          * \param[in]  v                    Velocities before constraining.
 184          * \param[in]  shakeTolerance       Target tolerance for SHAKE.
 185          * \param[in]  shakeUseSOR          Use successive over-relaxation method for SHAKE iterations.
 186          *                                  The general formula is:
 187          *                                     x_n+1 = (1-omega)*x_n + omega*f(x_n),
 188          *                                  where omega = 1 if SOR is off and may be < 1 if SOR is on.
 189          * \param[in]  lincsNumIterations   Number of iterations used to compute the inverse matrix.
 190          * \param[in]  lincsExpansionOrder  The order for algorithm that adjusts the direction of the
 191          *                                  bond after constraints are applied.
 192          * \param[in]  lincsWarnAngle       The threshold value for the change in bond angle. When
 193          *                                  exceeded the program will issue a warning.
 194          *
 195          */
 196         ConstraintsTestData(const std::string &title,
 197                             int numAtoms, std::vector<real> masses,
 198                             std::vector<int> constraints, std::vector<real> constraintsR0,
 199                             bool computeVirial, tensor virialScaledRef,
 200                             bool compute_dHdLambda, float dHdLambdaRef,
 201                             real initialTime, real timestep,
 202                             const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
 203                             real shakeTolerance, gmx_bool shakeUseSOR,
 204                             int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle);
 205
 206         /*! \brief
 207          * Reset the data structure so it can be reused.
 208          *
 209          * Set the coordinates and velocities back to their values before
 210          * constraining. The scaled virial tensor and dHdLambda are zeroed.
 211          *
 212          */
 213         void reset();
 214
 215         /*! \brief
 216          * Cleaning up the memory.
 217          */
 218         ~ConstraintsTestData();
 219 };
 220
 221 /*! \brief Apply SHAKE constraints to the test data.
 222  */
 223 void applyShake(ConstraintsTestData *testData, t_pbc pbc);
 224 /*! \brief Apply LINCS constraints to the test data.
 225  */
 226 void applyLincs(ConstraintsTestData *testData, t_pbc pbc);
 227 /*! \brief Apply CUDA version of LINCS constraints to the test data.
 228  *
 229  * All the data is copied to the GPU device, then LINCS is applied and
 230  * the resulting coordinates are copied back.
 231  */
 232 void applyLincsCuda(ConstraintsTestData *testData, t_pbc pbc);
 233
 234 }      // namespace test
 235 }      // namespace gmx
 236
 237 #endif // GMX_MDLIB_TESTS_CONSTR_IMPL_H