2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2017,2018,2019, 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 eawhtargetCONSTANT: numStepsUpdateTarget = 0; break;
87 case eawhtargetCUTOFF:
88 case eawhtargetBOLTZMANN:
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 eawhtargetLOCALBOLTZMANN:
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 Grid 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++)
127 GMX_RELEASE_ASSERT(dimParams[d].betak > 0,
128 "Inverse temperature (beta) and force constant (k) should be positive.");
129 double sigma = 1.0 / std::sqrt(dimParams[d].betak);
131 /* The additional sample is here because to cover a discretized
132 axis of length sigma one needs two samples, one for each
134 GMX_RELEASE_ASSERT(gridAxis[d].length() / sigma < std::numeric_limits<int>::max(),
135 "The axis length in units of sigma should fit in an int");
136 int numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length() / sigma)) + 1;
138 /* The minimum number of samples needed for simultaneously
139 covering all axes is limited by the axis requiring most
141 minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
144 /* Convert to number of steps using the sampling frequency. The
145 check interval should be a multiple of the update step
147 int numStepsUpdate = awhParams.numSamplesUpdateFreeEnergy * awhParams.nstSampleCoord;
148 GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy > 0,
149 "When checking for AWH coverings, the number of samples per AWH update need "
151 int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy);
152 int numStepsCheck = numUpdatesCheck * numStepsUpdate;
154 GMX_RELEASE_ASSERT(numStepsCheck % numStepsUpdate == 0,
155 "Only check covering at free energy update steps");
157 return numStepsCheck;
161 * Returns an approximation of the geometry factor used for initializing the AWH update size.
163 * The geometry factor is defined as the following sum of Gaussians:
164 * sum_{k!=0} exp(-0.5*(k*pi*x)^2)/(pi*k)^2,
165 * where k is a xArray.size()-dimensional integer vector with k_i in {0,1,..}.
167 * \param[in] xArray Array to evaluate.
168 * \returns the geometry factor.
170 double gaussianGeometryFactor(gmx::ArrayRef<const double> xArray)
172 /* For convenience we give the geometry factor function a name: zeta(x) */
173 constexpr size_t tableSize = 5;
174 std::array<const double, tableSize> xTabulated = { { 1e-5, 1e-4, 1e-3, 1e-2, 1e-1 } };
175 std::array<const double, tableSize> zetaTable1d = { { 0.166536811948, 0.16653116886, 0.166250075882,
176 0.162701098306, 0.129272430287 } };
177 std::array<const double, tableSize> zetaTable2d = { { 2.31985974274, 1.86307292523, 1.38159772648,
178 0.897554759158, 0.405578211115 } };
180 gmx::ArrayRef<const double> zetaTable;
182 if (xArray.size() == 1)
184 zetaTable = zetaTable1d;
186 else if (xArray.size() == 2)
188 zetaTable = zetaTable2d;
192 /* TODO... but this is anyway a rough estimate and > 2 dimensions is not so popular. */
193 zetaTable = zetaTable2d;
196 /* TODO. Really zeta is a function of an ndim-dimensional vector x and we shoudl have a ndim-dimensional lookup-table.
197 Here we take the geometric average of the components of x which is ok if the x-components are not very different. */
199 for (const double& x : xArray)
204 GMX_ASSERT(!xArray.empty(), "We should have a non-empty input array");
205 xScalar = std::pow(xScalar, 1.0 / xArray.size());
207 /* Look up zeta(x) */
209 while ((xIndex < xTabulated.size()) && (xScalar > xTabulated[xIndex]))
215 if (xIndex == xTabulated.size())
217 /* Take last value */
218 zEstimate = zetaTable[xTabulated.size() - 1];
220 else if (xIndex == 0)
222 zEstimate = zetaTable[xIndex];
227 double x0 = xTabulated[xIndex - 1];
228 double x1 = xTabulated[xIndex];
229 double w = (xScalar - x0) / (x1 - x0);
230 zEstimate = w * zetaTable[xIndex - 1] + (1 - w) * zetaTable[xIndex];
237 * Estimate a reasonable initial reference weight histogram size.
239 * \param[in] dimParams Parameters for the dimensions of the coordinate.
240 * \param[in] awhBiasParams Bias parameters.
241 * \param[in] gridAxis The Grid axes.
242 * \param[in] beta 1/(k_B T).
243 * \param[in] samplingTimestep Sampling frequency of probability weights.
244 * \returns estimate of initial histogram size.
246 double getInitialHistogramSizeEstimate(const std::vector<DimParams>& dimParams,
247 const AwhBiasParams& awhBiasParams,
248 const std::vector<GridAxis>& gridAxis,
250 double samplingTimestep)
252 /* Get diffusion factor */
253 double crossingTime = 0.;
254 std::vector<double> x;
255 for (size_t d = 0; d < gridAxis.size(); d++)
257 double axisLength = gridAxis[d].length();
260 crossingTime += awhBiasParams.dimParams[d].diffusion / (axisLength * axisLength);
261 /* The sigma of the Gaussian distribution in the umbrella */
262 double sigma = 1. / std::sqrt(dimParams[d].betak);
263 x.push_back(sigma / axisLength);
266 GMX_RELEASE_ASSERT(crossingTime > 0, "We need at least one dimension with non-zero length");
267 double errorInitialInKT = beta * awhBiasParams.errorInitial;
268 double histogramSize = gaussianGeometryFactor(x)
269 / (crossingTime * gmx::square(errorInitialInKT) * samplingTimestep);
271 return histogramSize;
274 /*! \brief Returns the number of simulations sharing bias updates.
276 * \param[in] awhBiasParams Bias parameters.
277 * \param[in] numSharingSimulations The number of simulations to share the bias across.
278 * \returns the number of shared updates.
280 int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulations)
282 GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
286 if (awhBiasParams.shareGroup > 0)
288 /* We do not yet support sharing within a simulation */
289 int numSharedWithinThisSimulation = 1;
290 numShared = numSharingSimulations * numSharedWithinThisSimulation;
298 BiasParams::BiasParams(const AwhParams& awhParams,
299 const AwhBiasParams& awhBiasParams,
300 const std::vector<DimParams>& dimParams,
303 DisableUpdateSkips disableUpdateSkips,
304 int numSharingSimulations,
305 const std::vector<GridAxis>& gridAxis,
307 invBeta(beta > 0 ? 1 / beta : 0),
308 numStepsSampleCoord_(awhParams.nstSampleCoord),
309 numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy),
310 numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
311 numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
312 eTarget(awhBiasParams.eTarget),
313 freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff),
314 temperatureScaleFactor(awhBiasParams.targetBetaScaling),
315 idealWeighthistUpdate(eTarget != eawhtargetLOCALBOLTZMANN),
316 numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
317 updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
318 localWeightScaling(eTarget == eawhtargetLOCALBOLTZMANN ? temperatureScaleFactor : 1),
319 initialErrorInKT(beta * awhBiasParams.errorInitial),
320 initialHistogramSize(getInitialHistogramSizeEstimate(dimParams,
324 numStepsSampleCoord_ * mdTimeStep)),
325 convolveForce(awhParams.ePotential == eawhpotentialCONVOLVED),
326 biasIndex(biasIndex),
327 disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
331 GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
334 for (int d = 0; d < awhBiasParams.ndim; d++)
336 double coverRadiusInNm = 0.5 * awhBiasParams.dimParams[d].coverDiameter;
337 double spacing = gridAxis[d].spacing();
338 coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;