src/programs/mdrun/tests/pull.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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
  38  * Tests for mdrun pull functionality.
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \ingroup module_mdrun_integration_tests
  42  */
  43 #include "gmxpre.h"
  44
  45 #include <map>
  46 #include <memory>
  47 #include <string>
  48 #include <unordered_map>
  49 #include <vector>
  50 #include <regex>
  51
  52 #include <gtest/gtest.h>
  53
  54 #include "gromacs/options/filenameoption.h"
  55 #include "gromacs/topology/idef.h"
  56 #include "gromacs/topology/ifunc.h"
  57 #include "gromacs/utility/stringutil.h"
  58
  59 #include "testutils/mpitest.h"
  60 #include "testutils/refdata.h"
  61 #include "testutils/simulationdatabase.h"
  62 #include "testutils/testasserts.h"
  63
  64 #include "energycomparison.h"
  65 #include "moduletest.h"
  66
  67 namespace gmx
  68 {
  69 namespace test
  70 {
  71 namespace
  72 {
  73
  74 /*! \brief Database of energy tolerances on the various systems. */
  75 std::unordered_map<std::string, double> energyToleranceForSystem_g = { { { "spc216", 1e-4 } } };
  76
  77 /*! \brief Database of pressure tolerances on the various systems. */
  78 std::unordered_map<std::string, double> pressureToleranceForSystem_g = { { { "spc216", 2e-4 } } };
  79
  80 const std::unordered_map<std::string, std::vector<std::pair<std::string, std::string>>> c_mdpPullParams = {
  81     { { "umbrella-3D",
  82         { { "pull-ngroups", "2" },
  83           { "pull-ncoords", "1" },
  84           { "pull-nstxout", "0" },
  85           { "pull-nstfout", "0" },
  86           { "pull-group1-name", "r_1" },
  87           { "pull-group2-name", "r_2" },
  88           { "pull-coord1-groups", "1 2" },
  89           { "pull-coord1-type", "umbrella" },
  90           { "pull-coord1-geometry", "distance" },
  91           { "pull-coord1-dim", "Y Y Y" },
  92           { "pull-coord1-init", "0.6" },
  93           { "pull-coord1-k", "100" } } },
  94       { "umbrella-2D",
  95         { { "pull-ngroups", "2" },
  96           { "pull-ncoords", "1" },
  97           { "pull-nstxout", "0" },
  98           { "pull-nstfout", "0" },
  99           { "pull-group1-name", "r_1" },
 100           { "pull-group2-name", "r_2" },
 101           { "pull-coord1-groups", "1 2" },
 102           { "pull-coord1-type", "umbrella" },
 103           { "pull-coord1-geometry", "distance" },
 104           { "pull-coord1-dim", "Y Y N" },
 105           { "pull-coord1-init", "0.4" },
 106           { "pull-coord1-k", "100" } } },
 107       { "constraint-flatbottom",
 108         { { "pull-ngroups", "3" },
 109           { "pull-ncoords", "2" },
 110           { "pull-nstxout", "0" },
 111           { "pull-nstfout", "0" },
 112           { "pull-group1-name", "r_1" },
 113           { "pull-group2-name", "r_2" },
 114           { "pull-group3-name", "r_3" },
 115           { "pull-coord1-groups", "1 2" },
 116           { "pull-coord1-type", "constraint" },
 117           { "pull-coord1-geometry", "distance" },
 118           { "pull-coord1-dim", "Y Y Y" },
 119           { "pull-coord1-init", "0.5" },
 120           { "pull-coord2-groups", "1 3" },
 121           { "pull-coord2-type", "flat-bottom" },
 122           { "pull-coord2-geometry", "distance" },
 123           { "pull-coord2-dim", "Y Y Y" },
 124           { "pull-coord2-init", "0.4" },
 125           { "pull-coord2-k", "100" } } },
 126       { "transformation-coord-umbrella-3D", // should yield the same force on the first pull coordinate as in the umbrella-3D test above
 127         {
 128                 { "pull-ngroups", "2" },
 129                 { "pull-ncoords", "2" },
 130                 { "pull-nstxout", "0" },
 131                 { "pull-nstfout", "0" },
 132                 { "pull-group1-name", "r_1" },
 133                 { "pull-group2-name", "r_2" },
 134                 { "pull-coord1-groups", "1 2" },
 135                 { "pull-coord1-type", "umbrella" },
 136                 { "pull-coord1-geometry", "distance" },
 137                 { "pull-coord1-dim", "Y Y Y" },
 138                 { "pull-coord1-k", "0" }, // set to zero since the force comes from the transformation coordinate
 139                 { "pull-coord2-geometry", "transformation" },
 140                 { "pull-coord2-type", "umbrella" },
 141                 { "pull-coord2-expression", "2*(0.6 - x1)" },
 142                 { "pull-coord2-init", "0" }, // -> pull-coord1-init = 0.6
 143                 { "pull-coord2-k", "25" }    // -> pull-coord1-k = 2^2*25 = 100
 144         } } }
 145 };
 146
 147 //! Helper type
 148 using MdpField = MdpFieldValues::value_type;
 149
 150 /*! \brief Test fixture base for simple mdrun systems
 151  *
 152  * This test ensures mdrun can run a simulation, reaching
 153  * reproducible energies.
 154  *
 155  * The choices for tolerance are arbitrary but sufficient. */
 156 class PullIntegrationTest :
 157     public MdrunTestFixture,
 158     public ::testing::WithParamInterface<std::tuple<std::string, std::string>>
 159 {
 160 };
 161
 162 //! Adds integrator and nonbonded parameter setup
 163 void addBasicMdpValues(MdpFieldValues* mdpFieldValues)
 164 {
 165     (*mdpFieldValues)["nsteps"]        = "20";
 166     (*mdpFieldValues)["nstcomm"]       = "10";
 167     (*mdpFieldValues)["nstlist"]       = "10";
 168     (*mdpFieldValues)["nstcalcenergy"] = "5";
 169     (*mdpFieldValues)["nstenergy"]     = "5";
 170     (*mdpFieldValues)["coulombtype"]   = "Reaction-field";
 171     (*mdpFieldValues)["vdwtype"]       = "Cut-off";
 172 }
 173
 174 bool isTransformationPullSetup(const std::string& pullSetupName)
 175 {
 176     const auto* pattern{ "transformation" };
 177     // std::regex_constants::icase - TO IGNORE CASE.
 178     auto rx = std::regex{ pattern, std::regex_constants::icase };
 179     return std::regex_search(pullSetupName, rx);
 180 }
 181
 182 TEST_P(PullIntegrationTest, WithinTolerances)
 183 {
 184     auto params         = GetParam();
 185     auto simulationName = std::get<0>(params);
 186     auto pullSetup      = std::get<1>(params);
 187     SCOPED_TRACE(formatString("Comparing simple mdrun for '%s'", simulationName.c_str()));
 188
 189     // TODO At some point we should also test PME-only ranks.
 190     int numRanksAvailable = getNumberOfTestMpiRanks();
 191     if (!isNumberOfPpRanksSupported(simulationName, numRanksAvailable))
 192     {
 193         fprintf(stdout,
 194                 "Test system '%s' cannot run with %d ranks.\n"
 195                 "The supported numbers are: %s\n",
 196                 simulationName.c_str(),
 197                 numRanksAvailable,
 198                 reportNumbersOfPpRanksSupported(simulationName).c_str());
 199         return;
 200     }
 201     auto mdpFieldValues = prepareMdpFieldValues(simulationName.c_str(), "md", "no", "no");
 202     addBasicMdpValues(&mdpFieldValues);
 203
 204     // Add the pull parameters
 205     mdpFieldValues["pull"]    = "yes";
 206     const auto& mdpPullParams = c_mdpPullParams.at(pullSetup);
 207     for (const auto& param : mdpPullParams)
 208     {
 209         mdpFieldValues[param.first] = param.second;
 210     }
 211
 212     // Prepare the .tpr file
 213     {
 214         CommandLine caller;
 215         runner_.useTopGroAndNdxFromDatabase(simulationName);
 216         runner_.useNdxFromDatabase(simulationName + "_pull");
 217         runner_.useStringAsMdpFile(prepareMdpFileContents(mdpFieldValues));
 218         EXPECT_EQ(0, runner_.callGrompp(caller));
 219     }
 220     // Do mdrun
 221     {
 222         // conver the tolerance to relative floating point tolerance
 223         auto energyTolerance   = energyToleranceForSystem_g.at(simulationName);
 224         auto pressureTolerance = pressureToleranceForSystem_g.at(simulationName);
 225         if (isTransformationPullSetup(pullSetup)) // need to increase the tolerance a bit due to numerical evaluations
 226         {
 227             energyTolerance *= 10;
 228             pressureTolerance *= 10;
 229         }
 230         auto relativeEnergyTolerance   = relativeToleranceAsFloatingPoint(1, energyTolerance);
 231         auto relativePressureTolerance = relativeToleranceAsFloatingPoint(1, pressureTolerance);
 232         CommandLine mdrunCaller;
 233         ASSERT_EQ(0, runner_.callMdrun(mdrunCaller));
 234         EnergyTermsToCompare energyTermsToCompare{ {
 235                 { interaction_function[F_COM_PULL].longname, relativeEnergyTolerance },
 236                 { interaction_function[F_EPOT].longname, relativeEnergyTolerance },
 237                 { interaction_function[F_EKIN].longname, relativeEnergyTolerance },
 238                 { interaction_function[F_PRES].longname, relativePressureTolerance },
 239         } };
 240         TestReferenceData    refData;
 241         auto                 checker = refData.rootChecker()
 242                                .checkCompound("Simulation", simulationName)
 243                                .checkCompound("PullSetup", pullSetup);
 244         checkEnergiesAgainstReferenceData(runner_.edrFileName_, energyTermsToCompare, &checker);
 245     }
 246 }
 247
 248 INSTANTIATE_TEST_SUITE_P(PullTest,
 249                          PullIntegrationTest,
 250                          ::testing::Combine(::testing::Values("spc216"),
 251                                             ::testing::Values("umbrella-3D",
 252                                                               "umbrella-2D",
 253                                                               "constraint-flatbottom",
 254                                                               "transformation-coord-umbrella-3D")));
 255
 256 } // namespace
 257 } // namespace test
 258 } // namespace gmx