2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
39 * Declares the BiasParams class.
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.
45 * \author Viveca Lindahl
46 * \author Berk Hess <hess@kth.se>
50 #ifndef GMX_AWH_BIASPARAMS_H
51 #define GMX_AWH_BIASPARAMS_H
55 #include "gromacs/math/vectypes.h"
56 #include "gromacs/utility/basedefinitions.h"
58 #include "dimparams.h"
68 /*! \internal \brief Constant parameters for the bias.
73 /*! \brief Switch to turn off update skips, useful for testing.
75 enum class DisableUpdateSkips
77 no
, /**< Allow update skips (when supported by the method) */
78 yes
/**< Disable update skips */
82 * Check if the parameters permit skipping updates.
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.
99 * \returns true when we can skip updates.
101 inline bool skipUpdates() const
103 return (!disableUpdateSkips_
&& localWeightScaling
== 1);
107 * Returns the the radius that needs to be sampled around a point before it is considered covered.
109 inline const awh_ivec
&coverRadius() const
115 * Returns whether we should sample the coordinate.
117 * \param[in] step The MD step number.
119 inline bool isSampleCoordStep(int64_t step
) const
121 return (step
> 0 && step
% numStepsSampleCoord_
== 0);
125 * Returns whether we should update the free energy.
127 * \param[in] step The MD step number.
129 inline bool isUpdateFreeEnergyStep(int64_t step
) const
131 int stepIntervalUpdateFreeEnergy
= numSamplesUpdateFreeEnergy_
*numStepsSampleCoord_
;
132 return (step
> 0 && step
% stepIntervalUpdateFreeEnergy
== 0);
136 * Returns whether we should update the target distribution.
138 * \param[in] step The MD step number.
140 inline bool isUpdateTargetStep(int64_t step
) const
142 return step
% numStepsUpdateTarget_
== 0;
146 * Returns if to do checks for covering in the initial stage.
148 * To avoid overhead due to expensive checks, we do not check
149 * at every free energy update. However, if checks are
150 * performed too rarely the detection of coverings will be
151 * delayed, ultimately affecting free energy convergence.
153 * \param[in] step Time step.
154 * \returns true at steps where checks should be performed.
155 * \note Only returns true at free energy update steps.
157 bool isCheckCoveringStep(int64_t step
) const
159 return step
% numStepsCheckCovering_
== 0;
163 * Returns if to perform checks for anomalies in the histogram.
165 * To avoid overhead due to expensive checks, we do not check
166 * at every free energy update. These checks are only used for
167 * warning the user and can be made as infrequently as
168 * neccessary without affecting the algorithm itself.
170 * \param[in] step Time step.
171 * \returns true at steps where checks should be performed.
172 * \note Only returns true at free energy update steps.
173 * \todo Currently this function just calls isCheckCoveringStep but the checks could be done less frequently.
175 bool isCheckHistogramForAnomaliesStep(int64_t step
) const
177 return isCheckCoveringStep(step
);
180 /*! \brief Constructor.
182 * The local Boltzmann target distibution is defined by
183 * 1) Adding the sampled weights instead of the target weights to the reference weight histogram.
184 * 2) Scaling the weights of these samples by the beta scaling factor.
185 * 3) Setting the target distribution equal the reference weight histogram.
186 * This requires the following special update settings:
187 * localWeightScaling = targetParam
188 * idealWeighthistUpdate = false
189 * Note: these variables could in principle be set to something else also for other target distribution types.
190 * However, localWeightScaling < 1 is in general expected to give lower efficiency and, except for local Boltzmann,
191 * idealWeightHistUpdate = false gives (in my experience) unstable, non-converging results.
193 * \param[in] awhParams AWH parameters.
194 * \param[in] awhBiasParams Bias parameters.
195 * \param[in] dimParams Bias dimension parameters.
196 * \param[in] beta 1/(k_B T) in units of 1/(kJ/mol), should be > 0.
197 * \param[in] mdTimeStep The MD time step.
198 * \param[in] numSharingSimulations The number of simulations to share the bias across.
199 * \param[in] gridAxis The grid axes.
200 * \param[in] disableUpdateSkips If to disable update skips, useful for testing.
201 * \param[in] biasIndex Index of the bias.
203 BiasParams(const AwhParams
&awhParams
,
204 const AwhBiasParams
&awhBiasParams
,
205 const std::vector
<DimParams
> &dimParams
,
208 DisableUpdateSkips disableUpdateSkips
,
209 int numSharingSimulations
,
210 const std::vector
<GridAxis
> &gridAxis
,
214 const double invBeta
; /**< 1/beta = kT in kJ/mol */
216 const int64_t numStepsSampleCoord_
; /**< Number of steps per coordinate value sample. */
218 const int numSamplesUpdateFreeEnergy_
; /**< Number of samples per free energy update. */
220 const int64_t numStepsUpdateTarget_
; /**< Number of steps per updating the target distribution. */
221 const int64_t numStepsCheckCovering_
; /**< Number of steps per checking for covering. */
223 const int eTarget
; /**< Type of target distribution. */
224 const double freeEnergyCutoffInKT
; /**< Free energy cut-off in kT for cut-off target distribution. */
225 const double temperatureScaleFactor
; /**< Temperature scaling factor for temperature scaled targed distributions. */
226 const bool idealWeighthistUpdate
; /**< Update reference weighthistogram using the target distribution? Otherwise use the realized distribution. */
227 const int numSharedUpdate
; /**< The number of (multi-)simulations sharing the bias update */
228 const double updateWeight
; /**< The probability weight accumulated for each update. */
229 const double localWeightScaling
; /**< Scaling factor applied to a sample before adding it to the reference weight histogram (= 1, usually). */
230 const double initialErrorInKT
; /**< Estimated initial free energy error in kT. */
231 const double initialHistogramSize
; /**< Initial reference weight histogram size. */
233 awh_ivec coverRadius_
; /**< The radius (in points) that needs to be sampled around a point before it is considered covered. */
235 const bool convolveForce
; /**< True if we convolve the force, false means use MC between umbrellas. */
236 const int biasIndex
; /**< Index of the bias, used as a second random seed and for priting. */
238 const bool disableUpdateSkips_
; /**< If true, we disallow update skips, even when the method supports it. */
243 #endif /* GMX_AWH_BIASPARAMS_H */