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 * 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 */ skipUpdates()101 inline bool skipUpdates() const { return (!disableUpdateSkips_ && localWeightScaling == 1); } 102 103 /*! \brief 104 * Returns the radius that needs to be sampled around a point before it is considered covered. 105 */ coverRadius()106 inline const awh_ivec& coverRadius() const { return coverRadius_; } 107 108 /*! \brief 109 * Returns whether we should sample the coordinate. 110 * 111 * \param[in] step The MD step number. 112 */ isSampleCoordStep(int64_t step)113 inline bool isSampleCoordStep(int64_t step) const 114 { 115 return (step > 0 && step % numStepsSampleCoord_ == 0); 116 } 117 118 /*! \brief 119 * Returns whether we should update the free energy. 120 * 121 * \param[in] step The MD step number. 122 */ isUpdateFreeEnergyStep(int64_t step)123 inline bool isUpdateFreeEnergyStep(int64_t step) const 124 { 125 int stepIntervalUpdateFreeEnergy = numSamplesUpdateFreeEnergy_ * numStepsSampleCoord_; 126 return (step > 0 && step % stepIntervalUpdateFreeEnergy == 0); 127 } 128 129 /*! \brief 130 * Returns whether we should update the target distribution. 131 * 132 * \param[in] step The MD step number. 133 */ isUpdateTargetStep(int64_t step)134 inline bool isUpdateTargetStep(int64_t step) const { return step % numStepsUpdateTarget_ == 0; } 135 136 /*! \brief 137 * Returns if to do checks for covering in the initial stage. 138 * 139 * To avoid overhead due to expensive checks, we do not check 140 * at every free energy update. However, if checks are 141 * performed too rarely the detection of coverings will be 142 * delayed, ultimately affecting free energy convergence. 143 * 144 * \param[in] step Time step. 145 * \returns true at steps where checks should be performed. 146 * \note Only returns true at free energy update steps. 147 */ isCheckCoveringStep(int64_t step)148 bool isCheckCoveringStep(int64_t step) const 149 { 150 return step > 0 && (step % numStepsCheckCovering_ == 0); 151 } 152 153 /*! \brief 154 * Returns if to perform checks for anomalies in the histogram. 155 * 156 * To avoid overhead due to expensive checks, we do not check 157 * at every free energy update. These checks are only used for 158 * warning the user and can be made as infrequently as 159 * neccessary without affecting the algorithm itself. 160 * 161 * \param[in] step Time step. 162 * \returns true at steps where checks should be performed. 163 * \note Only returns true at free energy update steps. 164 * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently. 165 */ isCheckHistogramForAnomaliesStep(int64_t step)166 bool isCheckHistogramForAnomaliesStep(int64_t step) const { return isCheckCoveringStep(step); } 167 168 /*! \brief Constructor. 169 * 170 * The local Boltzmann target distibution is defined by 171 * 1) Adding the sampled weights instead of the target weights to the reference weight histogram. 172 * 2) Scaling the weights of these samples by the beta scaling factor. 173 * 3) Setting the target distribution equal the reference weight histogram. 174 * This requires the following special update settings: 175 * localWeightScaling = targetParam 176 * idealWeighthistUpdate = false 177 * Note: these variables could in principle be set to something else also for other target distribution types. 178 * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann, 179 * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results. 180 * 181 * \param[in] awhParams AWH parameters. 182 * \param[in] awhBiasParams Bias parameters. 183 * \param[in] dimParams Bias dimension parameters. 184 * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0. 185 * \param[in] mdTimeStep The MD time step. 186 * \param[in] numSharingSimulations The number of simulations to share the bias across. 187 * \param[in] gridAxis The grid axes. 188 * \param[in] disableUpdateSkips If to disable update skips, useful for testing. 189 * \param[in] biasIndex Index of the bias. 190 */ 191 BiasParams(const AwhParams& awhParams, 192 const AwhBiasParams& awhBiasParams, 193 const std::vector<DimParams>& dimParams, 194 double beta, 195 double mdTimeStep, 196 DisableUpdateSkips disableUpdateSkips, 197 int numSharingSimulations, 198 const std::vector<GridAxis>& gridAxis, 199 int biasIndex); 200 201 /* Data members */ 202 const double invBeta; /**< 1/beta = kT in kJ/mol */ 203 private: 204 const int64_t numStepsSampleCoord_; /**< Number of steps per coordinate value sample. */ 205 public: 206 const int numSamplesUpdateFreeEnergy_; /**< Number of samples per free energy update. */ 207 private: 208 const int64_t numStepsUpdateTarget_; /**< Number of steps per updating the target distribution. */ 209 const int64_t numStepsCheckCovering_; /**< Number of steps per checking for covering. */ 210 public: 211 const int eTarget; /**< Type of target distribution. */ 212 const double freeEnergyCutoffInKT; /**< Free energy cut-off in kT for cut-off target distribution. */ 213 const double temperatureScaleFactor; /**< Temperature scaling factor for temperature scaled targed distributions. */ 214 const bool idealWeighthistUpdate; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */ 215 const int numSharedUpdate; /**< The number of (multi-)simulations sharing the bias update */ 216 const double updateWeight; /**< The probability weight accumulated for each update. */ 217 const double localWeightScaling; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */ 218 const double initialErrorInKT; /**< Estimated initial free energy error in kT. */ 219 const double initialHistogramSize; /**< Initial reference weight histogram size. */ 220 private: 221 awh_ivec coverRadius_; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */ 222 public: 223 const bool convolveForce; /**< True if we convolve the force, false means use MC between umbrellas. */ 224 const int biasIndex; /**< Index of the bias, used as a second random seed and for priting. */ 225 private: 226 const bool disableUpdateSkips_; /**< If true, we disallow update skips, even when the method supports it. */ 227 }; 228 229 } // namespace gmx 230 231 #endif /* GMX_AWH_BIASPARAMS_H */ 232