src/gromacs/awh/biasparams.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017, 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  *
  38  * \brief
  39  * Declares the BiasParams class.
  40  *
  41  * This class holds the parameters for the bias. Most are direct copies
  42  * of the input that the user provided. Some are a combination of user
  43  * input and properties of the simulated system.
  44  *
  45  * \author Viveca Lindahl
  46  * \author Berk Hess <hess@kth.se>
  47  * \ingroup module_awh
  48  */
  49
  50 #ifndef GMX_AWH_BIASPARAMS_H
  51 #define GMX_AWH_BIASPARAMS_H
  52
  53 #include <vector>
  54
  55 #include "gromacs/math/vectypes.h"
  56 #include "gromacs/utility/basedefinitions.h"
  57
  58 #include "dimparams.h"
  59
  60 namespace gmx
  61 {
  62
  63 struct AwhBiasParams;
  64 struct AwhParams;
  65 struct DimParams;
  66 class GridAxis;
  67
  68 /*! \internal \brief Constant parameters for the bias.
  69  */
  70 class BiasParams
  71 {
  72     public:
  73         /*! \brief Switch to turn off update skips, useful for testing.
  74          */
  75         enum class DisableUpdateSkips
  76         {
  77             no,  /**< Allow update skips (when supported by the method) */
  78             yes  /**< Disable update skips */
  79         };
  80
  81         /*! \brief
  82          * Check if the parameters permit skipping updates.
  83          *
  84          * Generally, we can skip updates of points that are non-local
  85          * at the time of the update if we for later times, when the points
  86          * with skipped updates have become local, know exactly how to apply
  87          * the previous updates. The free energy updates only depend
  88          * on local sampling, but the histogram rescaling factors
  89          * generally depend on the histogram size (all samples).
  90          * If the histogram size is kept constant or the scaling factors
  91          * are trivial, this is not a problem. However, if the histogram growth
  92          * is scaled down by some factor the size at the time of the update
  93          * needs to be known. It would be fairly simple to, for a deterministically
  94          * growing histogram, backtrack and calculate this value, but currently
  95          * we just disallow this case. This is not a restriction because it
  96          * only affects the local Boltzmann target type for which every update
  97          * is currently anyway global because the target is always updated globally.
  98          *
  99          * \returns true when we can skip updates.
 100          */
 101         inline bool skipUpdates() const
 102         {
 103             return (!disableUpdateSkips_ && localWeightScaling == 1);
 104         }
 105
 106         /*! \brief
 107          * Returns the the radius that needs to be sampled around a point before it is considered covered.
 108          */
 109         inline const awh_ivec &coverRadius() const
 110         {
 111             return coverRadius_;
 112         }
 113
 114         /*! \brief
 115          * Returns whether we should sample the coordinate.
 116          *
 117          * \param[in] step  The MD step number.
 118          */
 119         inline bool isSampleCoordStep(gmx_int64_t step) const
 120         {
 121             return (step > 0 && step % numStepsSampleCoord_ == 0);
 122         }
 123
 124         /*! \brief
 125          * Returns whether we should update the free energy.
 126          *
 127          * \param[in] step  The MD step number.
 128          */
 129         inline bool isUpdateFreeEnergyStep(gmx_int64_t step) const
 130         {
 131             int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_*numStepsSampleCoord_;
 132             return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0);
 133         }
 134
 135         /*! \brief
 136          * Returns whether we should update the target distribution.
 137          *
 138          * \param[in] step  The MD step number.
 139          */
 140         inline bool isUpdateTargetStep(gmx_int64_t step) const
 141         {
 142             return step % numStepsUpdateTarget_ == 0;
 143         }
 144
 145         /*! \brief
 146          * Returns if to do checks, only returns true at free-energy update steps.
 147          *
 148          * To avoid overhead due to expensive checks, we only do checks when we
 149          * have taken at least as many samples as we have points.
 150          *
 151          * \param[in] numPointsInHistogram  The total number of points in the bias histogram.
 152          * \param[in] step                  Time step.
 153          * \returns true at steps where checks should be performed.
 154          */
 155         bool isCheckStep(std::size_t numPointsInHistogram,
 156                          gmx_int64_t step) const;
 157
 158         /*! \brief Constructor.
 159          *
 160          * The local Boltzmann target distibution is defined by
 161          * 1) Adding the sampled weights instead of the target weights to the reference weight histogram.
 162          * 2) Scaling the weights of these samples by the beta scaling factor.
 163          * 3) Setting the target distribution equal the reference weight histogram.
 164          * This requires the following special update settings:
 165          *   localWeightScaling      = targetParam
 166          *   idealWeighthistUpdate   = false
 167          * Note: these variables could in principle be set to something else also for other target distribution types.
 168          * However, localWeightScaling < 1  is in general expected to give lower efficiency and, except for local Boltzmann,
 169          * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results.
 170          *
 171          * \param[in] awhParams              AWH parameters.
 172          * \param[in] awhBiasParams          Bias parameters.
 173          * \param[in] dimParams              Bias dimension parameters.
 174          * \param[in] beta                   1/(k_B T) in units of 1/(kJ/mol), should be > 0.
 175          * \param[in] mdTimeStep             The MD time step.
 176          * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 177          * \param[in] gridAxis               The grid axes.
 178          * \param[in] disableUpdateSkips     If to disable update skips, useful for testing.
 179          * \param[in] biasIndex              Index of the bias.
 180          */
 181         BiasParams(const AwhParams              &awhParams,
 182                    const AwhBiasParams          &awhBiasParams,
 183                    const std::vector<DimParams> &dimParams,
 184                    double                        beta,
 185                    double                        mdTimeStep,
 186                    DisableUpdateSkips            disableUpdateSkips,
 187                    int                           numSharingSimulations,
 188                    const std::vector<GridAxis>  &gridAxis,
 189                    int                           biasIndex);
 190
 191         /* Data members */
 192         const double      invBeta;                     /**< 1/beta = kT in kJ/mol */
 193     private:
 194         const gmx_int64_t numStepsSampleCoord_;        /**< Number of steps per coordinate value sample. */
 195         const int         numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */
 196         const gmx_int64_t numStepsUpdateTarget_;       /**< Number of steps per updating the target distribution. */
 197     public:
 198         const int         eTarget;                     /**< Type of target distribution. */
 199         const double      freeEnergyCutoffInKT;        /**< Free energy cut-off in kT for cut-off target distribution. */
 200         const double      temperatureScaleFactor;      /**< Temperature scaling factor for temperature scaled targed distributions. */
 201         const bool        idealWeighthistUpdate;       /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */
 202         const int         numSharedUpdate;             /**< The number of (multi-)simulations sharing the bias update */
 203         const double      updateWeight;                /**< The probability weight accumulated for each update. */
 204         const double      localWeightScaling;          /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
 205         const double      initialHistogramSize;        /**< Initial reference weight histogram size. */
 206     private:
 207         awh_ivec          coverRadius_;                /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
 208     public:
 209         const bool        convolveForce;               /**< True if we convolve the force, false means use MC between umbrellas. */
 210         const int         biasIndex;                   /**< Index of the bias, used as a second random seed and for priting. */
 211     private:
 212         const bool        disableUpdateSkips_;         /**< If true, we disallow update skips, even when the method supports it. */
 213 };
 214
 215 }      // namespace gmx
 216
 217 #endif /* GMX_AWH_BIASPARAMS_H */