2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,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.
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 Tests for SETTLE constraints
38 * The test runs on several small systems, containing 1 to 17 water molecules,
39 * with and without periodic boundary conditions, with and without velocity
40 * and virial updates. The CPU and GPU versions are tested, if the code was
41 * compiled with CUDA support and there is a CUDA-capable GPU in the system.
44 * 1. If the final distances between constrained atoms are within tolerance
45 * from the target distance.
46 * 2. If the velocities were updated when needed.
47 * 3. If the virial was computed.
49 * The test also compares the results from the CPU and GPU versions of the
50 * algorithm: final coordinates, velocities and virial should be within
51 * tolerance to one another.
53 * \todo This also tests that if the calling code requires velocities
54 * and virial updates, that those outputs do change, but does not
55 * test that those changes are correct.
57 * \todo Only no-PBC and cubic-PBC are tested here, but the correct
58 * function of the SIMD version of set_pbx_auic in all cases
59 * should be tested elsewhere.
61 * \todo The CPU and GPU versions are tested against each other. This
62 * should be changed to a proper test against pre-computed
63 * reference values. Also, these test will dry-run on a CUDA
64 * build if no CUDA-capable GPU is available.
66 * \author Mark Abraham <mark.j.abraham@gmail.com>
67 * \author Artem Zhmurov <zhmurov@gmail.com>
69 * \ingroup module_mdlib
73 #include "gromacs/mdlib/settle.h"
78 #include <unordered_map>
81 #include <gtest/gtest.h>
83 #include "gromacs/hardware/device_management.h"
84 #include "gromacs/math/vec.h"
85 #include "gromacs/math/vectypes.h"
86 #include "gromacs/mdtypes/mdatom.h"
87 #include "gromacs/pbcutil/pbc.h"
88 #include "gromacs/topology/idef.h"
89 #include "gromacs/topology/ifunc.h"
90 #include "gromacs/topology/topology.h"
91 #include "gromacs/utility/stringutil.h"
92 #include "gromacs/utility/unique_cptr.h"
94 #include "gromacs/mdlib/tests/watersystem.h"
95 #include "testutils/refdata.h"
96 #include "testutils/test_hardware_environment.h"
97 #include "testutils/testasserts.h"
99 #include "settletestdata.h"
100 #include "settletestrunners.h"
109 /*! \brief Parameters that will vary from test to test.
111 struct SettleTestParameters
113 //! Number of water molecules (SETTLEs) [1, 2, 4, 5, 7, 10, 12, 15, 17]
115 //! If the velocities should be updated while constraining [true/false]
116 bool updateVelocities;
117 //! If the virial should be computed [true/false]
119 //! Periodic boundary conditions [PBCXYZ/PBCNone]
123 /*! \brief Sets of parameters on which to run the tests.
125 const SettleTestParameters parametersSets[] = {
126 { 1, false, false, "PBCXYZ" }, // 1 water molecule
127 { 2, false, false, "PBCXYZ" }, // 2 water molecules
128 { 4, false, false, "PBCXYZ" }, // 4 water molecules
129 { 5, false, false, "PBCXYZ" }, // 5 water molecules
130 { 6, false, false, "PBCXYZ" }, // 6 water molecules
131 { 10, false, false, "PBCXYZ" }, // 10 water molecules
132 { 12, false, false, "PBCXYZ" }, // 12 water molecules
133 { 15, false, false, "PBCXYZ" }, // 15 water molecules
134 { 17, true, false, "PBCXYZ" }, // Update velocities
135 { 17, false, true, "PBCXYZ" }, // Compute virial
136 { 17, false, false, "PBCNone" }, // No periodic boundary
137 { 17, true, true, "PBCNone" }, // Update velocities, compute virial, without PBC
138 { 17, true, true, "PBCXYZ" }
139 }; // Update velocities, compute virial, with PBC
141 /*! \brief Test fixture for testing SETTLE.
143 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
147 std::unordered_map<std::string, t_pbc> pbcs_;
149 TestReferenceData refData_;
150 //! Checker for reference data
151 TestReferenceChecker checker_;
153 /*! \brief Test setup function.
155 * Setting up the PBCs and algorithms. Note, that corresponding string keywords
156 * have to be explicitly specified when parameters are initialied.
159 SettleTest() : checker_(refData_.rootChecker())
163 // PBC initialization
167 // Infinitely small box
168 matrix boxNone = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
169 set_pbc(&pbc, PbcType::No, boxNone);
170 pbcs_["PBCNone"] = pbc;
173 matrix boxXyz = { { real(1.86206), 0, 0 }, { 0, real(1.86206), 0 }, { 0, 0, real(1.86206) } };
174 set_pbc(&pbc, PbcType::Xyz, boxXyz);
175 pbcs_["PBCXYZ"] = pbc;
178 /*! \brief Check if the final interatomic distances are equal to target set by constraints.
180 * \param[in] numSettles Number of water molecules in the tested system.
181 * \param[in] tolerance Tolerance to compare floating point numbers.
182 * \param[in] testData An object, containing all the data structures needed by SETTLE.
184 static void checkConstrainsSatisfied(const int numSettles,
185 const FloatingPointTolerance tolerance,
186 const SettleTestData& testData)
188 for (int i = 0; i < numSettles; ++i)
190 const gmx::RVec& positionO = testData.xPrime_[i * testData.atomsPerSettle_ + 0];
191 const gmx::RVec& positionH1 = testData.xPrime_[i * testData.atomsPerSettle_ + 1];
192 const gmx::RVec& positionH2 = testData.xPrime_[i * testData.atomsPerSettle_ + 2];
194 real dOH = testData.dOH_;
195 real dHH = testData.dHH_;
197 EXPECT_REAL_EQ_TOL(dOH * dOH, distance2(positionO, positionH1), tolerance)
198 << formatString("for water %d. ", i);
199 EXPECT_REAL_EQ_TOL(dOH * dOH, distance2(positionO, positionH2), tolerance)
200 << formatString("for water %d. ", i);
201 EXPECT_REAL_EQ_TOL(dHH * dHH, distance2(positionH1, positionH2), tolerance)
202 << formatString("for water %d. ", i);
206 /*! \brief Check if the virial was updated and symmetric.
208 * The two tests on virial are:
209 * 1. If it was updated in case calcVirial is true.
210 * 2. If it is symmetrical.
212 * \param[in] calcVirial If the virial is computed.
213 * \param[in] tolerance Tolerance to compare floating point numbers.
214 * \param[in] testData An object, containing all the data structures needed by SETTLE.
216 static void checkVirialSymmetric(const bool calcVirial,
217 const FloatingPointTolerance tolerance,
218 const SettleTestData& testData)
220 for (int d = 0; d < DIM; ++d)
222 for (int dd = 0; dd < DIM; ++dd)
225 EXPECT_TRUE(calcVirial == (0. != testData.virial_[d][dd]))
226 << formatString("for virial component[%d][%d]. ", d, dd);
230 EXPECT_REAL_EQ_TOL(testData.virial_[d][dd], testData.virial_[dd][d], tolerance)
231 << formatString("Virial is not symmetrical for [%d][%d]. ", d, dd);
237 /*! \brief Check if the final positions correspond to reference values.
239 * \param[in] numSettles Number of water molecules in the tested system.
240 * \param[in] testData An object, containing all the data structures needed by SETTLE.
242 void checkFinalPositions(const int numSettles, const SettleTestData& testData)
244 TestReferenceChecker finalCoordinatesRef(
245 checker_.checkSequenceCompound("FinalCoordinates", numSettles));
246 for (int i = 0; i < numSettles; ++i)
248 TestReferenceChecker settlerRef(finalCoordinatesRef.checkCompound("Settler", nullptr));
249 TestReferenceChecker atomsRef(
250 settlerRef.checkSequenceCompound("Atoms", testData.atomsPerSettle_));
251 for (int j = 0; j < testData.atomsPerSettle_; ++j)
253 const gmx::RVec& xPrime = testData.xPrime_[testData.atomsPerSettle_ * i + j];
254 TestReferenceChecker xPrimeRef(atomsRef.checkCompound("Atom", nullptr));
255 xPrimeRef.checkReal(xPrime[XX], "XX");
256 xPrimeRef.checkReal(xPrime[YY], "YY");
257 xPrimeRef.checkReal(xPrime[ZZ], "ZZ");
262 /*! \brief Check if the final velocities correspond to reference values.
264 * \param[in] numSettles Number of water molecules in the tested system.
265 * \param[in] testData An object, containing all the data structures needed by SETTLE.
267 void checkFinalVelocities(const int numSettles, const SettleTestData& testData)
269 TestReferenceChecker finalCoordinatesRef(
270 checker_.checkSequenceCompound("FinalVelocities", numSettles));
271 for (int i = 0; i < numSettles; ++i)
273 TestReferenceChecker settlerRef(finalCoordinatesRef.checkCompound("Settler", nullptr));
274 TestReferenceChecker atomsRef(
275 settlerRef.checkSequenceCompound("Atoms", testData.atomsPerSettle_));
276 for (int j = 0; j < testData.atomsPerSettle_; ++j)
278 const gmx::RVec& v = testData.v_[testData.atomsPerSettle_ * i + j];
279 TestReferenceChecker vRef(atomsRef.checkCompound("Atom", nullptr));
280 vRef.checkReal(v[XX], "XX");
281 vRef.checkReal(v[YY], "YY");
282 vRef.checkReal(v[ZZ], "ZZ");
287 /*! \brief Check if the computed virial correspond to reference values.
289 * \param[in] testData An object, containing all the data structures needed by SETTLE.
291 void checkVirial(const SettleTestData& testData)
293 const tensor& virial = testData.virial_;
294 TestReferenceChecker virialRef(checker_.checkCompound("Virial", nullptr));
296 // TODO: Is it worth it to make this in a loop??
297 virialRef.checkReal(virial[XX][XX], "XX");
298 virialRef.checkReal(virial[XX][YY], "XY");
299 virialRef.checkReal(virial[XX][ZZ], "XZ");
300 virialRef.checkReal(virial[YY][XX], "YX");
301 virialRef.checkReal(virial[YY][YY], "YY");
302 virialRef.checkReal(virial[YY][ZZ], "YZ");
303 virialRef.checkReal(virial[ZZ][XX], "ZX");
304 virialRef.checkReal(virial[ZZ][YY], "ZY");
305 virialRef.checkReal(virial[ZZ][ZZ], "ZZ");
309 TEST_P(SettleTest, SatisfiesConstraints)
311 // Construct the list of runners
312 std::vector<std::unique_ptr<ISettleTestRunner>> runners;
313 // Add runners for CPU version
314 runners.emplace_back(std::make_unique<SettleHostTestRunner>());
315 // If supported, add runners for the GPU version for each available GPU
316 const bool addGpuRunners = GPU_SETTLE_SUPPORTED;
319 for (const auto& testDevice : getTestHardwareEnvironment()->getTestDeviceList())
321 runners.emplace_back(std::make_unique<SettleDeviceTestRunner>(*testDevice));
324 for (const auto& runner : runners)
326 // Make some symbolic names for the parameter combination.
327 SettleTestParameters params = GetParam();
329 int numSettles = params.numSettles;
330 bool updateVelocities = params.updateVelocities;
331 bool calcVirial = params.calcVirial;
332 std::string pbcName = params.pbcName;
335 // Make a string that describes which parameter combination is
336 // being tested, to help make failing tests comprehensible.
337 std::string testDescription = formatString(
338 "Testing %s with %d SETTLEs, %s, %svelocities and %scalculating the virial.",
339 runner->hardwareDescription().c_str(),
342 updateVelocities ? "with " : "without ",
343 calcVirial ? "" : "not ");
345 SCOPED_TRACE(testDescription);
347 auto testData = std::make_unique<SettleTestData>(numSettles);
349 ASSERT_LE(numSettles, testData->xPrime_.size() / testData->atomsPerSettle_)
350 << "cannot test that many SETTLEs. " << testDescription;
352 t_pbc pbc = pbcs_.at(pbcName);
355 runner->applySettle(testData.get(), pbc, updateVelocities, calcVirial, testDescription);
357 // The necessary tolerances for the test to pass were determined
358 // empirically. This isn't nice, but the required behavior that
359 // SETTLE produces constrained coordinates consistent with
360 // sensible sampling needs to be tested at a much higher level.
361 // TODO: Re-evaluate the tolerances.
362 real dOH = testData->dOH_;
363 FloatingPointTolerance tolerance = relativeToleranceAsPrecisionDependentUlp(dOH * dOH, 80, 380);
364 FloatingPointTolerance toleranceVirial = absoluteTolerance(0.000001);
366 FloatingPointTolerance tolerancePositions = absoluteTolerance(0.000001);
367 FloatingPointTolerance toleranceVelocities = absoluteTolerance(0.0001);
369 checkConstrainsSatisfied(numSettles, tolerance, *testData);
370 checkVirialSymmetric(calcVirial, toleranceVirial, *testData);
372 checker_.setDefaultTolerance(tolerancePositions);
373 checkFinalPositions(numSettles, *testData);
375 if (updateVelocities)
377 checker_.setDefaultTolerance(toleranceVelocities);
378 checkFinalVelocities(numSettles, *testData);
383 checker_.setDefaultTolerance(toleranceVirial);
384 checkVirial(*testData);
389 // Run test on pre-determined set of combinations for test parameters, which include the numbers of SETTLEs (water
390 // molecules), whether or not velocities are updated and virial contribution is computed, was the PBC enabled.
391 // The test will cycle through all available runners, including CPU and, if applicable, GPU implementations of SETTLE.
392 INSTANTIATE_TEST_SUITE_P(WithParameters, SettleTest, ::testing::ValuesIn(parametersSets));