estimate used by AWH to initialize :math:`N` in terms of more meaningful
quantities
-.. math:: \frac{1}{N_0} = \frac{1}{N_0(\varepsilon_0, D)} \sim D\varepsilon_0^2.
+.. math:: \frac{1}{N_0} = \frac{1}{N_0(\varepsilon_0, D)} = \frac{1}{\Delta
+ t_\mathrm{sample}} \max_d \frac{L_d^2}{2D_d} \varepsilon_0^2
:label: eqawhn0
-Essentially, this tells us that a slower system (small :math:`D`)
+where :math:`L_d` is the length of the interval and :math:`D_d` is
+the diffusion constant along dimension :math:`d` of the AWH bias.
+For one dimension, :math:`L^2/2D` is the average time to diffuse
+over a distance of :math:`L`. We then takes the maximum crossing
+time over all dimensions involved in the bias.
+Essentially, this formula tells us that a slower system (small :math:`D`)
requires more samples (larger :math:`N^0`) to attain the same level of
accuracy (:math:`\varepsilon_0`) at a given sampling rate. Conversely,
for a system of given diffusion, how to choose the initial biasing rate
maximum free energy difference of the PMF estimate. If this is much
larger than the expected magnitude of the free energy barriers that
should be crossed, then the system is probably being pulled too hard and
-:math:`D` should be decreased. :math:`\varepsilon_0` on the other hand,
-would only be tweaked when starting an AWH simulation using a fairly
-accurate guess of the PMF as input.
+:math:`D` should be decreased. An accurate estimate of the diffusion
+can be obtaining from an AWH simulation with the :ref:`gmx awh` tool.
+:math:`\varepsilon_0` on the other hand, should be a rough estimate
+of the initial error.
Tips for efficient sampling
^^^^^^^^^^^^^^^^^^^^^^^^^^^
return numStepsCheck;
}
-/*! \brief
- * Returns an approximation of the geometry factor used for initializing the AWH update size.
- *
- * The geometry factor is defined as the following sum of Gaussians:
- * sum_{k!=0} exp(-0.5*(k*pi*x)^2)/(pi*k)^2,
- * where k is a xArray.size()-dimensional integer vector with k_i in {0,1,..}.
- *
- * \param[in] xArray Array to evaluate.
- * \returns the geometry factor.
- */
-double gaussianGeometryFactor(gmx::ArrayRef<const double> xArray)
-{
- /* For convenience we give the geometry factor function a name: zeta(x) */
- constexpr size_t tableSize = 5;
- std::array<const double, tableSize> xTabulated = { { 1e-5, 1e-4, 1e-3, 1e-2, 1e-1 } };
- std::array<const double, tableSize> zetaTable1d = { { 0.166536811948, 0.16653116886, 0.166250075882,
- 0.162701098306, 0.129272430287 } };
- std::array<const double, tableSize> zetaTable2d = { { 2.31985974274, 1.86307292523, 1.38159772648,
- 0.897554759158, 0.405578211115 } };
-
- gmx::ArrayRef<const double> zetaTable;
-
- if (xArray.size() == 1)
- {
- zetaTable = zetaTable1d;
- }
- else if (xArray.size() == 2)
- { // NOLINT bugprone-branch-clone
- zetaTable = zetaTable2d;
- }
- else
- {
- /* TODO... but this is anyway a rough estimate and > 2 dimensions is not so popular.
- * Remove the above NOLINT when addressing this */
- zetaTable = zetaTable2d;
- }
-
- /* TODO. Really zeta is a function of an ndim-dimensional vector x and we shoudl have a ndim-dimensional lookup-table.
- Here we take the geometric average of the components of x which is ok if the x-components are not very different. */
- double xScalar = 1;
- for (const double& x : xArray)
- {
- xScalar *= x;
- }
-
- GMX_ASSERT(!xArray.empty(), "We should have a non-empty input array");
- xScalar = std::pow(xScalar, 1.0 / xArray.size());
-
- /* Look up zeta(x) */
- size_t xIndex = 0;
- while ((xIndex < xTabulated.size()) && (xScalar > xTabulated[xIndex]))
- {
- xIndex++;
- }
-
- double zEstimate;
- if (xIndex == xTabulated.size())
- {
- /* Take last value */
- zEstimate = zetaTable[xTabulated.size() - 1];
- }
- else if (xIndex == 0)
- {
- zEstimate = zetaTable[xIndex];
- }
- else
- {
- /* Interpolate */
- double x0 = xTabulated[xIndex - 1];
- double x1 = xTabulated[xIndex];
- double w = (xScalar - x0) / (x1 - x0);
- zEstimate = w * zetaTable[xIndex - 1] + (1 - w) * zetaTable[xIndex];
- }
-
- return zEstimate;
-}
-
/*! \brief
* Estimate a reasonable initial reference weight histogram size.
*
- * \param[in] dimParams Parameters for the dimensions of the coordinate.
* \param[in] awhBiasParams Bias parameters.
* \param[in] gridAxis The BiasGrid axes.
* \param[in] beta 1/(k_B T).
* \param[in] samplingTimestep Sampling frequency of probability weights.
* \returns estimate of initial histogram size.
*/
-double getInitialHistogramSizeEstimate(const std::vector<DimParams>& dimParams,
- const AwhBiasParams& awhBiasParams,
- const std::vector<GridAxis>& gridAxis,
- double beta,
- double samplingTimestep)
+double getInitialHistogramSizeEstimate(const AwhBiasParams& awhBiasParams,
+ const std::vector<GridAxis>& gridAxis,
+ double beta,
+ double samplingTimestep)
{
/* Get diffusion factor */
- double crossingTime = 0.;
+ double maxCrossingTime = 0.;
std::vector<double> x;
for (size_t d = 0; d < gridAxis.size(); d++)
{
- double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
- if (axisLength > 0)
- {
- crossingTime += awhBiasParams.dimParams[d].diffusion / (axisLength * axisLength);
- /* The sigma of the Gaussian distribution in the umbrella */
- double sigma = 1.;
- if (dimParams[d].isPullDimension())
- {
- GMX_RELEASE_ASSERT(dimParams[d].pullDimParams().betak != 0,
- "beta*k cannot be zero");
- sigma /= std::sqrt(dimParams[d].pullDimParams().betak);
- }
- x.push_back(sigma / axisLength);
- }
+ GMX_RELEASE_ASSERT(awhBiasParams.dimParams[d].diffusion > 0, "We need positive diffusion");
+ // With diffusion it takes on average T = L^2/2D time to cross length L
+ double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
+ double crossingTime = (axisLength * axisLength) / (2 * awhBiasParams.dimParams[d].diffusion);
+ maxCrossingTime = std::max(maxCrossingTime, crossingTime);
}
- GMX_RELEASE_ASSERT(crossingTime > 0, "We need at least one dimension with non-zero length");
+ GMX_RELEASE_ASSERT(maxCrossingTime > 0, "We need at least one dimension with non-zero length");
double errorInitialInKT = beta * awhBiasParams.errorInitial;
- double histogramSize = gaussianGeometryFactor(x)
- / (crossingTime * gmx::square(errorInitialInKT) * samplingTimestep);
+ double histogramSize = maxCrossingTime / (gmx::square(errorInitialInKT) * samplingTimestep);
return histogramSize;
}
updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
localWeightScaling(eTarget == eawhtargetLOCALBOLTZMANN ? temperatureScaleFactor : 1),
initialErrorInKT(beta * awhBiasParams.errorInitial),
- initialHistogramSize(getInitialHistogramSizeEstimate(dimParams,
- awhBiasParams,
- gridAxis,
- beta,
- numStepsSampleCoord_ * mdTimeStep)),
+ initialHistogramSize(
+ getInitialHistogramSizeEstimate(awhBiasParams, gridAxis, beta, numStepsSampleCoord_ * mdTimeStep)),
convolveForce(awhParams.ePotential == eawhpotentialCONVOLVED),
biasIndex(biasIndex),
disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff