2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 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.
37 * \brief Tests that mdrun restarts are exact, that is that two
38 * successive runs without appending reproduce a single-part run.
40 * \todo Extend the coverage to the appending case.
42 * \author Mark Abraham <mark.j.abraham@gmail.com>
43 * \ingroup module_mdrun_integration_tests
52 #include <gtest/gtest.h>
54 #include "gromacs/topology/idef.h"
55 #include "gromacs/trajectory/energyframe.h"
56 #include "gromacs/utility/stringutil.h"
58 #include "testutils/mpitest.h"
59 #include "testutils/simulationdatabase.h"
60 #include "testutils/testasserts.h"
62 #include "energycomparison.h"
63 #include "energyreader.h"
64 #include "moduletest.h"
74 /*! \brief Manages comparing a pair of matching energy frames from a
75 * normal run and the same run in two parts.
77 * The passed frame readers must contain logically equivalent
78 * contents, with no extra or missing frames from either the
79 * one- or two-part run.
81 * \todo This class has a similar interface with one in
82 * mdcomparison.h, but not enough to warrant extracting an interface
83 * class. Perhaps parts of this should be cast as a class that returns
84 * iterators to successive frames, so that the fancy footwork for
85 * pretending a two-part run is one sequence is separate from the
86 * comparison of that sequene with that of the one-part run?
88 * \tparam FrameReader Has readNextFrame() and frame() methods
89 * useful for returning successive Frame objects */
90 template<class FrameReader>
91 class ContinuationFramePairManager
94 //! Convenience typedef
95 typedef std::unique_ptr<FrameReader> FrameReaderPtr;
97 ContinuationFramePairManager(FrameReaderPtr full, FrameReaderPtr firstPart, FrameReaderPtr secondPart) :
98 full_(std::move(full)),
99 firstPart_(std::move(firstPart)),
100 secondPart_(std::move(secondPart)),
104 /*! \brief Probe for a pair of valid frames, and return true if both are found.
106 * Gives a test failure if exactly one frame is found, because
107 * it is an error for either run to have missing or extra
108 * frames. Note that the frame where the two-part run ends
109 * and begins is duplicated between the two output files by
110 * mdrun, and the test accommodates this.
112 * \todo This would be straightforward if velocity Verlet
113 * behaved like other integrators. */
114 bool shouldContinueComparing()
116 if (full_->readNextFrame())
120 if (firstPart_->readNextFrame())
122 // Two valid next frames exist, so we should continue comparing.
127 // First part ran out of frames, move on to the second part
128 isFirstPart_ = false;
129 if (secondPart_->readNextFrame())
131 // Skip a second-part frame so the one we will
132 // read can compare with the next full-run
134 secondPart_->frame();
135 if (secondPart_->readNextFrame())
137 // Two valid next frames exist, so we should continue comparing.
142 ADD_FAILURE() << "Second-part energy file had no (new) frames";
147 ADD_FAILURE() << "Second-part energy file had no frames";
153 if (secondPart_->readNextFrame())
155 // Two valid next frames exist, so we should continue comparing.
160 ADD_FAILURE() << "Full run energy file had at least one more frame than "
161 "two-part run energy file";
169 ADD_FAILURE() << "Full-run energy file ran out of frames before the first part of "
170 "the two-part run completed";
174 if (secondPart_->readNextFrame())
176 ADD_FAILURE() << "Two-part run energy file had at least one more frame than "
177 "full-run energy file";
181 // Both files ran out of frames at the same time, which is the expected behaviour.
185 // At least one file is out of frames, so should not continue comparing.
188 /*! \brief Compare all possible pairs of frames using \c compareTwoFrames.
190 * \tparam Frame The type of frame used in the comparison (returned
191 * by FrameReader and used by compareTwoFrames). */
192 template<class Frame>
193 void compareAllFramePairs(std::function<void(const Frame&, const Frame&)> compareTwoFrames)
195 while (shouldContinueComparing())
197 EnergyFrame firstFrame = full_->frame();
198 EnergyFrame secondFrame = isFirstPart_ ? firstPart_->frame() : secondPart_->frame();
199 SCOPED_TRACE("Comparing frames from two runs '" + firstFrame.frameName() + "' and '"
200 + secondFrame.frameName() + "'");
201 compareTwoFrames(firstFrame, secondFrame);
206 EnergyFrameReaderPtr full_;
207 EnergyFrameReaderPtr firstPart_;
208 EnergyFrameReaderPtr secondPart_;
212 /*! \brief Run grompp for a normal mdrun, the same mdrun stopping part
213 * way, doing a continuation, and compare the results. */
214 void runTest(TestFileManager* fileManager,
215 SimulationRunner* runner,
216 const std::string& simulationName,
217 int maxWarningsTolerated,
218 const MdpFieldValues& mdpFieldValues,
219 const EnergyTermsToCompare& energyTermsToCompare)
221 int numRanksAvailable = getNumberOfTestMpiRanks();
222 if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
225 "Test system '%s' cannot run with %d ranks.\n"
226 "The supported numbers are: %s\n",
227 simulationName.c_str(), numRanksAvailable,
228 reportNumbersOfPpRanksSupported(simulationName).c_str());
232 // prepare some names for files to use with the two mdrun calls
233 std::string fullRunTprFileName = fileManager->getTemporaryFilePath("full.tpr");
234 std::string firstPartRunTprFileName = fileManager->getTemporaryFilePath("firstpart.tpr");
235 std::string fullRunEdrFileName = fileManager->getTemporaryFilePath("full.edr");
236 std::string firstPartRunEdrFileName = fileManager->getTemporaryFilePath("firstpart.edr");
237 std::string firstPartRunCheckpointFileName = fileManager->getTemporaryFilePath("firstpart.cpt");
238 std::string secondPartRunEdrFileName = fileManager->getTemporaryFilePath("secondpart");
240 // prepare the full run .tpr file, which will be used for the full
241 // run, and for the second part of the two-part run.
243 // TODO evolve grompp to report the number of warnings issued, so
244 // tests always expect the right number.
246 caller.append("grompp");
247 caller.addOption("-maxwarn", maxWarningsTolerated);
248 runner->useTopGroAndNdxFromDatabase(simulationName);
249 runner->useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
250 runner->tprFileName_ = fullRunTprFileName;
251 EXPECT_EQ(0, runner->callGrompp(caller));
254 // prepare the .tpr file for the first part of the two-part run
256 // TODO evolve grompp to report the number of warnings issued, so
257 // tests always expect the right number.
259 caller.append("grompp");
260 caller.addOption("-maxwarn", maxWarningsTolerated);
261 runner->useTopGroAndNdxFromDatabase(simulationName);
262 auto firstPartMdpFieldValues = mdpFieldValues;
263 firstPartMdpFieldValues["nsteps"] = std::to_string(std::stoi(mdpFieldValues.at("nsteps")) / 2);
264 runner->useStringAsMdpFile(prepareMdpFileContents(firstPartMdpFieldValues));
265 runner->tprFileName_ = firstPartRunTprFileName;
266 EXPECT_EQ(0, runner->callGrompp(caller));
271 runner->tprFileName_ = fullRunTprFileName;
272 runner->edrFileName_ = fullRunEdrFileName;
273 CommandLine fullRunCaller;
274 fullRunCaller.append("mdrun");
275 ASSERT_EQ(0, runner->callMdrun(fullRunCaller));
278 // do a repeat of the first part of the same mdrun
280 runner->tprFileName_ = firstPartRunTprFileName;
281 runner->edrFileName_ = firstPartRunEdrFileName;
282 CommandLine firstPartRunCaller;
283 firstPartRunCaller.append("mdrun");
284 firstPartRunCaller.addOption("-cpo", firstPartRunCheckpointFileName);
285 ASSERT_EQ(0, runner->callMdrun(firstPartRunCaller));
288 // do a continuation (without appending) from the first part of
291 runner->tprFileName_ = fullRunTprFileName;
292 runner->edrFileName_ = secondPartRunEdrFileName;
293 CommandLine secondPartRunCaller;
294 secondPartRunCaller.append("mdrun");
295 // TODO We could test with appending but it would need a
296 // different implementation.
297 secondPartRunCaller.append("-noappend");
298 secondPartRunCaller.addOption("-cpi", firstPartRunCheckpointFileName);
299 ASSERT_EQ(0, runner->callMdrun(secondPartRunCaller));
300 // Cope with how -noappend works
301 secondPartRunEdrFileName += ".part0002.edr";
304 // Build the functor that will compare energy frames on the chosen
306 EnergyComparison energyComparison(energyTermsToCompare);
308 // Build the manager that will present matching pairs of frames to compare.
310 // TODO Here is an unnecessary copy of keys (ie. the energy term
311 // names), for convenience. In the future, use a range.
312 auto namesOfEnergiesToMatch = energyComparison.getEnergyNames();
313 ContinuationFramePairManager<EnergyFrameReader> energyManager(
314 openEnergyFileToReadTerms(fullRunEdrFileName, namesOfEnergiesToMatch),
315 openEnergyFileToReadTerms(firstPartRunEdrFileName, namesOfEnergiesToMatch),
316 openEnergyFileToReadTerms(secondPartRunEdrFileName, namesOfEnergiesToMatch));
317 // Compare the energy frames.
318 energyManager.compareAllFramePairs<EnergyFrame>(energyComparison);
321 /*! \brief Test fixture for mdrun exact continuations
323 * This test ensures mdrun can run a simulation, writing a trajectory
324 * and matching energies, and reproduce to within a tolerance the same
325 * energies from runs that stopped part of the way, and restarted from
328 * \todo Is there value in testing with mdrun -reprod? As well as
331 * \todo Add FEP case. */
332 class MdrunNoAppendContinuationIsExact :
333 public MdrunTestFixture,
334 public ::testing::WithParamInterface<std::tuple<std::string, std::string, std::string, std::string>>
338 MdrunNoAppendContinuationIsExact() {}
341 /* Listing all of these is tedious, but there's no other way to get a
342 * usefully descriptive string into the test-case name, so that when
343 * one breaks we can find out which one is broken without referring to
344 * this source code file.
346 * NOTE The choices for the tolerances are arbitrary but sufficient
347 * for comparison of runs to work on different hardware, and kinds and
348 * degrees parallelism. */
350 TEST_P(MdrunNoAppendContinuationIsExact, WithinTolerances)
352 auto params = GetParam();
353 auto simulationName = std::get<0>(params);
354 auto integrator = std::get<1>(params);
355 auto temperatureCoupling = std::get<2>(params);
356 auto pressureCoupling = std::get<3>(params);
358 // Check for unimplemented functionality
359 // TODO: Update this as modular simulator gains functionality
360 const bool isModularSimulatorExplicitlyDisabled = (getenv("GMX_DISABLE_MODULAR_SIMULATOR") != nullptr);
361 const bool isTCouplingCompatibleWithModularSimulator =
362 (temperatureCoupling == "no" || temperatureCoupling == "v-rescale");
363 if (integrator == "md-vv" && pressureCoupling == "parrinello-rahman"
364 && (isModularSimulatorExplicitlyDisabled || !isTCouplingCompatibleWithModularSimulator))
366 // Under md-vv, Parrinello-Rahman is only implemented for the modular simulator
369 if (integrator == "md-vv" && temperatureCoupling == "nose-hoover"
370 && pressureCoupling == "berendsen")
372 // This combination is not implemented in either legacy or modular simulator
377 formatString("Comparing normal and two-part run of simulation '%s' "
378 "with integrator '%s'",
379 simulationName.c_str(), integrator.c_str()));
381 auto mdpFieldValues = prepareMdpFieldValues(simulationName.c_str(), integrator.c_str(),
382 temperatureCoupling.c_str(), pressureCoupling.c_str());
383 // The exact lambda state choice is unimportant, so long as there
384 // is one when using an FEP input.
385 mdpFieldValues["other"] += formatString("\ninit-lambda-state = %d", 3);
386 mdpFieldValues["nsteps"] = "16";
388 // Forces on GPUs are generally not reproducible enough for a tight
389 // tolerance. Similarly, the propagation of sd and bd are not as
390 // reproducible as the others. So we use several ULP tolerance
391 // in all cases. This is looser than needed e.g. for md and md-vv
392 // with forces on CPUs, but there is no real risk of a bug with
393 // those propagators that would only be caught with a tighter
394 // tolerance in this particular test.
395 int ulpToleranceInMixed = 32;
396 int ulpToleranceInDouble = 64;
397 if (integrator == "bd")
399 // Somehow, the bd integrator has never been as reproducible
400 // as the others, either in continuations or reruns.
401 ulpToleranceInMixed = 200;
403 EnergyTermsToCompare energyTermsToCompare{
404 { { interaction_function[F_EPOT].longname,
405 relativeToleranceAsPrecisionDependentUlp(10.0, ulpToleranceInMixed, ulpToleranceInDouble) },
406 { interaction_function[F_EKIN].longname,
407 relativeToleranceAsPrecisionDependentUlp(10.0, ulpToleranceInMixed, ulpToleranceInDouble) } }
410 if (temperatureCoupling != "no" || pressureCoupling != "no")
412 energyTermsToCompare.insert({ interaction_function[F_ECONSERVED].longname,
413 relativeToleranceAsPrecisionDependentUlp(
414 10.0, ulpToleranceInMixed, ulpToleranceInDouble) });
417 if (pressureCoupling == "parrinello-rahman")
419 energyTermsToCompare.insert(
420 { "Box-Vel-XX", relativeToleranceAsPrecisionDependentUlp(1e-12, ulpToleranceInMixed,
421 ulpToleranceInDouble) });
422 energyTermsToCompare.insert(
423 { "Box-Vel-YY", relativeToleranceAsPrecisionDependentUlp(1e-12, ulpToleranceInMixed,
424 ulpToleranceInDouble) });
425 energyTermsToCompare.insert(
426 { "Box-Vel-ZZ", relativeToleranceAsPrecisionDependentUlp(1e-12, ulpToleranceInMixed,
427 ulpToleranceInDouble) });
430 int numWarningsToTolerate = 1;
431 runTest(&fileManager_, &runner_, simulationName, numWarningsToTolerate, mdpFieldValues,
432 energyTermsToCompare);
435 // TODO The time for OpenCL kernel compilation means these tests time
436 // out. Once that compilation is cached for the whole process, these
437 // tests can run in such configurations.
440 INSTANTIATE_TEST_CASE_P(
442 MdrunNoAppendContinuationIsExact,
443 ::testing::Combine(::testing::Values("argon12", "spc2", "alanine_vsite_vacuo"),
444 ::testing::Values("md", "md-vv", "bd", "sd"),
445 ::testing::Values("no"),
446 ::testing::Values("no")));
448 INSTANTIATE_TEST_CASE_P(NormalIntegratorsWithFEP,
449 MdrunNoAppendContinuationIsExact,
450 ::testing::Combine(::testing::Values("nonanol_vacuo"),
451 ::testing::Values("md", "md-vv", "bd", "sd"),
452 ::testing::Values("no"),
453 ::testing::Values("no")));
455 INSTANTIATE_TEST_CASE_P(
457 MdrunNoAppendContinuationIsExact,
458 ::testing::Combine(::testing::Values("argon12"),
459 ::testing::Values("md", "md-vv"),
460 ::testing::Values("berendsen", "v-rescale", "nose-hoover"),
461 ::testing::Values("no")));
463 INSTANTIATE_TEST_CASE_P(NPH,
464 MdrunNoAppendContinuationIsExact,
465 ::testing::Combine(::testing::Values("argon12"),
466 ::testing::Values("md", "md-vv"),
467 ::testing::Values("no"),
468 ::testing::Values("berendsen", "parrinello-rahman")));
470 INSTANTIATE_TEST_CASE_P(
472 MdrunNoAppendContinuationIsExact,
473 ::testing::Combine(::testing::Values("argon12"),
474 ::testing::Values("md", "md-vv"),
475 ::testing::Values("berendsen", "v-rescale", "nose-hoover"),
476 ::testing::Values("berendsen", "parrinello-rahman")));
478 INSTANTIATE_TEST_CASE_P(MTTK,
479 MdrunNoAppendContinuationIsExact,
480 ::testing::Combine(::testing::Values("argon12"),
481 ::testing::Values("md-vv"),
482 ::testing::Values("nose-hoover"),
483 ::testing::Values("mttk")));