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 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> calcForceAndUpdateBias(const awh_dvec        coordValue,
 206                                                        double*               awhPotential,
 207                                                        double*               potentialJump,
 208                                                        const t_commrec*      commRecord,
 209                                                        const gmx_multisim_t* ms,
 210                                                        double                t,
 211                                                        int64_t               step,
 212                                                        int64_t               seed,
 213                                                        FILE*                 fplog);
 214
 215     /*! \brief
 216      * Calculates the convolved bias for a given coordinate value.
 217      *
 218      * The convolved bias is the effective bias acting on the coordinate.
 219      * Since the bias here has arbitrary normalization, this only makes
 220      * sense as a relative, to other coordinate values, measure of the bias.
 221      *
 222      * \param[in] coordValue  The coordinate value.
 223      * \returns the convolved bias >= -GMX_FLOAT_MAX.
 224      */
 225     double calcConvolvedBias(const awh_dvec& coordValue) const
 226     {
 227         return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
 228     }
 229
 230     /*! \brief
 231      * Restore the bias state from history on the master rank and broadcast it.
 232      *
 233      * \param[in] biasHistory  Bias history struct, only allowed to be nullptr on non-master ranks.
 234      * \param[in] cr           The communication record.
 235      */
 236     void restoreStateFromHistory(const AwhBiasHistory* biasHistory, const t_commrec* cr);
 237
 238     /*! \brief
 239      * Allocate and initialize a bias history with the given bias state.
 240      *
 241      * This function will be called at the start of a new simulation.
 242      * Note that only constant data will be initialized here.
 243      * History data is set by \ref updateHistory.
 244      *
 245      * \param[in,out] biasHistory  AWH history to initialize.
 246      */
 247     void initHistoryFromState(AwhBiasHistory* biasHistory) const;
 248
 249     /*! \brief
 250      * Update the bias history with the current state.
 251      *
 252      * \param[out] biasHistory  Bias history struct.
 253      */
 254     void updateHistory(AwhBiasHistory* biasHistory) const;
 255
 256     /*! \brief
 257      * Do all previously skipped updates.
 258      * Public for use by tests.
 259      */
 260     void doSkippedUpdatesForAllPoints();
 261
 262     //! Returns the dimensionality of the bias.
 263     inline int ndim() const { return dimParams_.size(); }
 264
 265     /*! \brief Returns the dimension parameters.
 266      */
 267     inline const std::vector<DimParams>& dimParams() const { return dimParams_; }
 268
 269     //! Returns the bias parameters
 270     inline const BiasParams& params() const { return params_; }
 271
 272     //! Returns the global state of the bias.
 273     inline const BiasState& state() const { return state_; }
 274
 275     //! Returns the index of the bias.
 276     inline int biasIndex() const { return params_.biasIndex; }
 277
 278     /*! \brief Return the coordinate value for a grid point.
 279      *
 280      * \param[in] gridPointIndex  The index of the grid point.
 281      */
 282     inline const awh_dvec& getGridCoordValue(size_t gridPointIndex) const
 283     {
 284         GMX_ASSERT(gridPointIndex < grid_.numPoints(),
 285                    "gridPointIndex should be in the range of the grid");
 286
 287         return grid_.point(gridPointIndex).coordValue;
 288     }
 289
 290 private:
 291     /*! \brief
 292      * Performs statistical checks on the collected histograms and warns if issues are detected.
 293      *
 294      * \param[in]     t        Time.
 295      * \param[in]     step     Time step.
 296      * \param[in,out] fplog    Output file for warnings.
 297      */
 298     void warnForHistogramAnomalies(double t, int64_t step, FILE* fplog);
 299
 300     /*! \brief
 301      * Collect samples for the force correlation analysis on the grid.
 302      *
 303      * \param[in] probWeightNeighbor  Probability weight of the neighboring points.
 304      * \param[in] t                   The time.
 305      */
 306     void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor, double t);
 307
 308 public:
 309     /*! \brief Return a const reference to the force correlation grid.
 310      */
 311     const CorrelationGrid& forceCorrelationGrid() const
 312     {
 313         GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr,
 314                            "forceCorrelationGrid() should only be called with a valid force "
 315                            "correlation object");
 316
 317         return *forceCorrelationGrid_;
 318     }
 319
 320     /*! \brief Return the number of data blocks that have been prepared for writing.
 321      */
 322     int numEnergySubblocksToWrite() const;
 323
 324     /*! \brief Write bias data blocks to energy subblocks.
 325      *
 326      * \param[in,out] subblock  Energy subblocks to write to.
 327      * \returns the number of subblocks written.
 328      */
 329     int writeToEnergySubblocks(t_enxsubblock* subblock) const;
 330
 331     /* Data members. */
 332 private:
 333     const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
 334     const Grid grid_; /**< The multidimensional grid organizing the coordinate point locations. */
 335
 336     const BiasParams params_; /**< Constant parameters for the method. */
 337
 338     BiasState        state_; /**< The state, both global and of the grid points */
 339     std::vector<int> updateList_; /**< List of points for update for temporary use (could be made another tempWorkSpace) */
 340
 341     const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
 342
 343     std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
 344
 345     /* Force correlation grid */
 346     std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
 347
 348     /* I/O */
 349     std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
 350
 351     /* Temporary working vectors used during the update.
 352      * These are only here to avoid allocation at every MD step.
 353      */
 354     std::vector<double, AlignedAllocator<double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
 355     std::vector<double>                           tempForce_; /**< Bias force work buffer. */
 356
 357     /* Run-local counter to avoid flooding log with warnings. */
 358     int numWarningsIssued_; /**< The number of warning issued in the current run. */
 359 };
 360
 361 } // namespace gmx
 362
 363 #endif /* GMX_AWH_BIAS_H */