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"
68 struct gmx_multisim_t;
75 struct AwhBiasHistory;
79 struct AwhPointStateHistory;
80 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] numSharingSimulations The number of simulations to share the bias across.
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 int numSharingSimulations,
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.
210 * \param[in] commRecord Struct for intra-simulation communication.
211 * \param[in] ms Struct for multi-simulation communication.
213 * \param[in] step Time step.
214 * \param[in] seed Random seed.
215 * \param[in,out] fplog Log file.
216 * \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.
218 gmx::ArrayRef<const double> calcForceAndUpdateBias(const awh_dvec coordValue,
219 ArrayRef<const double> neighborLambdaEnergies,
220 ArrayRef<const double> neighborLambdaDhdl,
221 double* awhPotential,
222 double* potentialJump,
223 const t_commrec* commRecord,
224 const gmx_multisim_t* ms,
231 * Calculates the convolved bias for a given coordinate value.
233 * The convolved bias is the effective bias acting on the coordinate.
234 * Since the bias here has arbitrary normalization, this only makes
235 * sense as a relative, to other coordinate values, measure of the bias.
237 * \param[in] coordValue The coordinate value.
238 * \returns the convolved bias >= -GMX_FLOAT_MAX.
240 double calcConvolvedBias(const awh_dvec& coordValue) const
242 return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
246 * Restore the bias state from history on the master rank and broadcast it.
248 * \param[in] biasHistory Bias history struct, only allowed to be nullptr on non-master ranks.
249 * \param[in] cr The communication record.
251 void restoreStateFromHistory(const AwhBiasHistory* biasHistory, const t_commrec* cr);
254 * Allocate and initialize a bias history with the given bias state.
256 * This function will be called at the start of a new simulation.
257 * Note that only constant data will be initialized here.
258 * History data is set by \ref updateHistory.
260 * \param[in,out] biasHistory AWH history to initialize.
262 void initHistoryFromState(AwhBiasHistory* biasHistory) const;
265 * Update the bias history with the current state.
267 * \param[out] biasHistory Bias history struct.
269 void updateHistory(AwhBiasHistory* biasHistory) const;
272 * Do all previously skipped updates.
273 * Public for use by tests.
275 void doSkippedUpdatesForAllPoints();
277 //! Returns the dimensionality of the bias.
278 inline int ndim() const { return dimParams_.size(); }
280 /*! \brief Returns the dimension parameters.
282 inline ArrayRef<const DimParams> dimParams() const { return dimParams_; }
284 //! Returns the bias parameters
285 inline const BiasParams& params() const { return params_; }
287 //! Returns the global state of the bias.
288 inline const BiasState& state() const { return state_; }
290 //! Returns the index of the bias.
291 inline int biasIndex() const { return params_.biasIndex; }
293 /*! \brief Return the coordinate value for a grid point.
295 * \param[in] gridPointIndex The index of the grid point.
297 inline const awh_dvec& getGridCoordValue(size_t gridPointIndex) const
299 GMX_ASSERT(gridPointIndex < grid_.numPoints(),
300 "gridPointIndex should be in the range of the grid");
302 return grid_.point(gridPointIndex).coordValue;
307 * Performs statistical checks on the collected histograms and warns if issues are detected.
310 * \param[in] step Time step.
311 * \param[in,out] fplog Output file for warnings.
313 void warnForHistogramAnomalies(double t, int64_t step, FILE* fplog);
316 * Collect samples for the force correlation analysis on the grid.
318 * \param[in] probWeightNeighbor Probability weight of the neighboring points.
319 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
320 * points. The array is of length numLambdas+1, where numLambdas is the number of free
321 * energy lambda states. Element 0 in the array is the dHdL
322 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
323 * neighboring lambda states (also including the current state). When there are no free
324 * energy lambda state dimensions this can be empty.
325 * \param[in] t The time.
327 void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
328 ArrayRef<const double> neighborLambdaDhdl,
332 /*! \brief Return a const reference to the force correlation grid.
334 const CorrelationGrid& forceCorrelationGrid() const
336 GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr,
337 "forceCorrelationGrid() should only be called with a valid force "
338 "correlation object");
340 return *forceCorrelationGrid_;
343 /*! \brief Return the number of data blocks that have been prepared for writing.
345 int numEnergySubblocksToWrite() const;
347 /*! \brief Write bias data blocks to energy subblocks.
349 * \param[in,out] subblock Energy subblocks to write to.
350 * \returns the number of subblocks written.
352 int writeToEnergySubblocks(t_enxsubblock* subblock) const;
354 /*! \brief Returns true if the bias has a free energy lambda state dimension at all.
356 bool hasFepLambdaDimension() const
358 return std::any_of(std::begin(dimParams_), std::end(dimParams_), [](const auto& dimParam) {
359 return dimParam.isFepLambdaDimension();
364 * Returns whether we should sample the coordinate.
366 * \param[in] step The MD step number.
368 bool isSampleCoordStep(int64_t step) const;
372 const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
373 const BiasGrid grid_; /**< The multidimensional grid organizing the coordinate point locations. */
375 const BiasParams params_; /**< Constant parameters for the method. */
377 BiasState state_; /**< The state, both global and of the grid points */
378 std::vector<int> updateList_; /**< List of points for update for temporary use (could be made another tempWorkSpace) */
380 const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
382 std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
384 /* Force correlation grid */
385 std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
388 std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
390 /* Temporary working vectors used during the update.
391 * These are only here to avoid allocation at every MD step.
393 std::vector<double, AlignedAllocator<double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
394 std::vector<double> tempForce_; /**< Bias force work buffer. */
396 /* Run-local counter to avoid flooding log with warnings. */
397 int numWarningsIssued_; /**< The number of warning issued in the current run. */
402 #endif /* GMX_AWH_BIAS_H */