src/gromacs/mdlib/tests/settle.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2016,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 #include "gmxpre.h"
  36
  37 #include "gromacs/mdlib/settle.h"
  38
  39 #include "config.h"
  40
  41 #include <tuple>
  42 #include <vector>
  43
  44 #include <gtest/gtest.h>
  45
  46 #include "gromacs/gpu_utils/gpu_utils.h"
  47 #include "gromacs/math/paddedvector.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/math/vectypes.h"
  50 #include "gromacs/mdlib/settle_cuda.h"
  51 #include "gromacs/mdtypes/mdatom.h"
  52 #include "gromacs/pbcutil/pbc.h"
  53 #include "gromacs/topology/idef.h"
  54 #include "gromacs/topology/ifunc.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/smalloc.h"
  57 #include "gromacs/utility/stringutil.h"
  58 #include "gromacs/utility/unique_cptr.h"
  59
  60 #include "gromacs/mdlib/tests/watersystem.h"
  61 #include "testutils/testasserts.h"
  62
  63 namespace gmx
  64 {
  65
  66 namespace test
  67 {
  68
  69 //! Simple cubic simulation box to use in tests
  70 matrix g_box = {{real(1.86206), 0, 0}, {0, real(1.86206), 0}, {0, 0, real(1.86206)}};
  71
  72 //! Convenience typedef
  73 typedef std::tuple<int, bool, bool, bool> SettleTestParameters;
  74
  75 /*! \brief Test fixture for testing SETTLE position updates
  76  *
  77  * \todo This also tests that if the calling code requires velocities
  78  * and virial updates, that those outputs do change, but does not test
  79  * that those changes are correct.
  80  *
  81  * \todo Only no-PBC and cubic-PBC are tested here, but the correct
  82  * function of the SIMD version of set_pbx_auic in all cases should be
  83  * tested elsewhere.
  84  */
  85 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
  86 {
  87     public:
  88         //! Updated water atom positions to constrain (DIM reals per atom)
  89         PaddedVector<gmx::RVec> updatedPositions_;
  90         //! Water atom velocities to constrain (DIM reals per atom)
  91         PaddedVector<gmx::RVec> velocities_;
  92         //! No periodic boundary conditions
  93         t_pbc                   pbcNone_;
  94         //! Rectangular periodic box
  95         t_pbc                   pbcXyz_;
  96
  97         SettleTest()
  98         {
  99             // Moved these to the constructor body so that uncrustify will not confuse doxygen
 100             updatedPositions_.resizeWithPadding(c_waterPositions.size());
 101             velocities_.resizeWithPadding(c_waterPositions.size());
 102
 103             std::copy(c_waterPositions.begin(), c_waterPositions.end(), updatedPositions_.begin());
 104
 105             // No periodic boundary conditions
 106             set_pbc(&pbcNone_, epbcNONE, g_box);
 107
 108             // Rectangular periodic box
 109             set_pbc(&pbcXyz_, epbcXYZ, g_box);
 110
 111             // Perturb the atom positions, to appear like an
 112             // "update," and where there is definitely constraining
 113             // work to do.
 114             for (int i = 0; i < updatedPositions_.size()*DIM; ++i)
 115             {
 116                 if (i % 4 == 0)
 117                 {
 118                     updatedPositions_[i / DIM][i % DIM] += 0.01;
 119                 }
 120                 else if (i % 4 == 1)
 121                 {
 122                     updatedPositions_[i / DIM][i % DIM] -= 0.01;
 123                 }
 124                 else if (i % 4 == 2)
 125                 {
 126                     updatedPositions_[i / DIM][i % DIM] += 0.02;
 127                 }
 128                 else if (i % 4 == 3)
 129                 {
 130                     updatedPositions_[i / DIM][i % DIM] -= 0.02;
 131                 }
 132                 velocities_[i / DIM][i % DIM] = 0.0;
 133             }
 134         }
 135 };
 136
 137 TEST_P(SettleTest, SatisfiesConstraints)
 138 {
 139     int  numSettles;
 140     bool usePbc, useVelocities, calcVirial;
 141     // Make some symbolic names for the parameter combination.
 142     std::tie(numSettles, usePbc, useVelocities, calcVirial) = GetParam();
 143
 144     // Make a string that describes which parameter combination is
 145     // being tested, to help make failing tests comprehensible.
 146     std::string testDescription = formatString("while testing %d SETTLEs, %sPBC, %svelocities and %scalculating the virial",
 147                                                numSettles,
 148                                                usePbc ? "with " : "without ",
 149                                                useVelocities ? "with " : "without ",
 150                                                calcVirial ? "" : "not ");
 151
 152     const int settleType     = 0;
 153     const int atomsPerSettle = NRAL(F_SETTLE);
 154     ASSERT_LE(numSettles, updatedPositions_.size() / atomsPerSettle) << "cannot test that many SETTLEs " << testDescription;
 155
 156     // Set up the topology.
 157     gmx_mtop_t mtop;
 158     mtop.moltype.resize(1);
 159     mtop.molblock.resize(1);
 160     mtop.molblock[0].type = 0;
 161     std::vector<int> &iatoms = mtop.moltype[0].ilist[F_SETTLE].iatoms;
 162     for (int i = 0; i < numSettles; ++i)
 163     {
 164         iatoms.push_back(settleType);
 165         iatoms.push_back(i*atomsPerSettle + 0);
 166         iatoms.push_back(i*atomsPerSettle + 1);
 167         iatoms.push_back(i*atomsPerSettle + 2);
 168     }
 169
 170     // Set up the SETTLE parameters.
 171     const real     dOH = 0.09572;
 172     const real     dHH = 0.15139;
 173     t_iparams      iparams;
 174     iparams.settle.doh = dOH;
 175     iparams.settle.dhh = dHH;
 176     mtop.ffparams.iparams.push_back(iparams);
 177
 178     // Set up the masses.
 179     t_mdatoms         mdatoms;
 180     std::vector<real> mass, massReciprocal;
 181     mtop.moltype[0].atoms.atom = static_cast<t_atom*>(calloc(numSettles*atomsPerSettle, sizeof(t_atom)));
 182     const real        oxygenMass = 15.9994, hydrogenMass = 1.008;
 183     for (int i = 0; i < numSettles; ++i)
 184     {
 185         mass.push_back(oxygenMass);
 186         mass.push_back(hydrogenMass);
 187         mass.push_back(hydrogenMass);
 188         massReciprocal.push_back(1./oxygenMass);
 189         massReciprocal.push_back(1./hydrogenMass);
 190         massReciprocal.push_back(1./hydrogenMass);
 191         mtop.moltype[0].atoms.atom[i*atomsPerSettle + 0].m = oxygenMass;
 192         mtop.moltype[0].atoms.atom[i*atomsPerSettle + 1].m = hydrogenMass;
 193         mtop.moltype[0].atoms.atom[i*atomsPerSettle + 2].m = hydrogenMass;
 194     }
 195     mdatoms.massT   = mass.data();
 196     mdatoms.invmass = massReciprocal.data();
 197     mdatoms.homenr  = numSettles * atomsPerSettle;
 198
 199     const t_ilist  ilist   = { mtop.moltype[0].ilist[F_SETTLE].size(), 0, mtop.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 };
 200
 201     // Copy the original positions from the array of doubles to a vector of reals
 202     PaddedVector<gmx::RVec> startingPositions;
 203     startingPositions.resizeWithPadding(c_waterPositions.size());
 204     std::copy(c_waterPositions.begin(), c_waterPositions.end(), startingPositions.begin());
 205
 206     const real              reciprocalTimeStep = 1.0/0.002;
 207     tensor                  virial             = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 208
 209 #if GMX_GPU == GMX_GPU_CUDA
 210     // Make a copy of all data-structures for GPU code testing
 211     PaddedVector<gmx::RVec> updatedPositionsGpu   = updatedPositions_;
 212     PaddedVector<gmx::RVec> velocitiesGpu         = velocities_;
 213     tensor                  virialGpu             = {{0, 0, 0}, {0, 0, 0}, {0, 0, 0}};
 214 #endif
 215
 216     // Finally make the settle data structures
 217     settledata    *settled = settle_init(mtop);
 218
 219     settle_set_constraints(settled, &ilist, mdatoms);
 220
 221     // Run the test
 222     bool       errorOccured;
 223     int        numThreads         = 1, threadIndex = 0;
 224     csettle(settled, numThreads, threadIndex,
 225             usePbc ? &pbcXyz_ : &pbcNone_,
 226             static_cast<real *>(startingPositions[0]), static_cast<real *>(updatedPositions_[0]), reciprocalTimeStep,
 227             useVelocities ? static_cast<real *>(velocities_[0]) : nullptr,
 228             calcVirial, virial, &errorOccured);
 229     settle_free(settled);
 230     EXPECT_FALSE(errorOccured) << testDescription;
 231
 232     // The necessary tolerances for the test to pass were determined
 233     // empirically. This isn't nice, but the required behavior that
 234     // SETTLE produces constrained coordinates consistent with
 235     // sensible sampling needs to be tested at a much higher level.
 236     FloatingPointTolerance tolerance =
 237         relativeToleranceAsPrecisionDependentUlp(dOH*dOH, 80, 380);
 238
 239     // Verify the updated coordinates match the requirements
 240     for (int i = 0; i < numSettles; ++i)
 241     {
 242         const gmx::RVec &positionO  = updatedPositions_[i*3 + 0];
 243         const gmx::RVec &positionH1 = updatedPositions_[i*3 + 1];
 244         const gmx::RVec &positionH2 = updatedPositions_[i*3 + 2];
 245
 246         EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH1), tolerance) << formatString("for water %d ", i) << testDescription;
 247         EXPECT_REAL_EQ_TOL(dOH*dOH, distance2(positionO, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
 248         EXPECT_REAL_EQ_TOL(dHH*dHH, distance2(positionH1, positionH2), tolerance) << formatString("for water %d ", i) << testDescription;
 249
 250         // This merely tests whether the velocities were
 251         // updated from the starting values of zero (or not),
 252         // but not whether the update was correct.
 253         for (int j = 0; j < atomsPerSettle; ++j)
 254         {
 255             for (int d = 0; d < DIM; ++d)
 256             {
 257                 EXPECT_TRUE(useVelocities == (0. != velocities_[i*3 + j][d]))
 258                 << formatString("for water %d velocity atom %d dim %d", i, j, d)
 259                 << testDescription;
 260             }
 261         }
 262     }
 263
 264     FloatingPointTolerance toleranceVirial = absoluteTolerance(0.000001);
 265     // This merely tests whether the viral was updated from
 266     // the starting values of zero (or not), but not whether
 267     // any update was correct.
 268     for (int d = 0; d < DIM; ++d)
 269     {
 270         for (int dd = 0; dd < DIM; ++dd)
 271         {
 272
 273             EXPECT_TRUE(calcVirial == (0. != virial[d][dd]))
 274             << formatString("for virial component[%d][%d] ", d, dd)
 275             << testDescription;
 276
 277             if (calcVirial)
 278             {
 279                 EXPECT_REAL_EQ_TOL(virial[d][dd], virial[dd][d], toleranceVirial)
 280                 << formatString("Virial is not symmetrical for [%d][%d] ", d, dd)
 281                 << testDescription;
 282             }
 283         }
 284     }
 285
 286     // CUDA version will be tested only if
 287     // 1. The code was compiled with cuda
 288     // 2. There is a CUDA-capable GPU in a system
 289 #if GMX_GPU == GMX_GPU_CUDA
 290     // TODO: Here we should check that at least 1 suitable GPU is available
 291     if (canPerformGpuDetection())
 292     {
 293         // Run the CUDA code and check if it gives identical results to CPU code
 294         t_idef                        idef;
 295         idef.il[F_SETTLE] = ilist;
 296
 297         std::unique_ptr<SettleCuda>   settleCuda = std::make_unique<SettleCuda>(mtop);
 298         settleCuda->setPbc(usePbc ? &pbcXyz_ : &pbcNone_);
 299         settleCuda->set(idef, mdatoms);
 300
 301         settleCuda->copyApplyCopy(mdatoms.homenr,
 302                                   as_rvec_array(startingPositions.data()),
 303                                   as_rvec_array(updatedPositionsGpu.data()),
 304                                   useVelocities,
 305                                   as_rvec_array(velocitiesGpu.data()),
 306                                   reciprocalTimeStep,
 307                                   calcVirial,
 308                                   virialGpu);
 309
 310         FloatingPointTolerance toleranceGpuCpu = absoluteTolerance(0.0001);
 311
 312         for (unsigned i = 0; i < updatedPositionsGpu.size(); ++i)
 313         {
 314             for (int d = 0; d < DIM; ++d)
 315             {
 316                 EXPECT_REAL_EQ_TOL(updatedPositionsGpu[i][d], updatedPositions_[i][d], toleranceGpuCpu)
 317                 << formatString("Different coordinates in GPU and CPU codes for particle %d component %d ", i, d)
 318                 << testDescription;
 319             }
 320         }
 321         for (unsigned i = 0; i < velocitiesGpu.size(); ++i)
 322         {
 323             for (int d = 0; d < DIM; ++d)
 324             {
 325                 EXPECT_REAL_EQ_TOL(velocitiesGpu[i][d], velocities_[i][d], toleranceGpuCpu)
 326                 << formatString("Different velocities in GPU and CPU codes for particle %d component %d ", i, d)
 327                 << testDescription;
 328             }
 329         }
 330         for (int d1 = 0; d1 < DIM; ++d1)
 331         {
 332             for (int d2 = 0; d2 < DIM; ++d2)
 333             {
 334                 EXPECT_REAL_EQ_TOL(virialGpu[d1][d2], virial[d1][d2], toleranceGpuCpu)
 335                 << formatString("Different virial components in GPU and CPU codes for components %d and %d ", d1, d2)
 336                 << testDescription;
 337             }
 338         }
 339     }
 340 #endif
 341 }
 342
 343 // Scan the full Cartesian product of numbers of SETTLE interactions
 344 // (4 and 17 are chosen to test cases that do and do not match
 345 // hardware SIMD widths), and whether or not we use PBC, velocities or
 346 // calculate the virial contribution.
 347 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest,
 348                             ::testing::Combine(::testing::Values(1, 2, 4, 5, 7, 10, 12, 15, 17),
 349                                                    ::testing::Bool(),
 350                                                    ::testing::Bool(),
 351                                                    ::testing::Bool()));
 352
 353 }  // namespace test
 354 }  // namespace gmx