2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 * Tests to compare free energy simulations to reference
40 * \author Pascal Merz <pascal.merz@me.com>
41 * \ingroup module_mdrun_integration_tests
47 #include "gromacs/topology/ifunc.h"
48 #include "gromacs/utility/filestream.h"
49 #include "gromacs/utility/path.h"
50 #include "gromacs/utility/stringutil.h"
52 #include "testutils/mpitest.h"
53 #include "testutils/refdata.h"
54 #include "testutils/setenv.h"
55 #include "testutils/simulationdatabase.h"
56 #include "testutils/xvgtest.h"
58 #include "moduletest.h"
59 #include "simulatorcomparison.h"
66 /*! \brief Test fixture base for free energy calculations
68 * This test ensures that selected free energy perturbation calculations produce
69 * results identical to an earlier version. The results of this earlier version
70 * have been verified manually to ensure physical correctness.
72 using MaxNumWarnings = int;
73 using ListOfInteractionsToTest = std::vector<int>;
74 using FreeEnergyReferenceTestParams = std::tuple<std::string, MaxNumWarnings, ListOfInteractionsToTest>;
75 class FreeEnergyReferenceTest :
76 public MdrunTestFixture,
77 public ::testing::WithParamInterface<FreeEnergyReferenceTestParams>
80 struct PrintParametersToString
82 template<class ParamType>
83 std::string operator()(const testing::TestParamInfo<ParamType>& parameter) const
85 auto simulationName = std::get<0>(parameter.param);
86 std::replace(simulationName.begin(), simulationName.end(), '-', '_');
87 return simulationName + (GMX_DOUBLE ? "_d" : "_s");
92 TEST_P(FreeEnergyReferenceTest, WithinTolerances)
94 const auto& simulationName = std::get<0>(GetParam());
95 const auto maxNumWarnings = std::get<1>(GetParam());
96 const auto& interactionsList = std::get<2>(GetParam());
98 // As these tests check reproducibility, we restrict the maximum number
99 // of ranks to allow us to keep the tolerances tight. See also #3741.
100 const int numRanksAvailable = getNumberOfTestMpiRanks();
101 constexpr int maxNumRanks = 8;
102 if (numRanksAvailable > maxNumRanks)
105 "The FEP tests cannot run with %d ranks.\n"
106 "The maximum number of ranks supported is %d.",
107 numRanksAvailable, maxNumRanks);
111 SCOPED_TRACE(formatString("Comparing FEP simulation '%s' to reference", simulationName.c_str()));
113 // Tolerance set to pass with identical code version and a range of different test setups for most tests
114 const auto defaultEnergyTolerance = relativeToleranceAsFloatingPoint(50.0, GMX_DOUBLE ? 1e-5 : 1e-4);
115 // Some simulations are significantly longer, so they need a larger tolerance
116 const auto longEnergyTolerance = relativeToleranceAsFloatingPoint(50.0, GMX_DOUBLE ? 1e-4 : 1e-3);
117 const bool isLongSimulation = (simulationName == "expanded");
118 const auto energyTolerance = isLongSimulation ? longEnergyTolerance : defaultEnergyTolerance;
120 EnergyTermsToCompare energyTermsToCompare{ { interaction_function[F_EPOT].longname, energyTolerance } };
121 for (const auto& interaction : interactionsList)
123 energyTermsToCompare.emplace(interaction_function[interaction].longname, defaultEnergyTolerance);
126 // Specify how trajectory frame matching must work (only testing forces).
127 TrajectoryFrameMatchSettings trajectoryMatchSettings{ false,
130 ComparisonConditions::NoComparison,
131 ComparisonConditions::NoComparison,
132 ComparisonConditions::MustCompare };
133 TrajectoryTolerances trajectoryTolerances = TrajectoryComparison::s_defaultTrajectoryTolerances;
134 trajectoryTolerances.forces = relativeToleranceAsFloatingPoint(100.0, GMX_DOUBLE ? 5.0e-5 : 5.0e-4);
136 // Build the functor that will compare reference and test
137 // trajectory frames in the chosen way.
138 TrajectoryComparison trajectoryComparison{ trajectoryMatchSettings, trajectoryTolerances };
140 // Set simulation file names
141 auto simulationTrajectoryFileName = fileManager_.getTemporaryFilePath("trajectory.trr");
142 auto simulationEdrFileName = fileManager_.getTemporaryFilePath("energy.edr");
143 auto simulationDhdlFileName = fileManager_.getTemporaryFilePath("dhdl.xvg");
146 runner_.tprFileName_ = fileManager_.getTemporaryFilePath("sim.tpr");
147 runner_.useTopGroAndMdpFromFepTestDatabase(simulationName);
148 runGrompp(&runner_, { SimulationOptionTuple("-maxwarn", std::to_string(maxNumWarnings)) });
151 runner_.fullPrecisionTrajectoryFileName_ = simulationTrajectoryFileName;
152 runner_.edrFileName_ = simulationEdrFileName;
153 runner_.dhdlFileName_ = simulationDhdlFileName;
156 /* Currently used tests write trajectory (x/v/f) frames every 20 steps.
157 * Except for the expanded ensemble test, all tests run for 20 steps total.
158 * As the tolerances are relatively strict, we need to restrict the number of
159 * force frames we can expect to match.
160 * Testing more than the first force frame is only feasible in double precision
161 * using a single rank.
162 * Testing one force frame is only feasible in double precision.
163 * Note that this only concerns trajectory frames, energy frames are checked
165 const bool testTwoTrajectoryFrames = (GMX_DOUBLE && (getNumberOfTestMpiRanks() == 1));
166 const bool testOneTrajectoryFrame = GMX_DOUBLE;
168 // Compare simulation results
169 TestReferenceData refData;
170 TestReferenceChecker rootChecker(refData.rootChecker());
171 // Check that the energies agree with the refdata within tolerance.
172 checkEnergiesAgainstReferenceData(simulationEdrFileName, energyTermsToCompare, &rootChecker);
173 // Check that the trajectories agree with the refdata within tolerance.
174 if (testTwoTrajectoryFrames)
176 checkTrajectoryAgainstReferenceData(simulationTrajectoryFileName, trajectoryComparison,
177 &rootChecker, MaxNumFrames(2));
179 else if (testOneTrajectoryFrame)
181 checkTrajectoryAgainstReferenceData(simulationTrajectoryFileName, trajectoryComparison,
182 &rootChecker, MaxNumFrames(1));
186 checkTrajectoryAgainstReferenceData(simulationTrajectoryFileName, trajectoryComparison,
187 &rootChecker, MaxNumFrames(0));
189 if (File::exists(simulationDhdlFileName, File::returnFalseOnError))
191 TextInputFile dhdlFile(simulationDhdlFileName);
192 auto settings = XvgMatchSettings();
193 settings.tolerance = defaultEnergyTolerance;
194 checkXvgFile(&dhdlFile, &rootChecker, settings);
198 // TODO: The time for OpenCL kernel compilation means these tests time
199 // out. Once that compilation is cached for the whole process, these
200 // tests can run in such configurations.
202 INSTANTIATE_TEST_CASE_P(
203 FreeEnergyCalculationsAreEquivalentToReference,
204 FreeEnergyReferenceTest,
206 FreeEnergyReferenceTestParams{ "coulandvdwsequential_coul",
208 { F_DVDL_COUL, F_DVDL_VDW } },
209 FreeEnergyReferenceTestParams{ "coulandvdwsequential_vdw",
211 { F_DVDL_COUL, F_DVDL_VDW } },
212 FreeEnergyReferenceTestParams{ "coulandvdwtogether", MaxNumWarnings(0), { F_DVDL } },
213 FreeEnergyReferenceTestParams{ "expanded", MaxNumWarnings(0), { F_DVDL_COUL, F_DVDL_VDW } },
214 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
215 FreeEnergyReferenceTestParams{ "relative",
217 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED } },
218 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
219 FreeEnergyReferenceTestParams{
220 "relative-position-restraints",
222 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT } },
223 FreeEnergyReferenceTestParams{ "restraints", MaxNumWarnings(0), { F_DVDL_RESTRAINT } },
224 FreeEnergyReferenceTestParams{ "simtemp", MaxNumWarnings(0), {} },
225 FreeEnergyReferenceTestParams{ "transformAtoB", MaxNumWarnings(0), { F_DVDL } },
226 FreeEnergyReferenceTestParams{ "vdwalone", MaxNumWarnings(0), { F_DVDL } }),
227 FreeEnergyReferenceTest::PrintParametersToString());
229 INSTANTIATE_TEST_CASE_P(
230 DISABLED_FreeEnergyCalculationsAreEquivalentToReference,
231 FreeEnergyReferenceTest,
233 FreeEnergyReferenceTestParams{ "coulandvdwsequential_coul",
235 { F_DVDL_COUL, F_DVDL_VDW } },
236 FreeEnergyReferenceTestParams{ "coulandvdwsequential_vdw",
238 { F_DVDL_COUL, F_DVDL_VDW } },
239 FreeEnergyReferenceTestParams{ "coulandvdwtogether", MaxNumWarnings(0), { F_DVDL } },
240 FreeEnergyReferenceTestParams{ "expanded", MaxNumWarnings(0), { F_DVDL_COUL, F_DVDL_VDW } },
241 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
242 FreeEnergyReferenceTestParams{ "relative",
244 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED } },
245 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
246 FreeEnergyReferenceTestParams{
247 "relative-position-restraints",
249 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT } },
250 FreeEnergyReferenceTestParams{ "restraints", MaxNumWarnings(0), { F_DVDL_RESTRAINT } },
251 FreeEnergyReferenceTestParams{ "simtemp", MaxNumWarnings(0), {} },
252 FreeEnergyReferenceTestParams{ "transformAtoB", MaxNumWarnings(0), { F_DVDL } },
253 FreeEnergyReferenceTestParams{ "vdwalone", MaxNumWarnings(0), { F_DVDL } }),
254 FreeEnergyReferenceTest::PrintParametersToString());
258 } // namespace gmx::test