src/gromacs/mdlib/tests/shake.cpp

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