src/gromacs/applied_forces/awh/biasparams.h

   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  *
  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 template<typename>
  64 class ArrayRef;
  65 struct AwhBiasParams;
  66 struct AwhParams;
  67 struct DimParams;
  68 class GridAxis;
  69 enum class AwhTargetType : int;
  70
  71 /*! \internal \brief Constant parameters for the bias.
  72  */
  73 class BiasParams
  74 {
  75 public:
  76     /*! \brief Switch to turn off update skips, useful for testing.
  77      */
  78     enum class DisableUpdateSkips
  79     {
  80         no, /**< Allow update skips (when supported by the method) */
  81         yes /**< Disable update skips */
  82     };
  83
  84     /*! \brief
  85      * Check if the parameters permit skipping updates.
  86      *
  87      * Generally, we can skip updates of points that are non-local
  88      * at the time of the update if we for later times, when the points
  89      * with skipped updates have become local, know exactly how to apply
  90      * the previous updates. The free energy updates only depend
  91      * on local sampling, but the histogram rescaling factors
  92      * generally depend on the histogram size (all samples).
  93      * If the histogram size is kept constant or the scaling factors
  94      * are trivial, this is not a problem. However, if the histogram growth
  95      * is scaled down by some factor the size at the time of the update
  96      * needs to be known. It would be fairly simple to, for a deterministically
  97      * growing histogram, backtrack and calculate this value, but currently
  98      * we just disallow this case. This is not a restriction because it
  99      * only affects the local Boltzmann target type for which every update
 100      * is currently anyway global because the target is always updated globally.
 101      *
 102      * \returns true when we can skip updates.
 103      */
 104     inline bool skipUpdates() const { return (!disableUpdateSkips_ && localWeightScaling == 1); }
 105
 106     /*! \brief
 107      * Returns the radius that needs to be sampled around a point before it is considered covered.
 108      */
 109     inline const awh_ivec& coverRadius() const { return coverRadius_; }
 110
 111     /*! \brief
 112      * Returns whether we should sample the coordinate.
 113      *
 114      * \param[in] step  The MD step number.
 115      */
 116     inline bool isSampleCoordStep(int64_t step) const
 117     {
 118         return (step > 0 && step % numStepsSampleCoord_ == 0);
 119     }
 120
 121     /*! \brief
 122      * Returns whether we should update the free energy.
 123      *
 124      * \param[in] step  The MD step number.
 125      */
 126     inline bool isUpdateFreeEnergyStep(int64_t step) const
 127     {
 128         int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_ * numStepsSampleCoord_;
 129         return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0);
 130     }
 131
 132     /*! \brief
 133      * Returns whether we should update the target distribution.
 134      *
 135      * \param[in] step  The MD step number.
 136      */
 137     inline bool isUpdateTargetStep(int64_t step) const { return step % numStepsUpdateTarget_ == 0; }
 138
 139     /*! \brief
 140      * Returns if to do checks for covering in the initial stage.
 141      *
 142      * To avoid overhead due to expensive checks, we do not check
 143      * at every free energy update. However, if checks are
 144      * performed too rarely the detection of coverings will be
 145      * delayed, ultimately affecting free energy convergence.
 146      *
 147      * \param[in] step                  Time step.
 148      * \returns true at steps where checks should be performed.
 149      * \note  Only returns true at free energy update steps.
 150      */
 151     bool isCheckCoveringStep(int64_t step) const
 152     {
 153         return step > 0 && (step % numStepsCheckCovering_ == 0);
 154     }
 155
 156     /*! \brief
 157      * Returns if to perform checks for anomalies in the histogram.
 158      *
 159      * To avoid overhead due to expensive checks, we do not check
 160      * at every free energy update. These checks are only used for
 161      * warning the user and can be made as infrequently as
 162      * neccessary without affecting the algorithm itself.
 163      *
 164      * \param[in] step                  Time step.
 165      * \returns true at steps where checks should be performed.
 166      * \note Only returns true at free energy update steps.
 167      * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently.
 168      */
 169     bool isCheckHistogramForAnomaliesStep(int64_t step) const { return isCheckCoveringStep(step); }
 170
 171     /*! \brief Constructor.
 172      *
 173      * The local Boltzmann target distibution is defined by
 174      * 1) Adding the sampled weights instead of the target weights to the reference weight histogram.
 175      * 2) Scaling the weights of these samples by the beta scaling factor.
 176      * 3) Setting the target distribution equal the reference weight histogram.
 177      * This requires the following special update settings:
 178      *   localWeightScaling      = targetParam
 179      *   idealWeighthistUpdate   = false
 180      * Note: these variables could in principle be set to something else also for other target distribution types.
 181      * However, localWeightScaling < 1  is in general expected to give lower efficiency and, except for local Boltzmann,
 182      * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results.
 183      *
 184      * \param[in] awhParams              AWH parameters.
 185      * \param[in] awhBiasParams          Bias parameters.
 186      * \param[in] dimParams              Bias dimension parameters.
 187      * \param[in] beta                   1/(k_B T) in units of 1/(kJ/mol), should be > 0.
 188      * \param[in] mdTimeStep             The MD time step.
 189      * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 190      * \param[in] gridAxis               The grid axes.
 191      * \param[in] disableUpdateSkips     If to disable update skips, useful for testing.
 192      * \param[in] biasIndex              Index of the bias.
 193      */
 194     BiasParams(const AwhParams&          awhParams,
 195                const AwhBiasParams&      awhBiasParams,
 196                ArrayRef<const DimParams> dimParams,
 197                double                    beta,
 198                double                    mdTimeStep,
 199                DisableUpdateSkips        disableUpdateSkips,
 200                int                       numSharingSimulations,
 201                ArrayRef<const GridAxis>  gridAxis,
 202                int                       biasIndex);
 203
 204     /* Data members */
 205     const double invBeta; /**< 1/beta = kT in kJ/mol */
 206 private:
 207     const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */
 208 public:
 209     const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */
 210 private:
 211     const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */
 212     const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */
 213 public:
 214     const AwhTargetType eTarget;              /**< Type of target distribution. */
 215     const double        freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */
 216     const double        temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */
 217     const bool          idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */
 218     const int    numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */
 219     const double updateWeight;    /**< The probability weight accumulated for each update. */
 220     const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
 221     const double initialErrorInKT;   /**< Estimated initial free energy error in kT. */
 222     const double initialHistogramSize; /**< Initial reference weight histogram size. */
 223 private:
 224     awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
 225 public:
 226     const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */
 227     const int  biasIndex; /**< Index of the bias, used as a second random seed and for priting. */
 228 private:
 229     const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */
 230 };
 231
 232 } // namespace gmx
 233
 234 #endif /* GMX_AWH_BIASPARAMS_H */