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