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