src/gromacs/awh/grid.h

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