src/gromacs/mdlib/tests/shake.cpp

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