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,2020, 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 "biasgrid.h"
  63 #include "biasparams.h"
  64 #include "biasstate.h"
  65 #include "biaswriter.h"
  66 #include "dimparams.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 BiasGrid;
  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),
 163      * normally (only) the master rank does I/O.
 164      * \param[in] disableUpdateSkips     If to disable update skips, useful for testing.
 165      */
 166     Bias(int                            biasIndexInCollection,
 167          const AwhParams&               awhParams,
 168          const AwhBiasParams&           awhBiasParams,
 169          const std::vector<DimParams>&  dimParams,
 170          double                         beta,
 171          double                         mdTimeStep,
 172          int                            numSharingSimulations,
 173          const std::string&             biasInitFilename,
 174          ThisRankWillDoIO               thisRankWillDoIO,
 175          BiasParams::DisableUpdateSkips disableUpdateSkips = BiasParams::DisableUpdateSkips::no);
 176
 177     /*! \brief
 178      * Print information about initialization to log file.
 179      *
 180      * Prints information about AWH variables that are set internally
 181      * but might be of interest to the user.
 182      *
 183      * \param[in,out] fplog  Log file, can be nullptr.
 184      */
 185     void printInitializationToLog(FILE* fplog) const;
 186
 187     /*! \brief
 188      * Evolves the bias at every step.
 189      *
 190      * At each step the bias step needs to:
 191      * - set the bias force and potential;
 192      * - update the free energy and bias if needed;
 193      * - reweight samples to extract the PMF.
 194      *
 195      * \param[in]     coordValue     The current coordinate value(s).
 196      * \param[in]     neighborLambdaEnergies An array containing the energy of the system
 197      * in neighboring lambdas. The array is of length numLambdas+1, where numLambdas is
 198      * the number of free energy lambda states. Element 0 in the array is the energy
 199      * of the current state and elements 1..numLambdas contain the energy of the system in the
 200      * neighboring lambda states (also including the current state). When there are no free
 201      * energy lambda state dimensions this can be empty.
 202      * \param[in]     neighborLambdaDhdl     An array containing the dHdL at the neighboring lambda
 203      * points. The array is of length numLambdas+1, where numLambdas is the number of free
 204      * energy lambda states. Element 0 in the array is the dHdL
 205      * of the current state and elements 1..numLambdas contain the dHdL of the system in the
 206      * neighboring lambda states (also including the current state). When there are no free
 207      * energy lambda state dimensions this can be empty.
 208      * \param[out]    awhPotential   Bias potential.
 209      * \param[out]    potentialJump  Change in bias potential for this bias.
 210      * \param[in]     commRecord     Struct for intra-simulation communication.
 211      * \param[in]     ms             Struct for multi-simulation communication.
 212      * \param[in]     t              Time.
 213      * \param[in]     step           Time step.
 214      * \param[in]     seed           Random seed.
 215      * \param[in,out] fplog          Log file.
 216      * \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.
 217      */
 218     gmx::ArrayRef<const double> calcForceAndUpdateBias(const awh_dvec         coordValue,
 219                                                        ArrayRef<const double> neighborLambdaEnergies,
 220                                                        ArrayRef<const double> neighborLambdaDhdl,
 221                                                        double*                awhPotential,
 222                                                        double*                potentialJump,
 223                                                        const t_commrec*       commRecord,
 224                                                        const gmx_multisim_t*  ms,
 225                                                        double                 t,
 226                                                        int64_t                step,
 227                                                        int64_t                seed,
 228                                                        FILE*                  fplog);
 229
 230     /*! \brief
 231      * Calculates the convolved bias for a given coordinate value.
 232      *
 233      * The convolved bias is the effective bias acting on the coordinate.
 234      * Since the bias here has arbitrary normalization, this only makes
 235      * sense as a relative, to other coordinate values, measure of the bias.
 236      *
 237      * \param[in] coordValue  The coordinate value.
 238      * \returns the convolved bias >= -GMX_FLOAT_MAX.
 239      */
 240     double calcConvolvedBias(const awh_dvec& coordValue) const
 241     {
 242         return state_.calcConvolvedBias(dimParams_, grid_, coordValue);
 243     }
 244
 245     /*! \brief
 246      * Restore the bias state from history on the master rank and broadcast it.
 247      *
 248      * \param[in] biasHistory  Bias history struct, only allowed to be nullptr on non-master ranks.
 249      * \param[in] cr           The communication record.
 250      */
 251     void restoreStateFromHistory(const AwhBiasHistory* biasHistory, const t_commrec* cr);
 252
 253     /*! \brief
 254      * Allocate and initialize a bias history with the given bias state.
 255      *
 256      * This function will be called at the start of a new simulation.
 257      * Note that only constant data will be initialized here.
 258      * History data is set by \ref updateHistory.
 259      *
 260      * \param[in,out] biasHistory  AWH history to initialize.
 261      */
 262     void initHistoryFromState(AwhBiasHistory* biasHistory) const;
 263
 264     /*! \brief
 265      * Update the bias history with the current state.
 266      *
 267      * \param[out] biasHistory  Bias history struct.
 268      */
 269     void updateHistory(AwhBiasHistory* biasHistory) const;
 270
 271     /*! \brief
 272      * Do all previously skipped updates.
 273      * Public for use by tests.
 274      */
 275     void doSkippedUpdatesForAllPoints();
 276
 277     //! Returns the dimensionality of the bias.
 278     inline int ndim() const { return dimParams_.size(); }
 279
 280     /*! \brief Returns the dimension parameters.
 281      */
 282     inline const std::vector<DimParams>& dimParams() const { return dimParams_; }
 283
 284     //! Returns the bias parameters
 285     inline const BiasParams& params() const { return params_; }
 286
 287     //! Returns the global state of the bias.
 288     inline const BiasState& state() const { return state_; }
 289
 290     //! Returns the index of the bias.
 291     inline int biasIndex() const { return params_.biasIndex; }
 292
 293     /*! \brief Return the coordinate value for a grid point.
 294      *
 295      * \param[in] gridPointIndex  The index of the grid point.
 296      */
 297     inline const awh_dvec& getGridCoordValue(size_t gridPointIndex) const
 298     {
 299         GMX_ASSERT(gridPointIndex < grid_.numPoints(),
 300                    "gridPointIndex should be in the range of the grid");
 301
 302         return grid_.point(gridPointIndex).coordValue;
 303     }
 304
 305 private:
 306     /*! \brief
 307      * Performs statistical checks on the collected histograms and warns if issues are detected.
 308      *
 309      * \param[in]     t        Time.
 310      * \param[in]     step     Time step.
 311      * \param[in,out] fplog    Output file for warnings.
 312      */
 313     void warnForHistogramAnomalies(double t, int64_t step, FILE* fplog);
 314
 315     /*! \brief
 316      * Collect samples for the force correlation analysis on the grid.
 317      *
 318      * \param[in] probWeightNeighbor  Probability weight of the neighboring points.
 319      * \param[in] neighborLambdaDhdl  An array containing the dHdL at the neighboring lambda
 320      * points. The array is of length numLambdas+1, where numLambdas is the number of free
 321      * energy lambda states. Element 0 in the array is the dHdL
 322      * of the current state and elements 1..numLambdas contain the dHdL of the system in the
 323      * neighboring lambda states (also including the current state). When there are no free
 324      * energy lambda state dimensions this can be empty.
 325      * \param[in] t                   The time.
 326      */
 327     void updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
 328                                     ArrayRef<const double>      neighborLambdaDhdl,
 329                                     double                      t);
 330
 331 public:
 332     /*! \brief Return a const reference to the force correlation grid.
 333      */
 334     const CorrelationGrid& forceCorrelationGrid() const
 335     {
 336         GMX_RELEASE_ASSERT(forceCorrelationGrid_ != nullptr,
 337                            "forceCorrelationGrid() should only be called with a valid force "
 338                            "correlation object");
 339
 340         return *forceCorrelationGrid_;
 341     }
 342
 343     /*! \brief Return the number of data blocks that have been prepared for writing.
 344      */
 345     int numEnergySubblocksToWrite() const;
 346
 347     /*! \brief Write bias data blocks to energy subblocks.
 348      *
 349      * \param[in,out] subblock  Energy subblocks to write to.
 350      * \returns the number of subblocks written.
 351      */
 352     int writeToEnergySubblocks(t_enxsubblock* subblock) const;
 353
 354     /*! \brief Returns true if the bias has a free energy lambda state dimension at all.
 355      */
 356     bool hasFepLambdaDimension() const
 357     {
 358         return std::any_of(std::begin(dimParams_), std::end(dimParams_),
 359                            [](const auto& dimParam) { return dimParam.isFepLambdaDimension(); });
 360     }
 361
 362     /*! \brief Returns whether the specified dimension is a free energy lambda
 363      * state dimension.
 364      *
 365      * \param[in] dim      The dimension to check.
 366      *
 367      * \returns true if the specified dimension is a free energy lambda state dimension.
 368      */
 369     bool isFepLambdaDimension(int dim) const;
 370
 371     /*! \brief
 372      * Returns whether we should sample the coordinate.
 373      *
 374      * \param[in] step  The MD step number.
 375      */
 376     bool isSampleCoordStep(int64_t step) const;
 377
 378     /* Data members. */
 379 private:
 380     const std::vector<DimParams> dimParams_; /**< Parameters for each dimension. */
 381     const BiasGrid grid_; /**< The multidimensional grid organizing the coordinate point locations. */
 382
 383     const BiasParams params_; /**< Constant parameters for the method. */
 384
 385     BiasState        state_; /**< The state, both global and of the grid points */
 386     std::vector<int> updateList_; /**< List of points for update for temporary use (could be made another tempWorkSpace) */
 387
 388     const bool thisRankDoesIO_; /**< Tells whether this MPI rank will do I/O (checkpointing, AWH output) */
 389
 390     std::vector<double> biasForce_; /**< Vector for returning the force to the caller. */
 391
 392     /* Force correlation grid */
 393     std::unique_ptr<CorrelationGrid> forceCorrelationGrid_; /**< Takes care of force correlation statistics for every grid point. */
 394
 395     /* I/O */
 396     std::unique_ptr<BiasWriter> writer_; /**< Takes care of AWH data output. */
 397
 398     /* Temporary working vectors used during the update.
 399      * These are only here to avoid allocation at every MD step.
 400      */
 401     std::vector<double, AlignedAllocator<double>> alignedTempWorkSpace_; /**< Working vector of doubles. */
 402     std::vector<double>                           tempForce_; /**< Bias force work buffer. */
 403
 404     /* Run-local counter to avoid flooding log with warnings. */
 405     int numWarningsIssued_; /**< The number of warning issued in the current run. */
 406 };
 407
 408 } // namespace gmx
 409
 410 #endif /* GMX_AWH_BIAS_H */