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