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 Utility functions for 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 "periodicactions.h"
  48
  49 namespace gmx
  50 {
  51 namespace test
  52 {
  53
  54 /*! \brief Mdp parameters that determine the manner of simulation
  55  * propagation. */
  56 using PropagationParameters = MdpFieldValues;
  57
  58 /*! \brief Mdp parameters that should only affect the observations,
  59  *  not the simulation propagation. */
  60 using PeriodicOutputParameters = MdpFieldValues;
  61
  62 //! Function type to produce sets of .mdp parameters for testing periodic output
  63 using OutputParameterGeneratorFunction = std::vector<PeriodicOutputParameters> (*)();
  64
  65 void PeriodicActionsTest::doMdrun(const PeriodicOutputParameters& output)
  66 {
  67     auto propagation = std::get<0>(GetParam());
  68     SCOPED_TRACE(
  69             formatString("Doing %s simulation with %s integrator, %s tcoupling and %s pcoupling\n",
  70                          propagation["simulationName"].c_str(),
  71                          propagation["integrator"].c_str(),
  72                          propagation["tcoupl"].c_str(),
  73                          propagation["pcoupl"].c_str()));
  74     auto mdpFieldValues = prepareMdpFieldValues(propagation["simulationName"],
  75                                                 propagation["integrator"],
  76                                                 propagation["tcoupl"],
  77                                                 propagation["pcoupl"]);
  78
  79     // This lambda writes all mdp options in `source` into `target`, overwriting options already
  80     // present in `target`. It also filters out non-mdp option entries in the source maps
  81     auto overWriteMdpMapValues = [](const MdpFieldValues& source, MdpFieldValues& target) {
  82         for (auto const& [key, value] : source)
  83         {
  84             if (key == "simulationName" || key == "maxGromppWarningsTolerated" || key == "description")
  85             {
  86                 // Remove non-mdp entries used in propagation and output
  87                 continue;
  88             }
  89             target[key] = value;
  90         }
  91     };
  92     // Add options in propagation and output to the mdp options
  93     overWriteMdpMapValues(propagation, mdpFieldValues);
  94     overWriteMdpMapValues(output, mdpFieldValues);
  95
  96     // prepare the tpr file
  97     {
  98         CommandLine caller;
  99         caller.append("grompp");
 100         caller.addOption("-maxwarn", propagation["maxGromppWarningsTolerated"]);
 101         runner_.useTopGroAndNdxFromDatabase(propagation["simulationName"]);
 102         runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 103         EXPECT_EQ(0, runner_.callGrompp(caller));
 104     }
 105     // do a normal mdrun
 106     {
 107         CommandLine mdrunCaller;
 108         mdrunCaller.append("mdrun");
 109         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
 110     }
 111 }
 112
 113 void PeriodicActionsTest::prepareReferenceData()
 114 {
 115     SCOPED_TRACE("Preparing reference data");
 116
 117     // Configure the SimulationRunner to write output to files we can
 118     // use in comparing each test run.
 119     std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
 120
 121     // Run the reference simulation with everything output every step
 122     PeriodicOutputParameters outputEveryStep = {
 123         { "nstenergy", "1" },
 124         { "nstlog", "1" },
 125         { "nstdhdl", "1" },
 126         { "description", "output everything every step" },
 127     };
 128     doMdrun(outputEveryStep);
 129
 130     // Restore the standard filenames we'll use for the runs whose behavior we are testing.
 131     std::swap(runner_.edrFileName_, referenceFileNames_.edrFileName_);
 132 }
 133
 134 /*! \brief Compare the next frame returned by \c testFrameReader with the
 135  * matching frame from \c referenceTestReader using \c comparator.
 136  *
 137  * \returns  True when a successful comparison was made
 138  */
 139 template<typename Reader, typename Comparator>
 140 bool compareFrames(Reader referenceFrameReader, Reader testFrameReader, const Comparator& comparator)
 141 {
 142     if (!testFrameReader->readNextFrame())
 143     {
 144         return false;
 145     }
 146
 147     auto        testFrame = testFrameReader->frame();
 148     std::string frameName = testFrame.frameName();
 149     SCOPED_TRACE("Found frame from test file named " + frameName);
 150     bool foundMatch = false;
 151     while (!foundMatch)
 152     {
 153         if (referenceFrameReader->readNextFrame())
 154         {
 155             auto referenceFrame = referenceFrameReader->frame();
 156             // Can the reference and test frames be compared, because they have the same name?
 157             if (frameName == referenceFrame.frameName())
 158             {
 159                 SCOPED_TRACE("Found frame from reference file named " + frameName);
 160                 comparator(referenceFrame, testFrame);
 161                 foundMatch = true;
 162             }
 163         }
 164         else
 165         {
 166             ADD_FAILURE() << "Ran out of reference frames to compare";
 167             return false;
 168         }
 169     }
 170     return foundMatch;
 171 }
 172
 173 TEST_P(PeriodicActionsTest, PeriodicActionsAgreeWithReference)
 174 {
 175     auto propagation = std::get<0>(GetParam());
 176     SCOPED_TRACE(formatString("Comparing two simulations of '%s' with integrator '%s'",
 177                               propagation["simulationName"].c_str(),
 178                               propagation["integrator"].c_str()));
 179
 180     prepareReferenceData();
 181
 182     auto outputParametersGenerator = std::get<1>(GetParam());
 183     for (const PeriodicOutputParameters& output : outputParametersGenerator())
 184     {
 185         SCOPED_TRACE("Comparing to observe " + output.at("description"));
 186
 187         // Run the test simulation
 188         doMdrun(output);
 189
 190         // Prepare readers for the reference and test output files
 191         auto referenceEnergyFrameReader =
 192                 openEnergyFileToReadTerms(referenceFileNames_.edrFileName_, namesOfEnergiesToMatch_);
 193         auto testEnergyFrameReader =
 194                 openEnergyFileToReadTerms(runner_.edrFileName_, namesOfEnergiesToMatch_);
 195
 196         const bool shouldCompareEnergies   = fromString<int>(output.at("nstenergy")) > 0;
 197         bool       shouldContinueComparing = shouldCompareEnergies;
 198         while (shouldContinueComparing)
 199         {
 200             if (shouldCompareEnergies)
 201             {
 202                 SCOPED_TRACE("Comparing energy frames from reference '" + referenceFileNames_.edrFileName_
 203                              + "' and test '" + runner_.edrFileName_ + "'");
 204                 shouldContinueComparing = shouldContinueComparing
 205                                           && compareFrames(referenceEnergyFrameReader.get(),
 206                                                            testEnergyFrameReader.get(),
 207                                                            energyComparison_);
 208             }
 209         }
 210     }
 211 }
 212
 213 /*! \brief Some common choices of periodic output mdp parameters to
 214  * simplify defining values for the combinations under test */
 215 static PeriodicOutputParameters g_basicPeriodicOutputParameters = {
 216     { "nstenergy", "0" }, { "nstlog", "0" },           { "nstxout", "0" },
 217     { "nstvout", "0" },   { "nstfout", "0" },          { "nstxout-compressed", "0" },
 218     { "nstdhdl", "0" },   { "description", "unknown" }
 219 };
 220
 221 // \todo Test nstlog, nstdhdl, nstxout-compressed
 222 std::vector<PeriodicOutputParameters> outputParameters()
 223 {
 224     std::vector<PeriodicOutputParameters> parameterSets;
 225
 226     parameterSets.push_back(g_basicPeriodicOutputParameters);
 227     parameterSets.back()["nstenergy"]   = "1";
 228     parameterSets.back()["description"] = "energies every step works";
 229
 230     parameterSets.push_back(g_basicPeriodicOutputParameters);
 231     parameterSets.back()["nstenergy"]   = "4";
 232     parameterSets.back()["description"] = "energies every fourth step works";
 233
 234     parameterSets.push_back(g_basicPeriodicOutputParameters);
 235     parameterSets.back()["nstxout"]     = "4";
 236     parameterSets.back()["description"] = "coordinates every fourth step works";
 237
 238     parameterSets.push_back(g_basicPeriodicOutputParameters);
 239     parameterSets.back()["nstvout"]     = "4";
 240     parameterSets.back()["description"] = "velocities every fourth step works";
 241
 242     parameterSets.push_back(g_basicPeriodicOutputParameters);
 243     parameterSets.back()["nstfout"]     = "4";
 244     parameterSets.back()["description"] = "forces every fourth step works";
 245
 246     return parameterSets;
 247 }
 248
 249 std::vector<PropagationParameters> simplePropagationParameters()
 250 {
 251     return {
 252         { { "simulationName", "argon12" },
 253           { "integrator", "md" },
 254           { "comm-mode", "linear" },
 255           { "nstcomm", "1" },
 256           { "tcoupl", "no" },
 257           { "nsttcouple", "0" },
 258           { "pcoupl", "no" },
 259           { "nstpcouple", "0" },
 260           { "maxGromppWarningsTolerated", "0" } },
 261     };
 262 }
 263
 264 std::vector<PropagationParameters> propagationParametersWithCoupling()
 265 {
 266     std::string nsttcouple = "2";
 267     std::string nstpcouple = "3";
 268     std::string comm_mode  = "linear";
 269     std::string nstcomm    = "5";
 270
 271     std::vector<PropagationParameters> parameterSets;
 272     std::vector<std::string>           simulations = { "argon12" };
 273     for (const std::string& simulationName : simulations)
 274     {
 275         std::vector<std::string> integrators{ "md", "sd", "md-vv" };
 276         for (const std::string& integrator : integrators)
 277         {
 278             std::vector<std::string> tcouplValues{ "no", "v-rescale", "Nose-Hoover" };
 279             for (const std::string& tcoupl : tcouplValues)
 280             {
 281                 // SD doesn't support temperature-coupling algorithms,
 282                 if (integrator == "sd" && tcoupl != "no")
 283                 {
 284                     continue;
 285                 }
 286                 for (std::string pcoupl : { "no", "Berendsen", "Parrinello-Rahman", "C-rescale" })
 287                 {
 288                     // VV supports few algorithm combinations
 289                     if (integrator == "md-vv")
 290                     {
 291                         // P-R with VV is called MTTK
 292                         if (pcoupl == "Parrinello-Rahman")
 293                         {
 294                             pcoupl = "MTTK";
 295                         }
 296                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 297                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 298                         {
 299                             continue;
 300                         }
 301                     }
 302                     if (pcoupl == "MTTK" && simulationName != "argon12")
 303                     {
 304                         // MTTK does not support constraints
 305                         continue;
 306                     }
 307
 308                     int maxGromppWarningsTolerated = 0;
 309                     if (pcoupl == "Berendsen")
 310                     {
 311                         ++maxGromppWarningsTolerated;
 312                     }
 313                     parameterSets.emplace_back(PropagationParameters{
 314                             { "simulationName", simulationName },
 315                             { "integrator", integrator },
 316                             { "comm-mode", comm_mode },
 317                             { "nstcomm", nstcomm },
 318                             { "tcoupl", tcoupl },
 319                             { "nsttcouple", nsttcouple },
 320                             { "pcoupl", pcoupl },
 321                             { "nstpcouple", nstpcouple },
 322                             { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
 323                 }
 324             }
 325         }
 326     }
 327     return parameterSets;
 328 }
 329
 330 std::vector<PropagationParameters> propagationParametersWithConstraints()
 331 {
 332     std::string nsttcouple = "2";
 333     std::string nstpcouple = "3";
 334     std::string comm_mode  = "linear";
 335     std::string nstcomm    = "5";
 336
 337     std::vector<PropagationParameters> parameterSets;
 338     std::vector<std::string>           simulations = { "tip3p5" };
 339     for (const std::string& simulationName : simulations)
 340     {
 341         std::vector<std::string> integrators{ "md", "sd", "md-vv" };
 342         for (const std::string& integrator : integrators)
 343         {
 344             std::vector<std::string> tcouplValues{ "no", "v-rescale" };
 345             for (const std::string& tcoupl : tcouplValues)
 346             {
 347                 // SD doesn't support temperature-coupling algorithms,
 348                 if (integrator == "sd" && tcoupl != "no")
 349                 {
 350                     continue;
 351                 }
 352                 for (std::string pcoupl : { "no", "Parrinello-Rahman", "C-rescale" })
 353                 {
 354                     // VV supports few algorithm combinations
 355                     if (integrator == "md-vv")
 356                     {
 357                         // P-R with VV is called MTTK
 358                         if (pcoupl == "Parrinello-Rahman")
 359                         {
 360                             pcoupl = "MTTK";
 361                         }
 362                         if ((tcoupl == "Nose-Hoover" && pcoupl == "Berendsen")
 363                             || (tcoupl != "Nose-Hoover" && pcoupl == "MTTK"))
 364                         {
 365                             continue;
 366                         }
 367                     }
 368                     if (pcoupl == "MTTK" && simulationName != "argon12")
 369                     {
 370                         // MTTK does not support constraints
 371                         continue;
 372                     }
 373
 374                     int maxGromppWarningsTolerated = 0;
 375                     parameterSets.emplace_back(PropagationParameters{
 376                             { "simulationName", simulationName },
 377                             { "integrator", integrator },
 378                             { "comm-mode", comm_mode },
 379                             { "nstcomm", nstcomm },
 380                             { "tcoupl", tcoupl },
 381                             { "nsttcouple", nsttcouple },
 382                             { "pcoupl", pcoupl },
 383                             { "nstpcouple", nstpcouple },
 384                             { "maxGromppWarningsTolerated", toString(maxGromppWarningsTolerated) } });
 385                 }
 386             }
 387         }
 388     }
 389     return parameterSets;
 390 }
 391
 392 } // namespace test
 393 } // namespace gmx