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.",
112 SCOPED_TRACE(formatString("Comparing FEP simulation '%s' to reference", simulationName.c_str()));
114 // Tolerance set to pass with identical code version and a range of different test setups for most tests
115 const auto defaultEnergyTolerance = relativeToleranceAsFloatingPoint(50.0, GMX_DOUBLE ? 1e-5 : 1e-4);
116 // Some simulations are significantly longer, so they need a larger tolerance
117 const auto longEnergyTolerance = relativeToleranceAsFloatingPoint(50.0, GMX_DOUBLE ? 1e-4 : 1e-3);
118 const bool isLongSimulation = (simulationName == "expanded");
119 const auto energyTolerance = isLongSimulation ? longEnergyTolerance : defaultEnergyTolerance;
121 EnergyTermsToCompare energyTermsToCompare{ { interaction_function[F_EPOT].longname, energyTolerance } };
122 for (const auto& interaction : interactionsList)
124 energyTermsToCompare.emplace(interaction_function[interaction].longname, defaultEnergyTolerance);
127 // Specify how trajectory frame matching must work (only testing forces).
128 TrajectoryFrameMatchSettings trajectoryMatchSettings{ false,
131 ComparisonConditions::NoComparison,
132 ComparisonConditions::NoComparison,
133 ComparisonConditions::MustCompare };
134 TrajectoryTolerances trajectoryTolerances = TrajectoryComparison::s_defaultTrajectoryTolerances;
135 trajectoryTolerances.forces = relativeToleranceAsFloatingPoint(100.0, GMX_DOUBLE ? 5.0e-5 : 5.0e-4);
137 // Build the functor that will compare reference and test
138 // trajectory frames in the chosen way.
139 TrajectoryComparison trajectoryComparison{ trajectoryMatchSettings, trajectoryTolerances };
141 // Set simulation file names
142 auto simulationTrajectoryFileName = fileManager_.getTemporaryFilePath("trajectory.trr");
143 auto simulationEdrFileName = fileManager_.getTemporaryFilePath("energy.edr");
144 auto simulationDhdlFileName = fileManager_.getTemporaryFilePath("dhdl.xvg");
147 runner_.tprFileName_ = fileManager_.getTemporaryFilePath("sim.tpr");
148 runner_.useTopGroAndMdpFromFepTestDatabase(simulationName);
149 runGrompp(&runner_, { SimulationOptionTuple("-maxwarn", std::to_string(maxNumWarnings)) });
152 runner_.fullPrecisionTrajectoryFileName_ = simulationTrajectoryFileName;
153 runner_.edrFileName_ = simulationEdrFileName;
154 runner_.dhdlFileName_ = simulationDhdlFileName;
157 /* Currently used tests write trajectory (x/v/f) frames every 20 steps.
158 * Except for the expanded ensemble test, all tests run for 20 steps total.
159 * As the tolerances are relatively strict, we need to restrict the number of
160 * force frames we can expect to match.
161 * Testing more than the first force frame is only feasible in double precision
162 * using a single rank.
163 * Testing one force frame is only feasible in double precision.
164 * Note that this only concerns trajectory frames, energy frames are checked
166 const bool testTwoTrajectoryFrames = (GMX_DOUBLE && (getNumberOfTestMpiRanks() == 1));
167 const bool testOneTrajectoryFrame = GMX_DOUBLE;
169 // Compare simulation results
170 TestReferenceData refData;
171 TestReferenceChecker rootChecker(refData.rootChecker());
172 // Check that the energies agree with the refdata within tolerance.
173 checkEnergiesAgainstReferenceData(simulationEdrFileName, energyTermsToCompare, &rootChecker);
174 // Check that the trajectories agree with the refdata within tolerance.
175 if (testTwoTrajectoryFrames)
177 checkTrajectoryAgainstReferenceData(
178 simulationTrajectoryFileName, trajectoryComparison, &rootChecker, MaxNumFrames(2));
180 else if (testOneTrajectoryFrame)
182 checkTrajectoryAgainstReferenceData(
183 simulationTrajectoryFileName, trajectoryComparison, &rootChecker, MaxNumFrames(1));
187 checkTrajectoryAgainstReferenceData(
188 simulationTrajectoryFileName, trajectoryComparison, &rootChecker, MaxNumFrames(0));
190 if (File::exists(simulationDhdlFileName, File::returnFalseOnError))
192 TextInputFile dhdlFile(simulationDhdlFileName);
193 auto settings = XvgMatchSettings();
194 settings.tolerance = defaultEnergyTolerance;
195 checkXvgFile(&dhdlFile, &rootChecker, settings);
199 // TODO: The time for OpenCL kernel compilation means these tests time
200 // out. Once that compilation is cached for the whole process, these
201 // tests can run in such configurations.
203 INSTANTIATE_TEST_CASE_P(
204 FreeEnergyCalculationsAreEquivalentToReference,
205 FreeEnergyReferenceTest,
207 FreeEnergyReferenceTestParams{ "coulandvdwsequential_coul",
209 { F_DVDL_COUL, F_DVDL_VDW } },
210 FreeEnergyReferenceTestParams{ "coulandvdwsequential_vdw",
212 { F_DVDL_COUL, F_DVDL_VDW } },
213 FreeEnergyReferenceTestParams{ "coulandvdwtogether", MaxNumWarnings(0), { F_DVDL } },
214 FreeEnergyReferenceTestParams{ "expanded", MaxNumWarnings(0), { F_DVDL_COUL, F_DVDL_VDW } },
215 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
216 FreeEnergyReferenceTestParams{ "relative",
218 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED } },
219 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
220 FreeEnergyReferenceTestParams{
221 "relative-position-restraints",
223 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT } },
224 FreeEnergyReferenceTestParams{ "restraints", MaxNumWarnings(0), { F_DVDL_RESTRAINT } },
225 FreeEnergyReferenceTestParams{ "simtemp", MaxNumWarnings(0), {} },
226 FreeEnergyReferenceTestParams{ "transformAtoB", MaxNumWarnings(0), { F_DVDL } },
227 FreeEnergyReferenceTestParams{ "vdwalone", MaxNumWarnings(0), { F_DVDL } }),
228 FreeEnergyReferenceTest::PrintParametersToString());
230 INSTANTIATE_TEST_CASE_P(
231 DISABLED_FreeEnergyCalculationsAreEquivalentToReference,
232 FreeEnergyReferenceTest,
234 FreeEnergyReferenceTestParams{ "coulandvdwsequential_coul",
236 { F_DVDL_COUL, F_DVDL_VDW } },
237 FreeEnergyReferenceTestParams{ "coulandvdwsequential_vdw",
239 { F_DVDL_COUL, F_DVDL_VDW } },
240 FreeEnergyReferenceTestParams{ "coulandvdwtogether", MaxNumWarnings(0), { F_DVDL } },
241 FreeEnergyReferenceTestParams{ "expanded", MaxNumWarnings(0), { F_DVDL_COUL, F_DVDL_VDW } },
242 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
243 FreeEnergyReferenceTestParams{ "relative",
245 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED } },
246 // Tolerated warnings: No default bonded interaction types for perturbed atoms (10x)
247 FreeEnergyReferenceTestParams{
248 "relative-position-restraints",
250 { F_DVDL, F_DVDL_COUL, F_DVDL_VDW, F_DVDL_BONDED, F_DVDL_RESTRAINT } },
251 FreeEnergyReferenceTestParams{ "restraints", MaxNumWarnings(0), { F_DVDL_RESTRAINT } },
252 FreeEnergyReferenceTestParams{ "simtemp", MaxNumWarnings(0), {} },
253 FreeEnergyReferenceTestParams{ "transformAtoB", MaxNumWarnings(0), { F_DVDL } },
254 FreeEnergyReferenceTestParams{ "vdwalone", MaxNumWarnings(0), { F_DVDL } }),
255 FreeEnergyReferenceTest::PrintParametersToString());
259 } // namespace gmx::test