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