2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 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 the Leap-Frog integrator
38 * The test creates a system of independent particles exerting constant
39 * external forces and makes several numerical integration timesteps.
40 * The results are compared with the analytical solution (for the systems
41 * without the temperature coupling) and with the pre-computed reference
42 * values. The tests use runners that are created for each available
43 * implementation of the tested algorithm.
45 * \todo Add tests for integrators with pressure control.
46 * \todo Add PBC handling test.
48 * \author Artem Zhmurov <zhmurov@gmail.com>
49 * \ingroup module_mdlib
58 #include <gtest/gtest.h>
60 #include "gromacs/hardware/device_management.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/math/vectypes.h"
63 #include "gromacs/mdtypes/mdatom.h"
64 #include "gromacs/utility/stringutil.h"
66 #include "testutils/refdata.h"
67 #include "testutils/test_hardware_environment.h"
68 #include "testutils/testasserts.h"
70 #include "leapfrogtestdata.h"
71 #include "leapfrogtestrunners.h"
80 /*! \brief The parameters for the test.
82 * The test will run for combinations of:
87 * 4. Velocity components
89 * 6. Number of temperature coupling groups
91 struct LeapFrogTestParameters
93 //! Total number of atoms
97 //! Number of integration steps
103 //! Number of temperature coupling group (zero for no temperature coupling)
104 int numTCoupleGroups
;
105 //! Number of steps between pressure coupling steps (zero for no pressure coupling).
109 //! The set of parameters combinations to run the test on
110 const LeapFrogTestParameters parametersSets
[] = {
111 { 1, 0.001, 1, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, 0, 0 }, // Zero velocity and force
112 { 1, 0.001, 1, { 0.0, 0.0, 0.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // Zero velocity
113 { 1, 0.001, 1, { 1.0, -2.0, 3.0 }, { 0.0, 0.0, 0.0 }, 0, 0 }, // Zero force
114 { 1, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 1 particle
115 { 10, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 10 particles
116 { 100, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 100 particles
117 { 300, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 300 particles
118 { 1, 0.0005, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 0.0005 ps timestep
119 { 1, 0.001, 10, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 10 step
120 { 1, 0.001, 100, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 0 }, // 100 steps
121 { 100, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 1, 0 }, // 1 temperature couple group
122 { 100, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 2, 0 }, // 2 temperature couple groups
123 { 100, 0.001, 1, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 10, 0 }, // 10 temperature couple groups
124 { 100, 0.001, 10, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 1 }, // With pressure coupling
125 { 100, 0.001, 10, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 2, 1 }, // With both temperature and pressure coupling
126 { 100, 0.001, 10, { 1.0, -2.0, 3.0 }, { -3.0, 2.0, -1.0 }, 0, 3 }
127 }; // Do pressure coupling not on every step
130 /*! \brief Test fixture for LeapFrog integrator.
132 class LeapFrogTest
: public ::testing::TestWithParam
<LeapFrogTestParameters
>
136 TestReferenceData refData_
;
137 //! Checker for reference data
138 TestReferenceChecker checker_
;
140 LeapFrogTest() : checker_(refData_
.rootChecker()) {}
142 /*! \brief Test the numerical integrator against analytical solution for simple constant force case.
144 * \param[in] tolerance Tolerance
145 * \param[in] testData Test data object
146 * \param[in] totalTime Total numerical integration time
148 static void testAgainstAnalyticalSolution(FloatingPointTolerance tolerance
,
149 const LeapFrogTestData
& testData
,
150 const real totalTime
)
152 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
156 for (int d
= 0; d
< DIM
; d
++)
158 // Analytical solution for constant-force particle movement
159 real x0
= testData
.x0_
[i
][d
];
160 real v0
= testData
.v0_
[i
][d
];
161 real f
= testData
.f_
[i
][d
];
162 real im
= testData
.inverseMasses_
[i
];
164 xAnalytical
[d
] = x0
+ v0
* totalTime
+ 0.5 * f
* totalTime
* totalTime
* im
;
165 vAnalytical
[d
] = v0
+ f
* totalTime
* im
;
167 EXPECT_REAL_EQ_TOL(xAnalytical
[d
], testData
.xPrime_
[i
][d
], tolerance
) << gmx::formatString(
168 "Coordinate %d of atom %d is different from analytical solution.", d
, i
);
170 EXPECT_REAL_EQ_TOL(vAnalytical
[d
], testData
.v_
[i
][d
], tolerance
) << gmx::formatString(
171 "Velocity component %d of atom %d is different from analytical solution.", d
, i
);
176 /*! \brief Test the numerical integrator against pre-computed reference values.
178 * \param[in] testData Test data object
180 void testAgainstReferenceData(const LeapFrogTestData
& testData
)
182 TestReferenceChecker
finalPositionsRef(
183 checker_
.checkSequenceCompound("FinalPositions", testData
.numAtoms_
));
184 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
186 const gmx::RVec
& xPrime
= testData
.xPrime_
[i
];
187 TestReferenceChecker
xPrimeRef(finalPositionsRef
.checkCompound("Atom", nullptr));
188 xPrimeRef
.checkReal(xPrime
[XX
], "XX");
189 xPrimeRef
.checkReal(xPrime
[YY
], "YY");
190 xPrimeRef
.checkReal(xPrime
[ZZ
], "ZZ");
193 TestReferenceChecker
finalVelocitiesRef(
194 checker_
.checkSequenceCompound("FinalVelocities", testData
.numAtoms_
));
195 for (int i
= 0; i
< testData
.numAtoms_
; i
++)
197 const gmx::RVec
& v
= testData
.v_
[i
];
198 TestReferenceChecker
vRef(finalVelocitiesRef
.checkCompound("Atom", nullptr));
199 vRef
.checkReal(v
[XX
], "XX");
200 vRef
.checkReal(v
[YY
], "YY");
201 vRef
.checkReal(v
[ZZ
], "ZZ");
206 TEST_P(LeapFrogTest
, SimpleIntegration
)
208 // Construct the list of runners
209 std::vector
<std::unique_ptr
<ILeapFrogTestRunner
>> runners
;
210 // Add runners for CPU version
211 runners
.emplace_back(std::make_unique
<LeapFrogHostTestRunner
>());
212 // If using CUDA, add runners for the GPU version for each available GPU
213 const bool addGpuRunners
= HAVE_GPU_LEAPFROG
;
216 for (const auto& testDevice
: getTestHardwareEnvironment()->getTestDeviceList())
218 runners
.emplace_back(std::make_unique
<LeapFrogDeviceTestRunner
>(*testDevice
));
222 for (const auto& runner
: runners
)
224 LeapFrogTestParameters parameters
= GetParam();
226 std::string testDescription
= formatString(
227 "Testing on %s with %d atoms for %d timesteps with %d temperature coupling "
229 "%s pressure coupling (dt = %f, v0=(%f, %f, %f), f0=(%f, %f, %f), nstpcouple = "
231 runner
->hardwareDescription().c_str(), parameters
.numAtoms
, parameters
.numSteps
,
232 parameters
.numTCoupleGroups
, parameters
.nstpcouple
== 0 ? "without" : "with",
233 parameters
.timestep
, parameters
.v
[XX
], parameters
.v
[YY
], parameters
.v
[ZZ
],
234 parameters
.f
[XX
], parameters
.f
[YY
], parameters
.f
[ZZ
], parameters
.nstpcouple
);
235 SCOPED_TRACE(testDescription
);
237 std::unique_ptr
<LeapFrogTestData
> testData
= std::make_unique
<LeapFrogTestData
>(
238 parameters
.numAtoms
, parameters
.timestep
, parameters
.v
, parameters
.f
,
239 parameters
.numTCoupleGroups
, parameters
.nstpcouple
);
241 runner
->integrate(testData
.get(), parameters
.numSteps
);
243 real totalTime
= parameters
.numSteps
* parameters
.timestep
;
244 // TODO For the case of constant force, the numerical scheme is exact and
245 // the only source of errors is floating point arithmetic. Hence,
246 // the tolerance can be calculated.
247 FloatingPointTolerance tolerance
= absoluteTolerance(parameters
.numSteps
* 0.000005);
249 // Test against the analytical solution (without temperature coupling)
250 if (parameters
.numTCoupleGroups
== 0 && parameters
.nstpcouple
== 0)
252 testAgainstAnalyticalSolution(tolerance
, *testData
, totalTime
);
255 checker_
.setDefaultTolerance(tolerance
);
256 testAgainstReferenceData(*testData
);
260 INSTANTIATE_TEST_CASE_P(WithParameters
, LeapFrogTest
, ::testing::ValuesIn(parametersSets
));