src/gromacs/awh/biasstate.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019, 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 Grid.
  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 <vector>
  58
  59 #include "gromacs/math/vectypes.h"
  60 #include "gromacs/utility/alignedallocator.h"
  61 #include "gromacs/utility/arrayref.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 struct AwhBiasHistory;
  76 struct AwhBiasParams;
  77 class BiasParams;
  78 class Grid;
  79 class GridAxis;
  80 class PointState;
  81
  82 /*! \internal
  83  * \brief The state of a bias.
  84  *
  85  * The bias state has the current coordinate state: its value and the grid point
  86  * it maps to (the grid point of the umbrella potential if needed). It contains
  87  * a vector with the state for each point on the grid. It also
  88  * counts the number of updates issued and tracks which points have been sampled
  89  * since last update. Finally, the convergence state is a global property set
  90  * ultimately by the histogram size histogramSize in the sub-class HistogramSize,
  91  * since the update sizes are ~ 1/histogramSize.
  92  */
  93 class BiasState
  94 {
  95 public:
  96     /*! \brief Constructor.
  97      *
  98      * Constructs the global state and the point states on a provided
  99      * geometric grid passed in \p grid.
 100      *
 101      * \param[in] awhBiasParams         The Bias parameters from inputrec.
 102      * \param[in] histogramSizeInitial  The estimated initial histogram size.
 103      *                                  This is floating-point, since histograms use weighted
 104      *                                  entries and grow by a floating-point scaling factor.
 105      * \param[in] dimParams             The dimension parameters.
 106      * \param[in] grid                  The bias grid.
 107      */
 108     BiasState(const AwhBiasParams&          awhBiasParams,
 109               double                        histogramSizeInitial,
 110               const std::vector<DimParams>& dimParams,
 111               const Grid&                   grid);
 112
 113     /*! \brief
 114      * Restore the bias state from history.
 115      *
 116      * \param[in] biasHistory  Bias history struct.
 117      * \param[in] grid         The bias grid.
 118      */
 119     void restoreFromHistory(const AwhBiasHistory& biasHistory, const Grid& grid);
 120
 121     /*! \brief
 122      * Broadcast the bias state over the MPI ranks in this simulation.
 123      *
 124      * \param[in] commRecord  Struct for communication.
 125      */
 126     void broadcast(const t_commrec* commRecord);
 127
 128     /*! \brief
 129      * Allocate and initialize a bias history with the given bias state.
 130      *
 131      * This function will be called at the start of a new simulation.
 132      * Note that this only sets the correct size and does produce
 133      * a valid history object, but with all data set to zero.
 134      * Actual history data is set by \ref updateHistory.
 135      *
 136      * \param[in,out] biasHistory  AWH history to initialize.
 137      */
 138     void initHistoryFromState(AwhBiasHistory* biasHistory) const;
 139
 140     /*! \brief
 141      * Update the bias state history with the current state.
 142      *
 143      * \param[out] biasHistory  Bias history struct.
 144      * \param[in]  grid         The bias grid.
 145      */
 146     void updateHistory(AwhBiasHistory* biasHistory, const Grid& grid) const;
 147
 148 private:
 149     /*! \brief Convolves the given PMF using the given AWH bias.
 150      *
 151      * \note: The PMF is in single precision, because it is a statistical
 152      *        quantity and therefore never reaches full float precision.
 153      *
 154      * \param[in] dimParams     The bias dimensions parameters
 155      * \param[in] grid          The grid.
 156      * \param[in,out] convolvedPmf  Array returned will be of the same length as the AWH grid to store the convolved PMF in.
 157      */
 158     void calcConvolvedPmf(const std::vector<DimParams>& dimParams,
 159                           const Grid&                   grid,
 160                           std::vector<float>*           convolvedPmf) const;
 161
 162     /*! \brief
 163      * Convolves the PMF and sets the initial free energy to its convolution.
 164      *
 165      * \param[in] dimParams  The bias dimensions parameters
 166      * \param[in] grid       The bias grid.
 167      */
 168     void setFreeEnergyToConvolvedPmf(const std::vector<DimParams>& dimParams, const Grid& grid);
 169
 170     /*! \brief
 171      * Normalize the PMF histogram.
 172      *
 173      * \param[in] numSharingSims  The number of simulations sharing the bias.
 174      */
 175     void normalizePmf(int numSharingSims);
 176
 177 public:
 178     /*! \brief
 179      * Initialize the state of grid coordinate points.
 180      *
 181      * \param[in] awhBiasParams   Bias parameters from inputrec.
 182      * \param[in] dimParams       The dimension parameters.
 183      * \param[in] grid            The grid.
 184      * \param[in] params          The bias parameters.
 185      * \param[in] filename        Name of file to read PMF and target from.
 186      * \param[in] numBias         The number of biases.
 187      */
 188     void initGridPointState(const AwhBiasParams&          awhBiasParams,
 189                             const std::vector<DimParams>& dimParams,
 190                             const Grid&                   grid,
 191                             const BiasParams&             params,
 192                             const std::string&            filename,
 193                             int                           numBias);
 194
 195     /*! \brief
 196      * Performs statistical checks on the collected histograms and warns if issues are detected.
 197      *
 198      * \param[in]     grid            The grid.
 199      * \param[in]     biasIndex       The index of the bias we are checking for.
 200      * \param[in]     t               Time.
 201      * \param[in,out] fplog           Output file for warnings.
 202      * \param[in]     maxNumWarnings  Don't issue more than this number of warnings.
 203      * \returns the number of warnings issued.
 204      */
 205     int warnForHistogramAnomalies(const Grid& grid, int biasIndex, double t, FILE* fplog, int maxNumWarnings) const;
 206
 207     /*! \brief
 208      * Calculates and sets the force the coordinate experiences from an umbrella centered at the given point.
 209      *
 210      * The umbrella potential is an harmonic potential given by 0.5k(coord value - point value)^2. This
 211      * value is also returned.
 212      *
 213      * \param[in]     dimParams  The bias dimensions parameters.
 214      * \param[in]     grid       The grid.
 215      * \param[in]     point      Point for umbrella center.
 216      * \param[in,out] force      Force vector to set.
 217      * Returns the umbrella potential.
 218      */
 219     double calcUmbrellaForceAndPotential(const std::vector<DimParams>& dimParams,
 220                                          const Grid&                   grid,
 221                                          int                           point,
 222                                          gmx::ArrayRef<double>         force) const;
 223
 224     /*! \brief
 225      * Calculates and sets the convolved force acting on the coordinate.
 226      *
 227      * The convolved force is the weighted sum of forces from umbrellas
 228      * located at each point in the grid.
 229      *
 230      * \param[in]     dimParams           The bias dimensions parameters.
 231      * \param[in]     grid                The grid.
 232      * \param[in]     probWeightNeighbor  Probability weights of the neighbors.
 233      * \param[in]     forceWorkBuffer     Force work buffer, values only used internally.
 234      * \param[in,out] force               Bias force vector to set.
 235      */
 236     void calcConvolvedForce(const std::vector<DimParams>& dimParams,
 237                             const Grid&                   grid,
 238                             gmx::ArrayRef<const double>   probWeightNeighbor,
 239                             gmx::ArrayRef<double>         forceWorkBuffer,
 240                             gmx::ArrayRef<double>         force) const;
 241
 242     /*! \brief
 243      * Move the center point of the umbrella potential.
 244      *
 245      * A new umbrella center is sampled from the biased distibution. Also, the bias
 246      * force is updated and the new potential is return.
 247      *
 248      * This function should only be called when the bias force is not being convolved.
 249      * It is assumed that the probability distribution has been updated.
 250      *
 251      * \param[in] dimParams           Bias dimension parameters.
 252      * \param[in] grid                The grid.
 253      * \param[in] probWeightNeighbor  Probability weights of the neighbors.
 254      * \param[in,out] biasForce       The AWH bias force.
 255      * \param[in] step                Step number, needed for the random number generator.
 256      * \param[in] seed                Random seed.
 257      * \param[in] indexSeed           Second random seed, should be the bias Index.
 258      * \returns the new potential value.
 259      */
 260     double moveUmbrella(const std::vector<DimParams>& dimParams,
 261                         const Grid&                   grid,
 262                         gmx::ArrayRef<const double>   probWeightNeighbor,
 263                         gmx::ArrayRef<double>         biasForce,
 264                         int64_t                       step,
 265                         int64_t                       seed,
 266                         int                           indexSeed);
 267
 268 private:
 269     /*! \brief
 270      * Gets the histogram rescaling factors needed for skipped updates.
 271      *
 272      * \param[in]  params             The bias parameters.
 273      * \param[out] weighthistScaling  Scaling factor for the reference weight histogram.
 274      * \param[out] logPmfsumScaling   Log of the scaling factor for the PMF histogram.
 275      */
 276     void getSkippedUpdateHistogramScaleFactors(const BiasParams& params,
 277                                                double*           weighthistScaling,
 278                                                double*           logPmfsumScaling) const;
 279
 280 public:
 281     /*! \brief
 282      * Do all previously skipped updates.
 283      * Public for use by tests.
 284      *
 285      * \param[in] params  The bias parameters.
 286      */
 287     void doSkippedUpdatesForAllPoints(const BiasParams& params);
 288
 289     /*! \brief
 290      * Do previously skipped updates in this neighborhood.
 291      *
 292      * \param[in] params  The bias parameters.
 293      * \param[in] grid    The grid.
 294      */
 295     void doSkippedUpdatesInNeighborhood(const BiasParams& params, const Grid& grid);
 296
 297 private:
 298     /*! \brief
 299      * Reset the range used to make the local update list.
 300      *
 301      * \param[in] grid  The grid.
 302      */
 303     void resetLocalUpdateRange(const Grid& grid);
 304
 305     /*! \brief
 306      * Returns the new size of the reference weight histogram in the initial stage.
 307      *
 308      * This function also takes care resetting the histogram used for covering checks
 309      * and for exiting the initial stage.
 310      *
 311      * \param[in]     params            The bias parameters.
 312      * \param[in]     t                 Time.
 313      * \param[in]     detectedCovering  True if we detected that the sampling interval has been sufficiently covered.
 314      * \param[in,out] fplog             Log file.
 315      * \returns the new histogram size.
 316      */
 317     double newHistogramSizeInitialStage(const BiasParams& params, double t, bool detectedCovering, FILE* fplog);
 318
 319     /*! \brief
 320      * Check if the sampling region has been covered "enough" or not.
 321      *
 322      * A one-dimensional interval is defined as covered if each point has
 323      * accumulated the same weight as is in the peak of a discretized normal
 324      * distribution. For multiple dimensions, the weights are simply projected
 325      * onto each dimension and the multidimensional space is covered if each
 326      * dimension is.
 327      *
 328      * \note The covering criterion for multiple dimensions could improved, e.g.
 329      * by using a path finding algorithm.
 330      *
 331      * \param[in] params        The bias parameters.
 332      * \param[in] dimParams     Bias dimension parameters.
 333      * \param[in] grid          The grid.
 334      * \param[in] commRecord    Struct for intra-simulation communication.
 335      * \param[in] multiSimComm  Struct for multi-simulation communication.
 336      * \returns true if covered.
 337      */
 338     bool isSamplingRegionCovered(const BiasParams&             params,
 339                                  const std::vector<DimParams>& dimParams,
 340                                  const Grid&                   grid,
 341                                  const t_commrec*              commRecord,
 342                                  const gmx_multisim_t*         multiSimComm) const;
 343
 344     /*! \brief
 345      * Return the new reference weight histogram size for the current update.
 346      *
 347      * This function also takes care of checking for covering in the initial stage.
 348      *
 349      * \param[in]     params   The bias parameters.
 350      * \param[in]     t        Time.
 351      * \param[in]     covered  True if the sampling interval has been covered enough.
 352      * \param[in,out] fplog    Log file.
 353      * \returns the new histogram size.
 354      */
 355     double newHistogramSize(const BiasParams& params, double t, bool covered, FILE* fplog);
 356
 357 public:
 358     /*! \brief
 359      * Update the reaction coordinate value.
 360      *
 361      * \param[in] grid        The bias grid.
 362      * \param[in] coordValue  The current reaction coordinate value (there are no limits on allowed values).
 363      */
 364     void setCoordValue(const Grid& grid, const awh_dvec coordValue)
 365     {
 366         coordState_.setCoordValue(grid, coordValue);
 367     }
 368
 369     /*! \brief
 370      * Performs an update of the bias.
 371      *
 372      * The objective of the update is to use collected samples (probability weights)
 373      * to improve the free energy estimate. For sake of efficiency, the update is
 374      * local whenever possible, meaning that only points that have actually been sampled
 375      * are accessed and updated here. For certain AWH settings or at certain steps
 376      * however, global need to be performed. Besides the actual free energy update, this
 377      * function takes care of ensuring future convergence of the free energy. Convergence
 378      * is obtained by increasing the size of the reference weight histogram in a controlled
 379      * (sometimes dynamic) manner. Also, there are AWH variables that are direct functions
 380      * of the free energy or sampling history that need to be updated here, namely the target
 381      * distribution and the bias function.
 382      *
 383      * \param[in]     dimParams   The dimension parameters.
 384      * \param[in]     grid        The grid.
 385      * \param[in]     params      The bias parameters.
 386      * \param[in]     commRecord  Struct for intra-simulation communication.
 387      * \param[in]     ms          Struct for multi-simulation communication.
 388      * \param[in]     t           Time.
 389      * \param[in]     step        Time step.
 390      * \param[in,out] fplog       Log file.
 391      * \param[in,out] updateList  Work space to store a temporary list.
 392      */
 393     void updateFreeEnergyAndAddSamplesToHistogram(const std::vector<DimParams>& dimParams,
 394                                                   const Grid&                   grid,
 395                                                   const BiasParams&             params,
 396                                                   const t_commrec*              commRecord,
 397                                                   const gmx_multisim_t*         ms,
 398                                                   double                        t,
 399                                                   int64_t                       step,
 400                                                   FILE*                         fplog,
 401                                                   std::vector<int>*             updateList);
 402
 403     /*! \brief
 404      * Update the probability weights and the convolved bias.
 405      *
 406      * Given a coordinate value, each grid point is assigned a probability
 407      * weight, w(point|value), that depends on the current bias function. The sum
 408      * of these weights is needed for normalizing the probability sum to 1 but
 409      * also equals the effective, or convolved, biasing weight for this coordinate
 410      * value. The convolved bias is needed e.g. for extracting the PMF, so we save
 411      * it here since this saves us from doing extra exponential function evaluations
 412      * later on.
 413      *
 414      * \param[in]  dimParams  The bias dimensions parameters
 415      * \param[in]  grid       The grid.
 416      * \param[out] weight     Probability weights of the neighbors, SIMD aligned.
 417      * \returns the convolved bias.
 418      */
 419
 420     double updateProbabilityWeightsAndConvolvedBias(const std::vector<DimParams>& dimParams,
 421                                                     const Grid&                   grid,
 422                                                     std::vector<double, AlignedAllocator<double>>* weight) const;
 423
 424     /*! \brief
 425      * Take samples of the current probability weights for future updates and analysis.
 426      *
 427      * Points in the current neighborhood will now have data meaning they
 428      * need to be included in the local update list of the next update.
 429      * Therefore, the local update range is also update here.
 430      *
 431      * \param[in] grid                The grid.
 432      * \param[in] probWeightNeighbor  Probability weights of the neighbors.
 433      */
 434     void sampleProbabilityWeights(const Grid& grid, gmx::ArrayRef<const double> probWeightNeighbor);
 435
 436     /*! \brief
 437      * Sample the reaction coordinate and PMF for future updates or analysis.
 438      *
 439      * These samples do not affect the (future) sampling and are thus
 440      * pure observables. Statisics of these are stored in the energy file.
 441      *
 442      * \param[in] grid                The grid.
 443      * \param[in] probWeightNeighbor  Probability weights of the neighbors.
 444      * \param[in] convolvedBias       The convolved bias.
 445      */
 446     void sampleCoordAndPmf(const Grid& grid, gmx::ArrayRef<const double> probWeightNeighbor, double convolvedBias);
 447     /*! \brief
 448      * Calculates the convolved bias for a given coordinate value.
 449      *
 450      * The convolved bias is the effective bias acting on the coordinate.
 451      * Since the bias here has arbitrary normalization, this only makes
 452      * sense as a relative, to other coordinate values, measure of the bias.
 453      *
 454      * \note If it turns out to be costly to calculate this pointwise
 455      * the convolved bias for the whole grid could be returned instead.
 456      *
 457      * \param[in] dimParams   The bias dimensions parameters
 458      * \param[in] grid        The grid.
 459      * \param[in] coordValue  Coordinate value.
 460      * \returns the convolved bias >= -GMX_FLOAT_MAX.
 461      */
 462     double calcConvolvedBias(const std::vector<DimParams>& dimParams,
 463                              const Grid&                   grid,
 464                              const awh_dvec&               coordValue) const;
 465
 466     /*! \brief
 467      * Fills the given array with PMF values.
 468      *
 469      * Points outside of the biasing target region will get PMF = GMX_FLOAT_MAX.
 470      * \note: The PMF is in single precision, because it is a statistical
 471      *        quantity and therefore never reaches full float precision.
 472      *
 473      * \param[out] pmf  Array(ref) to be filled with the PMF values, should have the same size as the bias grid.
 474      */
 475     void getPmf(gmx::ArrayRef<float> /*pmf*/) const;
 476
 477     /*! \brief Returns the current coordinate state.
 478      */
 479     const CoordState& coordState() const { return coordState_; }
 480
 481     /*! \brief Returns a const reference to the point state.
 482      */
 483     const std::vector<PointState>& points() const { return points_; }
 484
 485     /*! \brief Returns true if we are in the initial stage.
 486      */
 487     bool inInitialStage() const { return histogramSize_.inInitialStage(); }
 488
 489     /*! \brief Returns the current histogram size.
 490      */
 491     inline HistogramSize histogramSize() const { return histogramSize_; }
 492
 493     /* Data members */
 494 private:
 495     CoordState coordState_; /**< The Current coordinate state */
 496
 497     /* The grid point state */
 498     std::vector<PointState> points_; /**< Vector of state of the grid points */
 499
 500     /* Covering values for each point on the grid */
 501     std::vector<double> weightSumCovering_; /**< Accumulated weights for covering checks */
 502
 503     HistogramSize histogramSize_; /**< Global histogram size related values. */
 504
 505     /* Track the part of the grid sampled since the last update. */
 506     awh_ivec originUpdatelist_; /**< The origin of the rectangular region that has been sampled since last update. */
 507     awh_ivec endUpdatelist_; /**< The end of the rectangular region that has been sampled since last update. */
 508 };
 509
 510 //! Linewidth used for warning output
 511 static const int c_linewidth = 80 - 2;
 512
 513 //! Indent used for warning output
 514 static const int c_indent = 0;
 515
 516 } // namespace gmx
 517
 518 #endif /* GMX_AWH_BIASSTATE_H */