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 both, but for the FEP lambda state dimension MC will in practice be used,
91 * except that eawhpotentialCONVOLVED also gives a potential output):
92 * eawhpotentialUMBRELLA: MC on lambda state
93 * eawhpotentialCONVOLVED: MD on a convolved potential landscape (falling back to MC on lambda state)
94 * disableUpdateSkips (should not affect the results):
95 * BiasParams::DisableUpdateSkips::yes: update the point state for every sample
96 * BiasParams::DisableUpdateSkips::no: update the point state at an interval > 1 sample
98 * Note: It would be nice to explicitly check that eawhpotential
99 * and disableUpdateSkips do not affect the point state.
100 * But the reference data will also ensure this.
102 AwhHistogramGrowthType eawhgrowth;
103 AwhPotentialType eawhpotential;
104 BiasParams::DisableUpdateSkips disableUpdateSkips;
105 std::tie(eawhgrowth, eawhpotential, disableUpdateSkips) = GetParam();
107 /* Set up a basic AWH setup with a single, 1D bias with parameters
108 * such that we can measure the effects of different parameters.
110 constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
111 constexpr int coordIndex = 0;
112 constexpr double origin = 0;
113 constexpr double end = c_numLambdaStates - 1;
114 constexpr double period = 0;
115 // Correction for removal of GaussianGeometryFactor/2 in histogram size
116 constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
117 const auto awhDimBuffer =
118 awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
119 const auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
120 params_ = std::make_unique<AwhTestParameters>(getAwhTestParameters(
121 eawhgrowth, eawhpotential, awhDimArrayRef, false, 0.4, true, 1.0, c_numLambdaStates));
123 seed_ = params_->awhParams.seed();
125 double mdTimeStep = 0.1;
127 bias_ = std::make_unique<Bias>(-1,
129 params_->awhParams.awhBiasParams()[0],
135 Bias::ThisRankWillDoIO::No,
140 TEST_P(BiasFepLambdaStateTest, ForcesBiasPmf)
142 gmx::test::TestReferenceData data;
143 gmx::test::TestReferenceChecker checker(data.rootChecker());
147 /* Make strings with the properties we expect to be different in the tests.
148 * These also helps to interpret the reference data.
150 std::vector<std::string> props;
151 props.push_back(formatString("stage: %s", bias.state().inInitialStage() ? "initial" : "final"));
152 props.push_back(formatString("convolve forces: %s", bias.params().convolveForce ? "yes" : "no"));
153 props.push_back(formatString("skip updates: %s", bias.params().skipUpdates() ? "yes" : "no"));
155 SCOPED_TRACE(gmx::formatString("%s, %s, %s", props[0].c_str(), props[1].c_str(), props[2].c_str()));
157 std::vector<double> force, pot;
159 double potentialJump = 0;
160 double mdTimeStep = 0.1;
163 /* Some energies to use as base values (to which some noise is added later on). */
164 std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
165 std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
166 const double magnitude = 12.0;
167 for (int i = 0; i < c_numLambdaStates; i++)
169 neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
170 neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
173 for (int step = 0; step < nSteps; step++)
175 int umbrellaGridpointIndex = bias.state().coordState().umbrellaGridpoint();
176 awh_dvec coordValue = { bias.getGridCoordValue(umbrellaGridpointIndex)[0], 0, 0, 0 };
177 double potential = 0;
178 gmx::ArrayRef<const double> biasForce = bias.calcForceAndUpdateBias(coordValue,
179 neighborLambdaEnergies,
190 force.push_back(biasForce[0]);
191 pot.push_back(potential);
194 /* When skipping updates, ensure all skipped updates are performed here.
195 * This should result in the same bias state as at output in a normal run.
197 if (bias.params().skipUpdates())
199 bias.doSkippedUpdatesForAllPoints();
202 std::vector<double> pointBias, logPmfsum;
203 for (const auto& point : bias.state().points())
205 pointBias.push_back(point.bias());
206 logPmfsum.push_back(point.logPmfSum());
209 constexpr int ulpTol = 10;
211 checker.checkSequence(props.begin(), props.end(), "Properties");
212 checker.setDefaultTolerance(absoluteTolerance(magnitude * GMX_DOUBLE_EPS * ulpTol));
213 checker.checkSequence(force.begin(), force.end(), "Force");
214 checker.checkSequence(pot.begin(), pot.end(), "Potential");
215 checker.setDefaultTolerance(relativeToleranceAsUlp(1.0, ulpTol));
216 checker.checkSequence(pointBias.begin(), pointBias.end(), "PointBias");
217 checker.checkSequence(logPmfsum.begin(), logPmfsum.end(), "PointLogPmfsum");
220 /* Scan initial/final phase, MC/convolved force and update skip (not) allowed
221 * Both the convolving and skipping should not affect the bias and PMF.
222 * It would be nice if the test would explicitly check for this.
223 * Currently this is tested through identical reference data.
225 INSTANTIATE_TEST_CASE_P(WithParameters,
226 BiasFepLambdaStateTest,
227 ::testing::Combine(::testing::Values(AwhHistogramGrowthType::Linear,
228 AwhHistogramGrowthType::ExponentialLinear),
229 ::testing::Values(AwhPotentialType::Umbrella,
230 AwhPotentialType::Convolved),
231 ::testing::Values(BiasParams::DisableUpdateSkips::yes,
232 BiasParams::DisableUpdateSkips::no)));
234 // Test that we detect coverings and exit the initial stage at the correct step
235 TEST(BiasFepLambdaStateTest, DetectsCovering)
237 constexpr AwhCoordinateProviderType coordinateProvider = AwhCoordinateProviderType::FreeEnergyLambda;
238 constexpr int coordIndex = 0;
239 constexpr double origin = 0;
240 constexpr double end = c_numLambdaStates - 1;
241 constexpr double period = 0;
242 constexpr double diffusion = 1e-4 / (0.12927243028700 * 2);
244 awhDimParamSerialized(coordinateProvider, coordIndex, origin, end, period, diffusion);
245 auto awhDimArrayRef = gmx::arrayRefFromArray(&awhDimBuffer, 1);
246 const AwhTestParameters params(getAwhTestParameters(AwhHistogramGrowthType::ExponentialLinear,
247 AwhPotentialType::Convolved,
255 const double mdTimeStep = 0.1;
259 params.awhParams.awhBiasParams()[0],
265 Bias::ThisRankWillDoIO::No);
267 const int64_t exitStepRef = 320;
269 bool inInitialStage = bias.state().inInitialStage();
271 /* Some energies to use as base values (to which some noise is added later on). */
272 std::vector<double> neighborLambdaEnergies(c_numLambdaStates);
273 std::vector<double> neighborLambdaDhdl(c_numLambdaStates);
274 const double magnitude = 12.0;
275 for (int i = 0; i < c_numLambdaStates; i++)
277 neighborLambdaEnergies[i] = magnitude * std::sin(i * 0.1);
278 neighborLambdaDhdl[i] = magnitude * std::cos(i * 0.1);
282 /* Normally this loop exits at exitStepRef, but we extend with failure */
283 for (step = 0; step <= 2 * exitStepRef; step++)
285 int umbrellaGridpointIndex = bias.state().coordState().umbrellaGridpoint();
286 awh_dvec coordValue = { bias.getGridCoordValue(umbrellaGridpointIndex)[0], 0, 0, 0 };
288 double potential = 0;
289 double potentialJump = 0;
290 bias.calcForceAndUpdateBias(coordValue,
291 neighborLambdaEnergies,
299 params.awhParams.seed(),
302 inInitialStage = bias.state().inInitialStage();
309 EXPECT_EQ(false, inInitialStage);
312 EXPECT_EQ(exitStepRef, step);