src/programs/mdrun/tests/multiple_time_stepping.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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
  38  * Tests to compare that multiple time stepping is (nearly) identical to normal integration.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \ingroup module_mdrun_integration_tests
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "gromacs/topology/ifunc.h"
  46 #include "gromacs/utility/stringutil.h"
  47
  48 #include "testutils/mpitest.h"
  49 #include "testutils/setenv.h"
  50 #include "testutils/simulationdatabase.h"
  51
  52 #include "moduletest.h"
  53 #include "simulatorcomparison.h"
  54
  55 namespace gmx
  56 {
  57 namespace test
  58 {
  59 namespace
  60 {
  61
  62 /*! \brief Test fixture base for two integration schemes
  63  *
  64  * This test ensures that integration with(out) different multiple time stepping
  65  * scheems (called via different mdp options) yield near identical energies,
  66  * forces and virial at step 0 and similar energies and virial after 4 steps.
  67  */
  68 using MtsComparisonTestParams = std::tuple<std::string, std::string>;
  69 class MtsComparisonTest : public MdrunTestFixture, public ::testing::WithParamInterface<MtsComparisonTestParams>
  70 {
  71 };
  72
  73 //! Returns set of energy terms to compare with associated tolerances
  74 EnergyTermsToCompare energyTermsToCompare(const real energyTol, const real virialTol)
  75 {
  76     return EnergyTermsToCompare{ { { interaction_function[F_EPOT].longname,
  77                                      relativeToleranceAsFloatingPoint(100.0, energyTol) },
  78                                    { "Vir-XX", relativeToleranceAsFloatingPoint(30.0, virialTol) },
  79                                    { "Vir-YY", relativeToleranceAsFloatingPoint(30.0, virialTol) },
  80                                    { "Vir-ZZ", relativeToleranceAsFloatingPoint(30.0, virialTol) } } };
  81 }
  82
  83 TEST_P(MtsComparisonTest, WithinTolerances)
  84 {
  85     auto params         = GetParam();
  86     auto simulationName = std::get<0>(params);
  87     auto mtsScheme      = std::get<1>(params);
  88
  89     // Note that there should be no relevant limitation on MPI ranks and OpenMP threads
  90     SCOPED_TRACE(formatString("Comparing for '%s' no MTS with MTS scheme '%s'",
  91                               simulationName.c_str(), mtsScheme.c_str()));
  92
  93     const int numSteps         = 4;
  94     auto      sharedMdpOptions = gmx::formatString(
  95             "integrator   = md\n"
  96             "dt           = 0.001\n"
  97             "nsteps       = %d\n"
  98             "verlet-buffer-tolerance = -1\n"
  99             "rlist        = 1.0\n"
 100             "coulomb-type = PME\n"
 101             "vdw-type     = cut-off\n"
 102             "rcoulomb     = 0.9\n"
 103             "rvdw         = 0.9\n"
 104             "constraints  = h-bonds\n",
 105             numSteps);
 106
 107     // set nstfout to > numSteps so we only write forces at step 0
 108     const int nstfout       = 2 * numSteps;
 109     auto      refMdpOptions = sharedMdpOptions
 110                          + gmx::formatString(
 111                                    "mts       = no\n"
 112                                    "nstcalcenergy = %d\n"
 113                                    "nstenergy = %d\n"
 114                                    "nstxout   = 0\n"
 115                                    "nstvout   = 0\n"
 116                                    "nstfout   = %d\n",
 117                                    numSteps, numSteps, nstfout);
 118
 119     auto mtsMdpOptions = sharedMdpOptions
 120                          + gmx::formatString(
 121                                    "mts        = yes\n"
 122                                    "mts-levels = 2\n"
 123                                    "mts-level2-forces = %s\n"
 124                                    "mts-level2-factor = 2\n"
 125                                    "nstcalcenergy = %d\n"
 126                                    "nstenergy  = %d\n"
 127                                    "nstxout    = 0\n"
 128                                    "nstvout    = 0\n"
 129                                    "nstfout    = %d\n",
 130                                    mtsScheme.c_str(), numSteps, numSteps, nstfout);
 131
 132     // At step 0 the energy and virial should only differ due to rounding errors
 133     EnergyTermsToCompare energyTermsToCompareStep0 = energyTermsToCompare(0.001, 0.01);
 134     EnergyTermsToCompare energyTermsToCompareAllSteps =
 135             energyTermsToCompare(mtsScheme == "pme" ? 0.015 : 0.04, mtsScheme == "pme" ? 0.1 : 0.2);
 136
 137     // Specify how trajectory frame matching must work.
 138     TrajectoryFrameMatchSettings trajectoryMatchSettings{ true,
 139                                                           true,
 140                                                           true,
 141                                                           ComparisonConditions::NoComparison,
 142                                                           ComparisonConditions::NoComparison,
 143                                                           ComparisonConditions::MustCompare };
 144     TrajectoryTolerances trajectoryTolerances = TrajectoryComparison::s_defaultTrajectoryTolerances;
 145
 146     // Build the functor that will compare reference and test
 147     // trajectory frames in the chosen way.
 148     TrajectoryComparison trajectoryComparison{ trajectoryMatchSettings, trajectoryTolerances };
 149
 150     // Set file names
 151     auto simulator1TrajectoryFileName = fileManager_.getTemporaryFilePath("sim1.trr");
 152     auto simulator1EdrFileName        = fileManager_.getTemporaryFilePath("sim1.edr");
 153     auto simulator2TrajectoryFileName = fileManager_.getTemporaryFilePath("sim2.trr");
 154     auto simulator2EdrFileName        = fileManager_.getTemporaryFilePath("sim2.edr");
 155
 156     // Run grompp
 157     runner_.tprFileName_ = fileManager_.getTemporaryFilePath("sim.tpr");
 158     runner_.useTopGroAndNdxFromDatabase(simulationName);
 159     runner_.useStringAsMdpFile(refMdpOptions);
 160     runGrompp(&runner_);
 161
 162     // Do first mdrun
 163     runner_.fullPrecisionTrajectoryFileName_ = simulator1TrajectoryFileName;
 164     runner_.edrFileName_                     = simulator1EdrFileName;
 165     runMdrun(&runner_);
 166
 167     runner_.useStringAsMdpFile(mtsMdpOptions);
 168     runGrompp(&runner_);
 169
 170     // Do second mdrun
 171     runner_.fullPrecisionTrajectoryFileName_ = simulator2TrajectoryFileName;
 172     runner_.edrFileName_                     = simulator2EdrFileName;
 173     runMdrun(&runner_);
 174
 175     // Compare simulation results at step 0, which should be indentical
 176     compareEnergies(simulator1EdrFileName, simulator2EdrFileName, energyTermsToCompareStep0,
 177                     MaxNumFrames(1));
 178     compareTrajectories(simulator1TrajectoryFileName, simulator2TrajectoryFileName, trajectoryComparison);
 179
 180     // Compare energies at the last step (and step 0 again) with lower tolerance
 181     compareEnergies(simulator1EdrFileName, simulator2EdrFileName, energyTermsToCompareAllSteps,
 182                     MaxNumFrames::compareAllFrames());
 183 }
 184
 185 INSTANTIATE_TEST_CASE_P(
 186         MultipleTimeSteppingIsNearSingleTimeStepping,
 187         MtsComparisonTest,
 188         ::testing::Combine(::testing::Values("ala"),
 189                            ::testing::Values("longrange-nonbonded",
 190                                              "longrange-nonbonded nonbonded pair dihedral")));
 191
 192 } // namespace
 193 } // namespace test
 194 } // namespace gmx