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