2 * This file is part of the GROMACS molecular simulation package.
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.
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 BiasState class.
41 * The data members of this class are the state variables of the bias.
42 * All interaction from the outside happens through the Bias class, which
43 * holds important helper classes such as DimParams and BiasGrid.
44 * This class holds many methods, but more are const methods that compute
45 * properties of the state.
47 * \author Viveca Lindahl
48 * \author Berk Hess <hess@kth.se>
52 #ifndef GMX_AWH_BIASSTATE_H
53 #define GMX_AWH_BIASSTATE_H
60 #include "gromacs/math/vectypes.h"
61 #include "gromacs/utility/alignedallocator.h"
62 #include "gromacs/utility/basedefinitions.h"
63 #include "gromacs/utility/gmxassert.h"
65 #include "coordstate.h"
66 #include "dimparams.h"
67 #include "histogramsize.h"
69 struct gmx_multisim_t;
77 struct AwhBiasHistory;
85 * \brief The state of a bias.
87 * The bias state has the current coordinate state: its value and the grid point
88 * it maps to (the grid point of the umbrella potential if needed). It contains
89 * a vector with the state for each point on the grid. It also
90 * counts the number of updates issued and tracks which points have been sampled
91 * since last update. Finally, the convergence state is a global property set
92 * ultimately by the histogram size histogramSize in the sub-class HistogramSize,
93 * since the update sizes are ~ 1/histogramSize.
98 /*! \brief Constructor.
100 * Constructs the global state and the point states on a provided
101 * geometric grid passed in \p grid.
103 * \param[in] awhBiasParams The Bias parameters from inputrec.
104 * \param[in] histogramSizeInitial The estimated initial histogram size.
105 * This is floating-point, since histograms use weighted
106 * entries and grow by a floating-point scaling factor.
107 * \param[in] dimParams The dimension parameters.
108 * \param[in] grid The bias grid.
110 BiasState(const AwhBiasParams& awhBiasParams,
111 double histogramSizeInitial,
112 const std::vector<DimParams>& dimParams,
113 const BiasGrid& grid);
116 * Restore the bias state from history.
118 * \param[in] biasHistory Bias history struct.
119 * \param[in] grid The bias grid.
121 void restoreFromHistory(const AwhBiasHistory& biasHistory, const BiasGrid& grid);
124 * Broadcast the bias state over the MPI ranks in this simulation.
126 * \param[in] commRecord Struct for communication.
128 void broadcast(const t_commrec* commRecord);
131 * Allocate and initialize a bias history with the given bias state.
133 * This function will be called at the start of a new simulation.
134 * Note that this only sets the correct size and does produce
135 * a valid history object, but with all data set to zero.
136 * Actual history data is set by \ref updateHistory.
138 * \param[in,out] biasHistory AWH history to initialize.
140 void initHistoryFromState(AwhBiasHistory* biasHistory) const;
143 * Update the bias state history with the current state.
145 * \param[out] biasHistory Bias history struct.
146 * \param[in] grid The bias grid.
148 void updateHistory(AwhBiasHistory* biasHistory, const BiasGrid& grid) const;
151 /*! \brief Convolves the given PMF using the given AWH bias.
153 * \note: The PMF is in single precision, because it is a statistical
154 * quantity and therefore never reaches full float precision.
156 * \param[in] dimParams The bias dimensions parameters
157 * \param[in] grid The grid.
158 * \param[in,out] convolvedPmf Array returned will be of the same length as the AWH grid to store the convolved PMF in.
160 void calcConvolvedPmf(const std::vector<DimParams>& dimParams,
161 const BiasGrid& grid,
162 std::vector<float>* convolvedPmf) const;
165 * Convolves the PMF and sets the initial free energy to its convolution.
167 * \param[in] dimParams The bias dimensions parameters
168 * \param[in] grid The bias grid.
170 void setFreeEnergyToConvolvedPmf(const std::vector<DimParams>& dimParams, const BiasGrid& grid);
173 * Normalize the PMF histogram.
175 * \param[in] numSharingSims The number of simulations sharing the bias.
177 void normalizePmf(int numSharingSims);
181 * Initialize the state of grid coordinate points.
183 * \param[in] awhBiasParams Bias parameters from inputrec.
184 * \param[in] dimParams The dimension parameters.
185 * \param[in] grid The grid.
186 * \param[in] params The bias parameters.
187 * \param[in] filename Name of file to read PMF and target from.
188 * \param[in] numBias The number of biases.
190 void initGridPointState(const AwhBiasParams& awhBiasParams,
191 const std::vector<DimParams>& dimParams,
192 const BiasGrid& grid,
193 const BiasParams& params,
194 const std::string& filename,
198 * Performs statistical checks on the collected histograms and warns if issues are detected.
200 * \param[in] grid The grid.
201 * \param[in] biasIndex The index of the bias we are checking for.
203 * \param[in,out] fplog Output file for warnings.
204 * \param[in] maxNumWarnings Don't issue more than this number of warnings.
205 * \returns the number of warnings issued.
207 int warnForHistogramAnomalies(const BiasGrid& grid, int biasIndex, double t, FILE* fplog, int maxNumWarnings) const;
210 * Calculates and sets the force the coordinate experiences from an umbrella centered at the given point.
212 * The umbrella potential is an harmonic potential given by 0.5k(coord value - point value)^2. This
213 * value is also returned.
215 * \param[in] dimParams The bias dimensions parameters.
216 * \param[in] grid The grid.
217 * \param[in] point Point for umbrella center.
218 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
219 * points. The array is of length numLambdas+1, where numLambdas is the number of free
220 * energy lambda states. Element 0 in the array is the dHdL
221 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
222 * neighboring lambda states (also including the current state). When there are no free
223 * energy lambda state dimensions this can be empty.
224 * \param[in,out] force Force vector to set.
225 * Returns the umbrella potential.
227 double calcUmbrellaForceAndPotential(const std::vector<DimParams>& dimParams,
228 const BiasGrid& grid,
230 ArrayRef<const double> neighborLambdaDhdl,
231 gmx::ArrayRef<double> force) const;
234 * Calculates and sets the convolved force acting on the coordinate.
236 * The convolved force is the weighted sum of forces from umbrellas
237 * located at each point in the grid.
239 * \param[in] dimParams The bias dimensions parameters.
240 * \param[in] grid The grid.
241 * \param[in] probWeightNeighbor Probability weights of the neighbors.
242 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
243 * points. The array is of length numLambdas+1, where numLambdas is the number of free
244 * energy lambda states. Element 0 in the array is the dHdL
245 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
246 * neighboring lambda states (also including the current state). When there are no free
247 * energy lambda state dimensions this can be empty.
248 * \param[in] forceWorkBuffer Force work buffer, values only used internally.
249 * \param[in,out] force Bias force vector to set.
251 void calcConvolvedForce(const std::vector<DimParams>& dimParams,
252 const BiasGrid& grid,
253 gmx::ArrayRef<const double> probWeightNeighbor,
254 ArrayRef<const double> neighborLambdaDhdl,
255 gmx::ArrayRef<double> forceWorkBuffer,
256 gmx::ArrayRef<double> force) const;
259 * Move the center point of the umbrella potential.
261 * A new umbrella center is sampled from the biased distibution. Also, the bias
262 * force is updated and the new potential is return.
264 * This function should only be called when the bias force is not being convolved.
265 * It is assumed that the probability distribution has been updated.
267 * \param[in] dimParams Bias dimension parameters.
268 * \param[in] grid The grid.
269 * \param[in] probWeightNeighbor Probability weights of the neighbors.
270 * \param[in] neighborLambdaDhdl An array containing the dHdL at the neighboring lambda
271 * points. The array is of length numLambdas+1, where numLambdas is the number of free
272 * energy lambda states. Element 0 in the array is the dHdL
273 * of the current state and elements 1..numLambdas contain the dHdL of the system in the
274 * neighboring lambda states (also including the current state). When there are no free
275 * energy lambda state dimensions this can be empty.
276 * \param[in,out] biasForce The AWH bias force.
277 * \param[in] step Step number, needed for the random number generator.
278 * \param[in] seed Random seed.
279 * \param[in] indexSeed Second random seed, should be the bias Index.
280 * \param[in] onlySampleUmbrellaGridpoint Only sample the umbrella gridpoint without calculating
281 * force and potential.
282 * \returns the new potential value.
284 double moveUmbrella(const std::vector<DimParams>& dimParams,
285 const BiasGrid& grid,
286 gmx::ArrayRef<const double> probWeightNeighbor,
287 ArrayRef<const double> neighborLambdaDhdl,
288 gmx::ArrayRef<double> biasForce,
292 bool onlySampleUmbrellaGridpoint);
296 * Gets the histogram rescaling factors needed for skipped updates.
298 * \param[in] params The bias parameters.
299 * \param[out] weighthistScaling Scaling factor for the reference weight histogram.
300 * \param[out] logPmfsumScaling Log of the scaling factor for the PMF histogram.
302 void getSkippedUpdateHistogramScaleFactors(const BiasParams& params,
303 double* weighthistScaling,
304 double* logPmfsumScaling) const;
308 * Do all previously skipped updates.
309 * Public for use by tests.
311 * \param[in] params The bias parameters.
313 void doSkippedUpdatesForAllPoints(const BiasParams& params);
316 * Do previously skipped updates in this neighborhood.
318 * \param[in] params The bias parameters.
319 * \param[in] grid The grid.
321 void doSkippedUpdatesInNeighborhood(const BiasParams& params, const BiasGrid& grid);
325 * Reset the range used to make the local update list.
327 * \param[in] grid The grid.
329 void resetLocalUpdateRange(const BiasGrid& grid);
332 * Returns the new size of the reference weight histogram in the initial stage.
334 * This function also takes care resetting the histogram used for covering checks
335 * and for exiting the initial stage.
337 * \param[in] params The bias parameters.
339 * \param[in] detectedCovering True if we detected that the sampling interval has been sufficiently covered.
340 * \param[in,out] fplog Log file.
341 * \returns the new histogram size.
343 double newHistogramSizeInitialStage(const BiasParams& params, double t, bool detectedCovering, FILE* fplog);
346 * Check if the sampling region has been covered "enough" or not.
348 * A one-dimensional interval is defined as covered if each point has
349 * accumulated the same weight as is in the peak of a discretized normal
350 * distribution. For multiple dimensions, the weights are simply projected
351 * onto each dimension and the multidimensional space is covered if each
354 * \note The covering criterion for multiple dimensions could improved, e.g.
355 * by using a path finding algorithm.
357 * \param[in] params The bias parameters.
358 * \param[in] dimParams Bias dimension parameters.
359 * \param[in] grid The grid.
360 * \param[in] commRecord Struct for intra-simulation communication.
361 * \param[in] multiSimComm Struct for multi-simulation communication.
362 * \returns true if covered.
364 bool isSamplingRegionCovered(const BiasParams& params,
365 const std::vector<DimParams>& dimParams,
366 const BiasGrid& grid,
367 const t_commrec* commRecord,
368 const gmx_multisim_t* multiSimComm) const;
371 * Return the new reference weight histogram size for the current update.
373 * This function also takes care of checking for covering in the initial stage.
375 * \param[in] params The bias parameters.
377 * \param[in] covered True if the sampling interval has been covered enough.
378 * \param[in,out] fplog Log file.
379 * \returns the new histogram size.
381 double newHistogramSize(const BiasParams& params, double t, bool covered, FILE* fplog);
385 * Update the reaction coordinate value.
387 * \param[in] grid The bias grid.
388 * \param[in] coordValue The current reaction coordinate value (there are no limits on allowed values).
390 void setCoordValue(const BiasGrid& grid, const awh_dvec coordValue)
392 coordState_.setCoordValue(grid, coordValue);
396 * Performs an update of the bias.
398 * The objective of the update is to use collected samples (probability weights)
399 * to improve the free energy estimate. For sake of efficiency, the update is
400 * local whenever possible, meaning that only points that have actually been sampled
401 * are accessed and updated here. For certain AWH settings or at certain steps
402 * however, global need to be performed. Besides the actual free energy update, this
403 * function takes care of ensuring future convergence of the free energy. Convergence
404 * is obtained by increasing the size of the reference weight histogram in a controlled
405 * (sometimes dynamic) manner. Also, there are AWH variables that are direct functions
406 * of the free energy or sampling history that need to be updated here, namely the target
407 * distribution and the bias function.
409 * \param[in] dimParams The dimension parameters.
410 * \param[in] grid The grid.
411 * \param[in] params The bias parameters.
412 * \param[in] commRecord Struct for intra-simulation communication.
413 * \param[in] ms Struct for multi-simulation communication.
415 * \param[in] step Time step.
416 * \param[in,out] fplog Log file.
417 * \param[in,out] updateList Work space to store a temporary list.
419 void updateFreeEnergyAndAddSamplesToHistogram(const std::vector<DimParams>& dimParams,
420 const BiasGrid& grid,
421 const BiasParams& params,
422 const t_commrec* commRecord,
423 const gmx_multisim_t* ms,
427 std::vector<int>* updateList);
430 * Update the probability weights and the convolved bias.
432 * Given a coordinate value, each grid point is assigned a probability
433 * weight, w(point|value), that depends on the current bias function. The sum
434 * of these weights is needed for normalizing the probability sum to 1 but
435 * also equals the effective, or convolved, biasing weight for this coordinate
436 * value. The convolved bias is needed e.g. for extracting the PMF, so we save
437 * it here since this saves us from doing extra exponential function evaluations
440 * \param[in] dimParams The bias dimensions parameters
441 * \param[in] grid The grid.
442 * \param[in] neighborLambdaEnergies An array containing the energy of the system
443 * in neighboring lambdas. The array is of length numLambdas+1, where numLambdas is
444 * the number of free energy lambda states. Element 0 in the array is the energy
445 * of the current state and elements 1..numLambdas contain the energy of the system in the
446 * neighboring lambda states (also including the current state). When there are no free
447 * energy lambda state dimensions this can be empty.
448 * \param[out] weight Probability weights of the neighbors, SIMD aligned.
449 * \returns the convolved bias.
452 double updateProbabilityWeightsAndConvolvedBias(const std::vector<DimParams>& dimParams,
453 const BiasGrid& grid,
454 ArrayRef<const double> neighborLambdaEnergies,
455 std::vector<double, AlignedAllocator<double>>* weight) const;
458 * Take samples of the current probability weights for future updates and analysis.
460 * Points in the current neighborhood will now have data meaning they
461 * need to be included in the local update list of the next update.
462 * Therefore, the local update range is also update here.
464 * \param[in] grid The grid.
465 * \param[in] probWeightNeighbor Probability weights of the neighbors.
467 void sampleProbabilityWeights(const BiasGrid& grid, gmx::ArrayRef<const double> probWeightNeighbor);
470 * Sample the reaction coordinate and PMF for future updates or analysis.
472 * These samples do not affect the (future) sampling and are thus
473 * pure observables. Statisics of these are stored in the energy file.
475 * \param[in] dimParams The bias dimensions parameters
476 * \param[in] grid The grid.
477 * \param[in] probWeightNeighbor Probability weights of the neighbors.
478 * \param[in] convolvedBias The convolved bias.
480 void sampleCoordAndPmf(const std::vector<DimParams>& dimParams,
481 const BiasGrid& grid,
482 gmx::ArrayRef<const double> probWeightNeighbor,
483 double convolvedBias);
485 * Calculates the convolved bias for a given coordinate value.
487 * The convolved bias is the effective bias acting on the coordinate.
488 * Since the bias here has arbitrary normalization, this only makes
489 * sense as a relative, to other coordinate values, measure of the bias.
491 * \note If it turns out to be costly to calculate this pointwise
492 * the convolved bias for the whole grid could be returned instead.
494 * \param[in] dimParams The bias dimensions parameters
495 * \param[in] grid The grid.
496 * \param[in] coordValue Coordinate value.
497 * \returns the convolved bias >= -GMX_FLOAT_MAX.
499 double calcConvolvedBias(const std::vector<DimParams>& dimParams,
500 const BiasGrid& grid,
501 const awh_dvec& coordValue) const;
504 * Fills the given array with PMF values.
506 * Points outside of the biasing target region will get PMF = GMX_FLOAT_MAX.
507 * \note: The PMF is in single precision, because it is a statistical
508 * quantity and therefore never reaches full float precision.
510 * \param[out] pmf Array(ref) to be filled with the PMF values, should have the same size as the bias grid.
512 void getPmf(gmx::ArrayRef<float> /*pmf*/) const;
514 /*! \brief Returns the current coordinate state.
516 const CoordState& coordState() const { return coordState_; }
518 /*! \brief Returns a const reference to the point state.
520 const std::vector<PointState>& points() const { return points_; }
522 /*! \brief Returns true if we are in the initial stage.
524 bool inInitialStage() const { return histogramSize_.inInitialStage(); }
526 /*! \brief Returns the current histogram size.
528 inline HistogramSize histogramSize() const { return histogramSize_; }
530 /*! \brief Sets the umbrella grid point to the current grid point
532 void setUmbrellaGridpointToGridpoint() { coordState_.setUmbrellaGridpointToGridpoint(); }
536 CoordState coordState_; /**< The Current coordinate state */
538 /* The grid point state */
539 std::vector<PointState> points_; /**< Vector of state of the grid points */
541 /* Covering values for each point on the grid */
542 std::vector<double> weightSumCovering_; /**< Accumulated weights for covering checks */
544 HistogramSize histogramSize_; /**< Global histogram size related values. */
546 /* Track the part of the grid sampled since the last update. */
547 awh_ivec originUpdatelist_; /**< The origin of the rectangular region that has been sampled since last update. */
548 awh_ivec endUpdatelist_; /**< The end of the rectangular region that has been sampled since last update. */
551 //! Linewidth used for warning output
552 static const int c_linewidth = 80 - 2;
554 //! Indent used for warning output
555 static const int c_indent = 0;
559 #endif /* GMX_AWH_BIASSTATE_H */