src/programs/mdrun/tests/periodicactions.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   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
  36 /*! \internal \file
  37  * \brief Tests to verify that a simulator that only does some actions
  38  * periodically produces the expected results.
  39  *
  40  * \author Mark Abraham <mark.j.abraham@gmail.com>
  41  * \ingroup module_mdrun_integration_tests
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "config.h"
  46
  47 #include <tuple>
  48
  49 #include "gromacs/trajectory/energyframe.h"
  50 #include "gromacs/utility/strconvert.h"
  51 #include "gromacs/utility/stringutil.h"
  52
  53 #include "testutils/simulationdatabase.h"
  54 #include "testutils/testasserts.h"
  55
  56 #include "energycomparison.h"
  57 #include "energyreader.h"
  58 #include "moduletest.h"
  59
  60 namespace gmx
  61 {
  62 namespace test
  63 {
  64 namespace
  65 {
  66
  67 /*! \brief Mdp parameters that determine the manner of simulation
  68  * propagation. */
  69 using PropagationParameters = MdpFieldValues;
  70
  71 /*! \brief Mdp parameters that should only affect the observations,
  72  *  not the simulation propagation. */
  73 using PeriodicOutputParameters = MdpFieldValues;
  74
  75 //! Helper type of output file names for the reference mdrun call
  76 struct ReferenceFileNames
  77 {
  78     //! Name of energy file
  79     std::string edrFileName_;
  80 };
  81
  82 //! Function type to produce sets of .mdp parameters for testing periodic output
  83 using OutputParameterGeneratorFunction = std::vector<PeriodicOutputParameters> (*)();
  84
  85 /*! \brief Test fixture base for comparing a simulator with one that
  86  * does everything every step
  87  *
  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.
  91  *
  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.
  95  */
  96 class PeriodicActionsTest :
  97     public MdrunTestFixture,
  98     public ::testing::WithParamInterface<std::tuple<PropagationParameters, OutputParameterGeneratorFunction>>
  99 {
 100 public:
 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();
 113 };
 114
 115 void PeriodicActionsTest::doMdrun(const PeriodicOutputParameters& output)
 116 {
 117     auto propagation = std::get<0>(GetParam());
 118     SCOPED_TRACE(
 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());
 126
 127     // prepare the tpr file
 128     {
 129         CommandLine caller;
 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));
 135     }
 136     // do a normal mdrun
 137     {
 138         CommandLine mdrunCaller;
 139         mdrunCaller.append("mdrun");
 140         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
 141     }
 142 }
 143
 144 void PeriodicActionsTest::prepareReferenceData()
 145 {
 146     SCOPED_TRACE("Preparing reference data");
 147
 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_);
 151
 152     // Run the reference simulation with everything output every step
 153     PeriodicOutputParameters outputEveryStep = {
 154         { "nstenergy", "1" },
 155         { "nstlog", "1" },
 156         { "nstdhdl", "1" },
 157         { "description", "output everything every step" },
 158     };
 159     doMdrun(outputEveryStep);
 160
 161     // Restore the standard filenames we'll use for the runs whose behavior we are testing.
 162     std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
 163 }
 164
 165 /*! \brief Compare the next frame returned by \c testFrameReader with the
 166  * matching frame from \c referenceTestReader using \c comparator.
 167  *
 168  * \returns  True when a successful comparison was made
 169  */
 170 template<typename Reader, typename Comparator>
 171 bool compareFrames(Reader referenceFrameReader, Reader testFrameReader, const Comparator& comparator)
 172 {
 173     if (!testFrameReader->readNextFrame())
 174     {
 175         return false;
 176     }
 177
 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;
 182     while (!foundMatch)
 183     {
 184         if (referenceFrameReader->readNextFrame())
 185         {
 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())
 189             {
 190                 SCOPED_TRACE("Found frame from reference file named " + frameName);
 191                 comparator(referenceFrame, testFrame);
 192                 foundMatch = true;
 193             }
 194         }
 195         else
 196         {
 197             ADD_FAILURE() << "Ran out of reference frames to compare";
 198             return false;
 199         }
 200     }
 201     return foundMatch;
 202 }
 203
 204 TEST_P(PeriodicActionsTest, PeriodicActionsAgreeWithReference)
 205 {
 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()));
 209
 210     prepareReferenceData();
 211
 212     auto outputParametersGenerator = std::get<1>(GetParam());
 213     for (const PeriodicOutputParameters& output : outputParametersGenerator())
 214     {
 215         SCOPED_TRACE("Comparing to observe " + output.at("description"));
 216
 217         // Run the test simulation
 218         doMdrun(output);
 219
 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_);
 225
 226         const bool shouldCompareEnergies   = fromString<int>(output.at("nstenergy")) > 0;
 227         bool       shouldContinueComparing = shouldCompareEnergies;
 228         while (shouldContinueComparing)
 229         {
 230             if (shouldCompareEnergies)
 231             {
 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_);
 237             }
 238         }
 239     }
 240 }
 241
 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" }
 248 };
 249
 250 /*! \brief Return vector of mdp parameter sets to test
 251  *
 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().
 257  *
 258  * \todo Test nstlog, nstdhdl, nstxout-compressed */
 259 std::vector<PeriodicOutputParameters> outputParameters()
 260 {
 261     std::vector<PeriodicOutputParameters> parameterSets;
 262
 263     parameterSets.push_back(g_basicPeriodicOutputParameters);
 264     parameterSets.back()["nstenergy"]   = "1";
 265     parameterSets.back()["description"] = "energies every step works";
 266
 267     parameterSets.push_back(g_basicPeriodicOutputParameters);
 268     parameterSets.back()["nstenergy"]   = "4";
 269     parameterSets.back()["description"] = "energies every fourth step works";
 270
 271     parameterSets.push_back(g_basicPeriodicOutputParameters);
 272     parameterSets.back()["nstxout"]     = "4";
 273     parameterSets.back()["description"] = "coordinates every fourth step works";
 274
 275     parameterSets.push_back(g_basicPeriodicOutputParameters);
 276     parameterSets.back()["nstvout"]     = "4";
 277     parameterSets.back()["description"] = "velocities every fourth step works";
 278
 279     parameterSets.push_back(g_basicPeriodicOutputParameters);
 280     parameterSets.back()["nstfout"]     = "4";
 281     parameterSets.back()["description"] = "forces every fourth step works";
 282
 283     return parameterSets;
 284 }
 285
 286 //! Returns sets of simple simulation propagators
 287 std::vector<PropagationParameters> simplePropagationParameters()
 288 {
 289     return {
 290         { { "simulationName", "argon12" },
 291           { "integrator", "md" },
 292           { "comm-mode", "linear" },
 293           { "nstcomm", "1" },
 294           { "tcoupl", "no" },
 295           { "nsttcouple", "0" },
 296           { "pcoupl", "no" },
 297           { "nstpcouple", "0" },
 298           { "maxGromppWarningsTolerated", "0" } },
 299     };
 300 }
 301
 302 /*! \brief Returns sets of simulation propagators including coupling
 303  *
 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()
 311 {
 312     std::string nsttcouple = "2";
 313     std::string nstpcouple = "3";
 314     std::string comm_mode  = "linear";
 315     std::string nstcomm    = "5";
 316
 317     std::vector<PropagationParameters> parameterSets;
 318     for (const std::string& simulationName : { "argon12" })
 319     {
 320         for (const std::string& integrator : { "md", "sd", "md-vv" })
 321         {
 322             for (const std::string& tcoupl : { "no", "v-rescale", "Nose-Hoover" })
 323             {
 324                 // SD doesn't support temperature-coupling algorithms,
 325                 if (integrator == "sd" && tcoupl != "no")
 326                 {
 327                     continue;
 328                 }
 329                 for (std::string pcoupl : { "no", "Berendsen", "Parrinello-Rahman" })
 330                 {
 331                     // VV supports few algorithm combinations
 332                     if (integrator == "md-vv")
 333                     {
 334                         // P-R with VV is called MTTK
 335                         if (pcoupl == "Parrinello-Rahman")
 336                         {
 337                             pcoupl = "MTTK";
 338                         }
 339                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 340                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 341                         {
 342                             continue;
 343                         }
 344                     }
 345                     if (pcoupl == "MTTK" && simulationName != "argon12")
 346                     {
 347                         // MTTK does not support constraints
 348                         continue;
 349                     }
 350
 351                     int maxGromppWarningsTolerated = 0;
 352                     if (pcoupl == "Berendsen")
 353                     {
 354                         ++maxGromppWarningsTolerated;
 355                     }
 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) } });
 366                 }
 367             }
 368         }
 369     }
 370     return parameterSets;
 371 }
 372
 373 /*! \brief Returns sets of simulation propagators including coupling
 374  *
 375  * These are chosen to cover the commonly used space of propagation
 376  * algorithms on systems with constraints. */
 377 std::vector<PropagationParameters> propagationParametersWithConstraints()
 378 {
 379     std::string nsttcouple = "2";
 380     std::string nstpcouple = "3";
 381     std::string comm_mode  = "linear";
 382     std::string nstcomm    = "5";
 383
 384     std::vector<PropagationParameters> parameterSets;
 385     for (const std::string& simulationName : { "tip3p5" })
 386     {
 387         for (const std::string& integrator : { "md", "sd", "md-vv" })
 388         {
 389             for (const std::string& tcoupl : { "no", "v-rescale" })
 390             {
 391                 // SD doesn't support temperature-coupling algorithms,
 392                 if (integrator == "sd" && tcoupl != "no")
 393                 {
 394                     continue;
 395                 }
 396                 for (std::string pcoupl : { "no", "Parrinello-Rahman" })
 397                 {
 398                     // VV supports few algorithm combinations
 399                     if (integrator == "md-vv")
 400                     {
 401                         // P-R with VV is called MTTK
 402                         if (pcoupl == "Parrinello-Rahman")
 403                         {
 404                             pcoupl = "MTTK";
 405                         }
 406                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 407                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 408                         {
 409                             continue;
 410                         }
 411                     }
 412                     if (pcoupl == "MTTK" && simulationName != "argon12")
 413                     {
 414                         // MTTK does not support constraints
 415                         continue;
 416                     }
 417
 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) } });
 429                 }
 430             }
 431         }
 432     }
 433     return parameterSets;
 434 }
 435
 436 using ::testing::Combine;
 437 using ::testing::Values;
 438 using ::testing::ValuesIn;
 439
 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.
 443 #if !GMX_GPU_OPENCL
 444 INSTANTIATE_TEST_CASE_P(BasicPropagators,
 445                         PeriodicActionsTest,
 446                         Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
 447 INSTANTIATE_TEST_CASE_P(PropagatorsWithCoupling,
 448                         PeriodicActionsTest,
 449                         Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
 450 INSTANTIATE_TEST_CASE_P(PropagatorsWithConstraints,
 451                         PeriodicActionsTest,
 452                         Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));
 453 #else
 454 INSTANTIATE_TEST_CASE_P(DISABLED_BasicPropagators,
 455                         PeriodicActionsTest,
 456                         Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
 457 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithCoupling,
 458                         PeriodicActionsTest,
 459                         Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
 460 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithConstraints,
 461                         PeriodicActionsTest,
 462                         Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));
 463 #endif
 464
 465 } // namespace
 466 } // namespace test
 467 } // namespace gmx