src/gromacs/mdlib/tests/leapfrog.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 Integrator tests
  37  *
  38  * Test for CUDA implementation of the Leap-Frog integrator.
  39  *
  40  * \todo Connect to the CPU-based version.
  41  * \todo Prepare for temperature and pressure controlled integrators.
  42  * \todo Add PBC handling test.
  43  * \todo Take over coordinates/velocities/forces handling - this infrastructure
  44  *       will be removed from integrator.
  45  *
  46  * \author Artem Zhmurov <zhmurov@gmail.com>
  47  * \ingroup module_mdlib
  48  */
  49
  50 #include "gmxpre.h"
  51
  52 #include "config.h"
  53
  54 #include <assert.h>
  55
  56 #include <cmath>
  57
  58 #include <algorithm>
  59 #include <unordered_map>
  60 #include <vector>
  61
  62 #include <gtest/gtest.h>
  63
  64 #include "gromacs/gpu_utils/gpu_utils.h"
  65 #include "gromacs/math/vec.h"
  66 #include "gromacs/mdlib/leapfrog_cuda.h"
  67 #include "gromacs/utility/stringutil.h"
  68
  69 #include "testutils/refdata.h"
  70 #include "testutils/testasserts.h"
  71
  72 namespace gmx
  73 {
  74 namespace test
  75 {
  76
  77 #if GMX_GPU == GMX_GPU_CUDA
  78
  79 /*! \brief The parameter space for the test.
  80  *
  81  * The test will run for all possible combinations of accessible
  82  * values of the:
  83  * 1. Number of atoms
  84  * 2. Timestep
  85  * 3-5. Velocity components
  86  * 6-8. Force components
  87  * 9. Number of steps
  88  */
  89 typedef std::tuple<int, real, real, real, real, real, real, real, int> IntegratorTestParameters;
  90
  91 /*! \brief Test fixture for integrator.
  92  *
  93  *  Creates a system of independent particles exerting constant external forces,
  94  *  makes several numerical integration timesteps and compares the result
  95  *  with analytical solution.
  96  *
  97  */
  98 class IntegratorTest : public ::testing::TestWithParam<IntegratorTestParameters>
  99 {
 100     public:
 101         //! Number of atoms in the system
 102         int  numAtoms_;
 103         //! Integration timestep
 104         real timestep_;
 105
 106         //! Initial coordinates
 107         std::vector<RVec> x0_;
 108         //! Current coordinates
 109         std::vector<RVec> x_;
 110         //! Coordinates after integrator update
 111         std::vector<RVec> xPrime_;
 112         //! Initial velocities
 113         std::vector<RVec> v0_;
 114         //! Current velocities
 115         std::vector<RVec> v_;
 116         //! External forces
 117         std::vector<RVec> f_;
 118         //! Inverse masses of the particles
 119         std::vector<real> inverseMasses_;
 120         //! MD atoms structure in which inverse masses will be passed to the integrator
 121         t_mdatoms         mdAtoms_;
 122
 123         /*! \brief Test initialization function.
 124          *
 125          * \param[in]  numAtoms  Number of atoms in the system
 126          * \param[in]  timestep  Integration timestep
 127          * \param[in]  v0        Initial velocity (same for all particles)
 128          * \param[in]  f0        External constant force, acting on all particles
 129          */
 130         void init(int numAtoms, real timestep, rvec v0, rvec f0)
 131         {
 132             numAtoms_    = numAtoms;
 133             timestep_    = timestep;
 134
 135             x0_.resize(numAtoms_);
 136             x_.resize(numAtoms_);
 137             xPrime_.resize(numAtoms_);
 138             v0_.resize(numAtoms_);
 139             v_.resize(numAtoms_);
 140             f_.resize(numAtoms_);
 141             inverseMasses_.resize(numAtoms_);
 142
 143             mdAtoms_.nr = numAtoms_;
 144
 145             for (unsigned i = 0; i < x_.size(); i++)
 146             {
 147                 // Typical PBC box size is tens of nanometers
 148                 x_[i][XX] = (i%21)*1.0;
 149                 x_[i][YY] = 6.5 + (i%13)*(-1.0);
 150                 x_[i][ZZ] = (i%32)*(0.0);
 151
 152                 for (int d = 0; d < DIM; d++)
 153                 {
 154                     xPrime_[i][d] = 0.0;
 155                     // Thermal velocity is ~1 nm/ps (|v0| = 1-2 nm/ps)
 156                     v_[i][d] = v0[d];
 157                     // TODO Check what value typical MD forces have (now ~ 1 kJ/mol/nm)
 158                     f_[i][d] = f0[d];
 159
 160                     x0_[i][d] = x_[i][d];
 161                     v0_[i][d] = v_[i][d];
 162                 }
 163                 // Atom masses are ~1-100 g/mol
 164                 inverseMasses_[i] = 1.0/(1.0 + i%100);
 165             }
 166             mdAtoms_.invmass = inverseMasses_.data();
 167         }
 168
 169 };
 170
 171 TEST_P(IntegratorTest, SimpleIntegration)
 172 {
 173     // TODO: Here we should check that at least 1 suitable GPU is available
 174     if (!canPerformGpuDetection())
 175     {
 176         return;
 177     }
 178     int  numAtoms;    // 1. Number of atoms
 179     real timestep;    // 2. Timestep
 180     rvec v0;          // 3. Velocity
 181     rvec f0;          // 4. Force
 182     int  nStep;       // 5. Number of steps
 183     std::tie(numAtoms, timestep, v0[XX], v0[YY], v0[ZZ], f0[XX], f0[YY], f0[ZZ], nStep) = GetParam();
 184
 185     std::string testDescription = formatString("while testing %d atoms for %d timestep (dt = %f, v0=(%f, %f, %f), f0=(%f, %f, %f))",
 186                                                numAtoms, nStep, timestep,
 187                                                v0[XX], v0[YY], v0[ZZ],
 188                                                f0[XX], f0[YY], f0[ZZ]);
 189
 190
 191     init(numAtoms, timestep, v0, f0);
 192
 193     std::unique_ptr<LeapFrogCuda> integrator = std::make_unique<LeapFrogCuda>();
 194
 195     integrator->set(mdAtoms_);
 196
 197     for (int step = 0; step < nStep; step++)
 198     {
 199         integrator->copyIntegrateCopy(numAtoms,
 200                                       (rvec*)(x_.data()),
 201                                       (rvec*)(xPrime_.data()),
 202                                       (rvec*)(v_.data()),
 203                                       (rvec*)(f_.data()),
 204                                       timestep_);
 205
 206         for (int i = 0; i < numAtoms; i++)
 207         {
 208             for (int d = 0; d < DIM; d++)
 209             {
 210                 x_.at(i)[d] = xPrime_.at(i)[d];
 211             }
 212         }
 213     }
 214
 215     real totalTime = nStep*timestep;
 216     // TODO For the case of constant force, the numerical scheme is exact and
 217     //      the only source of errors is floating point arithmetic. Hence,
 218     //      the tolerance can be calculated.
 219     FloatingPointTolerance tolerance = absoluteTolerance(nStep*0.000005);
 220
 221     for (unsigned i = 0; i < x_.size(); i++)
 222     {
 223         rvec xAnalytical;
 224         rvec vAnalytical;
 225         for (int d = 0; d < DIM; d++)
 226         {
 227             // Analytical solution for constant-force particle movement
 228             xAnalytical[d] = x0_[i][d] + v0_[i][d]*totalTime + 0.5*f_[i][d]*totalTime*totalTime*inverseMasses_[i];
 229             vAnalytical[d] = v0_[i][d] + f_[i][d]*totalTime*inverseMasses_[i];
 230
 231             EXPECT_REAL_EQ_TOL(xPrime_[i][d], xAnalytical[d], tolerance)
 232             << gmx::formatString("Coordinate %d of atom %d is different from analytical solution", d, i)
 233             << testDescription;
 234
 235             EXPECT_REAL_EQ_TOL(v_[i][d], vAnalytical[d], tolerance)
 236             << gmx::formatString("Velocity component %d of atom %d is different from analytical solution", d, i)
 237             << testDescription;
 238         }
 239     }
 240 }
 241
 242 INSTANTIATE_TEST_CASE_P(WithParameters, IntegratorTest,
 243                             ::testing::Combine(
 244                                     ::testing::Values(1, 10, 300),    // Number of atoms
 245                                     ::testing::Values(0.001, 0.0005), // Timestep
 246                                     ::testing::Values(-2.0, 0.0),     // vx
 247                                     ::testing::Values( 0.0, 2.0),     // vy
 248                                     ::testing::Values( 0.0),          // vz
 249                                     ::testing::Values(-1.0, 0.0),     // fx
 250                                     ::testing::Values( 0.0, 1.0),     // fy
 251                                     ::testing::Values( 2.0),          // fz
 252                                     ::testing::Values(1, 10)          // Number of steps
 253                                 ));
 254 #endif                                                                // GMX_GPU == GMX_GPU_CUDA
 255
 256 }                                                                     // namespace test
 257 }                                                                     // namespace gmx