src/programs/mdrun/tests/rerun.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, 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
  38  * Tests for the mdrun -rerun functionality
  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 <map>
  48 #include <memory>
  49 #include <string>
  50 #include <tuple>
  51 #include <vector>
  52
  53 #include <gtest/gtest.h>
  54
  55 #include "gromacs/options/filenameoption.h"
  56 #include "gromacs/topology/idef.h"
  57 #include "gromacs/topology/ifunc.h"
  58 #include "gromacs/trajectory/energyframe.h"
  59 #include "gromacs/trajectory/trajectoryframe.h"
  60 #include "gromacs/utility/stringutil.h"
  61
  62 #include "testutils/mpitest.h"
  63 #include "testutils/simulationdatabase.h"
  64 #include "testutils/testasserts.h"
  65
  66 #include "energycomparison.h"
  67 #include "energyreader.h"
  68 #include "mdruncomparison.h"
  69 #include "moduletest.h"
  70 #include "trajectorycomparison.h"
  71 #include "trajectoryreader.h"
  72
  73 namespace gmx
  74 {
  75 namespace test
  76 {
  77 namespace
  78 {
  79
  80 //! Functor for comparing reference and test frames on particular energies to match.
  81 class EnergyComparator
  82 {
  83     public:
  84         //! Constructor
  85         EnergyComparator(const EnergyTolerances &energiesToMatch)
  86             : energiesToMatch_(energiesToMatch) {}
  87         //! The functor method.
  88         void operator()(const EnergyFrame &reference, const EnergyFrame &test)
  89         {
  90             compareEnergyFrames(reference, test, energiesToMatch_);
  91         }
  92         //! Container of the energies to match and the tolerance required.
  93         const EnergyTolerances &energiesToMatch_;
  94 };
  95
  96 //! Run grompp for a normal mdrun, the run, and its rerun.
  97 void executeRerunTest(TestFileManager        *fileManager,
  98                       SimulationRunner       *runner,
  99                       const std::string      &simulationName,
 100                       int                     maxWarningsTolerated,
 101                       const MdpFieldValues   &mdpFieldValues,
 102                       const EnergyTolerances &energiesToMatch)
 103 {
 104     // TODO At some point we should also test PME-only ranks.
 105     int numRanksAvailable = getNumberOfTestMpiRanks();
 106     if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
 107     {
 108         fprintf(stdout, "Test system '%s' cannot run with %d ranks.\n"
 109                 "The supported numbers are: %s\n",
 110                 simulationName.c_str(), numRanksAvailable,
 111                 reportNumbersOfPpRanksSupported(simulationName).c_str());
 112         return;
 113     }
 114
 115     auto normalRunTrajectoryFileName = fileManager->getTemporaryFilePath("normal.trr");
 116     auto normalRunEdrFileName        = fileManager->getTemporaryFilePath("normal.edr");
 117     auto rerunTrajectoryFileName     = fileManager->getTemporaryFilePath("rerun.trr");
 118     auto rerunEdrFileName            = fileManager->getTemporaryFilePath("rerun.edr");
 119
 120     // prepare the .tpr file
 121     {
 122         // TODO evolve grompp to report the number of warnings issued, so
 123         // tests always expect the right number.
 124         CommandLine caller;
 125         caller.append("grompp");
 126         caller.addOption("-maxwarn", maxWarningsTolerated);
 127         runner->useTopGroAndNdxFromDatabase(simulationName);
 128         runner->useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 129         EXPECT_EQ(0, runner->callGrompp(caller));
 130     }
 131
 132     // do the normal mdrun
 133     {
 134         runner->fullPrecisionTrajectoryFileName_ = normalRunTrajectoryFileName;
 135         runner->edrFileName_                     = normalRunEdrFileName;
 136         CommandLine normalRunCaller;
 137         normalRunCaller.append("mdrun");
 138         ASSERT_EQ(0, runner->callMdrun(normalRunCaller));
 139     }
 140
 141     // do a rerun on the .trr just produced
 142     {
 143         runner->fullPrecisionTrajectoryFileName_ = rerunTrajectoryFileName;
 144         runner->edrFileName_                     = rerunEdrFileName;
 145         CommandLine rerunCaller;
 146         rerunCaller.append("mdrun");
 147         rerunCaller.addOption("-rerun", normalRunTrajectoryFileName);
 148         ASSERT_EQ(0, runner->callMdrun(rerunCaller));
 149     }
 150
 151     // Build the functor that will compare reference and test
 152     // energy frames on the chosen energy fields.
 153     //
 154     // TODO It would be less code if we used a lambda for this, but either
 155     // clang 3.4 or libstdc++ 5.2.1 have an issue with capturing a
 156     // std::unordered_map
 157     EnergyComparator energyComparator(energiesToMatch);
 158     // Build the manager that will present matching pairs of frames to compare.
 159     //
 160     // TODO Here is an unnecessary copy of keys (ie. the energy field
 161     // names), for convenience. In the future, use a range.
 162     auto namesOfEnergiesToMatch = getKeys(energiesToMatch);
 163     FramePairManager<EnergyFrameReader, EnergyFrame>
 164          energyManager(openEnergyFileToReadFields(normalRunEdrFileName, namesOfEnergiesToMatch),
 165                   openEnergyFileToReadFields(rerunEdrFileName, namesOfEnergiesToMatch));
 166     // Compare the energy frames.
 167     energyManager.compareAllFramePairs(energyComparator);
 168
 169     // Specify how trajectory frame matching must work.
 170     TrajectoryFrameMatchSettings trajectoryMatchSettings {
 171         true, true, true, true, false, true
 172     };
 173     /* Specify the default expected tolerances for trajectory
 174      * components for all simulation systems. */
 175     TrajectoryTolerances trajectoryTolerances {
 176         defaultRealTolerance(),                                               // box
 177         relativeToleranceAsFloatingPoint(1.0, 1.0e-3),                        // positions
 178         defaultRealTolerance(),                                               // velocities - unused for rerun
 179         relativeToleranceAsFloatingPoint(100.0, GMX_DOUBLE ? 1.0e-7 : 1.0e-5) // forces
 180     };
 181
 182     // Build the functor that will compare reference and test
 183     // trajectory frames in the chosen way.
 184     auto trajectoryComparator = [&trajectoryMatchSettings, &trajectoryTolerances](const TrajectoryFrame &reference, const TrajectoryFrame &test)
 185         {
 186             compareTrajectoryFrames(reference, test, trajectoryMatchSettings, trajectoryTolerances);
 187         };
 188     // Build the manager that will present matching pairs of frames to compare
 189     FramePairManager<TrajectoryFrameReader, TrajectoryFrame>
 190     trajectoryManager(std::make_unique<TrajectoryFrameReader>(normalRunTrajectoryFileName),
 191                       std::make_unique<TrajectoryFrameReader>(rerunTrajectoryFileName));
 192     // Compare the trajectory frames.
 193     trajectoryManager.compareAllFramePairs(trajectoryComparator);
 194 }
 195
 196 /*! \brief Test fixture base for mdrun -rerun
 197  *
 198  * This test ensures mdrun can run a simulation, writing a trajectory
 199  * and matching energies, and reproduce the same energies from a rerun
 200  * to within a tight tolerance. It says nothing about whether a rerun
 201  * can reproduce energies from a trajectory generated with older code,
 202  * since that is not a useful property. Whether mdrun produced correct
 203  * energies then and now needs different kinds of testing, but if
 204  * true, this test ensures the rerun has the expected property.
 205  *
 206  * The limitations of mdrun and its output means that reproducing the
 207  * same energies is currently only meaningful for integration without
 208  * thermostats or barostats, however the present form of the test
 209  * infrastructure has in-principle support for such, if that is ever
 210  * needed/useful.
 211  *
 212  * We should also not compare pressure, because with constraints the
 213  * non-search steps need a much larger tolerance, and per Redmine 1868
 214  * we should stop computing pressure in reruns anyway.
 215  *
 216  * Similarly, per 1868, in the present implementation the kinetic
 217  * energy quantities are not generally reproducible, either.
 218  *
 219  * The choices for tolerance are arbitrary but sufficient.  Rerun does
 220  * pair search every frame, so it cannot in general exactly reproduce
 221  * quantities from a normal run, because the accumulation order
 222  * differs. (Nor does it reproduce pair-search frames exactly,
 223  * either). */
 224 class MdrunRerunTest : public MdrunTestFixture,
 225                        public ::testing::WithParamInterface <
 226                        std::tuple < std::string, std::string>>
 227 {
 228 };
 229
 230 TEST_P(MdrunRerunTest, WithinTolerances)
 231 {
 232     auto params         = GetParam();
 233     auto simulationName = std::get<0>(params);
 234     auto integrator     = std::get<1>(params);
 235     SCOPED_TRACE(formatString("Comparing normal and rerun of simulation '%s' "
 236                               "with integrator '%s'",
 237                               simulationName.c_str(), integrator.c_str()));
 238
 239     auto mdpFieldValues = prepareMdpFieldValues(simulationName.c_str(),
 240                                                 integrator.c_str(),
 241                                                 "no", "no");
 242
 243     // bd is much less reproducible in a rerun than the other integrators
 244     const int        toleranceScaleFactor = (integrator == "bd") ? 2 : 1;
 245     EnergyTolerances energiesToMatch
 246     {{
 247          {
 248              interaction_function[F_EPOT].longname,
 249              relativeToleranceAsPrecisionDependentUlp(10.0, 24 * toleranceScaleFactor, 40 * toleranceScaleFactor)
 250          },
 251      }};
 252
 253     int numWarningsToTolerate = 0;
 254     executeRerunTest(&fileManager_, &runner_,
 255                      simulationName, numWarningsToTolerate, mdpFieldValues,
 256                      energiesToMatch);
 257 }
 258
 259 // TODO The time for OpenCL kernel compilation means these tests time
 260 // out. Once that compilation is cached for the whole process, these
 261 // tests can run in such configurations.
 262 #if GMX_GPU != GMX_GPU_OPENCL
 263 INSTANTIATE_TEST_CASE_P(NormalMdrunIsReproduced, MdrunRerunTest,
 264                             ::testing::Combine(::testing::Values("argon12", "tip3p5", "alanine_vsite_vacuo"),
 265                                                    ::testing::Values("md", "md-vv", "bd", "sd")));
 266 #else
 267 INSTANTIATE_TEST_CASE_P(DISABLED_NormalMdrunIsReproduced, MdrunRerunTest,
 268                             ::testing::Combine(::testing::Values("argon12", "tip3p5", "alanine_vsite_vacuo"),
 269                                                    ::testing::Values("md", "md-vv", "bd", "sd")));
 270 #endif
 271
 272 class MdrunRerunFreeEnergyTest : public MdrunTestFixture,
 273                                  public ::testing::WithParamInterface <
 274                                  std::tuple < std::string, std::string, int>>
 275 {
 276 };
 277
 278 TEST_P(MdrunRerunFreeEnergyTest, WithinTolerances)
 279 {
 280     auto params          = GetParam();
 281     auto simulationName  = std::get<0>(params);
 282     auto integrator      = std::get<1>(params);
 283     auto initLambdaState = std::get<2>(params);
 284     SCOPED_TRACE(formatString("Comparing normal and rerun of simulation '%s' "
 285                               "with integrator '%s' for initial lambda state %d",
 286                               simulationName.c_str(), integrator.c_str(), initLambdaState));
 287
 288     auto mdpFieldValues = prepareMdpFieldValues(simulationName.c_str(),
 289                                                 integrator.c_str(),
 290                                                 "no", "no");
 291     mdpFieldValues["other"] += formatString("\ninit-lambda-state %d", initLambdaState);
 292
 293     EnergyTolerances energiesToMatch
 294     {{
 295          {
 296              interaction_function[F_EPOT].longname, relativeToleranceAsPrecisionDependentUlp(10.0, 24, 32)
 297          },
 298          {
 299              interaction_function[F_DVDL_COUL].longname, relativeToleranceAsPrecisionDependentUlp(1.0, 8, 8)
 300          },
 301          {
 302              interaction_function[F_DVDL_VDW].longname, relativeToleranceAsPrecisionDependentUlp(1.0, 8, 8)
 303          },
 304          {
 305              interaction_function[F_DVDL_BONDED].longname, relativeToleranceAsPrecisionDependentUlp(1.0, 8, 8)
 306          },
 307          {
 308              interaction_function[F_DVDL_RESTRAINT].longname, relativeToleranceAsPrecisionDependentUlp(1.0, 8, 8)
 309          }
 310      }};
 311
 312     // The md integrator triggers a warning for nearly decoupled
 313     // states, which we need to suppress. TODO sometimes?
 314     int numWarningsToTolerate = (integrator == "md") ? 1 : 0;
 315     executeRerunTest(&fileManager_, &runner_,
 316                      simulationName, numWarningsToTolerate, mdpFieldValues,
 317                      energiesToMatch);
 318 }
 319
 320 // TODO The time for OpenCL kernel compilation means these tests time
 321 // out. Once that compilation is cached for the whole process, these
 322 // tests can run in such configurations.
 323 #if GMX_GPU != GMX_GPU_OPENCL
 324 INSTANTIATE_TEST_CASE_P(MdrunIsReproduced, MdrunRerunFreeEnergyTest,
 325                             ::testing::Combine(::testing::Values("nonanol_vacuo"),
 326                                                    ::testing::Values("md", "md-vv", "sd"),
 327                                                    ::testing::Range(0, 11)));
 328 #else
 329 INSTANTIATE_TEST_CASE_P(DISABLED_MdrunIsReproduced, MdrunRerunFreeEnergyTest,
 330                             ::testing::Combine(::testing::Values("nonanol_vacuo"),
 331                                                    ::testing::Values("md", "md-vv", "sd"),
 332                                                    ::testing::Range(0, 11)));
 333 #endif
 334
 335 } // namespace
 336 } // namespace test
 337 } // namespace gmx