2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2014,2015,2016,2018,2019,2020, 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/gpu_utils/gpu_testutils.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/testasserts.h"
98 #include "settletestdata.h"
99 #include "settletestrunners.h"
108 /*! \brief Parameters that will vary from test to test.
110 struct SettleTestParameters
112 //! Number of water molecules (SETTLEs) [1, 2, 4, 5, 7, 10, 12, 15, 17]
114 //! If the velocities should be updated while constraining [true/false]
115 bool updateVelocities;
116 //! If the virial should be computed [true/false]
118 //! Periodic boundary conditions [PBCXYZ/PBCNone]
122 /*! \brief Sets of parameters on which to run the tests.
124 const SettleTestParameters parametersSets[] = {
125 { 1, false, false, "PBCXYZ" }, // 1 water molecule
126 { 2, false, false, "PBCXYZ" }, // 2 water molecules
127 { 4, false, false, "PBCXYZ" }, // 4 water molecules
128 { 5, false, false, "PBCXYZ" }, // 5 water molecules
129 { 6, false, false, "PBCXYZ" }, // 6 water molecules
130 { 10, false, false, "PBCXYZ" }, // 10 water molecules
131 { 12, false, false, "PBCXYZ" }, // 12 water molecules
132 { 15, false, false, "PBCXYZ" }, // 15 water molecules
133 { 17, true, false, "PBCXYZ" }, // Update velocities
134 { 17, false, true, "PBCXYZ" }, // Compute virial
135 { 17, false, false, "PBCNone" }, // No periodic boundary
136 { 17, true, true, "PBCNone" }, // Update velocities, compute virial, without PBC
137 { 17, true, true, "PBCXYZ" }
138 }; // Update velocities, compute virial, with PBC
140 /*! \brief Test fixture for testing SETTLE.
142 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
146 std::unordered_map<std::string, t_pbc> pbcs_;
147 //! Runners (CPU and GPU versions of SETTLE)
148 std::unordered_map<std::string,
149 void (*)(SettleTestData* testData, const t_pbc pbc, const bool updateVelocities, const bool calcVirial, const std::string& testDescription)>
152 TestReferenceData refData_;
153 //! Checker for reference data
154 TestReferenceChecker checker_;
156 /*! \brief Test setup function.
158 * Setting up the PBCs and algorithms. Note, that corresponding string keywords
159 * have to be explicitly specified when parameters are initialied.
162 SettleTest() : checker_(refData_.rootChecker())
166 // PBC initialization
170 // Infinitely small box
171 matrix boxNone = { { 0, 0, 0 }, { 0, 0, 0 }, { 0, 0, 0 } };
172 set_pbc(&pbc, PbcType::No, boxNone);
173 pbcs_["PBCNone"] = pbc;
176 matrix boxXyz = { { real(1.86206), 0, 0 }, { 0, real(1.86206), 0 }, { 0, 0, real(1.86206) } };
177 set_pbc(&pbc, PbcType::Xyz, boxXyz);
178 pbcs_["PBCXYZ"] = pbc;
181 // All SETTLE runners should be registered here under appropriate conditions
183 runners_["SETTLE"] = applySettle;
185 // CUDA version will be tested only if:
186 // 1. The code was compiled with CUDA
187 // 2. There is a CUDA-capable GPU in a system
188 // 3. This GPU is detectable
189 // 4. GPU detection was not disabled by GMX_DISABLE_GPU_DETECTION environment variable
190 if (s_hasCompatibleGpus)
192 if (GMX_GPU == GMX_GPU_CUDA && s_hasCompatibleGpus)
194 runners_["SETTLE_GPU"] = applySettleGpu;
199 /*! \brief Check if the final interatomic distances are equal to target set by constraints.
201 * \param[in] numSettles Number of water molecules in the tested system.
202 * \param[in] tolerance Tolerance to compare floating point numbers.
203 * \param[in] testData An object, containing all the data structures needed by SETTLE.
205 static void checkConstrainsSatisfied(const int numSettles,
206 const FloatingPointTolerance tolerance,
207 const SettleTestData& testData)
209 for (int i = 0; i < numSettles; ++i)
211 const gmx::RVec& positionO = testData.xPrime_[i * testData.atomsPerSettle_ + 0];
212 const gmx::RVec& positionH1 = testData.xPrime_[i * testData.atomsPerSettle_ + 1];
213 const gmx::RVec& positionH2 = testData.xPrime_[i * testData.atomsPerSettle_ + 2];
215 real dOH = testData.dOH_;
216 real dHH = testData.dHH_;
218 EXPECT_REAL_EQ_TOL(dOH * dOH, distance2(positionO, positionH1), tolerance)
219 << formatString("for water %d. ", i);
220 EXPECT_REAL_EQ_TOL(dOH * dOH, distance2(positionO, positionH2), tolerance)
221 << formatString("for water %d. ", i);
222 EXPECT_REAL_EQ_TOL(dHH * dHH, distance2(positionH1, positionH2), tolerance)
223 << formatString("for water %d. ", i);
227 /*! \brief Check if the virial was updated and symmetric.
229 * The two tests on virial are:
230 * 1. If it was updated in case calcVirial is true.
231 * 2. If it is symmetrical.
233 * \param[in] calcVirial If the virial is computed.
234 * \param[in] tolerance Tolerance to compare floating point numbers.
235 * \param[in] testData An object, containing all the data structures needed by SETTLE.
237 static void checkVirialSymmetric(const bool calcVirial,
238 const FloatingPointTolerance tolerance,
239 const SettleTestData& testData)
241 for (int d = 0; d < DIM; ++d)
243 for (int dd = 0; dd < DIM; ++dd)
246 EXPECT_TRUE(calcVirial == (0. != testData.virial_[d][dd]))
247 << formatString("for virial component[%d][%d]. ", d, dd);
251 EXPECT_REAL_EQ_TOL(testData.virial_[d][dd], testData.virial_[dd][d], tolerance)
252 << formatString("Virial is not symmetrical for [%d][%d]. ", d, dd);
258 /*! \brief Check if the final positions correspond to reference values.
260 * \param[in] numSettles Number of water molecules in the tested system.
261 * \param[in] testData An object, containing all the data structures needed by SETTLE.
263 void checkFinalPositions(const int numSettles, const SettleTestData& testData)
265 TestReferenceChecker finalCoordinatesRef(
266 checker_.checkSequenceCompound("FinalCoordinates", numSettles));
267 for (int i = 0; i < numSettles; ++i)
269 TestReferenceChecker settlerRef(finalCoordinatesRef.checkCompound("Settler", nullptr));
270 TestReferenceChecker atomsRef(
271 settlerRef.checkSequenceCompound("Atoms", testData.atomsPerSettle_));
272 for (int j = 0; j < testData.atomsPerSettle_; ++j)
274 const gmx::RVec& xPrime = testData.xPrime_[testData.atomsPerSettle_ * i + j];
275 TestReferenceChecker xPrimeRef(atomsRef.checkCompound("Atom", nullptr));
276 xPrimeRef.checkReal(xPrime[XX], "XX");
277 xPrimeRef.checkReal(xPrime[YY], "YY");
278 xPrimeRef.checkReal(xPrime[ZZ], "ZZ");
283 /*! \brief Check if the final velocities correspond to reference values.
285 * \param[in] numSettles Number of water molecules in the tested system.
286 * \param[in] testData An object, containing all the data structures needed by SETTLE.
288 void checkFinalVelocities(const int numSettles, const SettleTestData& testData)
290 TestReferenceChecker finalCoordinatesRef(
291 checker_.checkSequenceCompound("FinalVelocities", numSettles));
292 for (int i = 0; i < numSettles; ++i)
294 TestReferenceChecker settlerRef(finalCoordinatesRef.checkCompound("Settler", nullptr));
295 TestReferenceChecker atomsRef(
296 settlerRef.checkSequenceCompound("Atoms", testData.atomsPerSettle_));
297 for (int j = 0; j < testData.atomsPerSettle_; ++j)
299 const gmx::RVec& v = testData.v_[testData.atomsPerSettle_ * i + j];
300 TestReferenceChecker vRef(atomsRef.checkCompound("Atom", nullptr));
301 vRef.checkReal(v[XX], "XX");
302 vRef.checkReal(v[YY], "YY");
303 vRef.checkReal(v[ZZ], "ZZ");
308 /*! \brief Check if the computed virial correspond to reference values.
310 * \param[in] testData An object, containing all the data structures needed by SETTLE.
312 void checkVirial(const SettleTestData& testData)
314 const tensor& virial = testData.virial_;
315 TestReferenceChecker virialRef(checker_.checkCompound("Virial", nullptr));
317 // TODO: Is it worth it to make this in a loop??
318 virialRef.checkReal(virial[XX][XX], "XX");
319 virialRef.checkReal(virial[XX][YY], "XY");
320 virialRef.checkReal(virial[XX][ZZ], "XZ");
321 virialRef.checkReal(virial[YY][XX], "YX");
322 virialRef.checkReal(virial[YY][YY], "YY");
323 virialRef.checkReal(virial[YY][ZZ], "YZ");
324 virialRef.checkReal(virial[ZZ][XX], "ZX");
325 virialRef.checkReal(virial[ZZ][YY], "ZY");
326 virialRef.checkReal(virial[ZZ][ZZ], "ZZ");
329 //! Store whether any compatible GPUs exist.
330 static bool s_hasCompatibleGpus;
331 //! Before any test is run, work out whether any compatible GPUs exist.
332 static void SetUpTestCase() { s_hasCompatibleGpus = canComputeOnGpu(); }
335 bool SettleTest::s_hasCompatibleGpus = false;
337 TEST_P(SettleTest, SatisfiesConstraints)
339 // Cycle through all available runners
340 for (const auto& runner : runners_)
342 std::string runnerName = runner.first;
344 // Make some symbolic names for the parameter combination.
345 SettleTestParameters params = GetParam();
347 int numSettles = params.numSettles;
348 bool updateVelocities = params.updateVelocities;
349 bool calcVirial = params.calcVirial;
350 std::string pbcName = params.pbcName;
353 // Make a string that describes which parameter combination is
354 // being tested, to help make failing tests comprehensible.
355 std::string testDescription = formatString(
356 "Testing %s with %d SETTLEs, %s, %svelocities and %scalculating the virial.",
357 runnerName.c_str(), numSettles, pbcName.c_str(),
358 updateVelocities ? "with " : "without ", calcVirial ? "" : "not ");
360 SCOPED_TRACE(testDescription);
362 auto testData = std::make_unique<SettleTestData>(numSettles);
364 ASSERT_LE(numSettles, testData->xPrime_.size() / testData->atomsPerSettle_)
365 << "cannot test that many SETTLEs. " << testDescription;
367 t_pbc pbc = pbcs_.at(pbcName);
370 runner.second(testData.get(), pbc, updateVelocities, calcVirial, testDescription);
372 // The necessary tolerances for the test to pass were determined
373 // empirically. This isn't nice, but the required behavior that
374 // SETTLE produces constrained coordinates consistent with
375 // sensible sampling needs to be tested at a much higher level.
376 // TODO: Re-evaluate the tolerances.
377 real dOH = testData->dOH_;
378 FloatingPointTolerance tolerance = relativeToleranceAsPrecisionDependentUlp(dOH * dOH, 80, 380);
379 FloatingPointTolerance toleranceVirial = absoluteTolerance(0.000001);
381 FloatingPointTolerance tolerancePositions = absoluteTolerance(0.000001);
382 FloatingPointTolerance toleranceVelocities = absoluteTolerance(0.0001);
384 checkConstrainsSatisfied(numSettles, tolerance, *testData);
385 checkVirialSymmetric(calcVirial, toleranceVirial, *testData);
387 checker_.setDefaultTolerance(tolerancePositions);
388 checkFinalPositions(numSettles, *testData);
390 if (updateVelocities)
392 checker_.setDefaultTolerance(toleranceVelocities);
393 checkFinalVelocities(numSettles, *testData);
398 checker_.setDefaultTolerance(toleranceVirial);
399 checkVirial(*testData);
404 // Run test on pre-determined set of combinations for test parameters, which include the numbers of SETTLEs (water
405 // molecules), whether or not velocities are updated and virial contribution is computed, was the PBC enabled.
406 // The test will cycle through all available runners, including CPU and, if applicable, GPU implementations of SETTLE.
407 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest, ::testing::ValuesIn(parametersSets));