src/gromacs/applied_forces/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,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.
   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 #include <limits>
  53
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/mdtypes/awh_params.h"
  56 #include "gromacs/utility/arrayref.h"
  57 #include "gromacs/utility/exceptions.h"
  58 #include "gromacs/utility/gmxassert.h"
  59
  60 #include "biasgrid.h"
  61
  62 namespace gmx
  63 {
  64
  65 namespace
  66 {
  67
  68 /*! \brief Determines the interval for updating the target distribution.
  69  *
  70  * The interval value is based on the target distrbution type
  71  * (this could be made a user-option but there is most likely
  72  * no big need for tweaking this for most users).
  73  *
  74  * \param[in] awhParams      AWH parameters.
  75  * \param[in] awhBiasParams  Bias parameters.
  76  * \returns the target update interval in steps.
  77  */
  78 int64_t calcTargetUpdateInterval(const AwhParams& awhParams, 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.targetDistribution())
  86     {
  87         case AwhTargetType::Constant: numStepsUpdateTarget = 0; break;
  88         case AwhTargetType::Cutoff:
  89         case AwhTargetType::Boltzmann:
  90             /* Updating the target generally requires updating the whole grid so to keep the cost
  91                down we generally update the target less often than the free energy (unless the free
  92                energy update step is set to > 100 samples). */
  93             numStepsUpdateTarget = std::max(100 % awhParams.numSamplesUpdateFreeEnergy(),
  94                                             awhParams.numSamplesUpdateFreeEnergy())
  95                                    * awhParams.nstSampleCoord();
  96             break;
  97         case AwhTargetType::LocalBoltzmann:
  98             /* The target distribution is set equal to the reference histogram which is updated every free energy update.
  99                So the target has to be updated at the same time. This leads to a global update each time because it is
 100                assumed that a target distribution update can take any form. This is a bit unfortunate for a "local"
 101                target distribution. One could avoid the global update by making a local target update function (and
 102                postponing target updates for non-local points as for the free energy update). We avoid such additions
 103                for now and accept that this target type always does global updates. */
 104             numStepsUpdateTarget = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
 105             break;
 106         default: GMX_RELEASE_ASSERT(false, "Unknown AWH target type"); break;
 107     }
 108
 109     return numStepsUpdateTarget;
 110 }
 111
 112 /*! \brief Determines the step interval for checking for covering.
 113  *
 114  * \param[in] awhParams      AWH parameters.
 115  * \param[in] dimParams         Parameters for the dimensions of the coordinate.
 116  * \param[in] gridAxis          The BiasGrid axes.
 117  * \returns the check interval in steps.
 118  */
 119 int64_t calcCheckCoveringInterval(const AwhParams&          awhParams,
 120                                   ArrayRef<const DimParams> dimParams,
 121                                   ArrayRef<const GridAxis>  gridAxis)
 122 {
 123     /* Each sample will have a width of sigma. To cover the axis a
 124        minimum number of samples of width sigma is required. */
 125     int minNumSamplesCover = 0;
 126     for (size_t d = 0; d < gridAxis.size(); d++)
 127     {
 128         int numSamplesCover;
 129         if (dimParams[d].isPullDimension())
 130         {
 131             GMX_RELEASE_ASSERT(
 132                     dimParams[d].pullDimParams().betak > 0,
 133                     "Inverse temperature (beta) and force constant (k) should be positive.");
 134             double sigma = 1.0 / std::sqrt(dimParams[d].pullDimParams().betak);
 135
 136             /* The additional sample is here because to cover a discretized
 137             axis of length sigma one needs two samples, one for each
 138             end point. */
 139             GMX_RELEASE_ASSERT(gridAxis[d].length() / sigma < std::numeric_limits<int>::max(),
 140                                "The axis length in units of sigma should fit in an int");
 141             numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length() / sigma)) + 1;
 142         }
 143         else
 144         {
 145             numSamplesCover = dimParams[d].fepDimParams().numFepLambdaStates;
 146         }
 147         /* The minimum number of samples needed for simultaneously
 148         covering all axes is limited by the axis requiring most
 149         samples. */
 150         minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
 151     }
 152
 153     /* Convert to number of steps using the sampling frequency. The
 154        check interval should be a multiple of the update step
 155        interval. */
 156     int numStepsUpdate = awhParams.numSamplesUpdateFreeEnergy() * awhParams.nstSampleCoord();
 157     GMX_RELEASE_ASSERT(awhParams.numSamplesUpdateFreeEnergy() > 0,
 158                        "When checking for AWH coverings, the number of samples per AWH update need "
 159                        "to be > 0.");
 160     int numUpdatesCheck = std::max(1, minNumSamplesCover / awhParams.numSamplesUpdateFreeEnergy());
 161     int numStepsCheck   = numUpdatesCheck * numStepsUpdate;
 162
 163     GMX_RELEASE_ASSERT(numStepsCheck % numStepsUpdate == 0,
 164                        "Only check covering at free energy update steps");
 165
 166     return numStepsCheck;
 167 }
 168
 169 /*! \brief
 170  * Estimate a reasonable initial reference weight histogram size.
 171  *
 172  * \param[in] awhBiasParams     Bias parameters.
 173  * \param[in] gridAxis          The BiasGrid axes.
 174  * \param[in] beta              1/(k_B T).
 175  * \param[in] samplingTimestep  Sampling frequency of probability weights.
 176  * \returns estimate of initial histogram size.
 177  */
 178 double getInitialHistogramSizeEstimate(const AwhBiasParams&     awhBiasParams,
 179                                        ArrayRef<const GridAxis> gridAxis,
 180                                        double                   beta,
 181                                        double                   samplingTimestep)
 182 {
 183     /* Get diffusion factor */
 184     double              maxCrossingTime = 0.;
 185     std::vector<double> x;
 186     const auto          awhDimParams = awhBiasParams.dimParams();
 187     for (size_t d = 0; d < gridAxis.size(); d++)
 188     {
 189         GMX_RELEASE_ASSERT(awhDimParams[d].diffusion() > 0, "We need positive diffusion");
 190         // With diffusion it takes on average T = L^2/2D time to cross length L
 191         double axisLength   = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
 192         double crossingTime = (axisLength * axisLength) / (2 * awhDimParams[d].diffusion());
 193         maxCrossingTime     = std::max(maxCrossingTime, crossingTime);
 194     }
 195     GMX_RELEASE_ASSERT(maxCrossingTime > 0, "We need at least one dimension with non-zero length");
 196     double errorInitialInKT = beta * awhBiasParams.initialErrorEstimate();
 197     double histogramSize    = maxCrossingTime / (gmx::square(errorInitialInKT) * samplingTimestep);
 198
 199     return histogramSize;
 200 }
 201
 202 /*! \brief Returns the number of simulations sharing bias updates.
 203  *
 204  * \param[in] awhBiasParams          Bias parameters.
 205  * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 206  * \returns the number of shared updates.
 207  */
 208 int getNumSharedUpdate(const AwhBiasParams& awhBiasParams, int numSharingSimulations)
 209 {
 210     GMX_RELEASE_ASSERT(numSharingSimulations >= 1, "We should ''share'' at least with ourselves");
 211
 212     int numShared = 1;
 213
 214     if (awhBiasParams.shareGroup() > 0)
 215     {
 216         /* We do not yet support sharing within a simulation */
 217         int numSharedWithinThisSimulation = 1;
 218         numShared                         = numSharingSimulations * numSharedWithinThisSimulation;
 219     }
 220
 221     return numShared;
 222 }
 223
 224 } // namespace
 225
 226 BiasParams::BiasParams(const AwhParams&          awhParams,
 227                        const AwhBiasParams&      awhBiasParams,
 228                        ArrayRef<const DimParams> dimParams,
 229                        double                    beta,
 230                        double                    mdTimeStep,
 231                        DisableUpdateSkips        disableUpdateSkips,
 232                        int                       numSharingSimulations,
 233                        ArrayRef<const GridAxis>  gridAxis,
 234                        int                       biasIndex) :
 235     invBeta(beta > 0 ? 1 / beta : 0),
 236     numStepsSampleCoord_(awhParams.nstSampleCoord()),
 237     numSamplesUpdateFreeEnergy_(awhParams.numSamplesUpdateFreeEnergy()),
 238     numStepsUpdateTarget_(calcTargetUpdateInterval(awhParams, awhBiasParams)),
 239     numStepsCheckCovering_(calcCheckCoveringInterval(awhParams, dimParams, gridAxis)),
 240     eTarget(awhBiasParams.targetDistribution()),
 241     freeEnergyCutoffInKT(beta * awhBiasParams.targetCutoff()),
 242     temperatureScaleFactor(awhBiasParams.targetBetaScaling()),
 243     idealWeighthistUpdate(eTarget != AwhTargetType::LocalBoltzmann),
 244     numSharedUpdate(getNumSharedUpdate(awhBiasParams, numSharingSimulations)),
 245     updateWeight(numSamplesUpdateFreeEnergy_ * numSharedUpdate),
 246     localWeightScaling(eTarget == AwhTargetType::LocalBoltzmann ? temperatureScaleFactor : 1),
 247     initialErrorInKT(beta * awhBiasParams.initialErrorEstimate()),
 248     initialHistogramSize(
 249             getInitialHistogramSizeEstimate(awhBiasParams, gridAxis, beta, numStepsSampleCoord_ * mdTimeStep)),
 250     convolveForce(awhParams.potential() == AwhPotentialType::Convolved),
 251     biasIndex(biasIndex),
 252     disableUpdateSkips_(disableUpdateSkips == DisableUpdateSkips::yes)
 253 {
 254     if (beta <= 0)
 255     {
 256         GMX_THROW(InvalidInputError("To use AWH, the beta=1/(k_B T) should be > 0"));
 257     }
 258
 259     const auto& awhDimParams = awhBiasParams.dimParams();
 260     for (int d = 0; d < gmx::ssize(awhDimParams); d++)
 261     {
 262         /* The spacing in FEP dimensions is 1. The calculated coverRadius will be in lambda states
 263          * (cf points in other dimensions). */
 264         double coverRadiusInNm = 0.5 * awhDimParams[d].coverDiameter();
 265         double spacing         = gridAxis[d].spacing();
 266         coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;
 267     }
 268 }
 269
 270 } // namespace gmx