src/gromacs/awh/biasparams.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019,2020, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  * \brief
  38  * Implements the initialization of the BiasParams class.
  39  *
  40  * \author Viveca Lindahl
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_awh
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "biasparams.h"
  48
  49 #include <cmath>
  50
  51 #include <algorithm>
  52
  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"
  58
  59 #include "biasgrid.h"
  60
  61 namespace gmx
  62 {
  63
  64 namespace
  65 {
  66
  67 /*! \brief Determines the interval for updating the target distribution.
  68  *
  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).
  72  *
  73  * \param[in] awhParams      AWH parameters.
  74  * \param[in] awhBiasParams  Bias parameters.
  75  * \returns the target update interval in steps.
  76  */
  77 int64_t calcTargetUpdateInterval(const AwhParams& awhParams, const AwhBiasParams& awhBiasParams)
  78 {
  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).
  83      */
  84     switch (awhBiasParams.eTarget)
  85     {
  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;
  95             break;
  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;
 104             break;
 105         default: GMX_RELEASE_ASSERT(false, "Unknown AWH target type"); break;
 106     }
 107
 108     return numStepsUpdateTarget;
 109 }
 110
 111 /*! \brief Determines the step interval for checking for covering.
 112  *
 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.
 117  */
 118 int64_t calcCheckCoveringInterval(const AwhParams&              awhParams,
 119                                   const std::vector<DimParams>& dimParams,
 120                                   const std::vector<GridAxis>&  gridAxis)
 121 {
 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++)
 126     {
 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);
 130
 131         /* The additional sample is here because to cover a discretized
 132            axis of length sigma one needs two samples, one for each
 133            end point. */
 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;
 137
 138         /* The minimum number of samples needed for simultaneously
 139            covering all axes is limited by the axis requiring most
 140            samples. */
 141         minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
 142     }
 143
 144     /* Convert to number of steps using the sampling frequency. The
 145        check interval should be a multiple of the update step
 146        interval. */
 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 "
 150                        "to be > 0.");
 151     int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy);
 152     int numStepsCheck   = numUpdatesCheck * numStepsUpdate;
 153
 154     GMX_RELEASE_ASSERT(numStepsCheck % numStepsUpdate == 0,
 155                        "Only check covering at free energy update steps");
 156
 157     return numStepsCheck;
 158 }
 159
 160 /*! \brief
 161  * Returns an approximation of the geometry factor used for initializing the AWH update size.
 162  *
 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,..}.
 166  *
 167  * \param[in] xArray  Array to evaluate.
 168  * \returns the geometry factor.
 169  */
 170 double gaussianGeometryFactor(gmx::ArrayRef<const double> xArray)
 171 {
 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 } };
 179
 180     gmx::ArrayRef<const double> zetaTable;
 181
 182     if (xArray.size() == 1)
 183     {
 184         zetaTable = zetaTable1d;
 185     }
 186     else if (xArray.size() == 2)
 187     { // NOLINT bugprone-branch-clone
 188         zetaTable = zetaTable2d;
 189     }
 190     else
 191     {
 192         /* TODO... but this is anyway a rough estimate and > 2 dimensions is not so popular.
 193          * Remove the above NOLINT when addressing this */
 194         zetaTable = zetaTable2d;
 195     }
 196
 197     /* TODO. Really zeta is a function of an ndim-dimensional vector x and we shoudl have a ndim-dimensional lookup-table.
 198        Here we take the geometric average of the components of x which is ok if the x-components are not very different. */
 199     double xScalar = 1;
 200     for (const double& x : xArray)
 201     {
 202         xScalar *= x;
 203     }
 204
 205     GMX_ASSERT(!xArray.empty(), "We should have a non-empty input array");
 206     xScalar = std::pow(xScalar, 1.0 / xArray.size());
 207
 208     /* Look up zeta(x) */
 209     size_t xIndex = 0;
 210     while ((xIndex < xTabulated.size()) && (xScalar > xTabulated[xIndex]))
 211     {
 212         xIndex++;
 213     }
 214
 215     double zEstimate;
 216     if (xIndex == xTabulated.size())
 217     {
 218         /* Take last value */
 219         zEstimate = zetaTable[xTabulated.size() - 1];
 220     }
 221     else if (xIndex == 0)
 222     {
 223         zEstimate = zetaTable[xIndex];
 224     }
 225     else
 226     {
 227         /* Interpolate */
 228         double x0 = xTabulated[xIndex - 1];
 229         double x1 = xTabulated[xIndex];
 230         double w  = (xScalar - x0) / (x1 - x0);
 231         zEstimate = w * zetaTable[xIndex - 1] + (1 - w) * zetaTable[xIndex];
 232     }
 233
 234     return zEstimate;
 235 }
 236
 237 /*! \brief
 238  * Estimate a reasonable initial reference weight histogram size.
 239  *
 240  * \param[in] dimParams         Parameters for the dimensions of the coordinate.
 241  * \param[in] awhBiasParams     Bias parameters.
 242  * \param[in] gridAxis          The BiasGrid axes.
 243  * \param[in] beta              1/(k_B T).
 244  * \param[in] samplingTimestep  Sampling frequency of probability weights.
 245  * \returns estimate of initial histogram size.
 246  */
 247 double getInitialHistogramSizeEstimate(const std::vector<DimParams>& dimParams,
 248                                        const AwhBiasParams&          awhBiasParams,
 249                                        const std::vector<GridAxis>&  gridAxis,
 250                                        double                        beta,
 251                                        double                        samplingTimestep)
 252 {
 253     /* Get diffusion factor */
 254     double              crossingTime = 0.;
 255     std::vector<double> x;
 256     for (size_t d = 0; d < gridAxis.size(); d++)
 257     {
 258         double axisLength = gridAxis[d].length();
 259         if (axisLength > 0)
 260         {
 261             crossingTime += awhBiasParams.dimParams[d].diffusion / (axisLength * axisLength);
 262             /* The sigma of the Gaussian distribution in the umbrella */
 263             double sigma = 1. / std::sqrt(dimParams[d].betak);
 264             x.push_back(sigma / axisLength);
 265         }
 266     }
 267     GMX_RELEASE_ASSERT(crossingTime > 0, "We need at least one dimension with non-zero length");
 268     double errorInitialInKT = beta * awhBiasParams.errorInitial;
 269     double histogramSize    = gaussianGeometryFactor(x)
 270                            / (crossingTime * gmx::square(errorInitialInKT) * samplingTimestep);
 271
 272     return histogramSize;
 273 }
 274
 275 /*! \brief Returns the number of simulations sharing bias updates.
 276  *
 277  * \param[in] awhBiasParams          Bias parameters.
 278  * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 279  * \returns the number of shared updates.
 280  */
 281 int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulations)
 282 {
 283     GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
 284
 285     int numShared = 1;
 286
 287     if (awhBiasParams.shareGroup > 0)
 288     {
 289         /* We do not yet support sharing within a simulation */
 290         int numSharedWithinThisSimulation = 1;
 291         numShared                         = numSharingSimulations * numSharedWithinThisSimulation;
 292     }
 293
 294     return numShared;
 295 }
 296
 297 } // namespace
 298
 299 BiasParams::BiasParams(const AwhParams&              awhParams,
 300                        const AwhBiasParams&          awhBiasParams,
 301                        const std::vector<DimParams>& dimParams,
 302                        double                        beta,
 303                        double                        mdTimeStep,
 304                        DisableUpdateSkips            disableUpdateSkips,
 305                        int                           numSharingSimulations,
 306                        const std::vector<GridAxis>&  gridAxis,
 307                        int                           biasIndex) :
 308     invBeta(beta > 0 ? 1 / beta : 0),
 309     numStepsSampleCoord_(awhParams.nstSampleCoord),
 310     numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy),
 311     numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
 312     numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
 313     eTarget(awhBiasParams.eTarget),
 314     freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff),
 315     temperatureScaleFactor(awhBiasParams.targetBetaScaling),
 316     idealWeighthistUpdate(eTarget != eawhtargetLOCALBOLTZMANN),
 317     numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
 318     updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
 319     localWeightScaling(eTarget == eawhtargetLOCALBOLTZMANN ? temperatureScaleFactor : 1),
 320     initialErrorInKT(beta * awhBiasParams.errorInitial),
 321     initialHistogramSize(getInitialHistogramSizeEstimate(dimParams,
 322                                                          awhBiasParams,
 323                                                          gridAxis,
 324                                                          beta,
 325                                                          numStepsSampleCoord_ * mdTimeStep)),
 326     convolveForce(awhParams.ePotential == eawhpotentialCONVOLVED),
 327     biasIndex(biasIndex),
 328     disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
 329 {
 330     if (beta <= 0)
 331     {
 332         GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
 333     }
 334
 335     for (int d = 0; d < awhBiasParams.ndim; d++)
 336     {
 337         double coverRadiusInNm = 0.5 * awhBiasParams.dimParams[d].coverDiameter;
 338         double spacing         = gridAxis[d].spacing();
 339         coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;
 340     }
 341 }
 342
 343 } // namespace gmx