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