src/programs/mdrun/tests/periodicactions.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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(),
 121                          propagation["integrator"].c_str(),
 122                          propagation["tcoupl"].c_str(),
 123                          propagation["pcoupl"].c_str()));
 124     auto mdpFieldValues = prepareMdpFieldValues(propagation["simulationName"],
 125                                                 propagation["integrator"],
 126                                                 propagation["tcoupl"],
 127                                                 propagation["pcoupl"]);
 128
 129     // This lambda writes all mdp options in `source` into `target`, overwriting options already
 130     // present in `target`. It also filters out non-mdp option entries in the source maps
 131     auto overWriteMdpMapValues = [](const MdpFieldValues& source, MdpFieldValues& target) {
 132         for (auto const& [key, value] : source)
 133         {
 134             if (key == "simulationName" || key == "maxGromppWarningsTolerated" || key == "description")
 135             {
 136                 // Remove non-mdp entries used in propagation and output
 137                 continue;
 138             }
 139             target[key] = value;
 140         }
 141     };
 142     // Add options in propagation and output to the mdp options
 143     overWriteMdpMapValues(propagation, mdpFieldValues);
 144     overWriteMdpMapValues(output, mdpFieldValues);
 145
 146     // prepare the tpr file
 147     {
 148         CommandLine caller;
 149         caller.append("grompp");
 150         caller.addOption("-maxwarn", propagation["maxGromppWarningsTolerated"]);
 151         runner_.useTopGroAndNdxFromDatabase(propagation["simulationName"]);
 152         runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 153         EXPECT_EQ(0, runner_.callGrompp(caller));
 154     }
 155     // do a normal mdrun
 156     {
 157         CommandLine mdrunCaller;
 158         mdrunCaller.append("mdrun");
 159         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
 160     }
 161 }
 162
 163 void PeriodicActionsTest::prepareReferenceData()
 164 {
 165     SCOPED_TRACE("Preparing reference data");
 166
 167     // Configure the SimulationRunner to write output to files we can
 168     // use in comparing each test run.
 169     std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
 170
 171     // Run the reference simulation with everything output every step
 172     PeriodicOutputParameters outputEveryStep = {
 173         { "nstenergy", "1" },
 174         { "nstlog", "1" },
 175         { "nstdhdl", "1" },
 176         { "description", "output everything every step" },
 177     };
 178     doMdrun(outputEveryStep);
 179
 180     // Restore the standard filenames we'll use for the runs whose behavior we are testing.
 181     std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
 182 }
 183
 184 /*! \brief Compare the next frame returned by \c testFrameReader with the
 185  * matching frame from \c referenceTestReader using \c comparator.
 186  *
 187  * \returns  True when a successful comparison was made
 188  */
 189 template<typename Reader, typename Comparator>
 190 bool compareFrames(Reader referenceFrameReader, Reader testFrameReader, const Comparator& comparator)
 191 {
 192     if (!testFrameReader->readNextFrame())
 193     {
 194         return false;
 195     }
 196
 197     auto        testFrame = testFrameReader->frame();
 198     std::string frameName = testFrame.frameName();
 199     SCOPED_TRACE("Found frame from test file named " + frameName);
 200     bool foundMatch = false;
 201     while (!foundMatch)
 202     {
 203         if (referenceFrameReader->readNextFrame())
 204         {
 205             auto referenceFrame = referenceFrameReader->frame();
 206             // Can the reference and test frames be compared, because they have the same name?
 207             if (frameName == referenceFrame.frameName())
 208             {
 209                 SCOPED_TRACE("Found frame from reference file named " + frameName);
 210                 comparator(referenceFrame, testFrame);
 211                 foundMatch = true;
 212             }
 213         }
 214         else
 215         {
 216             ADD_FAILURE() << "Ran out of reference frames to compare";
 217             return false;
 218         }
 219     }
 220     return foundMatch;
 221 }
 222
 223 TEST_P(PeriodicActionsTest, PeriodicActionsAgreeWithReference)
 224 {
 225     auto propagation = std::get<0>(GetParam());
 226     SCOPED_TRACE(formatString("Comparing two simulations of '%s' with integrator '%s'",
 227                               propagation["simulationName"].c_str(),
 228                               propagation["integrator"].c_str()));
 229
 230     prepareReferenceData();
 231
 232     auto outputParametersGenerator = std::get<1>(GetParam());
 233     for (const PeriodicOutputParameters& output : outputParametersGenerator())
 234     {
 235         SCOPED_TRACE("Comparing to observe " + output.at("description"));
 236
 237         // Run the test simulation
 238         doMdrun(output);
 239
 240         // Prepare readers for the reference and test output files
 241         auto referenceEnergyFrameReader =
 242                 openEnergyFileToReadTerms(referenceFileNames_.edrFileName_, namesOfEnergiesToMatch_);
 243         auto testEnergyFrameReader =
 244                 openEnergyFileToReadTerms(runner_.edrFileName_, namesOfEnergiesToMatch_);
 245
 246         const bool shouldCompareEnergies   = fromString<int>(output.at("nstenergy")) > 0;
 247         bool       shouldContinueComparing = shouldCompareEnergies;
 248         while (shouldContinueComparing)
 249         {
 250             if (shouldCompareEnergies)
 251             {
 252                 SCOPED_TRACE("Comparing energy frames from reference '" + referenceFileNames_.edrFileName_
 253                              + "' and test '" + runner_.edrFileName_ + "'");
 254                 shouldContinueComparing = shouldContinueComparing
 255                                           && compareFrames(referenceEnergyFrameReader.get(),
 256                                                            testEnergyFrameReader.get(),
 257                                                            energyComparison_);
 258             }
 259         }
 260     }
 261 }
 262
 263 /*! \brief Some common choices of periodic output mdp parameters to
 264  * simplify defining values for the combinations under test */
 265 PeriodicOutputParameters g_basicPeriodicOutputParameters = {
 266     { "nstenergy", "0" }, { "nstlog", "0" },           { "nstxout", "0" },
 267     { "nstvout", "0" },   { "nstfout", "0" },          { "nstxout-compressed", "0" },
 268     { "nstdhdl", "0" },   { "description", "unknown" }
 269 };
 270
 271 /*! \brief Return vector of mdp parameter sets to test
 272  *
 273  * These are constructed to observe the mdp parameter choices that
 274  * only affect when output is written agree with those that were
 275  * written from a reference run where output was done every step. The
 276  * numbers are chosen in the context of the defaults in
 277  * prepareDefaultMdpFieldValues().
 278  *
 279  * \todo Test nstlog, nstdhdl, nstxout-compressed */
 280 std::vector<PeriodicOutputParameters> outputParameters()
 281 {
 282     std::vector<PeriodicOutputParameters> parameterSets;
 283
 284     parameterSets.push_back(g_basicPeriodicOutputParameters);
 285     parameterSets.back()["nstenergy"]   = "1";
 286     parameterSets.back()["description"] = "energies every step works";
 287
 288     parameterSets.push_back(g_basicPeriodicOutputParameters);
 289     parameterSets.back()["nstenergy"]   = "4";
 290     parameterSets.back()["description"] = "energies every fourth step works";
 291
 292     parameterSets.push_back(g_basicPeriodicOutputParameters);
 293     parameterSets.back()["nstxout"]     = "4";
 294     parameterSets.back()["description"] = "coordinates every fourth step works";
 295
 296     parameterSets.push_back(g_basicPeriodicOutputParameters);
 297     parameterSets.back()["nstvout"]     = "4";
 298     parameterSets.back()["description"] = "velocities every fourth step works";
 299
 300     parameterSets.push_back(g_basicPeriodicOutputParameters);
 301     parameterSets.back()["nstfout"]     = "4";
 302     parameterSets.back()["description"] = "forces every fourth step works";
 303
 304     return parameterSets;
 305 }
 306
 307 //! Returns sets of simple simulation propagators
 308 std::vector<PropagationParameters> simplePropagationParameters()
 309 {
 310     return {
 311         { { "simulationName", "argon12" },
 312           { "integrator", "md" },
 313           { "comm-mode", "linear" },
 314           { "nstcomm", "1" },
 315           { "tcoupl", "no" },
 316           { "nsttcouple", "0" },
 317           { "pcoupl", "no" },
 318           { "nstpcouple", "0" },
 319           { "maxGromppWarningsTolerated", "0" } },
 320     };
 321 }
 322
 323 /*! \brief Returns sets of simulation propagators including coupling
 324  *
 325  * These are chosen to cover the commonly used space of propagation
 326  * algorithms togther with the perdiods between their actions. The
 327  * periods tested are chosen to be mutually co-prime and distinct from
 328  * the pair search and user output period (4), so that over the
 329  * duration of a short simulation many kinds of simulation step
 330  * behavior are tested. */
 331 std::vector<PropagationParameters> propagationParametersWithCoupling()
 332 {
 333     std::string nsttcouple = "2";
 334     std::string nstpcouple = "3";
 335     std::string comm_mode  = "linear";
 336     std::string nstcomm    = "5";
 337
 338     std::vector<PropagationParameters> parameterSets;
 339     for (const std::string& simulationName : { "argon12" })
 340     {
 341         for (const std::string& integrator : { "md", "sd", "md-vv" })
 342         {
 343             for (const std::string& tcoupl : { "no", "v-rescale", "Nose-Hoover" })
 344             {
 345                 // SD doesn't support temperature-coupling algorithms,
 346                 if (integrator == "sd" && tcoupl != "no")
 347                 {
 348                     continue;
 349                 }
 350                 for (std::string pcoupl : { "no", "Berendsen", "Parrinello-Rahman", "C-rescale" })
 351                 {
 352                     // VV supports few algorithm combinations
 353                     if (integrator == "md-vv")
 354                     {
 355                         // P-R with VV is called MTTK
 356                         if (pcoupl == "Parrinello-Rahman")
 357                         {
 358                             pcoupl = "MTTK";
 359                         }
 360                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 361                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 362                         {
 363                             continue;
 364                         }
 365                     }
 366                     if (pcoupl == "MTTK" && simulationName != "argon12")
 367                     {
 368                         // MTTK does not support constraints
 369                         continue;
 370                     }
 371
 372                     int maxGromppWarningsTolerated = 0;
 373                     if (pcoupl == "Berendsen")
 374                     {
 375                         ++maxGromppWarningsTolerated;
 376                     }
 377                     parameterSets.emplace_back(PropagationParameters{
 378                             { "simulationName", simulationName },
 379                             { "integrator", integrator },
 380                             { "comm-mode", comm_mode },
 381                             { "nstcomm", nstcomm },
 382                             { "tcoupl", tcoupl },
 383                             { "nsttcouple", nsttcouple },
 384                             { "pcoupl", pcoupl },
 385                             { "nstpcouple", nstpcouple },
 386                             { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
 387                 }
 388             }
 389         }
 390     }
 391     return parameterSets;
 392 }
 393
 394 /*! \brief Returns sets of simulation propagators including coupling
 395  *
 396  * These are chosen to cover the commonly used space of propagation
 397  * algorithms on systems with constraints. */
 398 std::vector<PropagationParameters> propagationParametersWithConstraints()
 399 {
 400     std::string nsttcouple = "2";
 401     std::string nstpcouple = "3";
 402     std::string comm_mode  = "linear";
 403     std::string nstcomm    = "5";
 404
 405     std::vector<PropagationParameters> parameterSets;
 406     for (const std::string& simulationName : { "tip3p5" })
 407     {
 408         for (const std::string& integrator : { "md", "sd", "md-vv" })
 409         {
 410             for (const std::string& tcoupl : { "no", "v-rescale" })
 411             {
 412                 // SD doesn't support temperature-coupling algorithms,
 413                 if (integrator == "sd" && tcoupl != "no")
 414                 {
 415                     continue;
 416                 }
 417                 for (std::string pcoupl : { "no", "Parrinello-Rahman", "C-rescale" })
 418                 {
 419                     // VV supports few algorithm combinations
 420                     if (integrator == "md-vv")
 421                     {
 422                         // P-R with VV is called MTTK
 423                         if (pcoupl == "Parrinello-Rahman")
 424                         {
 425                             pcoupl = "MTTK";
 426                         }
 427                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 428                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 429                         {
 430                             continue;
 431                         }
 432                     }
 433                     if (pcoupl == "MTTK" && simulationName != "argon12")
 434                     {
 435                         // MTTK does not support constraints
 436                         continue;
 437                     }
 438
 439                     int maxGromppWarningsTolerated = 0;
 440                     parameterSets.emplace_back(PropagationParameters{
 441                             { "simulationName", simulationName },
 442                             { "integrator", integrator },
 443                             { "comm-mode", comm_mode },
 444                             { "nstcomm", nstcomm },
 445                             { "tcoupl", tcoupl },
 446                             { "nsttcouple", nsttcouple },
 447                             { "pcoupl", pcoupl },
 448                             { "nstpcouple", nstpcouple },
 449                             { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
 450                 }
 451             }
 452         }
 453     }
 454     return parameterSets;
 455 }
 456
 457 using ::testing::Combine;
 458 using ::testing::Values;
 459 using ::testing::ValuesIn;
 460
 461 // TODO The time for OpenCL kernel compilation means these tests time
 462 // out. Once that compilation is cached for the whole process, these
 463 // tests can run in such configurations.
 464 #if !GMX_GPU_OPENCL
 465 INSTANTIATE_TEST_CASE_P(BasicPropagators,
 466                         PeriodicActionsTest,
 467                         Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
 468 INSTANTIATE_TEST_CASE_P(PropagatorsWithCoupling,
 469                         PeriodicActionsTest,
 470                         Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
 471 INSTANTIATE_TEST_CASE_P(PropagatorsWithConstraints,
 472                         PeriodicActionsTest,
 473                         Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));
 474 #else
 475 INSTANTIATE_TEST_CASE_P(DISABLED_BasicPropagators,
 476                         PeriodicActionsTest,
 477                         Combine(ValuesIn(simplePropagationParameters()), Values(outputParameters)));
 478 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithCoupling,
 479                         PeriodicActionsTest,
 480                         Combine(ValuesIn(propagationParametersWithCoupling()), Values(outputParameters)));
 481 INSTANTIATE_TEST_CASE_P(DISABLED_PropagatorsWithConstraints,
 482                         PeriodicActionsTest,
 483                         Combine(ValuesIn(propagationParametersWithConstraints()), Values(outputParameters)));
 484 #endif
 485
 486 } // namespace
 487 } // namespace test
 488 } // namespace gmx