2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2018,2019, 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.
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.
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.
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.
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.
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.
36 * \brief SHAKE and LINCS tests.
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
44 * \author Artem Zhmurov <zhmurov@gmail.com>
45 * \ingroup module_mdlib
50 #include "gromacs/mdlib/constr.h"
59 #include <unordered_map>
62 #include <gtest/gtest.h>
64 #include "gromacs/fileio/gmxfio.h"
65 #include "gromacs/gmxlib/nrnb.h"
66 #include "gromacs/math/paddedvector.h"
67 #include "gromacs/math/vec.h"
68 #include "gromacs/math/vectypes.h"
69 #include "gromacs/mdlib/lincs.h"
70 #include "gromacs/mdlib/shake.h"
71 #include "gromacs/mdrunutility/multisim.h"
72 #include "gromacs/mdtypes/commrec.h"
73 #include "gromacs/mdtypes/inputrec.h"
74 #include "gromacs/mdtypes/mdatom.h"
75 #include "gromacs/pbcutil/pbc.h"
76 #include "gromacs/topology/block.h"
77 #include "gromacs/topology/idef.h"
78 #include "gromacs/topology/ifunc.h"
79 #include "gromacs/topology/topology.h"
80 #include "gromacs/utility/smalloc.h"
81 #include "gromacs/utility/stringutil.h"
82 #include "gromacs/utility/unique_cptr.h"
84 #include "testutils/refdata.h"
85 #include "testutils/testasserts.h"
87 #include "constr_impl.h"
94 ConstraintsTestData::ConstraintsTestData(const std::string &title,
95 int numAtoms, std::vector<real> masses,
96 std::vector<int> constraints, std::vector<real> constraintsR0,
97 bool computeVirial, tensor virialScaledRef,
98 bool compute_dHdLambda, float dHdLambdaRef,
99 real initialTime, real timestep,
100 const std::vector<RVec> &x, const std::vector<RVec> &xPrime, const std::vector<RVec> &v,
101 real shakeTolerance, gmx_bool shakeUseSOR,
102 int lincsNumIterations, int lincsExpansionOrder, real lincsWarnAngle)
104 // This is to trick Gerrit
107 title_ = title; // Human-friendly name of the system
108 numAtoms_ = numAtoms; // Number of atoms
112 invmass_.resize(numAtoms); // Vector of inverse masses
114 for (int i = 0; i < numAtoms; i++)
116 invmass_[i] = 1.0/masses.at(i);
119 // Saving constraints to check if they are satisfied after algorithm was applied
120 constraints_ = constraints; // Constraints indices (in type-i-j format)
121 constraintsR0_ = constraintsR0; // Equilibrium distances for each type of constraint
123 invdt_ = 1.0/timestep; // Inverse timestep
125 // Communication record
133 // Input record - data that usually comes from configuration file (.mdp)
135 ir_.init_t = initialTime;
136 ir_.delta_t = timestep;
140 md_.nMassPerturbed = 0;
142 md_.invmass = invmass_.data();
144 md_.homenr = numAtoms;
147 computeVirial_ = computeVirial;
150 for (int i = 0; i < DIM; i++)
152 for (int j = 0; j < DIM; j++)
154 virialScaled_[i][j] = 0;
155 virialScaledRef_[i][j] = virialScaledRef[i][j];
161 // Free energy evaluation
162 compute_dHdLambda_ = compute_dHdLambda;
164 if (compute_dHdLambda_)
167 dHdLambdaRef_ = dHdLambdaRef;
175 // Constraints and their parameters (local topology)
176 for (int i = 0; i < F_NRE; i++)
180 idef_.il[F_CONSTR].nr = constraints.size();
182 snew(idef_.il[F_CONSTR].iatoms, constraints.size());
184 for (unsigned i = 0; i < constraints.size(); i++)
188 if (maxType < constraints.at(i))
190 maxType = constraints.at(i);
193 idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
195 snew(idef_.iparams, maxType + 1);
196 for (unsigned i = 0; i < constraints.size()/3; i++)
198 idef_.iparams[constraints.at(3*i)].constr.dA = constraintsR0.at(constraints.at(3*i));
199 idef_.iparams[constraints.at(3*i)].constr.dB = constraintsR0.at(constraints.at(3*i));
202 // Constraints and their parameters (global topology)
203 InteractionList interactionList;
204 interactionList.iatoms.resize(constraints.size());
205 for (unsigned i = 0; i < constraints.size(); i++)
207 interactionList.iatoms.at(i) = constraints.at(i);
209 InteractionList interactionListEmpty;
210 interactionListEmpty.iatoms.resize(0);
212 gmx_moltype_t molType;
213 molType.atoms.nr = numAtoms;
214 molType.ilist.at(F_CONSTR) = interactionList;
215 molType.ilist.at(F_CONSTRNC) = interactionListEmpty;
216 mtop_.moltype.push_back(molType);
218 gmx_molblock_t molBlock;
221 mtop_.molblock.push_back(molBlock);
223 mtop_.natoms = numAtoms;
224 mtop_.ffparams.iparams.resize(maxType + 1);
225 for (int i = 0; i <= maxType; i++)
227 mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
229 mtop_.bIntermolecularInteractions = false;
231 // Coordinates and velocities
232 x_.resizeWithPadding(numAtoms);
233 xPrime_.resizeWithPadding(numAtoms);
234 xPrime0_.resizeWithPadding(numAtoms);
235 xPrime2_.resizeWithPadding(numAtoms);
237 v_.resizeWithPadding(numAtoms);
238 v0_.resizeWithPadding(numAtoms);
240 std::copy(x.begin(), x.end(), x_.begin());
241 std::copy(xPrime.begin(), xPrime.end(), xPrime_.begin());
242 std::copy(xPrime.begin(), xPrime.end(), xPrime0_.begin());
243 std::copy(xPrime.begin(), xPrime.end(), xPrime2_.begin());
245 std::copy(v.begin(), v.end(), v_.begin());
246 std::copy(v.begin(), v.end(), v0_.begin());
248 // SHAKE-specific parameters
249 ir_.shake_tol = shakeTolerance;
250 ir_.bShakeSOR = shakeUseSOR;
252 // LINCS-specific parameters
253 ir_.nLincsIter = lincsNumIterations;
254 ir_.nProjOrder = lincsExpansionOrder;
255 ir_.LincsWarnAngle = lincsWarnAngle;
261 * Reset the data structure so it can be reused.
263 * Set the coordinates and velocities back to their values before
264 * constraining. The scaled virial tensor and dHdLambda are zeroed.
267 void ConstraintsTestData::reset()
275 for (int i = 0; i < DIM; i++)
277 for (int j = 0; j < DIM; j++)
279 virialScaled_[i][j] = 0;
287 * Cleaning up the memory.
289 ConstraintsTestData::~ConstraintsTestData()
291 sfree(idef_.il[F_CONSTR].iatoms);
292 sfree(idef_.iparams);
298 /*! \brief The two-dimensional parameter space for test.
300 * The test will run for all possible combinations of accessible
302 * 1. PBC setup ("PBCNONE" or "PBCXYZ")
303 * 2. The algorithm ("SHAKE", "LINCS" or "LINCS_GPU").
305 typedef std::tuple<std::string, std::string> ConstraintsTestParameters;
307 //! Names of all availible algorithms
308 std::vector<std::string> algorithmsNames;
310 //! Method that fills and returns algorithmNames to the test macros.
311 std::vector<std::string> getAlgorithmsNames()
313 algorithmsNames.emplace_back("SHAKE");
314 algorithmsNames.emplace_back("LINCS");
315 std::string errorMessage;
316 if (GMX_GPU == GMX_GPU_CUDA && canDetectGpus(&errorMessage))
318 algorithmsNames.emplace_back("LINCS_CUDA");
320 return algorithmsNames;
323 /*! \brief Test fixture for constraints.
325 * The fixture uses following test systems:
326 * 1. Two atoms, connected with one constraint (e.g. NH).
327 * 2. Three atoms, connected consequently with two constraints (e.g. CH2).
328 * 3. Three atoms, constrained to the fourth atom (e.g. CH3).
329 * 4. Four atoms, connected by two independent constraints.
330 * 5. Three atoms, connected by three constraints in a triangle
331 * (e.g. H2O with constrained H-O-H angle).
332 * 6. Four atoms, connected by three consequential constraints.
334 * For all systems, the final lengths of the constraints are tested against the
335 * reference values, the direction of each constraint is checked.
336 * Test also verifies that the center of mass has not been
337 * shifted by the constraints and that its velocity has not changed.
338 * For some systems, the value for scaled virial tensor is checked against
341 class ConstraintsTest : public ::testing::TestWithParam<ConstraintsTestParameters>
345 std::unordered_map <std::string, t_pbc> pbcs_;
346 //! Algorithms (SHAKE and LINCS)
347 std::unordered_map <std::string, void(*)(ConstraintsTestData *testData, t_pbc pbc)> algorithms_;
349 /*! \brief Test setup function.
351 * Setting up the pbcs and algorithms. Note, that corresponding string keywords
352 * have to be explicitly added at the end of this file when the tests are called.
355 void SetUp() override
359 // PBC initialization
363 // Infinitely small box
364 matrix boxNone = { {0, 0, 0}, {0, 0, 0}, {0, 0, 0} };
365 set_pbc(&pbc, epbcNONE, boxNone);
366 pbcs_["PBCNone"] = pbc;
369 matrix boxXyz = { {10.0, 0.0, 0.0}, {0.0, 20.0, 0.0}, {0.0, 0.0, 15.0} };
370 set_pbc(&pbc, epbcXYZ, boxXyz);
371 pbcs_["PBCXYZ"] = pbc;
377 algorithms_["SHAKE"] = applyShake;
379 algorithms_["LINCS"] = applyLincs;
380 // LINCS using CUDA (will only be called if CUDA is available)
381 algorithms_["LINCS_CUDA"] = applyLincsCuda;
385 * The test on the final length of constrained bonds.
387 * Goes through all the constraints and checks if the final length of all the constraints is equal
388 * to the target length with provided tolerance.
390 * \param[in] tolerance Allowed tolerance in final lengths.
391 * \param[in] testData Test data structure.
392 * \param[in] pbc Periodic boundary data.
394 void checkConstrainsLength(FloatingPointTolerance tolerance, const ConstraintsTestData &testData, t_pbc pbc)
397 // Test if all the constraints are satisfied
398 for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
400 real r0 = testData.constraintsR0_.at(testData.constraints_.at(3*c));
401 int i = testData.constraints_.at(3*c + 1);
402 int j = testData.constraints_.at(3*c + 2);
405 if (pbc.ePBC == epbcXYZ)
407 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
408 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
412 rvec_sub(testData.x_[i], testData.x_[j], xij0);
413 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
417 EXPECT_REAL_EQ_TOL(r0, d1, tolerance) << gmx::formatString(
418 "rij = %f, which is not equal to r0 = %f for constraint #%u, between atoms %d and %d"
419 " (before constraining rij was %f).", d1, r0, c, i, j, d0);
424 * The test on the final length of constrained bonds.
426 * Goes through all the constraints and checks if the direction of constraint has not changed
427 * by the algorithm (i.e. the constraints algorithm arrived to the solution that is closest
428 * to the initial system conformation).
430 * \param[in] testData Test data structure.
431 * \param[in] pbc Periodic boundary data.
433 void checkConstrainsDirection(const ConstraintsTestData &testData, t_pbc pbc)
436 for (unsigned c = 0; c < testData.constraints_.size()/3; c++)
438 int i = testData.constraints_.at(3*c + 1);
439 int j = testData.constraints_.at(3*c + 2);
441 if (pbc.ePBC == epbcXYZ)
443 pbc_dx_aiuc(&pbc, testData.x_[i], testData.x_[j], xij0);
444 pbc_dx_aiuc(&pbc, testData.xPrime_[i], testData.xPrime_[j], xij1);
448 rvec_sub(testData.x_[i], testData.x_[j], xij0);
449 rvec_sub(testData.xPrime_[i], testData.xPrime_[j], xij1);
452 real dot = xij0.dot(xij1);
454 EXPECT_GE(dot, 0.0) << gmx::formatString(
455 "The constraint %u changed direction. Constraining algorithm might have returned the wrong root "
456 "of the constraints equation.", c);
462 * The test on the coordinates of the center of the mass (COM) of the system.
464 * Checks if the center of mass has not been shifted by the constraints. Note,
465 * that this test does not take into account the periodic boundary conditions.
466 * Hence it will not work should the constraints decide to move atoms across
469 * \param[in] tolerance Allowed tolerance in COM coordinates.
470 * \param[in] testData Test data structure.
472 void checkCOMCoordinates(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
475 RVec comPrime0({0.0, 0.0, 0.0});
476 RVec comPrime({0.0, 0.0, 0.0});
477 for (int i = 0; i < testData.numAtoms_; i++)
479 comPrime0 += testData.masses_[i]*testData.xPrime0_[i];
480 comPrime += testData.masses_[i]*testData.xPrime_[i];
483 comPrime0 /= testData.numAtoms_;
484 comPrime /= testData.numAtoms_;
486 EXPECT_REAL_EQ_TOL(comPrime[XX], comPrime0[XX], tolerance)
487 << "Center of mass was shifted by constraints in x-direction.";
488 EXPECT_REAL_EQ_TOL(comPrime[YY], comPrime0[YY], tolerance)
489 << "Center of mass was shifted by constraints in y-direction.";
490 EXPECT_REAL_EQ_TOL(comPrime[ZZ], comPrime0[ZZ], tolerance)
491 << "Center of mass was shifted by constraints in z-direction.";
496 * The test on the velocity of the center of the mass (COM) of the system.
498 * Checks if the velocity of the center of mass has not changed.
500 * \param[in] tolerance Allowed tolerance in COM velocity components.
501 * \param[in] testData Test data structure.
503 void checkCOMVelocity(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
506 RVec comV0({0.0, 0.0, 0.0});
507 RVec comV({0.0, 0.0, 0.0});
508 for (int i = 0; i < testData.numAtoms_; i++)
510 comV0 += testData.masses_[i]*testData.v0_[i];
511 comV += testData.masses_[i]*testData.v_[i];
513 comV0 /= testData.numAtoms_;
514 comV /= testData.numAtoms_;
516 EXPECT_REAL_EQ_TOL(comV[XX], comV0[XX], tolerance)
517 << "Velocity of the center of mass in x-direction has been changed by constraints.";
518 EXPECT_REAL_EQ_TOL(comV[YY], comV0[YY], tolerance)
519 << "Velocity of the center of mass in y-direction has been changed by constraints.";
520 EXPECT_REAL_EQ_TOL(comV[ZZ], comV0[ZZ], tolerance)
521 << "Velocity of the center of mass in z-direction has been changed by constraints.";
525 * The test of virial tensor.
527 * Checks if the values in the scaled virial tensor are equal to pre-computed values.
529 * \param[in] tolerance Tolerance for the tensor values.
530 * \param[in] testData Test data structure.
532 void checkVirialTensor(FloatingPointTolerance tolerance, const ConstraintsTestData &testData)
534 for (int i = 0; i < DIM; i++)
536 for (int j = 0; j < DIM; j++)
538 EXPECT_REAL_EQ_TOL(testData.virialScaledRef_[i][j], testData.virialScaled_[i][j],
539 tolerance) << gmx::formatString(
540 "Values in virial tensor at [%d][%d] are not within the tolerance from reference value.", i, j);
546 TEST_P(ConstraintsTest, SingleConstraint){
547 std::string title = "one constraint (e.g. OH)";
550 std::vector<real> masses = {1.0, 12.0};
551 std::vector<int> constraints = {0, 0, 1};
552 std::vector<real> constraintsR0 = {0.1};
554 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
556 std::vector<RVec> x = {{ 0.0, oneTenthOverSqrtTwo, 0.0 },
557 { oneTenthOverSqrtTwo, 0.0, 0.0 }};
558 std::vector<RVec> xPrime = {{ 0.01, 0.08, 0.01 },
559 { 0.06, 0.01, -0.01 }};
560 std::vector<RVec> v = {{ 1.0, 2.0, 3.0 },
563 tensor virialScaledRef = {{-5.58e-04, 5.58e-04, 0.00e+00 },
564 { 5.58e-04, -5.58e-04, 0.00e+00 },
565 { 0.00e+00, 0.00e+00, 0.00e+00 }};
567 real shakeTolerance = 0.0001;
568 gmx_bool shakeUseSOR = false;
571 int lincslincsExpansionOrder = 4;
572 real lincsWarnAngle = 30.0;
574 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
575 (title, numAtoms, masses,
576 constraints, constraintsR0,
577 true, virialScaledRef,
579 real(0.0), real(0.001),
581 shakeTolerance, shakeUseSOR,
582 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
584 std::string algorithmName;
585 std::tie(pbcName, algorithmName) = GetParam();
586 t_pbc pbc = pbcs_.at(pbcName);
589 algorithms_.at(algorithmName)(testData.get(), pbc);
591 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
592 checkConstrainsDirection(*testData, pbc);
593 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
594 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
596 checkVirialTensor(absoluteTolerance(0.0001), *testData);
600 TEST_P(ConstraintsTest, TwoDisjointConstraints){
602 std::string title = "two disjoint constraints";
604 std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
605 std::vector<int> constraints = {0, 0, 1, 1, 2, 3};
606 std::vector<real> constraintsR0 = {2.0, 1.0};
609 std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
610 { 0.51, -3.02, 15.55 },
611 { -0.50, -3.00, 15.20 },
612 { -1.51, -2.95, 15.05 }};
614 std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
615 { 0.51, -3.02, 15.55 },
616 { -0.50, -3.00, 15.20 },
617 { -1.51, -2.95, 15.05 }};
619 std::vector<RVec> v = {{ 0.0, 1.0, 0.0 },
624 tensor virialScaledRef = {{ 3.3e-03, -1.7e-04, 5.6e-04 },
625 {-1.7e-04, 8.9e-06, -2.8e-05 },
626 { 5.6e-04, -2.8e-05, 8.9e-05 }};
628 real shakeTolerance = 0.0001;
629 gmx_bool shakeUseSOR = false;
632 int lincslincsExpansionOrder = 4;
633 real lincsWarnAngle = 30.0;
635 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
636 (title, numAtoms, masses,
637 constraints, constraintsR0,
638 true, virialScaledRef,
640 real(0.0), real(0.001),
642 shakeTolerance, shakeUseSOR,
643 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
646 std::string algorithmName;
647 std::tie(pbcName, algorithmName) = GetParam();
648 t_pbc pbc = pbcs_.at(pbcName);
651 algorithms_.at(algorithmName)(testData.get(), pbc);
653 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
654 checkConstrainsDirection(*testData, pbc);
655 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
656 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
658 checkVirialTensor(absoluteTolerance(0.0001), *testData);
662 TEST_P(ConstraintsTest, ThreeSequentialConstraints){
664 std::string title = "three atoms, connected longitudinally (e.g. CH2)";
666 std::vector<real> masses = {1.0, 12.0, 16.0 };
667 std::vector<int> constraints = {0, 0, 1, 1, 1, 2};
668 std::vector<real> constraintsR0 = {0.1, 0.2};
670 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
671 real twoTenthsOverSqrtThree = 0.2_real / std::sqrt(3.0_real);
673 std::vector<RVec> x = {{ oneTenthOverSqrtTwo, oneTenthOverSqrtTwo, 0.0 },
675 { twoTenthsOverSqrtThree, twoTenthsOverSqrtThree, twoTenthsOverSqrtThree }};
677 std::vector<RVec> xPrime = {{ 0.08, 0.07, 0.01 },
678 { -0.02, 0.01, -0.02 },
679 { 0.10, 0.12, 0.11 }};
681 std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
685 tensor virialScaledRef = {{ 4.14e-03, 4.14e-03, 3.31e-03},
686 { 4.14e-03, 4.14e-03, 3.31e-03},
687 { 3.31e-03, 3.31e-03, 3.31e-03}};
689 real shakeTolerance = 0.0001;
690 gmx_bool shakeUseSOR = false;
693 int lincslincsExpansionOrder = 4;
694 real lincsWarnAngle = 30.0;
696 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
697 (title, numAtoms, masses,
698 constraints, constraintsR0,
699 true, virialScaledRef,
701 real(0.0), real(0.001),
703 shakeTolerance, shakeUseSOR,
704 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
707 std::string algorithmName;
708 std::tie(pbcName, algorithmName) = GetParam();
709 t_pbc pbc = pbcs_.at(pbcName);
712 algorithms_.at(algorithmName)(testData.get(), pbc);
714 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
715 checkConstrainsDirection(*testData, pbc);
716 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
717 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
719 checkVirialTensor(absoluteTolerance(0.0001), *testData);
723 TEST_P(ConstraintsTest, ThreeConstraintsWithCentralAtom){
725 std::string title = "three atoms, connected to the central atom (e.g. CH3)";
727 std::vector<real> masses = {12.0, 1.0, 1.0, 1.0 };
728 std::vector<int> constraints = {0, 0, 1, 0, 0, 2, 0, 0, 3};
729 std::vector<real> constraintsR0 = {0.1};
732 std::vector<RVec> x = {{ 0.00, 0.00, 0.00 },
733 { 0.10, 0.00, 0.00 },
734 { 0.00, -0.10, 0.00 },
735 { 0.00, 0.00, 0.10 }};
737 std::vector<RVec> xPrime = {{ 0.004, 0.009, -0.010 },
738 { 0.110, -0.006, 0.003 },
739 {-0.007, -0.102, -0.007 },
740 {-0.005, 0.011, 0.102 }};
742 std::vector<RVec> v = {{ 1.0, 0.0, 0.0 },
747 tensor virialScaledRef = {{7.14e-04, 0.00e+00, 0.00e+00},
748 {0.00e+00, 1.08e-03, 0.00e+00},
749 {0.00e+00, 0.00e+00, 1.15e-03}};
751 real shakeTolerance = 0.0001;
752 gmx_bool shakeUseSOR = false;
755 int lincslincsExpansionOrder = 4;
756 real lincsWarnAngle = 30.0;
758 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
759 (title, numAtoms, masses,
760 constraints, constraintsR0,
761 true, virialScaledRef,
763 real(0.0), real(0.001),
765 shakeTolerance, shakeUseSOR,
766 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
769 std::string algorithmName;
770 std::tie(pbcName, algorithmName) = GetParam();
771 t_pbc pbc = pbcs_.at(pbcName);
774 algorithms_.at(algorithmName)(testData.get(), pbc);
776 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
777 checkConstrainsDirection(*testData, pbc);
778 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
779 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
781 checkVirialTensor(absoluteTolerance(0.0001), *testData);
784 TEST_P(ConstraintsTest, FourSequentialConstraints){
786 std::string title = "four atoms, connected longitudinally";
788 std::vector<real> masses = {0.5, 1.0/3.0, 0.25, 1.0};
789 std::vector<int> constraints = {0, 0, 1, 1, 1, 2, 2, 2, 3};
790 std::vector<real> constraintsR0 = {2.0, 1.0, 1.0};
793 std::vector<RVec> x = {{ 2.50, -3.10, 15.70 },
794 { 0.51, -3.02, 15.55 },
795 { -0.50, -3.00, 15.20 },
796 { -1.51, -2.95, 15.05 }};
798 std::vector<RVec> xPrime = {{ 2.50, -3.10, 15.70 },
799 { 0.51, -3.02, 15.55 },
800 { -0.50, -3.00, 15.20 },
801 { -1.51, -2.95, 15.05 }};
803 std::vector<RVec> v = {{ 0.0, 0.0, 2.0 },
808 tensor virialScaledRef = {{ 1.15e-01, -4.20e-03, 2.12e-02},
809 {-4.20e-03, 1.70e-04, -6.41e-04},
810 { 2.12e-02, -6.41e-04, 5.45e-03}};
812 real shakeTolerance = 0.0001;
813 gmx_bool shakeUseSOR = false;
816 int lincslincsExpansionOrder = 8;
817 real lincsWarnAngle = 30.0;
819 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
820 (title, numAtoms, masses,
821 constraints, constraintsR0,
822 true, virialScaledRef,
824 real(0.0), real(0.001),
826 shakeTolerance, shakeUseSOR,
827 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
830 std::string algorithmName;
831 std::tie(pbcName, algorithmName) = GetParam();
832 t_pbc pbc = pbcs_.at(pbcName);
835 algorithms_.at(algorithmName)(testData.get(), pbc);
837 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
838 checkConstrainsDirection(*testData, pbc);
839 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
840 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
842 checkVirialTensor(absoluteTolerance(0.01), *testData);
846 TEST_P(ConstraintsTest, TriangleOfConstraints){
848 std::string title = "basic triangle (tree atoms, connected to each other)";
850 std::vector<real> masses = {1.0, 1.0, 1.0};
851 std::vector<int> constraints = {0, 0, 1, 2, 0, 2, 1, 1, 2};
852 std::vector<real> constraintsR0 = {0.1, 0.1, 0.1};
854 real oneTenthOverSqrtTwo = 0.1_real / std::sqrt(2.0_real);
856 std::vector<RVec> x = {{ oneTenthOverSqrtTwo, 0.0, 0.0 },
857 { 0.0, oneTenthOverSqrtTwo, 0.0 },
858 { 0.0, 0.0, oneTenthOverSqrtTwo }};
860 std::vector<RVec> xPrime = {{ 0.09, -0.02, 0.01 },
861 { -0.02, 0.10, -0.02 },
862 { 0.03, -0.01, 0.07 }};
864 std::vector<RVec> v = {{ 1.0, 1.0, 1.0 },
865 { -2.0, -2.0, -2.0 },
868 tensor virialScaledRef = {{ 6.00e-04, -1.61e-03, 1.01e-03},
869 {-1.61e-03, 2.53e-03, -9.25e-04},
870 { 1.01e-03, -9.25e-04, -8.05e-05}};
872 real shakeTolerance = 0.0001;
873 gmx_bool shakeUseSOR = false;
876 int lincslincsExpansionOrder = 4;
877 real lincsWarnAngle = 30.0;
879 std::unique_ptr<ConstraintsTestData> testData = std::make_unique<ConstraintsTestData>
880 (title, numAtoms, masses,
881 constraints, constraintsR0,
882 true, virialScaledRef,
884 real(0.0), real(0.001),
886 shakeTolerance, shakeUseSOR,
887 lincsNIter, lincslincsExpansionOrder, lincsWarnAngle);
890 std::string algorithmName;
891 std::tie(pbcName, algorithmName) = GetParam();
892 t_pbc pbc = pbcs_.at(pbcName);
895 algorithms_.at(algorithmName)(testData.get(), pbc);
897 checkConstrainsLength(absoluteTolerance(0.0002), *testData, pbc);
898 checkConstrainsDirection(*testData, pbc);
899 checkCOMCoordinates(absoluteTolerance(0.0001), *testData);
900 checkCOMVelocity(absoluteTolerance(0.0001), *testData);
902 checkVirialTensor(absoluteTolerance(0.00001), *testData);
907 INSTANTIATE_TEST_CASE_P(WithParameters, ConstraintsTest,
908 ::testing::Combine(::testing::Values("PBCNone", "PBCXYZ"),
909 ::testing::ValuesIn(getAlgorithmsNames())));