2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,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.
38 * Implements the initialization of the BiasParams class.
40 * \author Viveca Lindahl
41 * \author Berk Hess <hess@kth.se>
47 #include "biasparams.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/mdtypes/awh_params.h"
55 #include "gromacs/utility/arrayref.h"
56 #include "gromacs/utility/exceptions.h"
57 #include "gromacs/utility/gmxassert.h"
67 /*! \brief Determines the interval for updating the target distribution.
69 * The interval value is based on the target distrbution type
70 * (this could be made a user-option but there is most likely
71 * no big need for tweaking this for most users).
73 * \param[in] awhParams AWH parameters.
74 * \param[in] awhBiasParams Bias parameters.
75 * \returns the target update interval in steps.
77 int64_t calcTargetUpdateInterval(const AwhParams& awhParams, const AwhBiasParams& awhBiasParams)
79 int64_t numStepsUpdateTarget = 0;
80 /* Set the target update frequency based on the target distrbution type
81 * (this could be made a user-option but there is most likely no big need
82 * for tweaking this for most users).
84 switch (awhBiasParams.targetDistribution())
86 case AwhTargetType::Constant: numStepsUpdateTarget = 0; break;
87 case AwhTargetType::Cutoff:
88 case AwhTargetType::Boltzmann:
89 /* Updating the target generally requires updating the whole grid so to keep the cost
90 down we generally update the target less often than the free energy (unless the free
91 energy update step is set to > 100 samples). */
92 numStepsUpdateTarget = std::max(100 % awhParams.numSamplesUpdateFreeEnergy(),
93 awhParams.numSamplesUpdateFreeEnergy())
94 * awhParams.nstSampleCoord();
96 case AwhTargetType::LocalBoltzmann:
97 /* The target distribution is set equal to the reference histogram which is updated every free energy update.
98 So the target has to be updated at the same time. This leads to a global update each time because it is
99 assumed that a target distribution update can take any form. This is a bit unfortunate for a "local"
100 target distribution. One could avoid the global update by making a local target update function (and
101 postponing target updates for non-local points as for the free energy update). We avoid such additions
102 for now and accept that this target type always does global updates. */
103 numStepsUpdateTarget = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
105 default: GMX_RELEASE_ASSERT(false, "Unknown AWH target type"); break;
108 return numStepsUpdateTarget;
111 /*! \brief Determines the step interval for checking for covering.
113 * \param[in] awhParams AWH parameters.
114 * \param[in] dimParams Parameters for the dimensions of the coordinate.
115 * \param[in] gridAxis The BiasGrid axes.
116 * \returns the check interval in steps.
118 int64_t calcCheckCoveringInterval(const AwhParams& awhParams,
119 ArrayRef<const DimParams> dimParams,
120 ArrayRef<const GridAxis> gridAxis)
122 /* Each sample will have a width of sigma. To cover the axis a
123 minimum number of samples of width sigma is required. */
124 int minNumSamplesCover = 0;
125 for (size_t d = 0; d < gridAxis.size(); d++)
128 if (dimParams[d].isPullDimension())
131 dimParams[d].pullDimParams().betak > 0,
132 "Inverse temperature (beta) and force constant (k) should be positive.");
133 double sigma = 1.0 / std::sqrt(dimParams[d].pullDimParams().betak);
135 /* The additional sample is here because to cover a discretized
136 axis of length sigma one needs two samples, one for each
138 GMX_RELEASE_ASSERT(gridAxis[d].length() / sigma < std::numeric_limits<int>::max(),
139 "The axis length in units of sigma should fit in an int");
140 numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length() / sigma)) + 1;
144 numSamplesCover = dimParams[d].fepDimParams().numFepLambdaStates;
146 /* The minimum number of samples needed for simultaneously
147 covering all axes is limited by the axis requiring most
149 minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
152 /* Convert to number of steps using the sampling frequency. The
153 check interval should be a multiple of the update step
155 int numStepsUpdate = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
156 GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy() > 0,
157 "When checking for AWH coverings, the number of samples per AWH update need "
159 int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy());
160 int numStepsCheck = numUpdatesCheck * numStepsUpdate;
162 GMX_RELEASE_ASSERT(numStepsCheck % numStepsUpdate == 0,
163 "Only check covering at free energy update steps");
165 return numStepsCheck;
169 * Estimate a reasonable initial reference weight histogram size.
171 * \param[in] awhBiasParams Bias parameters.
172 * \param[in] gridAxis The BiasGrid axes.
173 * \param[in] beta 1/(k_B T).
174 * \param[in] samplingTimestep Sampling frequency of probability weights.
175 * \returns estimate of initial histogram size.
177 double getInitialHistogramSizeEstimate(const AwhBiasParams& awhBiasParams,
178 ArrayRef<const GridAxis> gridAxis,
180 double samplingTimestep)
182 /* Get diffusion factor */
183 double maxCrossingTime = 0.;
184 std::vector<double> x;
185 const auto awhDimParams = awhBiasParams.dimParams();
186 for (size_t d = 0; d < gridAxis.size(); d++)
188 GMX_RELEASE_ASSERT(awhDimParams[d].diffusion() > 0, "We need positive diffusion");
189 // With diffusion it takes on average T = L^2/2D time to cross length L
190 double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
191 double crossingTime = (axisLength * axisLength) / (2 * awhDimParams[d].diffusion());
192 maxCrossingTime = std::max(maxCrossingTime, crossingTime);
194 GMX_RELEASE_ASSERT(maxCrossingTime > 0, "We need at least one dimension with non-zero length");
195 double errorInitialInKT = beta * awhBiasParams.initialErrorEstimate();
196 double histogramSize = maxCrossingTime / (gmx::square(errorInitialInKT) * samplingTimestep);
198 return histogramSize;
201 /*! \brief Returns the number of simulations sharing bias updates.
203 * \param[in] awhBiasParams Bias parameters.
204 * \param[in] numSharingSimulations The number of simulations to share the bias across.
205 * \returns the number of shared updates.
207 int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulations)
209 GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
213 if (awhBiasParams.shareGroup() > 0)
215 /* We do not yet support sharing within a simulation */
216 int numSharedWithinThisSimulation = 1;
217 numShared = numSharingSimulations * numSharedWithinThisSimulation;
225 BiasParams::BiasParams(const AwhParams& awhParams,
226 const AwhBiasParams& awhBiasParams,
227 ArrayRef<const DimParams> dimParams,
230 DisableUpdateSkips disableUpdateSkips,
231 int numSharingSimulations,
232 ArrayRef<const GridAxis> gridAxis,
234 invBeta(beta > 0 ? 1 / beta : 0),
235 numStepsSampleCoord_(awhParams.nstSampleCoord()),
236 numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy()),
237 numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
238 numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
239 eTarget(awhBiasParams.targetDistribution()),
240 freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff()),
241 temperatureScaleFactor(awhBiasParams.targetBetaScaling()),
242 idealWeighthistUpdate(eTarget != AwhTargetType::LocalBoltzmann),
243 numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
244 updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
245 localWeightScaling(eTarget == AwhTargetType::LocalBoltzmann ? temperatureScaleFactor : 1),
246 initialErrorInKT(beta * awhBiasParams.initialErrorEstimate()),
247 initialHistogramSize(
248 getInitialHistogramSizeEstimate(awhBiasParams, gridAxis, beta, numStepsSampleCoord_ * mdTimeStep)),
249 convolveForce(awhParams.potential() == AwhPotentialType::Convolved),
250 biasIndex(biasIndex),
251 disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
255 GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
258 const auto& awhDimParams = awhBiasParams.dimParams();
259 for (int d = 0; d < gmx::ssize(awhDimParams); d++)
261 /* The spacing in FEP dimensions is 1. The calculated coverRadius will be in lambda states
262 * (cf points in other dimensions). */
263 double coverRadiusInNm = 0.5 * awhDimParams[d].coverDiameter();
264 double spacing = gridAxis[d].spacing();
265 coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;