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.
37 * \brief Tests to verify that a simulator that only does some actions
38 * periodically produces the expected results.
40 * \author Mark Abraham <mark.j.abraham@gmail.com>
41 * \ingroup module_mdrun_integration_tests
49 #include "gromacs/trajectory/energyframe.h"
50 #include "gromacs/utility/strconvert.h"
51 #include "gromacs/utility/stringutil.h"
53 #include "testutils/simulationdatabase.h"
54 #include "testutils/testasserts.h"
56 #include "energycomparison.h"
57 #include "energyreader.h"
58 #include "moduletest.h"
67 /*! \brief Mdp parameters that determine the manner of simulation
69 using PropagationParameters = MdpFieldValues;
71 /*! \brief Mdp parameters that should only affect the observations,
72 * not the simulation propagation. */
73 using PeriodicOutputParameters = MdpFieldValues;
75 //! Helper type of output file names for the reference mdrun call
76 struct ReferenceFileNames
78 //! Name of energy file
79 std::string edrFileName_;
82 //! Function type to produce sets of .mdp parameters for testing periodic output
83 using OutputParameterGeneratorFunction = std::vector<PeriodicOutputParameters> (*)();
85 /*! \brief Test fixture base for comparing a simulator with one that
86 * does everything every step
88 * This test ensures that two simulator code paths called via
89 * different mdp options yield identical energy trajectories,
90 * up to some (arbitrary) precision.
92 * These tests are useful to check that periodic actions implemented
93 * in simulators are correct, and that different code paths expected
94 * to yield identical results are equivalent.
96 class PeriodicActionsTest :
97 public MdrunTestFixture,
98 public ::testing::WithParamInterface<std::tuple<PropagationParameters, OutputParameterGeneratorFunction>>
101 // PeriodicActionsTest();
102 //! Run mdrun with given output parameters
103 void doMdrun(const PeriodicOutputParameters& output);
104 //! Generate reference data from mdrun writing everything every step.
105 void prepareReferenceData();
106 //! Names for the output files from the reference mdrun call
107 ReferenceFileNames referenceFileNames_ = { fileManager_.getTemporaryFilePath("reference.edr") };
108 //! Functor for energy comparison
109 EnergyComparison energyComparison_{ EnergyComparison::defaultEnergyTermsToCompare(),
110 MaxNumFrames::compareAllFrames() };
111 //! Names of energies compared by energyComparison_
112 std::vector<std::string> namesOfEnergiesToMatch_ = energyComparison_.getEnergyNames();
115 void PeriodicActionsTest::doMdrun(const PeriodicOutputParameters& output)
117 auto propagation = std::get<0>(GetParam());
119 formatString("Doing %s simulation with %s integrator, %s tcoupling and %s pcoupling\n",
120 propagation["simulationName"].c_str(), propagation["integrator"].c_str(),
121 propagation["tcoupl"].c_str(), propagation["pcoupl"].c_str()));
122 auto mdpFieldValues = prepareMdpFieldValues(propagation["simulationName"], propagation["integrator"],
123 propagation["tcoupl"], propagation["pcoupl"]);
124 mdpFieldValues.insert(propagation.begin(), propagation.end());
125 mdpFieldValues.insert(output.begin(), output.end());
127 // prepare the tpr file
130 caller.append("grompp");
131 caller.addOption("-maxwarn", propagation["maxGromppWarningsTolerated"]);
132 runner_.useTopGroAndNdxFromDatabase(propagation["simulationName"]);
133 runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
134 EXPECT_EQ(0, runner_.callGrompp(caller));
138 CommandLine mdrunCaller;
139 mdrunCaller.append("mdrun");
140 ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
144 void PeriodicActionsTest::prepareReferenceData()
146 SCOPED_TRACE("Preparing reference data");
148 // Configure the SimulationRunner to write output to files we can
149 // use in comparing each test run.
150 std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
152 // Run the reference simulation with everything output every step
153 PeriodicOutputParameters outputEveryStep = {
154 { "nstenergy", "1" },
157 { "description", "output everything every step" },
159 doMdrun(outputEveryStep);
161 // Restore the standard filenames we'll use for the runs whose behavior we are testing.
162 std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
165 /*! \brief Compare the next frame returned by \c testFrameReader with the
166 * matching frame from \c referenceTestReader using \c comparator.
168 * \returns True when a successful comparison was made
170 template<typename Reader, typename Comparator>
171 bool compareFrames(Reader referenceFrameReader, Reader testFrameReader, const Comparator& comparator)
173 if (!testFrameReader->readNextFrame())
178 auto testFrame = testFrameReader->frame();
179 std::string frameName = testFrame.frameName();
180 SCOPED_TRACE("Found frame from test file named " + frameName);
181 bool foundMatch = false;
184 if (referenceFrameReader->readNextFrame())
186 auto referenceFrame = referenceFrameReader->frame();
187 // Can the reference and test frames be compared, because they have the same name?
188 if (frameName == referenceFrame.frameName())
190 SCOPED_TRACE("Found frame from reference file named " + frameName);
191 comparator(referenceFrame, testFrame);
197 ADD_FAILURE() << "Ran out of reference frames to compare";
204 TEST_P(PeriodicActionsTest, PeriodicActionsAgreeWithReference)
206 auto propagation = std::get<0>(GetParam());
207 SCOPED_TRACE(formatString("Comparing two simulations of '%s' with integrator '%s'",
208 propagation["simulationName"].c_str(), propagation["integrator"].c_str()));
210 prepareReferenceData();
212 auto outputParametersGenerator = std::get<1>(GetParam());
213 for (const PeriodicOutputParameters& output : outputParametersGenerator())
215 SCOPED_TRACE("Comparing to observe " + output.at("description"));
217 // Run the test simulation
220 // Prepare readers for the reference and test output files
221 auto referenceEnergyFrameReader =
222 openEnergyFileToReadTerms(referenceFileNames_.edrFileName_, namesOfEnergiesToMatch_);
223 auto testEnergyFrameReader =
224 openEnergyFileToReadTerms(runner_.edrFileName_, namesOfEnergiesToMatch_);
226 const bool shouldCompareEnergies = fromString<int>(output.at("nstenergy")) > 0;
227 bool shouldContinueComparing = shouldCompareEnergies;
228 while (shouldContinueComparing)
230 if (shouldCompareEnergies)
232 SCOPED_TRACE("Comparing energy frames from reference '" + referenceFileNames_.edrFileName_
233 + "' and test '" + runner_.edrFileName_ + "'");
234 shouldContinueComparing = shouldContinueComparing
235 && compareFrames(referenceEnergyFrameReader.get(),
236 testEnergyFrameReader.get(), energyComparison_);
242 /*! \brief Some common choices of periodic output mdp parameters to
243 * simplify defining values for the combinations under test */
244 PeriodicOutputParameters g_basicPeriodicOutputParameters = {
245 { "nstenergy", "0" }, { "nstlog", "0" }, { "nstxout", "0" },
246 { "nstvout", "0" }, { "nstfout", "0" }, { "nstxout-compressed", "0" },
247 { "nstdhdl", "0" }, { "description", "unknown" }
250 /*! \brief Return vector of mdp parameter sets to test
252 * These are constructed to observe the mdp parameter choices that
253 * only affect when output is written agree with those that were
254 * written from a reference run where output was done every step. The
255 * numbers are chosen in the context of the defaults in
256 * prepareDefaultMdpFieldValues().
258 * \todo Test nstlog, nstdhdl, nstxout-compressed */
259 std::vector<PeriodicOutputParameters> outputParameters()
261 std::vector<PeriodicOutputParameters> parameterSets;
263 parameterSets.push_back(g_basicPeriodicOutputParameters);
264 parameterSets.back()["nstenergy"] = "1";
265 parameterSets.back()["description"] = "energies every step works";
267 parameterSets.push_back(g_basicPeriodicOutputParameters);
268 parameterSets.back()["nstenergy"] = "4";
269 parameterSets.back()["description"] = "energies every fourth step works";
271 parameterSets.push_back(g_basicPeriodicOutputParameters);
272 parameterSets.back()["nstxout"] = "4";
273 parameterSets.back()["description"] = "coordinates every fourth step works";
275 parameterSets.push_back(g_basicPeriodicOutputParameters);
276 parameterSets.back()["nstvout"] = "4";
277 parameterSets.back()["description"] = "velocities every fourth step works";
279 parameterSets.push_back(g_basicPeriodicOutputParameters);
280 parameterSets.back()["nstfout"] = "4";
281 parameterSets.back()["description"] = "forces every fourth step works";
283 return parameterSets;
286 //! Returns sets of simple simulation propagators
287 std::vector<PropagationParameters> simplePropagationParameters()
290 { { "simulationName", "argon12" },
291 { "integrator", "md" },
292 { "comm-mode", "linear" },
295 { "nsttcouple", "0" },
297 { "nstpcouple", "0" },
298 { "maxGromppWarningsTolerated", "0" } },
302 /*! \brief Returns sets of simulation propagators including coupling
304 * These are chosen to cover the commonly used space of propagation
305 * algorithms togther with the perdiods between their actions. The
306 * periods tested are chosen to be mutually co-prime and distinct from
307 * the pair search and user output period (4), so that over the
308 * duration of a short simulation many kinds of simulation step
309 * behavior are tested. */
310 std::vector<PropagationParameters> propagationParametersWithCoupling()
312 std::string nsttcouple = "2";
313 std::string nstpcouple = "3";
314 std::string comm_mode = "linear";
315 std::string nstcomm = "5";
317 std::vector<PropagationParameters> parameterSets;
318 for (const std::string& simulationName : { "argon12" })
320 for (const std::string& integrator : { "md", "sd", "md-vv" })
322 for (const std::string& tcoupl : { "no", "v-rescale", "Nose-Hoover" })
324 // SD doesn't support temperature-coupling algorithms,
325 if (integrator == "sd" && tcoupl != "no")
329 for (std::string pcoupl : { "no", "Berendsen", "Parrinello-Rahman", "C-rescale" })
331 // VV supports few algorithm combinations
332 if (integrator == "md-vv")
334 // P-R with VV is called MTTK
335 if (pcoupl == "Parrinello-Rahman")
339 if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
340 || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
345 if (pcoupl == "MTTK" && simulationName != "argon12")
347 // MTTK does not support constraints
351 int maxGromppWarningsTolerated = 0;
352 if (pcoupl == "Berendsen")
354 ++maxGromppWarningsTolerated;
356 parameterSets.emplace_back(PropagationParameters{
357 { "simulationName", simulationName },
358 { "integrator", integrator },
359 { "comm-mode", comm_mode },
360 { "nstcomm", nstcomm },
361 { "tcoupl", tcoupl },
362 { "nsttcouple", nsttcouple },
363 { "pcoupl", pcoupl },
364 { "nstpcouple", nstpcouple },
365 { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
370 return parameterSets;
373 /*! \brief Returns sets of simulation propagators including coupling
375 * These are chosen to cover the commonly used space of propagation
376 * algorithms on systems with constraints. */
377 std::vector<PropagationParameters> propagationParametersWithConstraints()
379 std::string nsttcouple = "2";
380 std::string nstpcouple = "3";
381 std::string comm_mode = "linear";
382 std::string nstcomm = "5";
384 std::vector<PropagationParameters> parameterSets;
385 for (const std::string& simulationName : { "tip3p5" })
387 for (const std::string& integrator : { "md", "sd", "md-vv" })
389 for (const std::string& tcoupl : { "no", "v-rescale" })
391 // SD doesn't support temperature-coupling algorithms,
392 if (integrator == "sd" && tcoupl != "no")
396 for (std::string pcoupl : { "no", "Parrinello-Rahman", "C-rescale" })
398 // VV supports few algorithm combinations
399 if (integrator == "md-vv")
401 // P-R with VV is called MTTK
402 if (pcoupl == "Parrinello-Rahman")
406 if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
407 || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
412 if (pcoupl == "MTTK" && simulationName != "argon12")
414 // MTTK does not support constraints
418 int maxGromppWarningsTolerated = 0;
419 parameterSets.emplace_back(PropagationParameters{
420 { "simulationName", simulationName },
421 { "integrator", integrator },
422 { "comm-mode", comm_mode },
423 { "nstcomm", nstcomm },
424 { "tcoupl", tcoupl },
425 { "nsttcouple", nsttcouple },
426 { "pcoupl", pcoupl },
427 { "nstpcouple", nstpcouple },
428 { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
433 return parameterSets;
436 using ::testing::Combine;
437 using ::testing::Values;
438 using ::testing::ValuesIn;
440 // TODO The time for OpenCL kernel compilation means these tests time
441 // out. Once that compilation is cached for the whole process, these
442 // tests can run in such configurations.
444 INSTANTIATE_TEST_CASE_P(BasicPropagators,
446 Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
447 INSTANTIATE_TEST_CASE_P(PropagatorsWithCoupling,
449 Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
450 INSTANTIATE_TEST_CASE_P(PropagatorsWithConstraints,
452 Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));
454 INSTANTIATE_TEST_CASE_P(DISABLED_BasicPropagators,
456 Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
457 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithCoupling,
459 Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
460 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithConstraints,
462 Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));