src/gromacs/mdlib/tests/shake.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2014,2015,2017,2018,2019,2020,2021, 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/shake.h"
  38
  39 #include <assert.h>
  40
  41 #include <cmath>
  42
  43 #include <algorithm>
  44 #include <vector>
  45
  46 #include <gtest/gtest.h>
  47
  48 #include "gromacs/utility/arrayref.h"
  49
  50 #include "testutils/refdata.h"
  51 #include "testutils/testasserts.h"
  52
  53 namespace gmx
  54 {
  55 namespace
  56 {
  57
  58 /*! \brief Stride of the vector of int used to describe each SHAKE
  59  * constraint
  60  *
  61  * Like other such code, SHAKE is hard-wired to use t_ilist.iatoms as
  62  * a flat vector of tuples of general data. Here, they are triples
  63  * containing the index of the constraint type, and then the indices
  64  * of the two atoms involved. So for each constraint, we must stride
  65  * this vector by three to get access to its information. */
  66 const int constraintStride = 3;
  67
  68 /*! \brief Compute the displacements between pairs of constrained
  69  * atoms described in the iatom "topology". */
  70 std::vector<RVec> computeDisplacements(ArrayRef<const int> iatom, const std::vector<RVec>& positions)
  71 {
  72     assert(0 == iatom.size() % constraintStride);
  73     int               numConstraints = iatom.size() / constraintStride;
  74     std::vector<RVec> displacements;
  75
  76     for (int ll = 0; ll != numConstraints; ++ll)
  77     {
  78         int atom_i = iatom[ll * constraintStride + 1];
  79         int atom_j = iatom[ll * constraintStride + 2];
  80
  81         displacements.push_back(positions[atom_i] - positions[atom_j]);
  82     }
  83
  84     return displacements;
  85 }
  86
  87 /*! \brief Compute half of the reduced mass of each pair of constrained
  88  * atoms in the iatom "topology".
  89  *
  90  * The reduced mass is m = 1/(1/m_i + 1/m_j)) */
  91 std::vector<real> computeHalfOfReducedMasses(const std::vector<int>&  iatom,
  92                                              const std::vector<real>& inverseMasses)
  93 {
  94     int               numConstraints = iatom.size() / constraintStride;
  95     std::vector<real> halfOfReducedMasses;
  96
  97     for (int ll = 0; ll != numConstraints; ++ll)
  98     {
  99         int atom_i = iatom[ll * constraintStride + 1];
 100         int atom_j = iatom[ll * constraintStride + 2];
 101
 102         halfOfReducedMasses.push_back(0.5 / (inverseMasses[atom_i] + inverseMasses[atom_j]));
 103     }
 104
 105     return halfOfReducedMasses;
 106 }
 107
 108 /*! \brief Compute the distances corresponding to the vector of displacements components */
 109 std::vector<real> computeDistancesSquared(ArrayRef<const RVec> displacements)
 110 {
 111     int               numDistancesSquared = displacements.size();
 112     std::vector<real> distanceSquared;
 113
 114     for (int i = 0; i != numDistancesSquared; ++i)
 115     {
 116         distanceSquared.push_back(0.0);
 117         for (int d = 0; d != DIM; ++d)
 118         {
 119             real displacement = displacements[i][d];
 120             distanceSquared.back() += displacement * displacement;
 121         }
 122     }
 123
 124     return distanceSquared;
 125 }
 126
 127 /*! \brief Test fixture for testing SHAKE */
 128 class ShakeTest : public ::testing::Test
 129 {
 130 public:
 131     /*! \brief Set up data for test cases to use when constructing
 132         their input */
 133     void SetUp() override
 134     {
 135         inverseMassesDatabase_.push_back(2.0);
 136         inverseMassesDatabase_.push_back(3.0);
 137         inverseMassesDatabase_.push_back(4.0);
 138         inverseMassesDatabase_.push_back(1.0);
 139
 140         positionsDatabase_.emplace_back(2.5, -3.1, 15.7);
 141
 142         positionsDatabase_.emplace_back(0.51, -3.02, 15.55);
 143
 144         positionsDatabase_.emplace_back(-0.5, -3.0, 15.2);
 145
 146         positionsDatabase_.emplace_back(-1.51, -2.95, 15.05);
 147     }
 148
 149     //! Run the test
 150     static void runTest(size_t gmx_unused        numAtoms,
 151                         size_t                   numConstraints,
 152                         const std::vector<int>&  iatom,
 153                         const std::vector<real>& constrainedDistances,
 154                         const std::vector<real>& inverseMasses,
 155                         const std::vector<RVec>& positions)
 156     {
 157         // Check the test input is consistent
 158         assert(numConstraints * constraintStride == iatom.size());
 159         assert(numConstraints == constrainedDistances.size());
 160         assert(numAtoms == inverseMasses.size());
 161         assert(numAtoms == positions.size());
 162         for (size_t i = 0; i != numConstraints; ++i)
 163         {
 164             for (size_t j = 1; j < 3; j++)
 165             {
 166                 // Check that the topology refers to atoms that have masses and positions
 167                 assert(iatom[i * constraintStride + j] >= 0);
 168                 assert(iatom[i * constraintStride + j] < static_cast<int>(numAtoms));
 169             }
 170         }
 171         std::vector<real> distanceSquaredTolerances;
 172         std::vector<real> lagrangianValues;
 173         std::vector<real> constrainedDistancesSquared;
 174
 175         real coordMax = 0;
 176         for (size_t i = 0; i != numConstraints; ++i)
 177         {
 178             constrainedDistancesSquared.push_back(constrainedDistances[i] * constrainedDistances[i]);
 179             distanceSquaredTolerances.push_back(
 180                     0.5 / (constrainedDistancesSquared.back() * ShakeTest::tolerance_));
 181             lagrangianValues.push_back(0.0);
 182
 183             for (size_t j = 1; j < constraintStride; j++)
 184             {
 185                 for (int d = 0; d < DIM; d++)
 186                 {
 187                     coordMax = std::max(coordMax, std::abs(positions[iatom[i * constraintStride + j]][d]));
 188                 }
 189             }
 190         }
 191         std::vector<real> halfOfReducedMasses  = computeHalfOfReducedMasses(iatom, inverseMasses);
 192         std::vector<RVec> initialDisplacements = computeDisplacements(iatom, positions);
 193
 194         std::vector<RVec> finalPositions = positions;
 195         int               numIterations  = 0;
 196         int               numErrors      = 0;
 197
 198         cshake(iatom.data(),
 199                numConstraints,
 200                &numIterations,
 201                ShakeTest::maxNumIterations_,
 202                constrainedDistancesSquared,
 203                finalPositions,
 204                nullptr,
 205                initialDisplacements,
 206                halfOfReducedMasses,
 207                omega_,
 208                inverseMasses,
 209                distanceSquaredTolerances,
 210                lagrangianValues,
 211                &numErrors);
 212
 213         std::vector<RVec> finalDisplacements    = computeDisplacements(iatom, finalPositions);
 214         std::vector<real> finalDistancesSquared = computeDistancesSquared(finalDisplacements);
 215         assert(numConstraints == finalDistancesSquared.size());
 216
 217         EXPECT_EQ(0, numErrors);
 218         EXPECT_GT(numIterations, 1);
 219         EXPECT_LT(numIterations, ShakeTest::maxNumIterations_);
 220         // TODO wrap this in a Google Mock matcher if there's
 221         // other tests like it some time?
 222         for (size_t i = 0; i != numConstraints; ++i)
 223         {
 224             // We need to allow for the requested tolerance plus rounding
 225             // errors due to the absolute size of the coordinate values
 226             test::FloatingPointTolerance constraintTolerance =
 227                     test::absoluteTolerance(std::sqrt(constrainedDistancesSquared[i]) * ShakeTest::tolerance_
 228                                             + coordMax * GMX_REAL_EPS);
 229             // Assert that the constrained distances are within the required tolerance
 230             EXPECT_FLOAT_EQ_TOL(std::sqrt(constrainedDistancesSquared[i]),
 231                                 std::sqrt(finalDistancesSquared[i]),
 232                                 constraintTolerance);
 233         }
 234     }
 235
 236     //! Tolerance for SHAKE conversion (ie. shake-tol .mdp setting)
 237     static const real tolerance_;
 238     //! Maximum number of iterations permitted in these tests
 239     static const int maxNumIterations_;
 240     //! SHAKE over-relaxation (SOR) factor
 241     static const real omega_;
 242     //! Database of inverse masses of atoms in the topology
 243     std::vector<real> inverseMassesDatabase_;
 244     //! Database of atom positions (three reals per atom)
 245     std::vector<RVec> positionsDatabase_;
 246 };
 247
 248 const real ShakeTest::tolerance_        = 1e-5;
 249 const int  ShakeTest::maxNumIterations_ = 30;
 250 const real ShakeTest::omega_            = 1.0;
 251
 252 TEST_F(ShakeTest, ConstrainsOneBond)
 253 {
 254     int numAtoms       = 2;
 255     int numConstraints = 1;
 256
 257     std::vector<int> iatom;
 258     iatom.push_back(-1); // unused
 259     iatom.push_back(0);  // i atom index
 260     iatom.push_back(1);  // j atom index
 261
 262     std::vector<real> constrainedDistances;
 263     constrainedDistances.push_back(2.0);
 264
 265     std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
 266                                     inverseMassesDatabase_.begin() + numAtoms);
 267     std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
 268
 269     runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
 270 }
 271
 272 TEST_F(ShakeTest, ConstrainsTwoDisjointBonds)
 273 {
 274     int numAtoms       = 4;
 275     int numConstraints = 2;
 276
 277     std::vector<int> iatom;
 278     iatom.push_back(-1); // unused
 279     iatom.push_back(0);  // i atom index
 280     iatom.push_back(1);  // j atom index
 281
 282     iatom.push_back(-1); // unused
 283     iatom.push_back(2);  // i atom index
 284     iatom.push_back(3);  // j atom index
 285
 286     std::vector<real> constrainedDistances;
 287     constrainedDistances.push_back(2.0);
 288     constrainedDistances.push_back(1.0);
 289
 290     std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
 291                                     inverseMassesDatabase_.begin() + numAtoms);
 292     std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
 293
 294     runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
 295 }
 296
 297 TEST_F(ShakeTest, ConstrainsTwoBondsWithACommonAtom)
 298 {
 299     int numAtoms       = 3;
 300     int numConstraints = 2;
 301
 302     std::vector<int> iatom;
 303     iatom.push_back(-1); // unused
 304     iatom.push_back(0);  // i atom index
 305     iatom.push_back(1);  // j atom index
 306
 307     iatom.push_back(-1); // unused
 308     iatom.push_back(1);  // i atom index
 309     iatom.push_back(2);  // j atom index
 310
 311     std::vector<real> constrainedDistances;
 312     constrainedDistances.push_back(2.0);
 313     constrainedDistances.push_back(1.0);
 314
 315     std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
 316                                     inverseMassesDatabase_.begin() + numAtoms);
 317     std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
 318
 319     runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
 320 }
 321
 322 TEST_F(ShakeTest, ConstrainsThreeBondsWithCommonAtoms)
 323 {
 324     int numAtoms       = 4;
 325     int numConstraints = 3;
 326
 327     std::vector<int> iatom;
 328     iatom.push_back(-1); // unused
 329     iatom.push_back(0);  // i atom index
 330     iatom.push_back(1);  // j atom index
 331
 332     iatom.push_back(-1); // unused
 333     iatom.push_back(1);  // i atom index
 334     iatom.push_back(2);  // j atom index
 335
 336     iatom.push_back(-1); // unused
 337     iatom.push_back(2);  // i atom index
 338     iatom.push_back(3);  // j atom index
 339
 340     std::vector<real> constrainedDistances;
 341     constrainedDistances.push_back(2.0);
 342     constrainedDistances.push_back(1.0);
 343     constrainedDistances.push_back(1.0);
 344
 345     std::vector<real> inverseMasses(inverseMassesDatabase_.begin(),
 346                                     inverseMassesDatabase_.begin() + numAtoms);
 347     std::vector<RVec> positions(positionsDatabase_.begin(), positionsDatabase_.begin() + numAtoms);
 348
 349     runTest(numAtoms, numConstraints, iatom, constrainedDistances, inverseMasses, positions);
 350 }
 351
 352 } // namespace
 353 } // namespace gmx