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, 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 "grid.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,
  78                                  const AwhBiasParams &awhBiasParams)
  79 {
  80     int64_t numStepsUpdateTarget = 0;
  81     /* Set the target update frequency based on the target distrbution type
  82      * (this could be made a user-option but there is most likely no big need
  83      * for tweaking this for most users).
  84      */
  85     switch (awhBiasParams.eTarget)
  86     {
  87         case eawhtargetCONSTANT:
  88             numStepsUpdateTarget = 0;
  89             break;
  90         case eawhtargetCUTOFF:
  91         case eawhtargetBOLTZMANN:
  92             /* Updating the target generally requires updating the whole grid so to keep the cost down
  93                we generally update the target less often than the free energy (unless the free energy
  94                update step is set to > 100 samples). */
  95             numStepsUpdateTarget = std::max(100 % awhParams.numSamplesUpdateFreeEnergy,
  96                                             awhParams.numSamplesUpdateFreeEnergy)*awhParams.nstSampleCoord;
  97             break;
  98         case eawhtargetLOCALBOLTZMANN:
  99             /* The target distribution is set equal to the reference histogram which is updated every free energy update.
 100                So the target has to be updated at the same time. This leads to a global update each time because it is
 101                assumed that a target distribution update can take any form. This is a bit unfortunate for a "local"
 102                target distribution. One could avoid the global update by making a local target update function (and
 103                postponing target updates for non-local points as for the free energy update). We avoid such additions
 104                for now and accept that this target type always does global updates. */
 105             numStepsUpdateTarget = awhParams.numSamplesUpdateFreeEnergy*awhParams.nstSampleCoord;
 106             break;
 107         default:
 108             GMX_RELEASE_ASSERT(false, "Unknown AWH target type");
 109             break;
 110     }
 111
 112     return numStepsUpdateTarget;
 113 }
 114
 115 /*! \brief Determines the step interval for checking for covering.
 116  *
 117  * \param[in] awhParams      AWH parameters.
 118  * \param[in] dimParams         Parameters for the dimensions of the coordinate.
 119  * \param[in] gridAxis          The Grid axes.
 120  * \returns the check interval in steps.
 121  */
 122 int64_t calcCheckCoveringInterval(const AwhParams              &awhParams,
 123                                   const std::vector<DimParams> &dimParams,
 124                                   const std::vector<GridAxis>  &gridAxis)
 125 {
 126     /* Each sample will have a width of sigma. To cover the axis a
 127        minimum number of samples of width sigma is required. */
 128     int minNumSamplesCover = 0;
 129     for (size_t d = 0; d < gridAxis.size(); d++)
 130     {
 131         GMX_RELEASE_ASSERT(dimParams[d].betak > 0, "Inverse temperature (beta) and force constant (k) should be positive.");
 132         double sigma           = 1.0/std::sqrt(dimParams[d].betak);
 133
 134         /* The additional sample is here because to cover a discretized
 135            axis of length sigma one needs two samples, one for each
 136            end point. */
 137         GMX_RELEASE_ASSERT(gridAxis[d].length()/sigma < std::numeric_limits<int>::max(), "The axis length in units of sigma should fit in an int");
 138         int    numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length()/sigma)) + 1;
 139
 140         /* The minimum number of samples needed for simultaneously
 141            covering all axes is limited by the axis requiring most
 142            samples. */
 143         minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
 144     }
 145
 146     /* Convert to number of steps using the sampling frequency. The
 147        check interval should be a multiple of the update step
 148        interval. */
 149     int       numStepsUpdate      = awhParams.numSamplesUpdateFreeEnergy*awhParams.nstSampleCoord;
 150     GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy > 0, "When checking for AWH coverings, the number of samples per AWH update need 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, "Only check covering at free energy update steps");
 155
 156     return numStepsCheck;
 157 }
 158
 159 /*! \brief
 160  * Returns an approximation of the geometry factor used for initializing the AWH update size.
 161  *
 162  * The geometry factor is defined as the following sum of Gaussians:
 163  * sum_{k!=0} exp(-0.5*(k*pi*x)^2)/(pi*k)^2,
 164  * where k is a xArray.size()-dimensional integer vector with k_i in {0,1,..}.
 165  *
 166  * \param[in] xArray  Array to evaluate.
 167  * \returns the geometry factor.
 168  */
 169 double gaussianGeometryFactor(gmx::ArrayRef<const double> xArray)
 170 {
 171     /* For convenience we give the geometry factor function a name: zeta(x) */
 172     constexpr size_t                    tableSize   = 5;
 173     std::array<const double, tableSize> xTabulated  =
 174     { {1e-5, 1e-4, 1e-3, 1e-2, 1e-1} };
 175     std::array<const double, tableSize> zetaTable1d =
 176     { { 0.166536811948, 0.16653116886, 0.166250075882, 0.162701098306, 0.129272430287 } };
 177     std::array<const double, tableSize> zetaTable2d =
 178     { { 2.31985974274, 1.86307292523, 1.38159772648, 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     {
 188         zetaTable = zetaTable2d;
 189     }
 190     else
 191     {
 192         /* TODO... but this is anyway a rough estimate and > 2 dimensions is not so popular. */
 193         zetaTable = zetaTable2d;
 194     }
 195
 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. */
 198     double xScalar = 1;
 199     for (const double &x : xArray)
 200     {
 201         xScalar *= x;
 202     }
 203
 204     GMX_ASSERT(!xArray.empty(), "We should have a non-empty input array");
 205     xScalar = std::pow(xScalar, 1.0/xArray.size());
 206
 207     /* Look up zeta(x) */
 208     size_t xIndex = 0;
 209     while ((xIndex < xTabulated.size()) && (xScalar > xTabulated[xIndex]))
 210     {
 211         xIndex++;
 212     }
 213
 214     double zEstimate;
 215     if (xIndex == xTabulated.size())
 216     {
 217         /* Take last value */
 218         zEstimate = zetaTable[xTabulated.size() - 1];
 219     }
 220     else if (xIndex == 0)
 221     {
 222         zEstimate = zetaTable[xIndex];
 223     }
 224     else
 225     {
 226         /* Interpolate */
 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];
 231     }
 232
 233     return zEstimate;
 234 }
 235
 236 /*! \brief
 237  * Estimate a reasonable initial reference weight histogram size.
 238  *
 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.
 245  */
 246 double getInitialHistogramSizeEstimate(const std::vector<DimParams> &dimParams,
 247                                        const AwhBiasParams          &awhBiasParams,
 248                                        const std::vector<GridAxis>  &gridAxis,
 249                                        double                        beta,
 250                                        double                        samplingTimestep)
 251 {
 252     /* Get diffusion factor */
 253     double              crossingTime = 0.;
 254     std::vector<double> x;
 255     for (size_t d = 0; d < gridAxis.size(); d++)
 256     {
 257         double axisLength = gridAxis[d].length();
 258         if (axisLength > 0)
 259         {
 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);
 264         }
 265     }
 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)/(crossingTime*gmx::square(errorInitialInKT)*samplingTimestep);
 269
 270     return histogramSize;
 271 }
 272
 273 /*! \brief Returns the number of simulations sharing bias updates.
 274  *
 275  * \param[in] awhBiasParams          Bias parameters.
 276  * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 277  * \returns the number of shared updates.
 278  */
 279 int getNumSharedUpdate(const AwhBiasParams &awhBiasParams,
 280                        int                  numSharingSimulations)
 281 {
 282     GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
 283
 284     int numShared = 1;
 285
 286     if (awhBiasParams.shareGroup > 0)
 287     {
 288         /* We do not yet support sharing within a simulation */
 289         int numSharedWithinThisSimulation = 1;
 290         numShared                         = numSharingSimulations*numSharedWithinThisSimulation;
 291     }
 292
 293     return numShared;
 294 }
 295
 296 }   // namespace
 297
 298 BiasParams::BiasParams(const AwhParams              &awhParams,
 299                        const AwhBiasParams          &awhBiasParams,
 300                        const std::vector<DimParams> &dimParams,
 301                        double                        beta,
 302                        double                        mdTimeStep,
 303                        DisableUpdateSkips            disableUpdateSkips,
 304                        int                           numSharingSimulations,
 305                        const std::vector<GridAxis>  &gridAxis,
 306                        int                           biasIndex) :
 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, awhBiasParams, gridAxis, beta, numStepsSampleCoord_*mdTimeStep)),
 321     convolveForce(awhParams.ePotential == eawhpotentialCONVOLVED),
 322     biasIndex(biasIndex),
 323     disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
 324 {
 325     if (beta <= 0)
 326     {
 327         GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
 328     }
 329
 330     for (int d = 0; d < awhBiasParams.ndim; d++)
 331     {
 332         double coverRadiusInNm = 0.5*awhBiasParams.dimParams[d].coverDiameter;
 333         double spacing         = gridAxis[d].spacing();
 334         coverRadius_[d]        = spacing > 0 ?  static_cast<int>(std::round(coverRadiusInNm/spacing)) : 0;
 335     }
 336 }
 337
 338 } // namespace gmx