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(
  91             "Comparing for '%s' no MTS with MTS scheme '%s'", simulationName.c_str(), mtsScheme.c_str()));
  92
  93     const bool isPullTest = (mtsScheme.find("pull") != std::string::npos);
  94
  95     const int numSteps         = 4;
  96     auto      sharedMdpOptions = gmx::formatString(
  97             "integrator   = md\n"
  98             "dt           = 0.001\n"
  99             "nsteps       = %d\n"
 100             "verlet-buffer-tolerance = -1\n"
 101             "rlist        = 1.0\n"
 102             "coulomb-type = %s\n"
 103             "vdw-type     = cut-off\n"
 104             "rcoulomb     = 0.9\n"
 105             "rvdw         = 0.9\n"
 106             "constraints  = h-bonds\n",
 107             numSteps,
 108             isPullTest ? "reaction-field" : "PME");
 109
 110     if (isPullTest)
 111     {
 112         sharedMdpOptions +=
 113                 "pull                 = yes\n"
 114                 "pull-ngroups         = 2\n"
 115                 "pull-group1-name     = FirstWaterMolecule\n"
 116                 "pull-group2-name     = SecondWaterMolecule\n"
 117                 "pull-ncoords         = 1\n"
 118                 "pull-coord1-type     = umbrella\n"
 119                 "pull-coord1-geometry = distance\n"
 120                 "pull-coord1-groups   = 1 2\n"
 121                 "pull-coord1-init     = 1\n"
 122                 "pull-coord1-k        = 10000\n";
 123     }
 124
 125     // set nstfout to > numSteps so we only write forces at step 0
 126     const int nstfout       = 2 * numSteps;
 127     auto      refMdpOptions = sharedMdpOptions
 128                          + gmx::formatString(
 129                                    "mts       = no\n"
 130                                    "nstcalcenergy = %d\n"
 131                                    "nstenergy = %d\n"
 132                                    "nstxout   = 0\n"
 133                                    "nstvout   = 0\n"
 134                                    "nstfout   = %d\n",
 135                                    numSteps,
 136                                    numSteps,
 137                                    nstfout);
 138
 139     auto mtsMdpOptions = sharedMdpOptions
 140                          + gmx::formatString(
 141                                    "mts        = yes\n"
 142                                    "mts-levels = 2\n"
 143                                    "mts-level2-forces = %s\n"
 144                                    "mts-level2-factor = 2\n"
 145                                    "nstcalcenergy = %d\n"
 146                                    "nstenergy  = %d\n"
 147                                    "nstxout    = 0\n"
 148                                    "nstvout    = 0\n"
 149                                    "nstfout    = %d\n",
 150                                    mtsScheme.c_str(),
 151                                    numSteps,
 152                                    numSteps,
 153                                    nstfout);
 154
 155     // At step 0 the energy and virial should only differ due to rounding errors
 156     EnergyTermsToCompare energyTermsToCompareStep0 = energyTermsToCompare(0.001, 0.01);
 157     EnergyTermsToCompare energyTermsToCompareAllSteps =
 158             energyTermsToCompare(mtsScheme == "pme" ? 0.015 : 0.04, mtsScheme == "pme" ? 0.1 : 0.2);
 159
 160     // Specify how trajectory frame matching must work.
 161     TrajectoryFrameMatchSettings trajectoryMatchSettings{ true,
 162                                                           true,
 163                                                           true,
 164                                                           ComparisonConditions::NoComparison,
 165                                                           ComparisonConditions::NoComparison,
 166                                                           ComparisonConditions::MustCompare };
 167     TrajectoryTolerances trajectoryTolerances = TrajectoryComparison::s_defaultTrajectoryTolerances;
 168
 169     // Build the functor that will compare reference and test
 170     // trajectory frames in the chosen way.
 171     TrajectoryComparison trajectoryComparison{ trajectoryMatchSettings, trajectoryTolerances };
 172
 173     // Set file names
 174     auto simulator1TrajectoryFileName = fileManager_.getTemporaryFilePath("sim1.trr");
 175     auto simulator1EdrFileName        = fileManager_.getTemporaryFilePath("sim1.edr");
 176     auto simulator2TrajectoryFileName = fileManager_.getTemporaryFilePath("sim2.trr");
 177     auto simulator2EdrFileName        = fileManager_.getTemporaryFilePath("sim2.edr");
 178
 179     // Run grompp
 180     runner_.tprFileName_ = fileManager_.getTemporaryFilePath("sim.tpr");
 181     runner_.useTopGroAndNdxFromDatabase(simulationName);
 182     runner_.useStringAsMdpFile(refMdpOptions);
 183     runGrompp(&runner_);
 184
 185     // Do first mdrun
 186     runner_.fullPrecisionTrajectoryFileName_ = simulator1TrajectoryFileName;
 187     runner_.edrFileName_                     = simulator1EdrFileName;
 188     runMdrun(&runner_);
 189
 190     runner_.useStringAsMdpFile(mtsMdpOptions);
 191     runGrompp(&runner_);
 192
 193     // Do second mdrun
 194     runner_.fullPrecisionTrajectoryFileName_ = simulator2TrajectoryFileName;
 195     runner_.edrFileName_                     = simulator2EdrFileName;
 196     runMdrun(&runner_);
 197
 198     // Compare simulation results at step 0, which should be indentical
 199     compareEnergies(
 200             simulator1EdrFileName, simulator2EdrFileName, energyTermsToCompareStep0, MaxNumFrames(1));
 201     compareTrajectories(simulator1TrajectoryFileName, simulator2TrajectoryFileName, trajectoryComparison);
 202
 203     // Compare energies at the last step (and step 0 again) with lower tolerance
 204     compareEnergies(simulator1EdrFileName,
 205                     simulator2EdrFileName,
 206                     energyTermsToCompareAllSteps,
 207                     MaxNumFrames::compareAllFrames());
 208 }
 209
 210 INSTANTIATE_TEST_CASE_P(
 211         MultipleTimeSteppingIsNearSingleTimeStepping,
 212         MtsComparisonTest,
 213         ::testing::Combine(::testing::Values("ala"),
 214                            ::testing::Values("longrange-nonbonded",
 215                                              "longrange-nonbonded nonbonded pair dihedral")));
 216
 217 INSTANTIATE_TEST_CASE_P(MultipleTimeSteppingIsNearSingleTimeSteppingPull,
 218                         MtsComparisonTest,
 219                         ::testing::Combine(::testing::Values("spc2"), ::testing::Values("pull")));
 220
 221 } // namespace
 222 } // namespace test
 223 } // namespace gmx