2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2017,2018,2019,2020,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.
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.
43 #include <gmock/gmock.h>
44 #include <gtest/gtest.h>
46 #include "gromacs/applied_forces/awh/bias.h"
47 #include "gromacs/applied_forces/awh/correlationgrid.h"
48 #include "gromacs/applied_forces/awh/pointstate.h"
49 #include "gromacs/mdtypes/awh_params.h"
50 #include "gromacs/utility/stringutil.h"
52 #include "gromacs/applied_forces/awh/tests/awh_setup.h"
53 #include "testutils/refdata.h"
54 #include "testutils/testasserts.h"
62 //! The number of lambda states to use in the tests.
63 static constexpr int c_numLambdaStates = 16;
66 //! Convenience typedef: growth type enum, potential type enum, disable update skips
67 typedef std::tuple<AwhHistogramGrowthType, AwhPotentialType, BiasParams::DisableUpdateSkips> BiasTestParameters;
69 /*! \brief Test fixture for testing Bias updates
71 class BiasFepLambdaStateTest : public ::testing::TestWithParam<BiasTestParameters>
74 //! Storage for test parameters.
75 std::unique_ptr<AwhTestParameters> params_;
78 //! Random seed for AWH MC sampling
82 std::unique_ptr<Bias> bias_;
84 BiasFepLambdaStateTest()
86 /* We test all combinations of:
88 * eawhgrowthLINEAR: final, normal update phase
89 * ewahgrowthEXP_LINEAR: intial phase, updated size is constant
90 * eawhpotential (test only eawhpotentialUMBRELLA (MC) for FEP lambda dimensions)
91 * disableUpdateSkips (should not affect the results):
92 * BiasParams::DisableUpdateSkips::yes: update the point state for every sample
93 * BiasParams::DisableUpdateSkips::no: update the point state at an interval > 1 sample
95 * Note: It would be nice to explicitly check that eawhpotential
96 * and disableUpdateSkips do not affect the point state.
97 * But the reference data will also ensure this.
99 AwhHistogramGrowthType eawhgrowth;
100 AwhPotentialType eawhpotential;
101 BiasParams::DisableUpdateSkips disableUpdateSkips;
102 std::tie(eawhgrowth, eawhpotential, disableUpdateSkips) = GetParam();
104 /* Set up a basic AWH setup with a single, 1D bias with parameters
105 * such that we can measure the effects of different parameters.
107 constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
108 constexpr int coordIndex = 0;
109 constexpr double origin = 0;
110 constexpr double end = c_numLambdaStates - 1;
111 constexpr double period = 0;
112 // Correction for removal of GaussianGeometryFactor/2 in histogram size
113 constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
114 const auto awhDimBuffer =
115 awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
116 const auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
117 params_ = std::make_unique<AwhTestParameters>(getAwhTestParameters(
118 eawhgrowth, eawhpotential, awhDimArrayRef, false, 0.4, true, 1.0, c_numLambdaStates));
120 seed_ = params_->awhParams.seed();
122 double mdTimeStep = 0.1;
124 bias_ = std::make_unique<Bias>(-1,
126 params_->awhParams.awhBiasParams()[0],
132 Bias::ThisRankWillDoIO::No,
137 TEST_P(BiasFepLambdaStateTest, ForcesBiasPmf)
139 gmx::test::TestReferenceData data;
140 gmx::test::TestReferenceChecker checker(data.rootChecker());
144 /* Make strings with the properties we expect to be different in the tests.
145 * These also helps to interpret the reference data.
147 std::vector<std::string> props;
148 props.push_back(formatString("stage: %s", bias.state().inInitialStage() ? "initial" : "final"));
149 props.push_back(formatString("convolve forces: %s", bias.params().convolveForce ? "yes" : "no"));
150 props.push_back(formatString("skip updates: %s", bias.params().skipUpdates() ? "yes" : "no"));
152 SCOPED_TRACE(gmx::formatString("%s, %s, %s", props[0].c_str(), props[1].c_str(), props[2].c_str()));
154 std::vector<double> force, pot;
156 double potentialJump = 0;
157 double mdTimeStep = 0.1;
160 /* Some energies to use as base values (to which some noise is added later on). */
161 std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
162 std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
163 const double magnitude = 12.0;
164 for (int i = 0; i < c_numLambdaStates; i++)
166 neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
167 neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
170 for (int step = 0; step < nSteps; step++)
172 int umbrellaGridpointIndex = bias.state().coordState().umbrellaGridpoint();
173 awh_dvec coordValue = { bias.getGridCoordValue(umbrellaGridpointIndex)[0], 0, 0, 0 };
174 double potential = 0;
175 gmx::ArrayRef<const double> biasForce = bias.calcForceAndUpdateBias(coordValue,
176 neighborLambdaEnergies,
185 force.push_back(biasForce[0]);
186 pot.push_back(potential);
189 /* When skipping updates, ensure all skipped updates are performed here.
190 * This should result in the same bias state as at output in a normal run.
192 if (bias.params().skipUpdates())
194 bias.doSkippedUpdatesForAllPoints();
197 std::vector<double> pointBias, logPmfsum;
198 for (const auto& point : bias.state().points())
200 pointBias.push_back(point.bias());
201 logPmfsum.push_back(point.logPmfSum());
204 constexpr int ulpTol = 10;
206 checker.checkSequence(props.begin(), props.end(), "Properties");
207 checker.setDefaultTolerance(absoluteTolerance(magnitude * GMX_DOUBLE_EPS * ulpTol));
208 checker.checkSequence(force.begin(), force.end(), "Force");
209 checker.checkSequence(pot.begin(), pot.end(), "Potential");
210 checker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpTol));
211 checker.checkSequence(pointBias.begin(), pointBias.end(), "PointBias");
212 checker.checkSequence(logPmfsum.begin(), logPmfsum.end(), "PointLogPmfsum");
215 /* Scan initial/final phase, MC/convolved force and update skip (not) allowed
216 * Both the convolving and skipping should not affect the bias and PMF.
217 * It would be nice if the test would explicitly check for this.
218 * Currently this is tested through identical reference data.
220 INSTANTIATE_TEST_SUITE_P(WithParameters,
221 BiasFepLambdaStateTest,
222 ::testing::Combine(::testing::Values(AwhHistogramGrowthType::Linear,
223 AwhHistogramGrowthType::ExponentialLinear),
224 ::testing::Values(AwhPotentialType::Umbrella),
225 ::testing::Values(BiasParams::DisableUpdateSkips::yes,
226 BiasParams::DisableUpdateSkips::no)));
228 // Test that we detect coverings and exit the initial stage at the correct step
229 TEST(BiasFepLambdaStateTest, DetectsCovering)
231 constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
232 constexpr int coordIndex = 0;
233 constexpr double origin = 0;
234 constexpr double end = c_numLambdaStates - 1;
235 constexpr double period = 0;
236 constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
238 awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
239 auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
240 const AwhTestParameters params(getAwhTestParameters(AwhHistogramGrowthType::ExponentialLinear,
241 AwhPotentialType::Umbrella,
249 const double mdTimeStep = 0.1;
253 params.awhParams.awhBiasParams()[0],
259 Bias::ThisRankWillDoIO::No);
261 const int64_t exitStepRef = 320;
263 bool inInitialStage = bias.state().inInitialStage();
265 /* Some energies to use as base values (to which some noise is added later on). */
266 std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
267 std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
268 const double magnitude = 12.0;
269 for (int i = 0; i < c_numLambdaStates; i++)
271 neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
272 neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
276 /* Normally this loop exits at exitStepRef, but we extend with failure */
277 for (step = 0; step <= 2 * exitStepRef; step++)
279 int umbrellaGridpointIndex = bias.state().coordState().umbrellaGridpoint();
280 awh_dvec coordValue = { bias.getGridCoordValue(umbrellaGridpointIndex)[0], 0, 0, 0 };
282 double potential = 0;
283 double potentialJump = 0;
284 bias.calcForceAndUpdateBias(coordValue,
285 neighborLambdaEnergies,
291 params.awhParams.seed(),
294 inInitialStage = bias.state().inInitialStage();
301 EXPECT_EQ(false, inInitialStage);
304 EXPECT_EQ(exitStepRef, step);