src/gromacs/awh/grid.cpp

   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  * \brief
  38  * Implements functions in grid.h.
  39  *
  40  * \author Viveca Lindahl
  41  * \author Berk Hess <hess@kth.se>
  42  * \ingroup module_awh
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "grid.h"
  48
  49 #include <cassert>
  50 #include <cmath>
  51 #include <cstring>
  52
  53 #include <algorithm>
  54
  55 #include "gromacs/math/functions.h"
  56 #include "gromacs/math/utilities.h"
  57 #include "gromacs/mdtypes/awh_params.h"
  58 #include "gromacs/utility/cstringutil.h"
  59 #include "gromacs/utility/exceptions.h"
  60 #include "gromacs/utility/gmxassert.h"
  61 #include "gromacs/utility/smalloc.h"
  62 #include "gromacs/utility/stringutil.h"
  63
  64 namespace gmx
  65 {
  66
  67 namespace
  68 {
  69
  70 /*! \brief
  71  * Modify x so that it is periodic in [-period/2, +period/2).
  72  *
  73  * x is modified by shifting its value by a +/- a period if
  74  * needed. Thus, it is assumed that x is at most one period
  75  * away from this interval. For period = 0, x is not modified.
  76  *
  77  * \param[in,out] x       Pointer to the value to modify.
  78  * \param[in]     period  The period, or 0 if not periodic.
  79  */
  80 void centerPeriodicValueAroundZero(double *x,
  81                                    double  period)
  82 {
  83     GMX_ASSERT(period >= 0, "Periodic should not be negative");
  84
  85     const double halfPeriod = period*0.5;
  86
  87     if (*x >= halfPeriod)
  88     {
  89         *x -= period;
  90     }
  91     else if (*x < -halfPeriod)
  92     {
  93         *x += period;
  94     }
  95 }
  96
  97 /*! \brief
  98  * If period>0, retrun x so that it is periodic in [0, period), else return x.
  99  *
 100  * Return x is shifted its value by a +/- a period, if
 101  * needed. Thus, it is assumed that x is at most one period
 102  * away from this interval. For this domain and period > 0
 103  * this is equivalent to x = x % period. For period = 0,
 104  * x is not modified.
 105  *
 106  * \param[in,out] x       Pointer to the value to modify, should be >= 0.
 107  * \param[in]     period  The period, or 0 if not periodic.
 108  * \returns for period>0: index value witin [0, period), otherwise: \p x.
 109  */
 110 int indexWithinPeriod(int x,
 111                       int period)
 112 {
 113     GMX_ASSERT(period >= 0, "Periodic should not be negative");
 114
 115     if (period == 0)
 116     {
 117         return x;
 118     }
 119
 120     GMX_ASSERT(x > -period && x < 2*period, "x should not be more shifted by more than one period");
 121
 122     if (x >= period)
 123     {
 124         return x - period;
 125     }
 126     else if (x < 0)
 127     {
 128         return x + period;
 129     }
 130     else
 131     {
 132         return x;
 133     }
 134 }
 135
 136 /*! \brief
 137  * Get the length of the interval (origin, end).
 138  *
 139  * This returns the distance obtained by connecting the origin point to
 140  * the end point in the positive direction. Note that this is generally
 141  * not the shortest distance. For period > 0, both origin and
 142  * end are expected to take values in the same periodic interval,
 143  * ie. |origin - end| < period.
 144  *
 145  * \param[in] origin    Start value of the interval.
 146  * \param[in] end       End value of the interval.
 147  * \param[in] period    The period, or 0 if not periodic.
 148  * \returns the interval length from origin to end.
 149  */
 150 double getIntervalLengthPeriodic(double origin,
 151                                  double end,
 152                                  double period)
 153 {
 154     double length = end - origin;
 155     if (length < 0)
 156     {
 157         /* The interval wraps around the +/- boundary which has a discontinuous jump of -period. */
 158         length += period;
 159     }
 160
 161     GMX_RELEASE_ASSERT(length >= 0, "Negative AWH grid axis length.");
 162     GMX_RELEASE_ASSERT(period == 0 || length <= period, "Interval length longer than period.");
 163
 164     return length;
 165 }
 166
 167 /*! \brief
 168  * Get the deviation x - x0.
 169  *
 170  * For period > 0, the deviation with minimum absolute value is returned,
 171  * i.e. with a value in the interval [-period/2, +period/2).
 172  * Also for period > 0, it is assumed that |x - x0| < period.
 173  *
 174  * \param[in] x        From value.
 175  * \param[in] x0       To value.
 176  * \param[in] period   The period, or 0 if not periodic.
 177  * \returns the deviation from x to x0.
 178  */
 179 double getDeviationPeriodic(double x,
 180                             double x0,
 181                             double period)
 182 {
 183     double dev = x - x0;
 184
 185     if (period > 0)
 186     {
 187         centerPeriodicValueAroundZero(&dev, period);
 188     }
 189
 190     return dev;
 191 }
 192
 193 }   // namespace
 194
 195 double getDeviationFromPointAlongGridAxis(const Grid &grid,
 196                                           int         dimIndex,
 197                                           int         pointIndex,
 198                                           double      value)
 199 {
 200     double coordValue = grid.point(pointIndex).coordValue[dimIndex];
 201
 202     return getDeviationPeriodic(value, coordValue, grid.axis(dimIndex).period());
 203 }
 204
 205 void linearArrayIndexToMultiDim(int indexLinear, int numDimensions, const awh_ivec numPointsDim, awh_ivec indexMulti)
 206 {
 207     for (int d = 0; d < numDimensions; d++)
 208     {
 209         int stride = 1;
 210
 211         for (int k = d + 1; k < numDimensions; k++)
 212         {
 213             stride *= numPointsDim[k];
 214         }
 215
 216         indexMulti[d] = indexLinear/stride;
 217         indexLinear  -= indexMulti[d]*stride;
 218     }
 219 }
 220
 221 void linearGridindexToMultiDim(const Grid &grid,
 222                                int         indexLinear,
 223                                awh_ivec    indexMulti)
 224 {
 225     awh_ivec  numPointsDim;
 226     const int numDimensions = grid.numDimensions();
 227     for (int d = 0; d < numDimensions; d++)
 228     {
 229         numPointsDim[d] = grid.axis(d).numPoints();
 230     }
 231
 232     linearArrayIndexToMultiDim(indexLinear, numDimensions, numPointsDim, indexMulti);
 233 }
 234
 235
 236 int multiDimArrayIndexToLinear(const awh_ivec indexMulti,
 237                                int            numDimensions,
 238                                const awh_ivec numPointsDim)
 239 {
 240     int stride      = 1;
 241     int indexLinear = 0;
 242     for (int d = numDimensions - 1; d >= 0; d--)
 243     {
 244         indexLinear += stride*indexMulti[d];
 245         stride      *= numPointsDim[d];
 246     }
 247
 248     return indexLinear;
 249 }
 250
 251 namespace
 252 {
 253
 254 /*! \brief Convert a multidimensional grid point index to a linear one.
 255  *
 256  * \param[in] axis       The grid axes.
 257  * \param[in] indexMulti Multidimensional grid point index to convert to a linear one.
 258  * \returns the linear index.
 259  */
 260 int multiDimGridIndexToLinear(const std::vector<GridAxis> &axis,
 261                               const awh_ivec               indexMulti)
 262 {
 263     awh_ivec numPointsDim = { 0 };
 264
 265     for (size_t d = 0; d < axis.size(); d++)
 266     {
 267         numPointsDim[d] = axis[d].numPoints();
 268     }
 269
 270     return multiDimArrayIndexToLinear(indexMulti, axis.size(), numPointsDim);
 271 }
 272
 273 }   // namespace
 274
 275 int multiDimGridIndexToLinear(const Grid     &grid,
 276                               const awh_ivec  indexMulti)
 277 {
 278     return multiDimGridIndexToLinear(grid.axis(), indexMulti);
 279 }
 280
 281 namespace
 282 {
 283
 284 /*! \brief
 285  * Take a step in a multidimensional array.
 286  *
 287  * The multidimensional index gives the starting point to step from. Dimensions are
 288  * stepped through in order of decreasing dimensional index such that the index is
 289  * incremented in the highest dimension possible. If the starting point is the end
 290  * of the array, a step cannot be taken and the index is not modified.
 291  *
 292  * \param[in] numDim        Number of dimensions of the array.
 293  * \param[in] numPoints     Vector with the number of points along each dimension.
 294  * \param[in,out] indexDim  Multidimensional index, each with values in [0, numPoints[d] - 1].
 295  * \returns true if a step was taken, false if not.
 296  */
 297 bool stepInMultiDimArray(int            numDim,
 298                          const awh_ivec numPoints,
 299                          awh_ivec       indexDim)
 300 {
 301     bool haveStepped = false;
 302
 303     for (int d = numDim - 1; d >= 0 && !haveStepped; d--)
 304     {
 305         if (indexDim[d] < numPoints[d] - 1)
 306         {
 307             /* Not at a boundary, just increase by 1. */
 308             indexDim[d]++;
 309             haveStepped = true;
 310         }
 311         else
 312         {
 313             /* At a boundary. If we are not at the end of the array,
 314                reset the index and check if we can step in higher dimensions */
 315             if (d > 0)
 316             {
 317                 indexDim[d] = 0;
 318             }
 319         }
 320     }
 321
 322     return haveStepped;
 323 }
 324
 325 /*! \brief
 326  * Transforms a grid point index to to the multidimensional index of a subgrid.
 327  *
 328  * The subgrid is defined by the location of its origin and the number of points
 329  * along each dimension. The index transformation thus consists of a projection
 330  * of the linear index onto each dimension, followed by a translation of the origin.
 331  * The subgrid may have parts that don't overlap with the grid. E.g. the origin
 332  * vector can have negative components meaning the origin lies outside of the grid.
 333  * However, the given point needs to be both a grid and subgrid point.
 334  *
 335  * Periodic boundaries are taken care of by wrapping the subgrid around the grid.
 336  * Thus, for periodic dimensions the number of subgrid points need to be less than
 337  * the number of points in a period to prevent problems of wrapping around.
 338  *
 339  * \param[in]     grid            The grid.
 340  * \param[in]     subgridOrigin   Vector locating the subgrid origin relative to the grid origin.
 341  * \param[in]     subgridNpoints  The number of subgrid points in each dimension.
 342  * \param[in]     point           Grid point to get subgrid index for.
 343  * \param[in,out] subgridIndex    Subgrid multidimensional index.
 344  */
 345 void gridToSubgridIndex(const Grid     &grid,
 346                         const awh_ivec  subgridOrigin,
 347                         const awh_ivec  subgridNpoints,
 348                         int             point,
 349                         awh_ivec        subgridIndex)
 350 {
 351     /* Get the subgrid index of the given grid point, for each dimension. */
 352     for (int d = 0; d < grid.numDimensions(); d++)
 353     {
 354         /* The multidimensional grid point index relative to the subgrid origin. */
 355         subgridIndex[d] =
 356             indexWithinPeriod(grid.point(point).index[d] - subgridOrigin[d],
 357                               grid.axis(d).numPointsInPeriod());
 358
 359         /* The given point should be in the subgrid. */
 360         GMX_RELEASE_ASSERT((subgridIndex[d] >= 0) && (subgridIndex[d] < subgridNpoints[d]),
 361                            "Attempted to convert an AWH grid point index not in subgrid to out of bounds subgrid index");
 362     }
 363 }
 364
 365 /*! \brief
 366  * Transform a multidimensional subgrid index to a grid point index.
 367  *
 368  * If the given subgrid point is not a grid point the transformation will not be successful
 369  * and the grid point index will not be set. Periodic boundaries are taken care of by
 370  * wrapping the subgrid around the grid.
 371  *
 372  * \param[in]     grid           The grid.
 373  * \param[in]     subgridOrigin  Vector locating the subgrid origin relative to the grid origin.
 374  * \param[in]     subgridIndex   Subgrid multidimensional index to get grid point index for.
 375  * \param[in,out] gridIndex      Grid point index.
 376  * \returns true if the transformation was successful.
 377  */
 378 bool subgridToGridIndex(const Grid     &grid,
 379                         const awh_ivec  subgridOrigin,
 380                         const awh_ivec  subgridIndex,
 381                         int            *gridIndex)
 382 {
 383     awh_ivec globalIndexDim;
 384
 385     /* Check and apply boundary conditions for each dimension */
 386     for (int d = 0; d < grid.numDimensions(); d++)
 387     {
 388         /* Transform to global multidimensional indexing by adding the origin */
 389         globalIndexDim[d] = subgridOrigin[d] + subgridIndex[d];
 390
 391         /* The local grid is allowed to stick out on the edges of the global grid. Here the boundary conditions are applied.*/
 392         if (globalIndexDim[d] < 0 ||  globalIndexDim[d] > grid.axis(d).numPoints() - 1)
 393         {
 394             /* Try to wrap around if periodic. Otherwise, the transformation failed so return. */
 395             if (!grid.axis(d).isPeriodic())
 396             {
 397                 return false;
 398             }
 399
 400             /* The grid might not contain a whole period. Can only wrap around if this gap is not too large. */
 401             int gap = grid.axis(d).numPointsInPeriod() - grid.axis(d).numPoints();
 402
 403             int bridge;
 404             int numWrapped;
 405             if (globalIndexDim[d] < 0)
 406             {
 407                 bridge     = -globalIndexDim[d];
 408                 numWrapped = bridge - gap;
 409                 if (numWrapped > 0)
 410                 {
 411                     globalIndexDim[d] = grid.axis(d).numPoints() - numWrapped;
 412                 }
 413             }
 414             else
 415             {
 416                 bridge     = globalIndexDim[d] - (grid.axis(d).numPoints() - 1);
 417                 numWrapped = bridge - gap;
 418                 if (numWrapped > 0)
 419                 {
 420                     globalIndexDim[d] = numWrapped - 1;
 421                 }
 422             }
 423
 424             if (numWrapped <= 0)
 425             {
 426                 return false;
 427             }
 428         }
 429     }
 430
 431     /* Translate from multidimensional to linear indexing and set the return value */
 432     (*gridIndex) = multiDimGridIndexToLinear(grid, globalIndexDim);
 433
 434     return true;
 435 }
 436
 437 }   // namespace
 438
 439 bool advancePointInSubgrid(const Grid     &grid,
 440                            const awh_ivec  subgridOrigin,
 441                            const awh_ivec  subgridNumPoints,
 442                            int            *gridPointIndex)
 443 {
 444     /* Initialize the subgrid index to the subgrid origin. */
 445     awh_ivec subgridIndex = { 0 };
 446
 447     /* Get the subgrid index of the given grid point index. */
 448     if (*gridPointIndex >= 0)
 449     {
 450         gridToSubgridIndex(grid, subgridOrigin, subgridNumPoints, *gridPointIndex, subgridIndex);
 451     }
 452     else
 453     {
 454         /* If no grid point is given we start at the subgrid origin (which subgridIndex is initialized to).
 455            If this is a valid grid point then we're done, otherwise keep looking below. */
 456         /* TODO: separate into a separate function (?) */
 457         if (subgridToGridIndex(grid, subgridOrigin, subgridIndex, gridPointIndex))
 458         {
 459             return true;
 460         }
 461     }
 462
 463     /* Traverse the subgrid and look for the first point that is also in the grid. */
 464     while (stepInMultiDimArray(grid.numDimensions(), subgridNumPoints, subgridIndex))
 465     {
 466         /* If this is a valid grid point, the grid point index is updated.*/
 467         if (subgridToGridIndex(grid, subgridOrigin, subgridIndex, gridPointIndex))
 468         {
 469             return true;
 470         }
 471     }
 472
 473     return false;
 474 }
 475
 476 /*! \brief
 477  * Returns the point distance between from value x to value x0 along the given axis.
 478  *
 479  * Note that the returned distance may be negative or larger than the
 480  * number of points in the axis. For a periodic axis, the distance is chosen
 481  * to be in [0, period), i.e. always positive but not the shortest one.
 482  *
 483  * \param[in]  axis   Grid axis.
 484  * \param[in]  x      From value.
 485  * \param[in]  x0     To value.
 486  * \returns (x - x0) in number of points.
 487  */
 488 static int pointDistanceAlongAxis(const GridAxis &axis, double x, double x0)
 489 {
 490     int distance = 0;
 491
 492     if (axis.spacing() > 0)
 493     {
 494         /* Get the real-valued distance. For a periodic axis, the shortest one. */
 495         double period = axis.period();
 496         double dx     = getDeviationPeriodic(x, x0, period);
 497
 498         /* Transform the distance into a point distance by rounding. */
 499         distance = gmx::roundToInt(dx/axis.spacing());
 500
 501         /* If periodic, shift the point distance to be in [0, period) */
 502         distance = indexWithinPeriod(distance, axis.numPointsInPeriod());
 503     }
 504
 505     return distance;
 506 }
 507
 508 /*! \brief
 509  * Query if a value is in range of the grid.
 510  *
 511  * \param[in] value   Value to check.
 512  * \param[in] axis    The grid axes.
 513  * \returns true if the value is in the grid.
 514  */
 515 static bool valueIsInGrid(const awh_dvec               value,
 516                           const std::vector<GridAxis> &axis)
 517 {
 518     /* For each dimension get the one-dimensional index and check if it is in range. */
 519     for (size_t d = 0; d < axis.size(); d++)
 520     {
 521         /* The index is computed as the point distance from the origin. */
 522         int index = pointDistanceAlongAxis(axis[d], value[d], axis[d].origin());
 523
 524         if (!(index >= 0 && index < axis[d].numPoints()))
 525         {
 526             return false;
 527         }
 528     }
 529
 530     return true;
 531 }
 532
 533 bool Grid::covers(const awh_dvec value) const
 534 {
 535     return valueIsInGrid(value, axis());
 536 }
 537
 538 int GridAxis::nearestIndex(double value) const
 539 {
 540     /* Get the point distance to the origin. This may by an out of index range for the axis. */
 541     int index = pointDistanceAlongAxis(*this, value, origin_);
 542
 543     if (index < 0 || index >= numPoints_)
 544     {
 545         if (isPeriodic())
 546         {
 547             GMX_RELEASE_ASSERT(index >= 0 && index < numPointsInPeriod_,
 548                                "Index not in periodic interval 0 for AWH periodic axis");
 549             int endDistance    = (index - (numPoints_ - 1));
 550             int originDistance = (numPointsInPeriod_ - index);
 551             index = originDistance < endDistance ? 0 : numPoints_ - 1;
 552         }
 553         else
 554         {
 555             index = (index < 0) ? 0 : (numPoints_ - 1);
 556         }
 557     }
 558
 559     return index;
 560 }
 561
 562 /*! \brief
 563  * Map a value to the nearest point in the grid.
 564  *
 565  * \param[in] value  Value.
 566  * \param[in] axis   The grid axes.
 567  * \returns the point index nearest to the value.
 568  */
 569 static int getNearestIndexInGrid(const awh_dvec               value,
 570                                  const std::vector<GridAxis> &axis)
 571 {
 572     awh_ivec indexMulti;
 573
 574     /* If the index is out of range, modify it so that it is in range by choosing the nearest point on the edge. */
 575     for (size_t d = 0; d < axis.size(); d++)
 576     {
 577         indexMulti[d] = axis[d].nearestIndex(value[d]);
 578     }
 579
 580     return multiDimGridIndexToLinear(axis, indexMulti);
 581 }
 582
 583 int Grid::nearestIndex(const awh_dvec value) const
 584 {
 585     return getNearestIndexInGrid(value, axis());
 586 }
 587
 588 namespace
 589 {
 590
 591 /*! \brief
 592  * Find and set the neighbors of a grid point.
 593  *
 594  * The search space for neighbors is a subgrid with size set by a scope cutoff.
 595  * In general not all point within scope will be valid grid points.
 596  *
 597  * \param[in]     pointIndex           Grid point index.
 598  * \param[in]     grid                 The grid.
 599  * \param[in,out] neighborIndexArray   Array to fill with neighbor indices.
 600  */
 601 void setNeighborsOfGridPoint(int               pointIndex,
 602                              const Grid       &grid,
 603                              std::vector<int> *neighborIndexArray)
 604 {
 605     const int c_maxNeighborsAlongAxis = 1 + 2*static_cast<int>(Grid::c_numPointsPerSigma*Grid::c_scopeCutoff);
 606
 607     awh_ivec  numCandidates           = {0};
 608     awh_ivec  subgridOrigin           = {0};
 609     for (int d = 0; d < grid.numDimensions(); d++)
 610     {
 611         /* The number of candidate points along this dimension is given by the scope cutoff. */
 612         numCandidates[d] = std::min(c_maxNeighborsAlongAxis,
 613                                     grid.axis(d).numPoints());
 614
 615         /* The origin of the subgrid to search */
 616         int centerIndex  = grid.point(pointIndex).index[d];
 617         subgridOrigin[d] = centerIndex - numCandidates[d]/2;
 618     }
 619
 620     /* Find and set the neighbors */
 621     int  neighborIndex = -1;
 622     bool aPointExists  = true;
 623
 624     /* Keep looking for grid points while traversing the subgrid. */
 625     while (aPointExists)
 626     {
 627         /* The point index is updated if a grid point was found. */
 628         aPointExists = advancePointInSubgrid(grid, subgridOrigin, numCandidates, &neighborIndex);
 629
 630         if (aPointExists)
 631         {
 632             neighborIndexArray->push_back(neighborIndex);
 633         }
 634     }
 635 }
 636
 637 }   // namespace
 638
 639 void Grid::initPoints()
 640 {
 641     awh_ivec     numPointsDimWork = { 0 };
 642     awh_ivec     indexWork        = { 0 };
 643
 644     for (size_t d = 0; d < axis_.size(); d++)
 645     {
 646         /* Temporarily gather the number of points in each dimension in one array */
 647         numPointsDimWork[d] = axis_[d].numPoints();
 648     }
 649
 650     for (auto &point : point_)
 651     {
 652         for (size_t d = 0; d < axis_.size(); d++)
 653         {
 654             point.coordValue[d] = axis_[d].origin() + indexWork[d]*axis_[d].spacing();
 655
 656             if (axis_[d].period() > 0)
 657             {
 658                 /* Do we always want the values to be centered around 0 ? */
 659                 centerPeriodicValueAroundZero(&point.coordValue[d], axis_[d].period());
 660             }
 661
 662             point.index[d] = indexWork[d];
 663         }
 664
 665         stepInMultiDimArray(axis_.size(), numPointsDimWork, indexWork);
 666     }
 667 }
 668
 669 GridAxis::GridAxis(double origin, double end,
 670                    double period, double pointDensity) :
 671     origin_(origin),
 672     period_(period)
 673 {
 674     length_ = getIntervalLengthPeriodic(origin_, end, period_);
 675
 676     /* Automatically determine number of points based on the user given endpoints
 677        and the expected fluctuations in the umbrella. */
 678     if (length_ == 0)
 679     {
 680         numPoints_ = 1;
 681     }
 682     else if (pointDensity == 0)
 683     {
 684         numPoints_ = 2;
 685     }
 686     else
 687     {
 688         /* An extra point is added here to account for the endpoints. The
 689            minimum number of points for a non-zero interval is 2. */
 690         numPoints_            = 1 + static_cast<int>(std::ceil(length_*pointDensity));
 691     }
 692
 693     /* Set point spacing based on the number of points */
 694     if (isPeriodic())
 695     {
 696         /* Set the grid spacing so that a period is matched exactly by an integer number of points.
 697            The number of points in a period is equal to the number of grid spacings in a period
 698            since the endpoints are connected.  */
 699         numPointsInPeriod_ = length_ > 0 ? static_cast<int>(std::ceil(period/length_*(numPoints_ - 1))) : 1;
 700         spacing_           = period_/numPointsInPeriod_;
 701
 702         /* Modify the number of grid axis points to be compatible with the period dependent spacing. */
 703         numPoints_         = std::min(static_cast<int>(round(length_/spacing_)) + 1,
 704                                       numPointsInPeriod_);
 705     }
 706     else
 707     {
 708         numPointsInPeriod_ = 0;
 709         spacing_           = numPoints_ > 1 ? length_/(numPoints_ - 1) : 0;
 710     }
 711 }
 712
 713 GridAxis::GridAxis(double origin, double end,
 714                    double period, int numPoints) :
 715     origin_(origin),
 716     period_(period),
 717     numPoints_(numPoints)
 718 {
 719     length_            = getIntervalLengthPeriodic(origin_, end, period_);
 720     spacing_           = numPoints_ > 1 ? length_/(numPoints_ - 1) : period_;
 721     numPointsInPeriod_ = static_cast<int>(std::round(period_/spacing_));
 722 }
 723
 724 Grid::Grid(const std::vector<DimParams> &dimParams,
 725            const AwhDimParams           *awhDimParams)
 726 {
 727     /* Define the discretization along each dimension */
 728     awh_dvec period;
 729     int      numPoints = 1;
 730     for (size_t d = 0; d < dimParams.size(); d++)
 731     {
 732         double origin = dimParams[d].scaleUserInputToInternal(awhDimParams[d].origin);
 733         double end    = dimParams[d].scaleUserInputToInternal(awhDimParams[d].end);
 734         period[d]     = dimParams[d].scaleUserInputToInternal(awhDimParams[d].period);
 735         static_assert(c_numPointsPerSigma >= 1.0, "The number of points per sigma should be at least 1.0 to get a uniformly covering the reaction using Gaussians");
 736         double pointDensity = std::sqrt(dimParams[d].betak)*c_numPointsPerSigma;
 737         axis_.emplace_back(origin, end, period[d], pointDensity);
 738         numPoints *= axis_[d].numPoints();
 739     }
 740
 741     point_.resize(numPoints);
 742
 743     /* Set their values */
 744     initPoints();
 745
 746     /* Keep a neighbor list for each point.
 747      * Note: could also generate neighbor list only when needed
 748      * instead of storing them for each point.
 749      */
 750     for (size_t m = 0; m < point_.size(); m++)
 751     {
 752         std::vector<int> *neighbor = &point_[m].neighbor;
 753
 754         setNeighborsOfGridPoint(m, *this, neighbor);
 755     }
 756 }
 757
 758 void mapGridToDataGrid(std::vector<int>    *gridpointToDatapoint,
 759                        const double* const *data,
 760                        int                  numDataPoints,
 761                        const std::string   &dataFilename,
 762                        const Grid          &grid,
 763                        const std::string   &correctFormatMessage)
 764 {
 765     /* Transform the data into a grid in order to map each grid point to a data point
 766        using the grid functions. */
 767
 768     /* Count the number of points for each dimension. Each dimension
 769        has its own stride. */
 770     int              stride           = 1;
 771     int              numPointsCounted = 0;
 772     std::vector<int> numPoints(grid.numDimensions());
 773     for (int d = grid.numDimensions() - 1; d >= 0; d--)
 774     {
 775         int    numPointsInDim = 0;
 776         int    pointIndex     = 0;
 777         double firstValue     = data[d][pointIndex];
 778         do
 779         {
 780             numPointsInDim++;
 781             pointIndex       += stride;
 782         }
 783         while (pointIndex < numDataPoints &&
 784                !gmx_within_tol(firstValue, data[d][pointIndex], GMX_REAL_EPS));
 785
 786         /* The stride in dimension dimension d - 1 equals the number of points
 787            dimension d. */
 788         stride = numPointsInDim;
 789
 790         numPointsCounted = (numPointsCounted == 0) ? numPointsInDim : numPointsCounted*numPointsInDim;
 791
 792         numPoints[d]     = numPointsInDim;
 793     }
 794
 795     if (numPointsCounted != numDataPoints)
 796     {
 797         std::string mesg = gmx::formatString("Could not extract data properly from %s. Wrong data format?"
 798                                              "\n\n%s",
 799                                              dataFilename.c_str(), correctFormatMessage.c_str());
 800         GMX_THROW(InvalidInputError(mesg));
 801     }
 802
 803     std::vector<GridAxis> axis_;
 804     axis_.reserve(grid.numDimensions());
 805     /* The data grid has the data that was read and the properties of the AWH grid */
 806     for (int d = 0; d < grid.numDimensions(); d++)
 807     {
 808         axis_.emplace_back(data[d][0], data[d][numDataPoints - 1],
 809                            grid.axis(d).period(), numPoints[d]);
 810     }
 811
 812     /* Map each grid point to a data point. No interpolation, just pick the nearest one.
 813      * It is assumed that the given data is uniformly spaced for each dimension.
 814      */
 815     for (size_t m = 0; m < grid.numPoints(); m++)
 816     {
 817         /* We only define what we need for the datagrid since it's not needed here which is a bit ugly */
 818
 819         if (!valueIsInGrid(grid.point(m).coordValue, axis_))
 820         {
 821             std::string mesg =
 822                 gmx::formatString("%s does not contain data for all coordinate values. "
 823                                   "Make sure your input data covers the whole sampling domain "
 824                                   "and is correctly formatted. \n\n%s",
 825                                   dataFilename.c_str(), correctFormatMessage.c_str());
 826             GMX_THROW(InvalidInputError(mesg));
 827         }
 828         (*gridpointToDatapoint)[m] = getNearestIndexInGrid(grid.point(m).coordValue, axis_);
 829     }
 830 }
 831
 832 } // namespace gmx