src/gromacs/awh/bias.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 Bias class.
  40  *
  41  * This class is essentially a wrapper around the BiasState class.
  42  * In addition to BiasState, it holds all data that BiasState needs
  43  * to update the bias. Interaction of the outside world, such as updating
  44  * BiasState or extracting bias data all happen through Bias.
  45  *
  46  * \author Viveca Lindahl
  47  * \author Berk Hess <hess@kth.se>
  48  * \ingroup module_awh
  49  */
  50
  51 #ifndef GMX_AWH_BIAS_H
  52 #define GMX_AWH_BIAS_H
  53
  54 #include <memory>
  55 #include <vector>
  56
  57 #include "gromacs/math/vectypes.h"
  58 #include "gromacs/utility/alignedallocator.h"
  59 #include "gromacs/utility/basedefinitions.h"
  60 #include "gromacs/utility/gmxassert.h"
  61
  62 #include "biasparams.h"
  63 #include "biasstate.h"
  64 #include "biaswriter.h"
  65 #include "dimparams.h"
  66 #include "grid.h"
  67
  68 struct gmx_multisim_t;
  69 struct t_commrec;
  70 struct t_enxsubblock;
  71
  72 namespace gmx
  73 {
  74
  75 struct AwhBiasHistory;
  76 struct AwhBiasParams;
  77 struct AwhHistory;
  78 struct AwhParams;
  79 struct AwhPointStateHistory;
  80 class CorrelationGrid;
  81 class Grid;
  82 class GridAxis;
  83 class PointState;
  84
  85 /*! \internal
  86  * \brief A bias acting on a multidimensional coordinate.
  87  *
  88  * At each step AWH should provide its biases with updated
  89  * values of their coordinates. Each bias provides AWH with an updated
  90  * bias forces and the corresponding potential.
  91  *
  92  * See the user manual for details on the algorithm and equations.
  93  *
  94  * The bias is responsible for keeping and updating a free energy estimate
  95  * along the coordinate. The bias potential is basically a function of the
  96  * free energy estimate and so also changes by the update.
  97  * The free energy update is based on information from coordinate samples
  98  * collected at a constant bias potential, between updates.
  99  *
 100  * The bias keeps a grid with coordinate points that organizes spatial
 101  * information about the coordinate. The grid has the the same geometry
 102  * as the coordinate, i.e. they have the same dimensionality and periodicity
 103  * (if any). The number of points in the grid sets the resolution of
 104  * the collected data and its extent defines the sampling region of interest.
 105  *
 106  * Each coordinate point has further statistical properties and function values
 107  * which a grid point does not know about. E.g., for the bias each coordinate point
 108  * is associated with values of the bias, free energy and target distribution,
 109  * accumulated sampling weight, etc. For this the bias attaches to each grid
 110  * point a state. The grid + vector of point states are the bias coordinate points.
 111  *
 112  * The bias has a fairly complex global state keeping track of where
 113  * the system (coordinate) currently is (CoordState), where it has
 114  * sampled since the last update (BiasState) and controlling the free energy
 115  * convergence rate (HistogramSize).
 116  *
 117  * Partly, the complexity comes from the bias having two convergence stages:
 118  * an initial stage which in an heuristic, non-deterministic way restricts
 119  * the early convergence rate for sake of robustness; and a final stage
 120  * where the convergence rate is constant. The length of the initial stage
 121  * depends on the sampling and is unknown beforehand.
 122  *
 123  * Another complexity comes from the fact that coordinate points,
 124  * for sake of efficiency in the case of many grid points, are typically
 125  * only accessed in recently sampled regions even though the free energy
 126  * update is inherently global and affects all points.
 127  * The bias allows points thay are non-local at the time the update
 128  * was issued to postpone ("skip", as it is called in the code) the update.
 129  * A non-local point is defined as a point which has not been sampled since
 130  * the last update. Local points are points that have been sampled since
 131  * the last update. The (current) set of local points are kept track of by
 132  * the bias state and reset after every update. An update is called local
 133  * if it only updates local points. Non-local points will temporarily "skip"
 134  * the update until next time they are local (or when a global update
 135  * is issued). For this to work, the bias keeps a global "clock"
 136  * (in HistogramSize) of the number of issued updates. Each PointState
 137  * also has its own local "clock" with the counting the number of updates
 138  * it has pulled through. When a point updates its state it asserts
 139  * that its local clock is synchronized with the global clock.
 140  */
 141 class Bias
 142 {
 143     public:
 144         //! Enum for requesting Bias set up with(out) I/O on this rank.
 145         enum class ThisRankWillDoIO
 146         {
 147             No,   //!< This rank will not do I/O.
 148             Yes   //!< This rank will do I/O.
 149         };
 150
 151         /*! \brief
 152          * Constructor.
 153          *
 154          * \param[in] biasIndexInCollection  Index of the bias in collection.
 155          * \param[in] awhParams              AWH parameters.
 156          * \param[in] awhBiasParams          Bias parameters.
 157          * \param[in] dimParams              Bias dimension parameters.
 158          * \param[in] beta                   1/(k_B T).
 159          * \param[in] mdTimeStep             The MD time step.
 160          * \param[in] numSharingSimulations  The number of simulations to share the bias across.
 161          * \param[in] biasInitFilename       Name of file to read PMF and target from.
 162          * \param[in] thisRankWillDoIO       Tells whether this MPI rank will do I/O (checkpointing, AWH output), normally (only) the master rank does I/O.
 163          * \param[in] disableUpdateSkips     If to disable update skips, useful for testing.
 164          */
 165         Bias(int                             biasIndexInCollection,
 166              const AwhParams                &awhParams,
 167              const AwhBiasParams            &awhBiasParams,
 168              const std::vector<DimParams>   &dimParams,
 169              double                          beta,
 170              double                          mdTimeStep,
 171              int                             numSharingSimulations,
 172              const std::string              &biasInitFilename,
 173              ThisRankWillDoIO                thisRankWillDoIO,
 174              BiasParams::DisableUpdateSkips  disableUpdateSkips = BiasParams::DisableUpdateSkips::no);
 175
 176         /*! \brief
 177          * Print information about initialization to log file.
 178          *
 179          * Prints information about AWH variables that are set internally
 180          * but might be of interest to the user.
 181          *
 182          * \param[in,out] fplog  Log file, can be nullptr.
 183          */
 184         void printInitializationToLog(FILE *fplog) const;
 185
 186         /*! \brief
 187          * Evolves the bias at every step.
 188          *
 189          * At each step the bias step needs to:
 190          * - set the bias force and potential;
 191          * - update the free energy and bias if needed;
 192          * - reweight samples to extract the PMF.
 193          *
 194          * \param[in]     coordValue     The current coordinate value(s).
 195          * \param[out]    awhPotential   Bias potential.
 196          * \param[out]    potentialJump  Change in bias potential for this bias.
 197          * \param[in]     commRecord     Struct for intra-simulation communication.
 198          * \param[in]     ms             Struct for multi-simulation communication.
 199          * \param[in]     t              Time.
 200          * \param[in]     step           Time step.
 201          * \param[in]     seed           Random seed.
 202          * \param[in,out] fplog          Log file.
 203          * \returns a reference to the bias force, size \ref ndim(), valid until the next call of this method or destruction of Bias, whichever comes first.
 204          */
 205         gmx::ArrayRef<const double>
 206         calcForceAndUpdateBias(const awh_dvec        coordValue,
 207                                double               *awhPotential,
 208                                double               *potentialJump,
 209                                const t_commrec      *commRecord,
 210                                const gmx_multisim_t *ms,
 211                                double                t,
 212                                int64_t               step,
 213                                int64_t               seed,
 214                                FILE                 *fplog);
 215
 216         /*! \brief
 217          * Calculates the convolved bias for a given coordinate value.
 218          *
 219          * The convolved bias is the effective bias acting on the coordinate.
 220          * Since the bias here has arbitrary normalization, this only makes
 221          * sense as a relative, to other coordinate values, measure of the bias.
 222          *
 223          * \param[in] coordValue  The coordinate value.
 224          * \returns the convolved bias >= -GMX_FLOAT_MAX.
 225          */
 226         double calcConvolvedBias(const awh_dvec &coordValue) const
 227         {
 228             return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
 229         }
 230
 231         /*! \brief
 232          * Restore the bias state from history on the master rank and broadcast it.
 233          *
 234          * \param[in] biasHistory  Bias history struct, only allowed to be nullptr on non-master ranks.
 235          * \param[in] cr           The communication record.
 236          */
 237         void restoreStateFromHistory(const AwhBiasHistory *biasHistory,
 238                                      const t_commrec      *cr);
 239
 240         /*! \brief
 241          * Allocate and initialize a bias history with the given bias state.
 242          *
 243          * This function will be called at the start of a new simulation.
 244          * Note that only constant data will be initialized here.
 245          * History data is set by \ref updateHistory.
 246          *
 247          * \param[in,out] biasHistory  AWH history to initialize.
 248          */
 249         void initHistoryFromState(AwhBiasHistory *biasHistory) const;
 250
 251         /*! \brief
 252          * Update the bias history with the current state.
 253          *
 254          * \param[out] biasHistory  Bias history struct.
 255          */
 256         void updateHistory(AwhBiasHistory *biasHistory) const;
 257
 258         /*! \brief
 259          * Do all previously skipped updates.
 260          * Public for use by tests.
 261          */
 262         void doSkippedUpdatesForAllPoints();
 263
 264         //! Returns the dimensionality of the bias.
 265         inline int ndim() const
 266         {
 267             return dimParams_.size();
 268         }
 269
 270         /*! \brief Returns the dimension parameters.
 271          */
 272         inline const std::vector<DimParams> &dimParams() const
 273         {
 274             return dimParams_;
 275         }
 276
 277         //! Returns the bias parameters
 278         inline const BiasParams &params() const
 279         {
 280             return params_;
 281         }
 282
 283         //! Returns the global state of the bias.
 284         inline const BiasState &state() const
 285         {
 286             return state_;
 287         }
 288
 289         //! Returns the index of the bias.
 290         inline int biasIndex() const
 291         {
 292             return params_.biasIndex;
 293         }
 294
 295         /*! \brief Return the coordinate value for a grid point.
 296          *
 297          * \param[in] gridPointIndex  The index of the grid point.
 298          */
 299         inline const awh_dvec &getGridCoordValue(size_t gridPointIndex) const
 300         {
 301             GMX_ASSERT(gridPointIndex < grid_.numPoints(), "gridPointIndex should be in the range of the grid");
 302
 303             return grid_.point(gridPointIndex).coordValue;
 304         }
 305
 306     private:
 307         /*! \brief
 308          * Performs statistical checks on the collected histograms and warns if issues are detected.
 309          *
 310          * \param[in]     t        Time.
 311          * \param[in]     step     Time step.
 312          * \param[in,out] fplog    Output file for warnings.
 313          */
 314         void warnForHistogramAnomalies(double       t,
 315                                        int64_t      step,
 316                                        FILE        *fplog);
 317
 318         /*! \brief
 319          * Collect samples for the force correlation analysis on the grid.
 320          *
 321          * \param[in] probWeightNeighbor  Probability weight of the neighboring points.
 322          * \param[in] t                   The time.
 323          */
 324         void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
 325                                         double                      t);
 326
 327     public:
 328         /*! \brief Return a const reference to the force correlation grid.
 329          */
 330         const CorrelationGrid &forceCorrelationGrid() const
 331         {
 332             GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr, "forceCorrelationGrid() should only be called with a valid force correlation object");
 333
 334             return *forceCorrelationGrid_;
 335         }
 336
 337         /*! \brief Return the number of data blocks that have been prepared for writing.
 338          */
 339         int numEnergySubblocksToWrite() const;
 340
 341         /*! \brief Write bias data blocks to energy subblocks.
 342          *
 343          * \param[in,out] subblock  Energy subblocks to write to.
 344          * \returns the number of subblocks written.
 345          */
 346         int writeToEnergySubblocks(t_enxsubblock *subblock) const;
 347
 348         /* Data members. */
 349     private:
 350         const std::vector<DimParams> dimParams_;         /**< Parameters for each dimension. */
 351         const Grid                   grid_;              /**< The multidimensional grid organizing the coordinate point locations. */
 352
 353         const BiasParams             params_;            /**< Constant parameters for the method. */
 354
 355         BiasState                    state_;             /**< The state, both global and of the grid points */
 356         std::vector<int>             updateList_;        /**< List of points for update for temporary use (could be made another tempWorkSpace) */
 357
 358         const bool                   thisRankDoesIO_;    /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
 359
 360         std::vector<double>          biasForce_;         /**< Vector for returning the force to the caller. */
 361
 362         /* Force correlation grid */
 363         std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
 364
 365         /* I/O */
 366         std::unique_ptr<BiasWriter>  writer_;      /**< Takes care of AWH data output. */
 367
 368         /* Temporary working vectors used during the update.
 369          * These are only here to avoid allocation at every MD step.
 370          */
 371         std::vector < double, AlignedAllocator < double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
 372         std::vector<double>   tempForce_;                                        /**< Bias force work buffer. */
 373
 374         /* Run-local counter to avoid flooding log with warnings. */
 375         int                          numWarningsIssued_; /**< The number of warning issued in the current run. */
 376 };
 377
 378 }      // namespace gmx
 379
 380 #endif /* GMX_AWH_BIAS_H */