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.eTarget)
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 const std::vector<DimParams>& dimParams,
120 const std::vector<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 const std::vector<GridAxis>& gridAxis,
180 double samplingTimestep)
182 /* Get diffusion factor */
183 double maxCrossingTime = 0.;
184 std::vector<double> x;
185 for (size_t d = 0; d < gridAxis.size(); d++)
187 GMX_RELEASE_ASSERT(awhBiasParams.dimParams[d].diffusion > 0, "We need positive diffusion");
188 // With diffusion it takes on average T = L^2/2D time to cross length L
189 double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
190 double crossingTime = (axisLength * axisLength) / (2 * awhBiasParams.dimParams[d].diffusion);
191 maxCrossingTime = std::max(maxCrossingTime, crossingTime);
193 GMX_RELEASE_ASSERT(maxCrossingTime > 0, "We need at least one dimension with non-zero length");
194 double errorInitialInKT = beta * awhBiasParams.errorInitial;
195 double histogramSize = maxCrossingTime / (gmx::square(errorInitialInKT) * samplingTimestep);
197 return histogramSize;
200 /*! \brief Returns the number of simulations sharing bias updates.
202 * \param[in] awhBiasParams Bias parameters.
203 * \param[in] numSharingSimulations The number of simulations to share the bias across.
204 * \returns the number of shared updates.
206 int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulations)
208 GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
212 if (awhBiasParams.shareGroup > 0)
214 /* We do not yet support sharing within a simulation */
215 int numSharedWithinThisSimulation = 1;
216 numShared = numSharingSimulations * numSharedWithinThisSimulation;
224 BiasParams::BiasParams(const AwhParams& awhParams,
225 const AwhBiasParams& awhBiasParams,
226 const std::vector<DimParams>& dimParams,
229 DisableUpdateSkips disableUpdateSkips,
230 int numSharingSimulations,
231 const std::vector<GridAxis>& gridAxis,
233 invBeta(beta > 0 ? 1 / beta : 0),
234 numStepsSampleCoord_(awhParams.nstSampleCoord),
235 numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy),
236 numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
237 numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
238 eTarget(awhBiasParams.eTarget),
239 freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff),
240 temperatureScaleFactor(awhBiasParams.targetBetaScaling),
241 idealWeighthistUpdate(eTarget != AwhTargetType::LocalBoltzmann),
242 numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
243 updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
244 localWeightScaling(eTarget == AwhTargetType::LocalBoltzmann ? temperatureScaleFactor : 1),
245 initialErrorInKT(beta * awhBiasParams.errorInitial),
246 initialHistogramSize(
247 getInitialHistogramSizeEstimate(awhBiasParams, gridAxis, beta, numStepsSampleCoord_ * mdTimeStep)),
248 convolveForce(awhParams.ePotential == AwhPotentialType::Convolved),
249 biasIndex(biasIndex),
250 disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
254 GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
257 for (int d = 0; d < awhBiasParams.ndim; d++)
259 /* The spacing in FEP dimensions is 1. The calculated coverRadius will be in lambda states
260 * (cf points in other dimensions). */
261 double coverRadiusInNm = 0.5 * awhBiasParams.dimParams[d].coverDiameter;
262 double spacing = gridAxis[d].spacing();
263 coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;