src/gromacs/mdlib/tests/constr.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 /*! \internal \file
  36  * \brief SHAKE and LINCS tests.
  37  *
  38  * \todo Better tests for virial are needed.
  39  * \todo Tests for bigger systems to test threads synchronization,
  40  *       reduction, etc. on the GPU.
  41  * \todo Tests for algorithms for derivatives.
  42  * \todo Free-energy perturbation tests
  43  *
  44  * \author Artem Zhmurov <zhmurov@gmail.com>
  45  * \ingroup module_mdlib
  46  */
  47 #include "gmxpre.h"
  48
  49 #include "config.h"
  50
  51 #include <assert.h>
  52
  53 #include <unordered_map>
  54 #include <vector>
  55
  56 #include <gtest/gtest.h>
  57
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/math/vectypes.h"
  60 #include "gromacs/pbcutil/pbc.h"
  61 #include "gromacs/utility/stringutil.h"
  62
  63 #include "testutils/refdata.h"
  64 #include "testutils/test_hardware_environment.h"
  65 #include "testutils/testasserts.h"
  66
  67 #include "constrtestdata.h"
  68 #include "constrtestrunners.h"
  69
  70 //! Helper function to convert t_pbc into string and make test failure messages readable
  71 static void PrintTo(const t_pbc& pbc, std::ostream* os)
  72 {
  73     *os << "PBC: " << c_pbcTypeNames[pbc.pbcType];
  74 }
  75
  76
  77 namespace gmx
  78 {
  79
  80 namespace test
  81 {
  82 namespace
  83 {
  84
  85 // Define the set of PBCs to run the test for
  86 const std::vector<t_pbc> c_pbcs = [] {
  87     std::vector<t_pbc> pbcs;
  88     t_pbc              pbc;
  89
  90     // Infinitely small box
  91     matrix boxNone = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
  92     set_pbc(&pbc, PbcType::No, boxNone);
  93     pbcs.emplace_back(pbc);
  94
  95     // Rectangular box
  96     matrix boxXyz = { { 10.0, 0.0, 0.0 }, { 0.0, 20.0, 0.0 }, { 0.0, 0.0, 15.0 } };
  97     set_pbc(&pbc, PbcType::Xyz, boxXyz);
  98     pbcs.emplace_back(pbc);
  99
 100     return pbcs;
 101 }();
 102
 103 struct ConstraintsTestSystem
 104 {
 105     //! Human-friendly name of the system.
 106     std::string title;
 107     //! Number of atoms in the system.
 108     int numAtoms;
 109     //! Atom masses. Size of this vector should be equal to numAtoms.
 110     std::vector<real> masses;
 111     /*! \brief List of constraints, organized in triples of integers.
 112      *
 113      * First integer is the index of type for a constraint, second
 114      * and third are the indices of constrained atoms. The types
 115      * of constraints should be sequential but not necessarily
 116      * start from zero (which is the way they normally are in
 117      * GROMACS).
 118      */
 119     std::vector<int> constraints;
 120     /*! \brief Target values for bond lengths for bonds of each type.
 121      *
 122      * The size of this vector should be equal to the total number of
 123      * unique types in constraints vector.
 124      */
 125     std::vector<real> constraintsR0;
 126     //! Coordinates before integration step.
 127     std::vector<RVec> x;
 128     //! Coordinates after integration step, but before constraining.
 129     std::vector<RVec> xPrime;
 130     //! Velocities before constraining.
 131     std::vector<RVec> v;
 132
 133     //! Target tolerance for SHAKE.
 134     real shakeTolerance = 0.0001;
 135     /*! \brief Use successive over-relaxation method for SHAKE iterations.
 136      *
 137      * The general formula is:
 138      * x_n+1 = (1-omega)*x_n + omega*f(x_n),
 139      * where omega = 1 if SOR is off and may be < 1 if SOR is on.
 140      */
 141     bool shakeUseSOR = false;
 142
 143     //! Number of iterations used to compute the inverse matrix.
 144     int lincsNIter = 1;
 145     //! The order for algorithm that adjusts the direction of the bond after constraints are applied.
 146     int lincslincsExpansionOrder = 4;
 147     //! The threshold value for the change in bond angle. When exceeded the program will issue a warning
 148     real lincsWarnAngle = 30.0;
 149 };
 150
 151 //! Helper function to convert ConstraintsTestSystem into string and make test failure messages readable
 152 void PrintTo(const ConstraintsTestSystem& constraintsTestSystem, std::ostream* os)
 153 {
 154     *os << constraintsTestSystem.title << " - " << constraintsTestSystem.numAtoms << " atoms";
 155 }
 156
 157 const std::vector<ConstraintsTestSystem> c_constraintsTestSystemList = [] {
 158     std::vector<ConstraintsTestSystem> constraintsTestSystemList;
 159     {
 160         ConstraintsTestSystem constraintsTestSystem;
 161
 162         constraintsTestSystem.title    = "one constraint (e.g. OH)";
 163         constraintsTestSystem.numAtoms = 2;
 164
 165         constraintsTestSystem.masses        = { 1.0, 12.0 };
 166         constraintsTestSystem.constraints   = { 0, 0, 1 };
 167         constraintsTestSystem.constraintsR0 = { 0.1 };
 168
 169         real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
 170
 171         constraintsTestSystem.x = { { 0.0, oneTenthOverSqrtTwo, 0.0 }, { oneTenthOverSqrtTwo, 0.0, 0.0 } };
 172         constraintsTestSystem.xPrime = { { 0.01, 0.08, 0.01 }, { 0.06, 0.01, -0.01 } };
 173         constraintsTestSystem.v      = { { 1.0, 2.0, 3.0 }, { 3.0, 2.0, 1.0 } };
 174
 175         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 176     }
 177     {
 178         ConstraintsTestSystem constraintsTestSystem;
 179
 180         constraintsTestSystem.title         = "two disjoint constraints";
 181         constraintsTestSystem.numAtoms      = 4;
 182         constraintsTestSystem.masses        = { 0.5, 1.0 / 3.0, 0.25, 1.0 };
 183         constraintsTestSystem.constraints   = { 0, 0, 1, 1, 2, 3 };
 184         constraintsTestSystem.constraintsR0 = { 2.0, 1.0 };
 185
 186
 187         constraintsTestSystem.x = { { 2.50, -3.10, 15.70 },
 188                                     { 0.51, -3.02, 15.55 },
 189                                     { -0.50, -3.00, 15.20 },
 190                                     { -1.51, -2.95, 15.05 } };
 191
 192         constraintsTestSystem.xPrime = { { 2.50, -3.10, 15.70 },
 193                                          { 0.51, -3.02, 15.55 },
 194                                          { -0.50, -3.00, 15.20 },
 195                                          { -1.51, -2.95, 15.05 } };
 196
 197         constraintsTestSystem.v = { { 0.0, 1.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 0.0, 0.0, 1.0 }, { 0.0, 0.0, 0.0 } };
 198
 199         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 200     }
 201     {
 202         ConstraintsTestSystem constraintsTestSystem;
 203
 204         constraintsTestSystem.title         = "three atoms, connected longitudinally (e.g. CH2)";
 205         constraintsTestSystem.numAtoms      = 3;
 206         constraintsTestSystem.masses        = { 1.0, 12.0, 16.0 };
 207         constraintsTestSystem.constraints   = { 0, 0, 1, 1, 1, 2 };
 208         constraintsTestSystem.constraintsR0 = { 0.1, 0.2 };
 209
 210         real oneTenthOverSqrtTwo    = 0.1_real / std::sqrt(2.0_real);
 211         real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
 212
 213         constraintsTestSystem.x = { { oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
 214                                     { 0.0, 0.0, 0.0 },
 215                                     { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree } };
 216
 217         constraintsTestSystem.xPrime = { { 0.08, 0.07, 0.01 }, { -0.02, 0.01, -0.02 }, { 0.10, 0.12, 0.11 } };
 218
 219         constraintsTestSystem.v = { { 1.0, 0.0, 0.0 }, { 0.0, 1.0, 0.0 }, { 0.0, 0.0, 1.0 } };
 220
 221         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 222     }
 223     {
 224         ConstraintsTestSystem constraintsTestSystem;
 225
 226         constraintsTestSystem.title         = "four atoms, connected longitudinally";
 227         constraintsTestSystem.numAtoms      = 4;
 228         constraintsTestSystem.masses        = { 0.5, 1.0 / 3.0, 0.25, 1.0 };
 229         constraintsTestSystem.constraints   = { 0, 0, 1, 1, 1, 2, 2, 2, 3 };
 230         constraintsTestSystem.constraintsR0 = { 2.0, 1.0, 1.0 };
 231
 232
 233         constraintsTestSystem.x = { { 2.50, -3.10, 15.70 },
 234                                     { 0.51, -3.02, 15.55 },
 235                                     { -0.50, -3.00, 15.20 },
 236                                     { -1.51, -2.95, 15.05 } };
 237
 238         constraintsTestSystem.xPrime = { { 2.50, -3.10, 15.70 },
 239                                          { 0.51, -3.02, 15.55 },
 240                                          { -0.50, -3.00, 15.20 },
 241                                          { -1.51, -2.95, 15.05 } };
 242
 243         constraintsTestSystem.v = {
 244             { 0.0, 0.0, 2.0 }, { 0.0, 0.0, 3.0 }, { 0.0, 0.0, -4.0 }, { 0.0, 0.0, -1.0 }
 245         };
 246
 247         // Overriding default values since LINCS converges slowly for this system.
 248         constraintsTestSystem.lincsNIter               = 4;
 249         constraintsTestSystem.lincslincsExpansionOrder = 8;
 250
 251         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 252     }
 253     {
 254         ConstraintsTestSystem constraintsTestSystem;
 255
 256         constraintsTestSystem.title       = "three atoms, connected to the central atom (e.g. CH3)";
 257         constraintsTestSystem.numAtoms    = 4;
 258         constraintsTestSystem.masses      = { 12.0, 1.0, 1.0, 1.0 };
 259         constraintsTestSystem.constraints = { 0, 0, 1, 0, 0, 2, 0, 0, 3 };
 260         constraintsTestSystem.constraintsR0 = { 0.1 };
 261
 262
 263         constraintsTestSystem.x = {
 264             { 0.00, 0.00, 0.00 }, { 0.10, 0.00, 0.00 }, { 0.00, -0.10, 0.00 }, { 0.00, 0.00, 0.10 }
 265         };
 266
 267         constraintsTestSystem.xPrime = { { 0.004, 0.009, -0.010 },
 268                                          { 0.110, -0.006, 0.003 },
 269                                          { -0.007, -0.102, -0.007 },
 270                                          { -0.005, 0.011, 0.102 } };
 271
 272         constraintsTestSystem.v = { { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 }, { 1.0, 0.0, 0.0 } };
 273
 274         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 275     }
 276     {
 277         ConstraintsTestSystem constraintsTestSystem;
 278
 279         constraintsTestSystem.title       = "basic triangle (three atoms, connected to each other)";
 280         constraintsTestSystem.numAtoms    = 3;
 281         constraintsTestSystem.masses      = { 1.0, 1.0, 1.0 };
 282         constraintsTestSystem.constraints = { 0, 0, 1, 2, 0, 2, 1, 1, 2 };
 283         constraintsTestSystem.constraintsR0 = { 0.1, 0.1, 0.1 };
 284
 285         real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
 286
 287         constraintsTestSystem.x = { { oneTenthOverSqrtTwo, 0.0, 0.0 },
 288                                     { 0.0, oneTenthOverSqrtTwo, 0.0 },
 289                                     { 0.0, 0.0, oneTenthOverSqrtTwo } };
 290
 291         constraintsTestSystem.xPrime = { { 0.09, -0.02, 0.01 }, { -0.02, 0.10, -0.02 }, { 0.03, -0.01, 0.07 } };
 292
 293         constraintsTestSystem.v = { { 1.0, 1.0, 1.0 }, { -2.0, -2.0, -2.0 }, { 1.0, 1.0, 1.0 } };
 294
 295         constraintsTestSystemList.emplace_back(constraintsTestSystem);
 296     }
 297
 298     return constraintsTestSystemList;
 299 }();
 300
 301
 302 /*! \brief Test fixture for constraints.
 303  *
 304  * The fixture uses following test systems:
 305  * 1. Two atoms, connected with one constraint (e.g. NH).
 306  * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
 307  * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
 308  * 4. Four atoms, connected by two independent constraints.
 309  * 5. Three atoms, connected by three constraints in a triangle
 310  *      (e.g. H2O with constrained H-O-H angle).
 311  * 6. Four atoms, connected by three consequential constraints.
 312  *
 313  * For all systems, the final lengths of the constraints are tested against the
 314  * reference values, the direction of each constraint is checked.
 315  * Test also verifies that the center of mass has not been
 316  * shifted by the constraints and that its velocity has not changed.
 317  * For some systems, the value for scaled virial tensor is checked against
 318  * pre-computed data.
 319  */
 320 class ConstraintsTest : public ::testing::TestWithParam<std::tuple<ConstraintsTestSystem, t_pbc>>
 321 {
 322 public:
 323     //! Reference data
 324     TestReferenceData refData_;
 325     //! Checker for reference data
 326     TestReferenceChecker checker_;
 327
 328     ConstraintsTest() : checker_(refData_.rootChecker()) {}
 329
 330     /*! \brief Test if the final position correspond to the reference data.
 331      *
 332      * \param[in] testData        Test data structure.
 333      */
 334     void checkFinalPositions(const ConstraintsTestData& testData)
 335     {
 336         TestReferenceChecker finalPositionsRef(
 337                 checker_.checkSequenceCompound("FinalPositions", testData.numAtoms_));
 338         for (int i = 0; i < testData.numAtoms_; i++)
 339         {
 340             TestReferenceChecker xPrimeRef(finalPositionsRef.checkCompound("Atom", nullptr));
 341             const gmx::RVec&     xPrime = testData.xPrime_[i];
 342             xPrimeRef.checkReal(xPrime[XX], "XX");
 343             xPrimeRef.checkReal(xPrime[YY], "YY");
 344             xPrimeRef.checkReal(xPrime[ZZ], "ZZ");
 345         }
 346     }
 347
 348     /*! \brief Test if the final velocities correspond to the reference data.
 349      *
 350      * \param[in] testData        Test data structure.
 351      */
 352     void checkFinalVelocities(const ConstraintsTestData& testData)
 353     {
 354         TestReferenceChecker finalVelocitiesRef(
 355                 checker_.checkSequenceCompound("FinalVelocities", testData.numAtoms_));
 356         for (int i = 0; i < testData.numAtoms_; i++)
 357         {
 358             TestReferenceChecker vRef(finalVelocitiesRef.checkCompound("Atom", nullptr));
 359             const gmx::RVec&     v = testData.v_[i];
 360             vRef.checkReal(v[XX], "XX");
 361             vRef.checkReal(v[YY], "YY");
 362             vRef.checkReal(v[ZZ], "ZZ");
 363         }
 364     }
 365
 366     /*! \brief
 367      * The test on the final length of constrained bonds.
 368      *
 369      * Goes through all the constraints and checks if the final length of all the constraints is
 370      * equal to the target length with provided tolerance.
 371      *
 372      * \param[in] tolerance       Allowed tolerance in final lengths.
 373      * \param[in] testData        Test data structure.
 374      * \param[in] pbc             Periodic boundary data.
 375      */
 376     static void checkConstrainsLength(FloatingPointTolerance     tolerance,
 377                                       const ConstraintsTestData& testData,
 378                                       t_pbc                      pbc)
 379     {
 380
 381         // Test if all the constraints are satisfied
 382         for (index c = 0; c < ssize(testData.constraints_) / 3; c++)
 383         {
 384             real r0 = testData.constraintsR0_.at(testData.constraints_.at(3 * c));
 385             int  i  = testData.constraints_.at(3 * c + 1);
 386             int  j  = testData.constraints_.at(3 * c + 2);
 387             RVec xij0, xij1;
 388             real d0, d1;
 389             if (pbc.pbcType == PbcType::Xyz)
 390             {
 391                 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
 392                 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
 393             }
 394             else
 395             {
 396                 rvec_sub(testData.x_[i], testData.x_[j], xij0);
 397                 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
 398             }
 399             d0 = norm(xij0);
 400             d1 = norm(xij1);
 401             EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
 402                     "rij = %f, which is not equal to r0 = %f for constraint #%zd, between atoms %d "
 403                     "and %d"
 404                     " (before constraining rij was %f).",
 405                     d1,
 406                     r0,
 407                     c,
 408                     i,
 409                     j,
 410                     d0);
 411         }
 412     }
 413
 414     /*! \brief
 415      * The test on the final length of constrained bonds.
 416      *
 417      * Goes through all the constraints and checks if the direction of constraint has not changed
 418      * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
 419      * to the initial system conformation).
 420      *
 421      * \param[in] testData        Test data structure.
 422      * \param[in] pbc             Periodic boundary data.
 423      */
 424     static void checkConstrainsDirection(const ConstraintsTestData& testData, t_pbc pbc)
 425     {
 426
 427         for (index c = 0; c < ssize(testData.constraints_) / 3; c++)
 428         {
 429             int  i = testData.constraints_.at(3 * c + 1);
 430             int  j = testData.constraints_.at(3 * c + 2);
 431             RVec xij0, xij1;
 432             if (pbc.pbcType == PbcType::Xyz)
 433             {
 434                 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
 435                 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
 436             }
 437             else
 438             {
 439                 rvec_sub(testData.x_[i], testData.x_[j], xij0);
 440                 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
 441             }
 442
 443             real dot = xij0.dot(xij1);
 444
 445             EXPECT_GE(dot, 0.0) << gmx::formatString(
 446                     "The constraint %zd changed direction. Constraining algorithm might have "
 447                     "returned the wrong root "
 448                     "of the constraints equation.",
 449                     c);
 450         }
 451     }
 452
 453     /*! \brief
 454      * The test on the coordinates of the center of the mass (COM) of the system.
 455      *
 456      * Checks if the center of mass has not been shifted by the constraints. Note,
 457      * that this test does not take into account the periodic boundary conditions.
 458      * Hence it will not work should the constraints decide to move atoms across
 459      * PBC borders.
 460      *
 461      * \param[in] tolerance       Allowed tolerance in COM coordinates.
 462      * \param[in] testData        Test data structure.
 463      */
 464     static void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData& testData)
 465     {
 466
 467         RVec comPrime0({ 0.0, 0.0, 0.0 });
 468         RVec comPrime({ 0.0, 0.0, 0.0 });
 469         for (int i = 0; i < testData.numAtoms_; i++)
 470         {
 471             comPrime0 += testData.masses_[i] * testData.xPrime0_[i];
 472             comPrime += testData.masses_[i] * testData.xPrime_[i];
 473         }
 474
 475         comPrime0 /= testData.numAtoms_;
 476         comPrime /= testData.numAtoms_;
 477
 478         EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
 479                 << "Center of mass was shifted by constraints in x-direction.";
 480         EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
 481                 << "Center of mass was shifted by constraints in y-direction.";
 482         EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
 483                 << "Center of mass was shifted by constraints in z-direction.";
 484     }
 485
 486     /*! \brief
 487      * The test on the velocity of the center of the mass (COM) of the system.
 488      *
 489      * Checks if the velocity of the center of mass has not changed.
 490      *
 491      * \param[in] tolerance       Allowed tolerance in COM velocity components.
 492      * \param[in] testData        Test data structure.
 493      */
 494     static void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData& testData)
 495     {
 496
 497         RVec comV0({ 0.0, 0.0, 0.0 });
 498         RVec comV({ 0.0, 0.0, 0.0 });
 499         for (int i = 0; i < testData.numAtoms_; i++)
 500         {
 501             comV0 += testData.masses_[i] * testData.v0_[i];
 502             comV += testData.masses_[i] * testData.v_[i];
 503         }
 504         comV0 /= testData.numAtoms_;
 505         comV /= testData.numAtoms_;
 506
 507         EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
 508                 << "Velocity of the center of mass in x-direction has been changed by constraints.";
 509         EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
 510                 << "Velocity of the center of mass in y-direction has been changed by constraints.";
 511         EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
 512                 << "Velocity of the center of mass in z-direction has been changed by constraints.";
 513     }
 514
 515     /*! \brief
 516      * The test of virial tensor.
 517      *
 518      * Checks if the values in the scaled virial tensor are equal to pre-computed values.
 519      *
 520      * \param[in] testData        Test data structure.
 521      */
 522     void checkVirialTensor(const ConstraintsTestData& testData)
 523     {
 524         const tensor&        virialScaled = testData.virialScaled_;
 525         TestReferenceChecker virialScaledRef(checker_.checkCompound("VirialScaled", nullptr));
 526
 527         virialScaledRef.checkReal(virialScaled[XX][XX], "XX");
 528         virialScaledRef.checkReal(virialScaled[XX][YY], "XY");
 529         virialScaledRef.checkReal(virialScaled[XX][ZZ], "XZ");
 530         virialScaledRef.checkReal(virialScaled[YY][XX], "YX");
 531         virialScaledRef.checkReal(virialScaled[YY][YY], "YY");
 532         virialScaledRef.checkReal(virialScaled[YY][ZZ], "YZ");
 533         virialScaledRef.checkReal(virialScaled[ZZ][XX], "ZX");
 534         virialScaledRef.checkReal(virialScaled[ZZ][YY], "ZY");
 535         virialScaledRef.checkReal(virialScaled[ZZ][ZZ], "ZZ");
 536     }
 537
 538     //! Before any test is run, work out whether any compatible GPUs exist.
 539     static std::vector<std::unique_ptr<IConstraintsTestRunner>> getRunners()
 540     {
 541         std::vector<std::unique_ptr<IConstraintsTestRunner>> runners;
 542         // Add runners for CPU versions of SHAKE and LINCS
 543         // runners.emplace_back(std::make_unique<ShakeConstraintsRunner>());
 544         runners.emplace_back(std::make_unique<LincsConstraintsRunner>());
 545         // If supported, add runners for the GPU version of LINCS for each available GPU
 546         const bool addGpuRunners = GPU_CONSTRAINTS_SUPPORTED;
 547         if (addGpuRunners)
 548         {
 549             for (const auto& testDevice : getTestHardwareEnvironment()->getTestDeviceList())
 550             {
 551                 runners.emplace_back(std::make_unique<LincsDeviceConstraintsRunner>(*testDevice));
 552             }
 553         }
 554         return runners;
 555     }
 556 };
 557
 558 TEST_P(ConstraintsTest, SatisfiesConstraints)
 559 {
 560     auto                  params                = GetParam();
 561     ConstraintsTestSystem constraintsTestSystem = std::get<0>(params);
 562     t_pbc                 pbc                   = std::get<1>(params);
 563
 564     ConstraintsTestData testData(constraintsTestSystem.title,
 565                                  constraintsTestSystem.numAtoms,
 566                                  constraintsTestSystem.masses,
 567                                  constraintsTestSystem.constraints,
 568                                  constraintsTestSystem.constraintsR0,
 569                                  true,
 570                                  false,
 571                                  0,
 572                                  real(0.0),
 573                                  real(0.001),
 574                                  constraintsTestSystem.x,
 575                                  constraintsTestSystem.xPrime,
 576                                  constraintsTestSystem.v,
 577                                  constraintsTestSystem.shakeTolerance,
 578                                  constraintsTestSystem.shakeUseSOR,
 579                                  constraintsTestSystem.lincsNIter,
 580                                  constraintsTestSystem.lincslincsExpansionOrder,
 581                                  constraintsTestSystem.lincsWarnAngle);
 582
 583     float maxX = 0.0F;
 584     float maxV = 0.0F;
 585     for (int i = 0; i < constraintsTestSystem.numAtoms; i++)
 586     {
 587         for (int d = 0; d < DIM; d++)
 588         {
 589             maxX = fmax(fabs(constraintsTestSystem.x[i][d]), maxX);
 590             maxV = fmax(fabs(constraintsTestSystem.v[i][d]), maxV);
 591         }
 592     }
 593
 594     float maxVirialScaled = 0.0F;
 595     for (int d1 = 0; d1 < DIM; d1++)
 596     {
 597         for (int d2 = 0; d2 < DIM; d2++)
 598         {
 599             maxVirialScaled = fmax(fabs(testData.virialScaled_[d1][d2]), maxVirialScaled);
 600         }
 601     }
 602
 603     FloatingPointTolerance positionsTolerance = relativeToleranceAsFloatingPoint(maxX, 0.002F);
 604     FloatingPointTolerance velocityTolerance  = relativeToleranceAsFloatingPoint(maxV, 0.002F);
 605     FloatingPointTolerance virialTolerance = relativeToleranceAsFloatingPoint(maxVirialScaled, 0.002F);
 606     FloatingPointTolerance lengthTolerance = relativeToleranceAsFloatingPoint(0.1, 0.002F);
 607
 608     // Cycle through all available runners
 609     for (const auto& runner : getRunners())
 610     {
 611         SCOPED_TRACE(formatString("Testing %s with %s PBC using %s.",
 612                                   testData.title_.c_str(),
 613                                   c_pbcTypeNames[pbc.pbcType].c_str(),
 614                                   runner->name().c_str()));
 615
 616         testData.reset();
 617
 618         // Apply constraints
 619         runner->applyConstraints(&testData, pbc);
 620
 621
 622         checker_.setDefaultTolerance(positionsTolerance);
 623         checkFinalPositions(testData);
 624         checker_.setDefaultTolerance(velocityTolerance);
 625         checkFinalVelocities(testData);
 626
 627         checkConstrainsLength(lengthTolerance, testData, pbc);
 628         checkConstrainsDirection(testData, pbc);
 629         checkCOMCoordinates(positionsTolerance, testData);
 630         checkCOMVelocity(velocityTolerance, testData);
 631
 632         checker_.setDefaultTolerance(virialTolerance);
 633         checkVirialTensor(testData);
 634     }
 635 }
 636
 637 INSTANTIATE_TEST_SUITE_P(WithParameters,
 638                          ConstraintsTest,
 639                          ::testing::Combine(::testing::ValuesIn(c_constraintsTestSystemList),
 640                                             ::testing::ValuesIn(c_pbcs)));
 641
 642 } // namespace
 643 } // namespace test
 644 } // namespace gmx