src/gromacs/applied_forces/awh/biasgrid.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2019,2020,2021, 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  * This file contains datatypes and function declarations necessary
  40  * for AWH to interface with the grid code.
  41  *
  42  * The grid organizes spatial properties of the AWH coordinate points.
  43  * This includes traversing points in a specific order, locating
  44  * neighboring points and calculating distances. Multiple dimensions
  45  * as well as periodic dimensions are supported.
  46  *
  47  * \todo: Replace this by a more generic grid class once that is available.
  48  *
  49  * \author Viveca Lindahl
  50  * \author Berk Hess <hess@kth.se>
  51  * \ingroup module_awh
  52  */
  53
  54 #ifndef GMX_AWH_BIASGRID_H
  55 #define GMX_AWH_BIASGRID_H
  56
  57 #include <memory>
  58 #include <optional>
  59 #include <string>
  60
  61 #include "gromacs/utility/arrayref.h"
  62
  63 #include "dimparams.h" /* This is needed for awh_dvec */
  64
  65 namespace gmx
  66 {
  67
  68 struct AwhDimParams;
  69
  70 /*! \internal
  71  * \brief An axis, i.e. dimension, of the grid.
  72  */
  73 class GridAxis
  74 {
  75 public:
  76     /*! \brief Constructor.
  77      *
  78      * The spacing and number of points are set such that we have
  79      * at least the requested point density.
  80      * Requesting 0 point density results in the minimum number
  81      * of points (2).
  82      *
  83      * \param[in] origin           Starting value.
  84      * \param[in] end              End value.
  85      * \param[in] period           Period, pass 0 if not periodic.
  86      * \param[in] pointDensity     Requested number of point per unit of axis length.
  87      */
  88     GridAxis(double origin, double end, double period, double pointDensity);
  89
  90     /*! \brief Constructor.
  91      *
  92      * \param[in] origin           Starting value.
  93      * \param[in] end              End value.
  94      * \param[in] period           Period, pass 0 if not periodic.
  95      * \param[in] numPoints        The number of points.
  96      * \param[in] isFepLambdaAxis     If this axis is controlling lambda.
  97      */
  98     GridAxis(double origin, double end, double period, int numPoints, bool isFepLambdaAxis);
  99
 100     /*! \brief Returns whether the axis has periodic boundaries.
 101      */
 102     bool isPeriodic() const { return period_ > 0; }
 103
 104     /*! \brief Returns the period of the grid along the axis.
 105      */
 106     double period() const { return period_; }
 107
 108     /*! \brief Returns the grid origin along the axis.
 109      */
 110     double origin() const { return origin_; }
 111
 112     /*! \brief Returns the grid point spacing along the axis.
 113      */
 114     double spacing() const { return spacing_; }
 115
 116     /*! \brief Returns the number of grid points along the axis.
 117      */
 118     int numPoints() const { return numPoints_; }
 119
 120     /*! \brief Returns the period of the grid in points along the axis.
 121      *
 122      * Returns 0 if the axis is not periodic.
 123      */
 124     int numPointsInPeriod() const { return numPointsInPeriod_; }
 125
 126     /*! \brief Returns the length of the interval.
 127      *
 128      * This returns the distance obtained by connecting the origin point to
 129      * the end point in the positive direction. Note that this is generally
 130      * not the shortest distance. For period > 0, both origin and
 131      * end are expected to take values in the same periodic interval.
 132      */
 133     double length() const { return length_; }
 134
 135     /*! \brief Map a value to the nearest point index along an axis.
 136      *
 137      * \param[in] value  Value along the axis.
 138      * \returns the index nearest to the value.
 139      */
 140     int nearestIndex(double value) const;
 141
 142     /*! \brief Returns whether this axis is coupled to the free energy lambda state.
 143      */
 144     bool isFepLambdaAxis() const { return isFepLambdaAxis_; }
 145
 146 private:
 147     double origin_;            /**< Interval start value */
 148     double length_;            /**< Interval length */
 149     double period_;            /**< The period, 0 if not periodic */
 150     double spacing_;           /**< Point spacing */
 151     int    numPoints_;         /**< Number of points in the interval */
 152     int    numPointsInPeriod_; /**< Number of points in a period (0 if no periodicity) */
 153     bool isFepLambdaAxis_; /**< If this axis is coupled to the system's free energy lambda state */
 154 };
 155
 156 /*! \internal
 157  * \brief A point in the grid.
 158  *
 159  * A grid point has a coordinate value and a coordinate index of the same dimensionality as the
 160  * grid. It knows the linear indices of its neighboring point (which are useful only when handed up
 161  * to the grid).
 162  */
 163 struct GridPoint
 164 {
 165     awh_dvec         coordValue; /**< Multidimensional coordinate value of this point */
 166     awh_ivec         index;      /**< Multidimensional point indices */
 167     std::vector<int> neighbor;   /**< Linear point indices of the neighboring points */
 168 };
 169
 170 /*! \internal
 171  * \brief The grid for a single bias, generally multidimensional and periodic.
 172  *
 173  * The grid discretizes a multidimensional space with some given resolution.
 174  * Each dimension is represented by an axis which sets the spatial extent,
 175  * point spacing and periodicity of the grid in that direction.
 176  */
 177 class BiasGrid
 178 {
 179 private:
 180     /*! \brief Initializes the grid points.
 181      */
 182     void initPoints();
 183
 184 public:
 185     /*! \brief
 186      * The point density per sigma of the Gaussian distribution in an umbrella.
 187      *
 188      * This value should be at least 1 to uniformly cover the reaction coordinate
 189      * range with density and having it larger than 1 does not add information.
 190      */
 191     static constexpr double c_numPointsPerSigma = 1.0;
 192
 193     //! Cut-off in sigma for considering points, neglects 4e-8 of the density.
 194     static constexpr double c_scopeCutoff = 5.5;
 195
 196     /*! \brief Construct a grid using AWH input parameters.
 197      *
 198      * \param[in] dimParams     Dimension parameters including the expected inverse variance of the
 199      * coordinate living on the grid (determines the grid spacing).
 200      * \param[in] awhDimParams  Dimension params from inputrec.
 201      */
 202     BiasGrid(ArrayRef<const DimParams> dimParams, const AwhDimParams* awhDimParams);
 203
 204     /*! \brief Returns the number of points in the grid.
 205      *
 206      * \returns the number of points in the grid.
 207      */
 208     size_t numPoints() const { return point_.size(); }
 209
 210     /*! \brief Returns a reference to a point on the grid.
 211      *
 212      * \returns a constant reference to a point on the grid.
 213      */
 214     const GridPoint& point(size_t pointIndex) const { return point_[pointIndex]; }
 215
 216     /*! \brief Returns the dimensionality of the grid.
 217      *
 218      * \returns the dimensionality of the grid.
 219      */
 220     int numDimensions() const { return axis_.size(); }
 221
 222     /*! \brief Returns the grid axes.
 223      *
 224      * \returns a constant reference to the grid axes.
 225      */
 226     ArrayRef<const GridAxis> axis() const { return axis_; }
 227
 228     /*! \brief Returns a grid axis.
 229      *
 230      * param[in] dim  Dimension to return the grid axis for.
 231      * \returns a constant reference to the grid axis.
 232      */
 233     const GridAxis& axis(int dim) const { return axis_[dim]; }
 234
 235     /*! \brief Find the grid point with value nearest to the given value.
 236      *
 237      * \param[in] value  Value vector.
 238      * \returns the grid point index.
 239      */
 240     int nearestIndex(const awh_dvec value) const;
 241
 242     /*! \brief Query if the value is in the grid.
 243      *
 244      * \param[in] value  Value vector.
 245      * \returns true if the value is in the grid.
 246      * \note It is assumed that any periodicity of value has already been taken care of.
 247      */
 248     bool covers(const awh_dvec value) const;
 249
 250     /*! \brief Returns true if the grid has a free energy lambda state axis at all.
 251      */
 252     bool hasLambdaAxis() const
 253     {
 254         return std::any_of(std::begin(axis_), std::end(axis_), [](const auto& axis) {
 255             return axis.isFepLambdaAxis();
 256         });
 257     }
 258
 259     /*! \brief
 260      * Returns the index of a free energy lambda state axis (there can be
 261      * no more than one) if there is one.
 262      */
 263     std::optional<int> lambdaAxisIndex() const;
 264
 265     /*! \brief
 266      * Returns the number of free energy lambda states in the grid (the number
 267      * of points along a free energy lambda state axis) or 0 if there are no free energy
 268      * lambda state axes.
 269      */
 270     int numFepLambdaStates() const;
 271
 272 private:
 273     std::vector<GridPoint> point_; /**< Points on the grid */
 274     std::vector<GridAxis>  axis_;  /**< Axes, one for each dimension. */
 275 };
 276
 277 /*! \endcond */
 278
 279 /*! \brief Convert a multidimensional grid point index to a linear one.
 280  *
 281  * \param[in] grid        The grid.
 282  * \param[in] indexMulti  Multidimensional grid point index to convert to a linear one.
 283  * \returns the linear index.
 284  */
 285 int multiDimGridIndexToLinear(const BiasGrid& grid, const awh_ivec indexMulti);
 286
 287 /*! \brief Convert multidimensional array index to a linear one.
 288  *
 289  * \param[in] indexMulti    Multidimensional index to convert to a linear one.
 290  * \param[in] numDim        Number of dimensions of the array.
 291  * \param[in] numPointsDim  Number of points of the array.
 292  * \returns the linear index.
 293  * \note This function can be used without having an initialized grid.
 294  */
 295 int multiDimArrayIndexToLinear(const awh_ivec indexMulti, int numDim, const awh_ivec numPointsDim);
 296
 297 /*! \brief Convert a linear grid point index to a multidimensional one.
 298  *
 299  * \param[in]  grid         The grid.
 300  * \param[in]  indexLinear  Linear grid point index to convert to a multidimensional one.
 301  * \param[out] indexMulti   The multidimensional index.
 302  */
 303 void linearGridindexToMultiDim(const BiasGrid& grid, int indexLinear, awh_ivec indexMulti);
 304
 305 /*! \brief Convert a linear array index to a multidimensional one.
 306  *
 307  * \param[in]  indexLinear   Linear array index
 308  * \param[in]  ndim          Number of dimensions of the array.
 309  * \param[in]  numPointsDim  Number of points for each dimension.
 310  * \param[out] indexMulti    The multidimensional index.
 311  */
 312 void linearArrayIndexToMultiDim(int indexLinear, int ndim, const awh_ivec numPointsDim, awh_ivec indexMulti);
 313
 314 /*! \brief
 315  * Find the next grid point in the sub-part of the grid given a starting point.
 316  *
 317  * The given grid point index is updated to the next valid grid point index
 318  * by traversing the sub-part of the grid, here termed the subgrid.
 319  * Since the subgrid range might extend beyond the actual size of the grid,
 320  * the subgrid is traversed until a point both in the subgrid and grid is
 321  * found. If no point is found, the function returns false and the index is
 322  * not modified. The starting point needs to be inside of the subgrid. However,
 323  * if this index is not given, meaning < 0, then the search is initialized at
 324  * the subgrid origin, i.e. in this case the "next" grid point index is
 325  * defined to be the first common grid/subgrid point.
 326  *
 327  * \param[in]     grid            The grid.
 328  * \param[in]     subgridOrigin   Vector locating the subgrid origin relative to the grid origin.
 329  * \param[in]     subgridNpoints  Number of points along each subgrid dimension.
 330  * \param[in,out] gridPointIndex  Pointer to the starting/next grid point index.
 331  * \returns true if the grid point was updated.
 332  */
 333 bool advancePointInSubgrid(const BiasGrid& grid,
 334                            const awh_ivec  subgridOrigin,
 335                            const awh_ivec  subgridNpoints,
 336                            int*            gridPointIndex);
 337
 338 /*! \brief Maps each point in the grid to a point in the data grid.
 339  *
 340  * This functions maps an AWH bias grid to a user provided input data grid.
 341  * The value of data grid point i along dimension d is given by data[d][i].
 342  * The number of dimensions of the data should equal that of the grid.
 343  * A fatal error is thrown if extracting the data fails or the data does not cover the whole grid.
 344  *
 345  * \param[out] gridpointToDatapoint  Array mapping each grid point to a data point index.
 346  * \param[in]  data                  2D array in format ndim x ndatapoints with data grid point values.
 347  * \param[in]  numDataPoints         Number of data points.
 348  * \param[in]  dataFilename          The data filename.
 349  * \param[in]  grid                  The grid.
 350  * \param[in]  correctFormatMessage  String to include in error message if extracting the data fails.
 351  */
 352 void mapGridToDataGrid(std::vector<int>*    gridpointToDatapoint,
 353                        const double* const* data,
 354                        int                  numDataPoints,
 355                        const std::string&   dataFilename,
 356                        const BiasGrid&      grid,
 357                        const std::string&   correctFormatMessage);
 358
 359 /*! \brief
 360  * Get the deviation along one dimension from the given value to a point in the grid.
 361  *
 362  * \param[in] grid        The grid.
 363  * \param[in] dimIndex    Dimensional index in [0, ndim -1].
 364  * \param[in] pointIndex  BiasGrid point index.
 365  * \param[in] value       Value along the given dimension.
 366  * \returns the deviation of the given value to the given point.
 367  */
 368 double getDeviationFromPointAlongGridAxis(const BiasGrid& grid, int dimIndex, int pointIndex, double value);
 369
 370 /*! \brief
 371  * Get the deviation from one point to another along one dimension in the grid.
 372  *
 373  * \param[in] grid        The grid.
 374  * \param[in] dimIndex    Dimensional index in [0, ndim -1].
 375  * \param[in] pointIndex1 Grid point index of the first point.
 376  * \param[in] pointIndex2 Grid point index of the second point.
 377  * \returns the deviation of the two points along the given axis.
 378  */
 379 double getDeviationFromPointAlongGridAxis(const BiasGrid& grid, int dimIndex, int pointIndex1, int pointIndex2);
 380
 381 /*! \brief
 382  * Checks whether two points are along a free energy lambda state axis.
 383  *
 384  * \param[in] grid        The grid.
 385  * \param[in] pointIndex1 Grid point index of the first point.
 386  * \param[in] pointIndex2 Grid point index of the second point.
 387  * \returns true if the two points are along a free energy lambda state axis.
 388  */
 389 bool pointsAlongLambdaAxis(const BiasGrid& grid, int pointIndex1, int pointIndex2);
 390
 391 /*! \brief
 392  * Checks whether two points are different in the free energy lambda state
 393  * dimension (if any).
 394  *
 395  * \param[in] grid        The grid.
 396  * \param[in] pointIndex1 Grid point index of the first point.
 397  * \param[in] pointIndex2 Grid point index of the second point.
 398  * \returns true if the two points have different lambda values.
 399  */
 400 bool pointsHaveDifferentLambda(const BiasGrid& grid, int pointIndex1, int pointIndex2);
 401
 402 } // namespace gmx
 403
 404 #endif