2 * This file is part of the GROMACS molecular simulation package.
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.
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 Bias class.
41 * This class is essentially a wrapper around the BiasState class.
42 * In addition to BiasState, it holds all data that BiasState needs
43 * to update the bias. Interaction of the outside world, such as updating
44 * BiasState or extracting bias data all happen through Bias.
46 * \author Viveca Lindahl
47 * \author Berk Hess <hess@kth.se>
51 #ifndef GMX_AWH_BIAS_H
52 #define GMX_AWH_BIAS_H
57 #include "gromacs/math/vectypes.h"
58 #include "gromacs/utility/alignedallocator.h"
59 #include "gromacs/utility/basedefinitions.h"
60 #include "gromacs/utility/gmxassert.h"
63 #include "biasparams.h"
64 #include "biasstate.h"
65 #include "biaswriter.h"
66 #include "dimparams.h"
74 struct AwhBiasHistory;
78 struct AwhPointStateHistory;
79 class CorrelationGrid;
86 * \brief A bias acting on a multidimensional coordinate.
88 * At each step AWH should provide its biases with updated
89 * values of their coordinates. Each bias provides AWH with an updated
90 * bias forces and the corresponding potential.
92 * See the user manual for details on the algorithm and equations.
94 * The bias is responsible for keeping and updating a free energy estimate
95 * along the coordinate. The bias potential is basically a function of the
96 * free energy estimate and so also changes by the update.
97 * The free energy update is based on information from coordinate samples
98 * collected at a constant bias potential, between updates.
100 * The bias keeps a grid with coordinate points that organizes spatial
101 * information about the coordinate. The grid has the same geometry
102 * as the coordinate, i.e. they have the same dimensionality and periodicity
103 * (if any). The number of points in the grid sets the resolution of
104 * the collected data and its extent defines the sampling region of interest.
106 * Each coordinate point has further statistical properties and function values
107 * which a grid point does not know about. E.g., for the bias each coordinate point
108 * is associated with values of the bias, free energy and target distribution,
109 * accumulated sampling weight, etc. For this the bias attaches to each grid
110 * point a state. The grid + vector of point states are the bias coordinate points.
112 * The bias has a fairly complex global state keeping track of where
113 * the system (coordinate) currently is (CoordState), where it has
114 * sampled since the last update (BiasState) and controlling the free energy
115 * convergence rate (HistogramSize).
117 * Partly, the complexity comes from the bias having two convergence stages:
118 * an initial stage which in an heuristic, non-deterministic way restricts
119 * the early convergence rate for sake of robustness; and a final stage
120 * where the convergence rate is constant. The length of the initial stage
121 * depends on the sampling and is unknown beforehand.
123 * Another complexity comes from the fact that coordinate points,
124 * for sake of efficiency in the case of many grid points, are typically
125 * only accessed in recently sampled regions even though the free energy
126 * update is inherently global and affects all points.
127 * The bias allows points thay are non-local at the time the update
128 * was issued to postpone ("skip", as it is called in the code) the update.
129 * A non-local point is defined as a point which has not been sampled since
130 * the last update. Local points are points that have been sampled since
131 * the last update. The (current) set of local points are kept track of by
132 * the bias state and reset after every update. An update is called local
133 * if it only updates local points. Non-local points will temporarily "skip"
134 * the update until next time they are local (or when a global update
135 * is issued). For this to work, the bias keeps a global "clock"
136 * (in HistogramSize) of the number of issued updates. Each PointState
137 * also has its own local "clock" with the counting the number of updates
138 * it has pulled through. When a point updates its state it asserts
139 * that its local clock is synchronized with the global clock.
144 //! Enum for requesting Bias set up with(out) I/O on this rank.
145 enum class ThisRankWillDoIO
147 No, //!< This rank will not do I/O.
148 Yes //!< This rank will do I/O.
154 * \param[in] biasIndexInCollection Index of the bias in collection.
155 * \param[in] awhParams AWH parameters.
156 * \param[in] awhBiasParams Bias parameters.
157 * \param[in] dimParams Bias dimension parameters.
158 * \param[in] beta 1/(k_B T).
159 * \param[in] mdTimeStep The MD time step.
160 * \param[in] biasSharing Pointer to object for sharing bias over simulations, can be nullptr
161 * \param[in] biasInitFilename Name of file to read PMF and target from.
162 * \param[in] thisRankWillDoIO Tells whether this MPI rank will do I/O (checkpointing, AWH output),
163 * normally (only) the master rank does I/O.
164 * \param[in] disableUpdateSkips If to disable update skips, useful for testing.
166 Bias(int biasIndexInCollection,
167 const AwhParams& awhParams,
168 const AwhBiasParams& awhBiasParams,
169 ArrayRef<const DimParams> dimParams,
172 const BiasSharing* biasSharing,
173 const std::string& biasInitFilename,
174 ThisRankWillDoIO thisRankWillDoIO,
175 BiasParams::DisableUpdateSkips disableUpdateSkips = BiasParams::DisableUpdateSkips::no);
178 * Print information about initialization to log file.
180 * Prints information about AWH variables that are set internally
181 * but might be of interest to the user.
183 * \param[in,out] fplog Log file, can be nullptr.
185 void printInitializationToLog(FILE* fplog) const;
188 * Evolves the bias at every step.
190 * At each step the bias step needs to:
191 * - set the bias force and potential;
192 * - update the free energy and bias if needed;
193 * - reweight samples to extract the PMF.
195 * \param[in] coordValue The current coordinate value(s).
196 * \param[in] neighborLambdaEnergies An array containing the energy of the system
197 * in neighboring lambdas. The array is of length numLambdas+1, where numLambdas is
198 * the number of free energy lambda states. Element 0 in the array is the energy
199 * of the current state and elements 1..numLambdas contain the energy of the system in the
200 * neighboring lambda states (also including the current state). When there are no free
201 * energy lambda state dimensions this can be empty.
202 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
203 * points. The array is of length numLambdas+1, where numLambdas is the number of free
204 * energy lambda states. Element 0 in the array is the dHdL
205 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
206 * neighboring lambda states (also including the current state). When there are no free
207 * energy lambda state dimensions this can be empty.
208 * \param[out] awhPotential Bias potential.
209 * \param[out] potentialJump Change in bias potential for this bias.
211 * \param[in] step Time step.
212 * \param[in] seed Random seed.
213 * \param[in,out] fplog Log file.
214 * \returns a reference to the bias force, size \ref ndim(), valid until the next call of this method or destruction of Bias, whichever comes first.
216 gmx::ArrayRef<const double> calcForceAndUpdateBias(const awh_dvec coordValue,
217 ArrayRef<const double> neighborLambdaEnergies,
218 ArrayRef<const double> neighborLambdaDhdl,
219 double* awhPotential,
220 double* potentialJump,
227 * Calculates the convolved bias for a given coordinate value.
229 * The convolved bias is the effective bias acting on the coordinate.
230 * Since the bias here has arbitrary normalization, this only makes
231 * sense as a relative, to other coordinate values, measure of the bias.
233 * \param[in] coordValue The coordinate value.
234 * \returns the convolved bias >= -GMX_FLOAT_MAX.
236 double calcConvolvedBias(const awh_dvec& coordValue) const
238 return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
242 * Restore the bias state from history on the master rank and broadcast it.
244 * \param[in] biasHistory Bias history struct, only allowed to be nullptr on non-master ranks.
245 * \param[in] cr The communication record.
247 void restoreStateFromHistory(const AwhBiasHistory* biasHistory, const t_commrec* cr);
250 * Allocate and initialize a bias history with the given bias state.
252 * This function will be called at the start of a new simulation.
253 * Note that only constant data will be initialized here.
254 * History data is set by \ref updateHistory.
256 * \param[in,out] biasHistory AWH history to initialize.
258 void initHistoryFromState(AwhBiasHistory* biasHistory) const;
261 * Update the bias history with the current state.
263 * \param[out] biasHistory Bias history struct.
265 void updateHistory(AwhBiasHistory* biasHistory) const;
268 * Do all previously skipped updates.
269 * Public for use by tests.
271 void doSkippedUpdatesForAllPoints();
273 //! Returns the dimensionality of the bias.
274 inline int ndim() const { return dimParams_.size(); }
276 /*! \brief Returns the dimension parameters.
278 inline ArrayRef<const DimParams> dimParams() const { return dimParams_; }
280 //! Returns the bias parameters
281 inline const BiasParams& params() const { return params_; }
283 //! Returns the global state of the bias.
284 inline const BiasState& state() const { return state_; }
286 //! Returns the index of the bias.
287 inline int biasIndex() const { return params_.biasIndex; }
289 /*! \brief Return the coordinate value for a grid point.
291 * \param[in] gridPointIndex The index of the grid point.
293 inline const awh_dvec& getGridCoordValue(size_t gridPointIndex) const
295 GMX_ASSERT(gridPointIndex < grid_.numPoints(),
296 "gridPointIndex should be in the range of the grid");
298 return grid_.point(gridPointIndex).coordValue;
303 * Performs statistical checks on the collected histograms and warns if issues are detected.
306 * \param[in] step Time step.
307 * \param[in,out] fplog Output file for warnings.
309 void warnForHistogramAnomalies(double t, int64_t step, FILE* fplog);
312 * Collect samples for the force correlation analysis on the grid.
314 * \param[in] probWeightNeighbor Probability weight of the neighboring points.
315 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
316 * points. The array is of length numLambdas+1, where numLambdas is the number of free
317 * energy lambda states. Element 0 in the array is the dHdL
318 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
319 * neighboring lambda states (also including the current state). When there are no free
320 * energy lambda state dimensions this can be empty.
321 * \param[in] t The time.
323 void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
324 ArrayRef<const double> neighborLambdaDhdl,
328 /*! \brief Return a const reference to the force correlation grid.
330 const CorrelationGrid& forceCorrelationGrid() const
332 GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr,
333 "forceCorrelationGrid() should only be called with a valid force "
334 "correlation object");
336 return *forceCorrelationGrid_;
339 /*! \brief Return the number of data blocks that have been prepared for writing.
341 int numEnergySubblocksToWrite() const;
343 /*! \brief Write bias data blocks to energy subblocks.
345 * \param[in,out] subblock Energy subblocks to write to.
346 * \returns the number of subblocks written.
348 int writeToEnergySubblocks(t_enxsubblock* subblock) const;
350 /*! \brief Returns true if the bias has a free energy lambda state dimension at all.
352 bool hasFepLambdaDimension() const
354 return std::any_of(std::begin(dimParams_), std::end(dimParams_), [](const auto& dimParam) {
355 return dimParam.isFepLambdaDimension();
360 * Returns whether we should sample the coordinate.
362 * \param[in] step The MD step number.
364 bool isSampleCoordStep(int64_t step) const;
368 const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
369 const BiasGrid grid_; /**< The multidimensional grid organizing the coordinate point locations. */
371 const BiasParams params_; /**< Constant parameters for the method. */
373 BiasState state_; /**< The state, both global and of the grid points */
374 std::vector<int> updateList_; /**< List of points for update for temporary use (could be made another tempWorkSpace) */
376 const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
378 std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
380 /* Force correlation grid */
381 std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
384 std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
386 /* Temporary working vectors used during the update.
387 * These are only here to avoid allocation at every MD step.
389 std::vector<double, AlignedAllocator<double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
390 std::vector<double> tempForce_; /**< Bias force work buffer. */
392 /* Run-local counter to avoid flooding log with warnings. */
393 int numWarningsIssued_; /**< The number of warning issued in the current run. */
398 #endif /* GMX_AWH_BIAS_H */