1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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.
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 Tests for SETTLE constraints
37 *
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.
42 *
43 * The tests check:
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.
48 *
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.
52 *
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.
56 *
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.
60 *
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.
65 *
66 * \author Mark Abraham <mark.j.abraham@gmail.com>
67 * \author Artem Zhmurov <zhmurov@gmail.com>
68 *
69 * \ingroup module_mdlib
70 */
71 #include "gmxpre.h"
72
73 #include "gromacs/mdlib/settle.h"
74
75 #include "config.h"
76
77 #include <algorithm>
78 #include <unordered_map>
79 #include <vector>
80
81 #include <gtest/gtest.h>
82
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"
93
94 #include "gromacs/mdlib/tests/watersystem.h"
95 #include "testutils/refdata.h"
96 #include "testutils/test_hardware_environment.h"
97 #include "testutils/testasserts.h"
98
99 #include "settletestdata.h"
100 #include "settletestrunners.h"
101
102 namespace gmx
103 {
104 namespace test
105 {
106 namespace
107 {
108
109 /*! \brief Parameters that will vary from test to test.
110 */
111 struct SettleTestParameters
112 {
113 //! Number of water molecules (SETTLEs) [1, 2, 4, 5, 7, 10, 12, 15, 17]
114 int numSettles;
115 //! If the velocities should be updated while constraining [true/false]
116 bool updateVelocities;
117 //! If the virial should be computed [true/false]
118 bool calcVirial;
119 //! Periodic boundary conditions [PBCXYZ/PBCNone]
120 std::string pbcName;
121 };
122
123 /*! \brief Sets of parameters on which to run the tests.
124 */
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
140
141 /*! \brief Test fixture for testing SETTLE.
142 */
143 class SettleTest : public ::testing::TestWithParam<SettleTestParameters>
144 {
145 public:
146 //! PBC setups
147 std::unordered_map<std::string, t_pbc> pbcs_;
148 //! Reference data
149 TestReferenceData refData_;
150 //! Checker for reference data
151 TestReferenceChecker checker_;
152
153 /*! \brief Test setup function.
154 *
155 * Setting up the PBCs and algorithms. Note, that corresponding string keywords
156 * have to be explicitly specified when parameters are initialied.
157 *
158 */
SettleTest()159 SettleTest() : checker_(refData_.rootChecker())
160 {
161
162 //
163 // PBC initialization
164 //
165 t_pbc pbc;
166
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;
171
172 // Rectangular box
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;
176 }
177
178 /*! \brief Check if the final interatomic distances are equal to target set by constraints.
179 *
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.
183 */
checkConstrainsSatisfied(const int numSettles,const FloatingPointTolerance tolerance,const SettleTestData & testData)184 static void checkConstrainsSatisfied(const int numSettles,
185 const FloatingPointTolerance tolerance,
186 const SettleTestData& testData)
187 {
188 for (int i = 0; i < numSettles; ++i)
189 {
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];
193
194 real dOH = testData.dOH_;
195 real dHH = testData.dHH_;
196
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);
203 }
204 }
205
206 /*! \brief Check if the virial was updated and symmetric.
207 *
208 * The two tests on virial are:
209 * 1. If it was updated in case calcVirial is true.
210 * 2. If it is symmetrical.
211 *
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.
215 */
checkVirialSymmetric(const bool calcVirial,const FloatingPointTolerance tolerance,const SettleTestData & testData)216 static void checkVirialSymmetric(const bool calcVirial,
217 const FloatingPointTolerance tolerance,
218 const SettleTestData& testData)
219 {
220 for (int d = 0; d < DIM; ++d)
221 {
222 for (int dd = 0; dd < DIM; ++dd)
223 {
224
225 EXPECT_TRUE(calcVirial == (0. != testData.virial_[d][dd]))
226 << formatString("for virial component[%d][%d]. ", d, dd);
227
228 if (calcVirial)
229 {
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);
232 }
233 }
234 }
235 }
236
237 /*! \brief Check if the final positions correspond to reference values.
238 *
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.
241 */
checkFinalPositions(const int numSettles,const SettleTestData & testData)242 void checkFinalPositions(const int numSettles, const SettleTestData& testData)
243 {
244 TestReferenceChecker finalCoordinatesRef(
245 checker_.checkSequenceCompound("FinalCoordinates", numSettles));
246 for (int i = 0; i < numSettles; ++i)
247 {
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)
252 {
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");
258 }
259 }
260 }
261
262 /*! \brief Check if the final velocities correspond to reference values.
263 *
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.
266 */
checkFinalVelocities(const int numSettles,const SettleTestData & testData)267 void checkFinalVelocities(const int numSettles, const SettleTestData& testData)
268 {
269 TestReferenceChecker finalCoordinatesRef(
270 checker_.checkSequenceCompound("FinalVelocities", numSettles));
271 for (int i = 0; i < numSettles; ++i)
272 {
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)
277 {
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");
283 }
284 }
285 }
286
287 /*! \brief Check if the computed virial correspond to reference values.
288 *
289 * \param[in] testData An object, containing all the data structures needed by SETTLE.
290 */
checkVirial(const SettleTestData & testData)291 void checkVirial(const SettleTestData& testData)
292 {
293 const tensor& virial = testData.virial_;
294 TestReferenceChecker virialRef(checker_.checkCompound("Virial", nullptr));
295
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");
306 }
307 };
308
TEST_P(SettleTest,SatisfiesConstraints)309 TEST_P(SettleTest, SatisfiesConstraints)
310 {
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 using CUDA, add runners for the GPU version for each available GPU
316 if (GMX_GPU_CUDA)
317 {
318 for (const auto& testDevice : getTestHardwareEnvironment()->getTestDeviceList())
319 {
320 runners.emplace_back(std::make_unique<SettleDeviceTestRunner>(*testDevice));
321 }
322 }
323 for (const auto& runner : runners)
324 {
325 // Make some symbolic names for the parameter combination.
326 SettleTestParameters params = GetParam();
327
328 int numSettles = params.numSettles;
329 bool updateVelocities = params.updateVelocities;
330 bool calcVirial = params.calcVirial;
331 std::string pbcName = params.pbcName;
332
333
334 // Make a string that describes which parameter combination is
335 // being tested, to help make failing tests comprehensible.
336 std::string testDescription = formatString(
337 "Testing %s with %d SETTLEs, %s, %svelocities and %scalculating the virial.",
338 runner->hardwareDescription().c_str(), numSettles, pbcName.c_str(),
339 updateVelocities ? "with " : "without ", calcVirial ? "" : "not ");
340
341 SCOPED_TRACE(testDescription);
342
343 auto testData = std::make_unique<SettleTestData>(numSettles);
344
345 ASSERT_LE(numSettles, testData->xPrime_.size() / testData->atomsPerSettle_)
346 << "cannot test that many SETTLEs. " << testDescription;
347
348 t_pbc pbc = pbcs_.at(pbcName);
349
350 // Apply SETTLE
351 runner->applySettle(testData.get(), pbc, updateVelocities, calcVirial, testDescription);
352
353 // The necessary tolerances for the test to pass were determined
354 // empirically. This isn't nice, but the required behavior that
355 // SETTLE produces constrained coordinates consistent with
356 // sensible sampling needs to be tested at a much higher level.
357 // TODO: Re-evaluate the tolerances.
358 real dOH = testData->dOH_;
359 FloatingPointTolerance tolerance = relativeToleranceAsPrecisionDependentUlp(dOH * dOH, 80, 380);
360 FloatingPointTolerance toleranceVirial = absoluteTolerance(0.000001);
361
362 FloatingPointTolerance tolerancePositions = absoluteTolerance(0.000001);
363 FloatingPointTolerance toleranceVelocities = absoluteTolerance(0.0001);
364
365 checkConstrainsSatisfied(numSettles, tolerance, *testData);
366 checkVirialSymmetric(calcVirial, toleranceVirial, *testData);
367
368 checker_.setDefaultTolerance(tolerancePositions);
369 checkFinalPositions(numSettles, *testData);
370
371 if (updateVelocities)
372 {
373 checker_.setDefaultTolerance(toleranceVelocities);
374 checkFinalVelocities(numSettles, *testData);
375 }
376
377 if (calcVirial)
378 {
379 checker_.setDefaultTolerance(toleranceVirial);
380 checkVirial(*testData);
381 }
382 }
383 }
384
385 // Run test on pre-determined set of combinations for test parameters, which include the numbers of SETTLEs (water
386 // molecules), whether or not velocities are updated and virial contribution is computed, was the PBC enabled.
387 // The test will cycle through all available runners, including CPU and, if applicable, GPU implementations of SETTLE.
388 INSTANTIATE_TEST_CASE_P(WithParameters, SettleTest, ::testing::ValuesIn(parametersSets));
389
390 } // namespace
391 } // namespace test
392 } // namespace gmx
393