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 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/compat/make_unique.h"
58 #include "gromacs/math/vectypes.h"
59 #include "gromacs/utility/alignedallocator.h"
60 #include "gromacs/utility/basedefinitions.h"
61 #include "gromacs/utility/gmxassert.h"
63 #include "biasparams.h"
64 #include "biasstate.h"
65 #include "biaswriter.h"
66 #include "dimparams.h"
69 struct gmx_multisim_t;
76 struct AwhBiasHistory;
80 struct AwhPointStateHistory;
81 class CorrelationGrid;
87 * \brief A bias acting on a multidimensional coordinate.
89 * At each step AWH should provide its biases with updated
90 * values of their coordinates. Each bias provides AWH with an updated
91 * bias forces and the corresponding potential.
93 * See the user manual for details on the algorithm and equations.
95 * The bias is responsible for keeping and updating a free energy estimate
96 * along the coordinate. The bias potential is basically a function of the
97 * free energy estimate and so also changes by the update.
98 * The free energy update is based on information from coordinate samples
99 * collected at a constant bias potential, between updates.
101 * The bias keeps a grid with coordinate points that organizes spatial
102 * information about the coordinate. The grid has the the same geometry
103 * as the coordinate, i.e. they have the same dimensionality and periodicity
104 * (if any). The number of points in the grid sets the resolution of
105 * the collected data and its extent defines the sampling region of interest.
107 * Each coordinate point has further statistical properties and function values
108 * which a grid point does not know about. E.g., for the bias each coordinate point
109 * is associated with values of the bias, free energy and target distribution,
110 * accumulated sampling weight, etc. For this the bias attaches to each grid
111 * point a state. The grid + vector of point states are the bias coordinate points.
113 * The bias has a fairly complex global state keeping track of where
114 * the system (coordinate) currently is (CoordState), where it has
115 * sampled since the last update (BiasState) and controlling the free energy
116 * convergence rate (HistogramSize).
118 * Partly, the complexity comes from the bias having two convergence stages:
119 * an initial stage which in an heuristic, non-deterministic way restricts
120 * the early convergence rate for sake of robustness; and a final stage
121 * where the convergence rate is constant. The length of the initial stage
122 * depends on the sampling and is unknown beforehand.
124 * Another complexity comes from the fact that coordinate points,
125 * for sake of efficiency in the case of many grid points, are typically
126 * only accessed in recently sampled regions even though the free energy
127 * update is inherently global and affects all points.
128 * The bias allows points thay are non-local at the time the update
129 * was issued to postpone ("skip", as it is called in the code) the update.
130 * A non-local point is defined as a point which has not been sampled since
131 * the last update. Local points are points that have been sampled since
132 * the last update. The (current) set of local points are kept track of by
133 * the bias state and reset after every update. An update is called local
134 * if it only updates local points. Non-local points will temporarily "skip"
135 * the update until next time they are local (or when a global update
136 * is issued). For this to work, the bias keeps a global "clock"
137 * (in HistogramSize) of the number of issued updates. Each PointState
138 * also has its own local "clock" with the counting the number of updates
139 * it has pulled through. When a point updates its state it asserts
140 * that its local clock is synchronized with the global clock.
145 //! Enum for requesting Bias set up with(out) I/O on this rank.
146 enum class ThisRankWillDoIO
148 No, //!< This rank will not do I/O.
149 Yes //!< This rank will do I/O.
155 * \param[in] biasIndexInCollection Index of the bias in collection.
156 * \param[in] awhParams AWH parameters.
157 * \param[in] awhBiasParams Bias parameters.
158 * \param[in] dimParams Bias dimension parameters.
159 * \param[in] beta 1/(k_B T).
160 * \param[in] mdTimeStep The MD time step.
161 * \param[in] numSharingSimulations The number of simulations to share the bias across.
162 * \param[in] biasInitFilename Name of file to read PMF and target from.
163 * \param[in] thisRankWillDoIO Tells whether this MPI rank will do I/O (checkpointing, AWH output), 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 const std::vector<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[out] awhPotential Bias potential.
197 * \param[out] potentialJump Change in bias potential for this bias.
198 * \param[in] commRecord Struct for intra-simulation communication.
199 * \param[in] ms Struct for multi-simulation communication.
201 * \param[in] step Time step.
202 * \param[in] seed Random seed.
203 * \param[in,out] fplog Log file.
204 * \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.
206 gmx::ArrayRef<const double>
207 calcForceAndUpdateBias(const awh_dvec coordValue,
208 double *awhPotential,
209 double *potentialJump,
210 const t_commrec *commRecord,
211 const gmx_multisim_t *ms,
218 * Calculates the convolved bias for a given coordinate value.
220 * The convolved bias is the effective bias acting on the coordinate.
221 * Since the bias here has arbitrary normalization, this only makes
222 * sense as a relative, to other coordinate values, measure of the bias.
224 * \param[in] coordValue The coordinate value.
225 * \returns the convolved bias >= -GMX_FLOAT_MAX.
227 double calcConvolvedBias(const awh_dvec &coordValue) const
229 return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
233 * Restore the bias state from history on the master rank and broadcast it.
235 * \param[in] biasHistory Bias history struct, only allowed to be nullptr on non-master ranks.
236 * \param[in] cr The communication record.
238 void restoreStateFromHistory(const AwhBiasHistory *biasHistory,
239 const t_commrec *cr);
242 * Allocate and initialize a bias history with the given bias state.
244 * This function will be called at the start of a new simulation.
245 * Note that only constant data will be initialized here.
246 * History data is set by \ref updateHistory.
248 * \param[in,out] biasHistory AWH history to initialize.
250 void initHistoryFromState(AwhBiasHistory *biasHistory) const;
253 * Update the bias history with the current state.
255 * \param[out] biasHistory Bias history struct.
257 void updateHistory(AwhBiasHistory *biasHistory) const;
260 * Do all previously skipped updates.
261 * Public for use by tests.
263 void doSkippedUpdatesForAllPoints();
265 //! Returns the dimensionality of the bias.
266 inline int ndim() const
268 return dimParams_.size();
271 /*! \brief Returns the dimension parameters.
273 inline const std::vector<DimParams> &dimParams() const
278 //! Returns the bias parameters
279 inline const BiasParams ¶ms() const
284 //! Returns the global state of the bias.
285 inline const BiasState &state() const
290 //! Returns the index of the bias.
291 inline int biasIndex() const
293 return params_.biasIndex;
296 /*! \brief Return the coordinate value for a grid point.
298 * \param[in] gridPointIndex The index of the grid point.
300 inline const awh_dvec &getGridCoordValue(size_t gridPointIndex) const
302 GMX_ASSERT(gridPointIndex < grid_.numPoints(), "gridPointIndex should be in the range of the grid");
304 return grid_.point(gridPointIndex).coordValue;
309 * Performs statistical checks on the collected histograms and warns if issues are detected.
312 * \param[in] step Time step.
313 * \param[in,out] fplog Output file for warnings.
315 void warnForHistogramAnomalies(double t,
320 * Collect samples for the force correlation analysis on the grid.
322 * \param[in] probWeightNeighbor Probability weight of the neighboring points.
323 * \param[in] t The time.
325 void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
329 /*! \brief Return a const reference to the force correlation grid.
331 const CorrelationGrid &forceCorrelationGrid() const
333 GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr, "forceCorrelationGrid() should only be called with a valid force correlation object");
335 return *forceCorrelationGrid_.get();
338 /*! \brief Return the number of data blocks that have been prepared for writing.
340 int numEnergySubblocksToWrite() const;
342 /*! \brief Write bias data blocks to energy subblocks.
344 * \param[in,out] subblock Energy subblocks to write to.
345 * \returns the number of subblocks written.
347 int writeToEnergySubblocks(t_enxsubblock *subblock) const;
351 const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
352 const Grid grid_; /**< The multidimensional grid organizing the coordinate point locations. */
354 const BiasParams params_; /**< Constant parameters for the method. */
356 BiasState state_; /**< The state, both global and of the grid points */
357 std::vector<int> updateList_; /**< List of points for update for temporary use (could be made another tempWorkSpace) */
359 const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
361 std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
363 /* Force correlation grid */
364 std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
367 std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
369 /* Temporary working vectors used during the update.
370 * These are only here to avoid allocation at every MD step.
372 std::vector < double, AlignedAllocator < double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
373 std::vector<double> tempForce_; /**< Bias force work buffer. */
375 /* Run-local counter to avoid flooding log with warnings. */
376 int numWarningsIssued_; /**< The number of warning issued in the current run. */
381 #endif /* GMX_AWH_BIAS_H */