src/gromacs/applied_forces/awh/biasstate.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  *
  38  * \brief
  39  * Declares the BiasState class.
  40  *
  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.
  46  *
  47  * \author Viveca Lindahl
  48  * \author Berk Hess <hess@kth.se>
  49  * \ingroup module_awh
  50  */
  51
  52 #ifndef GMX_AWH_BIASSTATE_H
  53 #define GMX_AWH_BIASSTATE_H
  54
  55 #include <cstdio>
  56
  57 #include <string>
  58 #include <vector>
  59
  60 #include "gromacs/math/vectypes.h"
  61 #include "gromacs/utility/alignedallocator.h"
  62 #include "gromacs/utility/basedefinitions.h"
  63 #include "gromacs/utility/gmxassert.h"
  64
  65 #include "coordstate.h"
  66 #include "dimparams.h"
  67 #include "histogramsize.h"
  68
  69 struct gmx_multisim_t;
  70 struct t_commrec;
  71
  72 namespace gmx
  73 {
  74
  75 template<typename>
  76 class ArrayRef;
  77 struct AwhBiasHistory;
  78 struct AwhBiasParams;
  79 class BiasParams;
  80 class BiasGrid;
  81 class GridAxis;
  82 class PointState;
  83
  84 /*! \internal
  85  * \brief The state of a bias.
  86  *
  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.
  94  */
  95 class BiasState
  96 {
  97 public:
  98     /*! \brief Constructor.
  99      *
 100      * Constructs the global state and the point states on a provided
 101      * geometric grid passed in \p grid.
 102      *
 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.
 109      */
 110     BiasState(const AwhBiasParams&      awhBiasParams,
 111               double                    histogramSizeInitial,
 112               ArrayRef<const DimParams> dimParams,
 113               const BiasGrid&           grid);
 114
 115     /*! \brief
 116      * Restore the bias state from history.
 117      *
 118      * \param[in] biasHistory  Bias history struct.
 119      * \param[in] grid         The bias grid.
 120      */
 121     void restoreFromHistory(const AwhBiasHistory& biasHistory, const BiasGrid& grid);
 122
 123     /*! \brief
 124      * Broadcast the bias state over the MPI ranks in this simulation.
 125      *
 126      * \param[in] commRecord  Struct for communication.
 127      */
 128     void broadcast(const t_commrec* commRecord);
 129
 130     /*! \brief
 131      * Allocate and initialize a bias history with the given bias state.
 132      *
 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.
 137      *
 138      * \param[in,out] biasHistory  AWH history to initialize.
 139      */
 140     void initHistoryFromState(AwhBiasHistory* biasHistory) const;
 141
 142     /*! \brief
 143      * Update the bias state history with the current state.
 144      *
 145      * \param[out] biasHistory  Bias history struct.
 146      * \param[in]  grid         The bias grid.
 147      */
 148     void updateHistory(AwhBiasHistory* biasHistory, const BiasGrid& grid) const;
 149
 150 private:
 151     /*! \brief Convolves the given PMF using the given AWH bias.
 152      *
 153      * \note: The PMF is in single precision, because it is a statistical
 154      *        quantity and therefore never reaches full float precision.
 155      *
 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.
 159      */
 160     void calcConvolvedPmf(ArrayRef<const DimParams> dimParams,
 161                           const BiasGrid&           grid,
 162                           std::vector<float>*       convolvedPmf) const;
 163
 164     /*! \brief
 165      * Convolves the PMF and sets the initial free energy to its convolution.
 166      *
 167      * \param[in] dimParams  The bias dimensions parameters
 168      * \param[in] grid       The bias grid.
 169      */
 170     void setFreeEnergyToConvolvedPmf(ArrayRef<const DimParams> dimParams, const BiasGrid& grid);
 171
 172     /*! \brief
 173      * Normalize the PMF histogram.
 174      *
 175      * \param[in] numSharingSims  The number of simulations sharing the bias.
 176      */
 177     void normalizePmf(int numSharingSims);
 178
 179 public:
 180     /*! \brief
 181      * Initialize the state of grid coordinate points.
 182      *
 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.
 189      */
 190     void initGridPointState(const AwhBiasParams&      awhBiasParams,
 191                             ArrayRef<const DimParams> dimParams,
 192                             const BiasGrid&           grid,
 193                             const BiasParams&         params,
 194                             const std::string&        filename,
 195                             int                       numBias);
 196
 197     /*! \brief
 198      * Performs statistical checks on the collected histograms and warns if issues are detected.
 199      *
 200      * \param[in]     grid            The grid.
 201      * \param[in]     biasIndex       The index of the bias we are checking for.
 202      * \param[in]     t               Time.
 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.
 206      */
 207     int warnForHistogramAnomalies(const BiasGrid& grid, int biasIndex, double t, FILE* fplog, int maxNumWarnings) const;
 208
 209     /*! \brief
 210      * Calculates and sets the force the coordinate experiences from an umbrella centered at the given point.
 211      *
 212      * The umbrella potential is an harmonic potential given by 0.5k(coord value - point value)^2. This
 213      * value is also returned.
 214      *
 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.
 226      */
 227     double calcUmbrellaForceAndPotential(ArrayRef<const DimParams> dimParams,
 228                                          const BiasGrid&           grid,
 229                                          int                       point,
 230                                          ArrayRef<const double>    neighborLambdaDhdl,
 231                                          ArrayRef<double>          force) const;
 232
 233     /*! \brief
 234      * Calculates and sets the convolved force acting on the coordinate.
 235      *
 236      * The convolved force is the weighted sum of forces from umbrellas
 237      * located at each point in the grid.
 238      *
 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.
 250      */
 251     void calcConvolvedForce(ArrayRef<const DimParams> dimParams,
 252                             const BiasGrid&           grid,
 253                             ArrayRef<const double>    probWeightNeighbor,
 254                             ArrayRef<const double>    neighborLambdaDhdl,
 255                             ArrayRef<double>          forceWorkBuffer,
 256                             ArrayRef<double>          force) const;
 257
 258     /*! \brief
 259      * Move the center point of the umbrella potential.
 260      *
 261      * A new umbrella center is sampled from the biased distibution. Also, the bias
 262      * force is updated and the new potential is return.
 263      *
 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.
 266      *
 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.
 283      */
 284     double moveUmbrella(ArrayRef<const DimParams> dimParams,
 285                         const BiasGrid&           grid,
 286                         ArrayRef<const double>    probWeightNeighbor,
 287                         ArrayRef<const double>    neighborLambdaDhdl,
 288                         ArrayRef<double>          biasForce,
 289                         int64_t                   step,
 290                         int64_t                   seed,
 291                         int                       indexSeed,
 292                         bool                      onlySampleUmbrellaGridpoint);
 293
 294 private:
 295     /*! \brief
 296      * Gets the histogram rescaling factors needed for skipped updates.
 297      *
 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.
 301      */
 302     void getSkippedUpdateHistogramScaleFactors(const BiasParams& params,
 303                                                double*           weighthistScaling,
 304                                                double*           logPmfsumScaling) const;
 305
 306 public:
 307     /*! \brief
 308      * Do all previously skipped updates.
 309      * Public for use by tests.
 310      *
 311      * \param[in] params  The bias parameters.
 312      */
 313     void doSkippedUpdatesForAllPoints(const BiasParams& params);
 314
 315     /*! \brief
 316      * Do previously skipped updates in this neighborhood.
 317      *
 318      * \param[in] params  The bias parameters.
 319      * \param[in] grid    The grid.
 320      */
 321     void doSkippedUpdatesInNeighborhood(const BiasParams& params, const BiasGrid& grid);
 322
 323 private:
 324     /*! \brief
 325      * Reset the range used to make the local update list.
 326      *
 327      * \param[in] grid  The grid.
 328      */
 329     void resetLocalUpdateRange(const BiasGrid& grid);
 330
 331     /*! \brief
 332      * Returns the new size of the reference weight histogram in the initial stage.
 333      *
 334      * This function also takes care resetting the histogram used for covering checks
 335      * and for exiting the initial stage.
 336      *
 337      * \param[in]     params            The bias parameters.
 338      * \param[in]     t                 Time.
 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.
 342      */
 343     double newHistogramSizeInitialStage(const BiasParams& params, double t, bool detectedCovering, FILE* fplog);
 344
 345     /*! \brief
 346      * Check if the sampling region has been covered "enough" or not.
 347      *
 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
 352      * dimension is.
 353      *
 354      * \note The covering criterion for multiple dimensions could improved, e.g.
 355      * by using a path finding algorithm.
 356      *
 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.
 363      */
 364     bool isSamplingRegionCovered(const BiasParams&         params,
 365                                  ArrayRef<const DimParams> dimParams,
 366                                  const BiasGrid&           grid,
 367                                  const t_commrec*          commRecord,
 368                                  const gmx_multisim_t*     multiSimComm) const;
 369
 370     /*! \brief
 371      * Return the new reference weight histogram size for the current update.
 372      *
 373      * This function also takes care of checking for covering in the initial stage.
 374      *
 375      * \param[in]     params   The bias parameters.
 376      * \param[in]     t        Time.
 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.
 380      */
 381     double newHistogramSize(const BiasParams& params, double t, bool covered, FILE* fplog);
 382
 383 public:
 384     /*! \brief
 385      * Update the reaction coordinate value.
 386      *
 387      * \param[in] grid        The bias grid.
 388      * \param[in] coordValue  The current reaction coordinate value (there are no limits on allowed values).
 389      */
 390     void setCoordValue(const BiasGrid& grid, const awh_dvec coordValue)
 391     {
 392         coordState_.setCoordValue(grid, coordValue);
 393     }
 394
 395     /*! \brief
 396      * Performs an update of the bias.
 397      *
 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.
 408      *
 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.
 414      * \param[in]     t           Time.
 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.
 418      */
 419     void updateFreeEnergyAndAddSamplesToHistogram(ArrayRef<const DimParams> dimParams,
 420                                                   const BiasGrid&           grid,
 421                                                   const BiasParams&         params,
 422                                                   const t_commrec*          commRecord,
 423                                                   const gmx_multisim_t*     ms,
 424                                                   double                    t,
 425                                                   int64_t                   step,
 426                                                   FILE*                     fplog,
 427                                                   std::vector<int>*         updateList);
 428
 429     /*! \brief
 430      * Update the probability weights and the convolved bias.
 431      *
 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
 438      * later on.
 439      *
 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.
 450      */
 451
 452     double updateProbabilityWeightsAndConvolvedBias(ArrayRef<const DimParams> dimParams,
 453                                                     const BiasGrid&           grid,
 454                                                     ArrayRef<const double> neighborLambdaEnergies,
 455                                                     std::vector<double, AlignedAllocator<double>>* weight) const;
 456
 457     /*! \brief
 458      * Take samples of the current probability weights for future updates and analysis.
 459      *
 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.
 463      *
 464      * \param[in] grid                The grid.
 465      * \param[in] probWeightNeighbor  Probability weights of the neighbors.
 466      */
 467     void sampleProbabilityWeights(const BiasGrid& grid, ArrayRef<const double> probWeightNeighbor);
 468
 469     /*! \brief
 470      * Sample the reaction coordinate and PMF for future updates or analysis.
 471      *
 472      * These samples do not affect the (future) sampling and are thus
 473      * pure observables. Statisics of these are stored in the energy file.
 474      *
 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.
 479      */
 480     void sampleCoordAndPmf(const std::vector<DimParams>& dimParams,
 481                            const BiasGrid&               grid,
 482                            ArrayRef<const double>        probWeightNeighbor,
 483                            double                        convolvedBias);
 484     /*! \brief
 485      * Calculates the convolved bias for a given coordinate value.
 486      *
 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.
 490      *
 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.
 493      *
 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.
 498      */
 499     double calcConvolvedBias(ArrayRef<const DimParams> dimParams,
 500                              const BiasGrid&           grid,
 501                              const awh_dvec&           coordValue) const;
 502
 503     /*! \brief
 504      * Fills the given array with PMF values.
 505      *
 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.
 509      *
 510      * \param[out] pmf  Array(ref) to be filled with the PMF values, should have the same size as the bias grid.
 511      */
 512     void getPmf(ArrayRef<float> /*pmf*/) const;
 513
 514     /*! \brief Returns the current coordinate state.
 515      */
 516     const CoordState& coordState() const { return coordState_; }
 517
 518     /*! \brief Returns a const reference to the point state.
 519      */
 520     const std::vector<PointState>& points() const { return points_; }
 521
 522     /*! \brief Returns true if we are in the initial stage.
 523      */
 524     bool inInitialStage() const { return histogramSize_.inInitialStage(); }
 525
 526     /*! \brief Returns the current histogram size.
 527      */
 528     inline HistogramSize histogramSize() const { return histogramSize_; }
 529
 530     /*! \brief Sets the umbrella grid point to the current grid point
 531      */
 532     void setUmbrellaGridpointToGridpoint() { coordState_.setUmbrellaGridpointToGridpoint(); }
 533
 534     /* Data members */
 535 private:
 536     CoordState coordState_; /**< The Current coordinate state */
 537
 538     /* The grid point state */
 539     std::vector<PointState> points_; /**< Vector of state of the grid points */
 540
 541     /* Covering values for each point on the grid */
 542     std::vector<double> weightSumCovering_; /**< Accumulated weights for covering checks */
 543
 544     HistogramSize histogramSize_; /**< Global histogram size related values. */
 545
 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. */
 549 };
 550
 551 //! Linewidth used for warning output
 552 static const int c_linewidth = 80 - 2;
 553
 554 //! Indent used for warning output
 555 static const int c_indent = 0;
 556
 557 } // namespace gmx
 558
 559 #endif /* GMX_AWH_BIASSTATE_H */