Make AWH parameters proper C++
[alexxy/gromacs.git] / src / gromacs / applied_forces / awh / coordstate.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
4  * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, by the GROMACS development team, led by
5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6  * and including many others, as listed in the AUTHORS file in the
7  * top-level source directory and at http://www.gromacs.org.
8  *
9  * GROMACS is free software; you can redistribute it and/or
10  * modify it under the terms of the GNU Lesser General Public License
11  * as published by the Free Software Foundation; either version 2.1
12  * of the License, or (at your option) any later version.
13  *
14  * GROMACS is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17  * Lesser General Public License for more details.
18  *
19  * You should have received a copy of the GNU Lesser General Public
20  * License along with GROMACS; if not, see
21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
23  *
24  * If you want to redistribute modifications to GROMACS, please
25  * consider that scientific software is very special. Version
26  * control is crucial - bugs must be traceable. We will be happy to
27  * consider code for inclusion in the official distribution, but
28  * derived work must not be called official GROMACS. Details are found
29  * in the README & COPYING files - if they are missing, get the
30  * official version at http://www.gromacs.org.
31  *
32  * To help us fund GROMACS development, we humbly ask that you cite
33  * the research papers on the package. Check out http://www.gromacs.org.
34  */
35
36 /*! \internal \file
37  * \brief
38  * Implements the CoordState class.
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 "coordstate.h"
48
49 #include <algorithm>
50
51 #include "gromacs/math/utilities.h"
52 #include "gromacs/mdtypes/awh_history.h"
53 #include "gromacs/mdtypes/awh_params.h"
54 #include "gromacs/random/threefry.h"
55 #include "gromacs/random/uniformrealdistribution.h"
56 #include "gromacs/utility/arrayref.h"
57 #include "gromacs/utility/exceptions.h"
58 #include "gromacs/utility/gmxassert.h"
59 #include "gromacs/utility/stringutil.h"
60
61 #include "biasgrid.h"
62
63 namespace gmx
64 {
65
66 CoordState::CoordState(const AwhBiasParams&      awhBiasParams,
67                        ArrayRef<const DimParams> dimParams,
68                        const BiasGrid&           grid)
69 {
70     GMX_RELEASE_ASSERT(awhBiasParams.ndim() == dimParams.ssize(),
71                        "Need to have identical size for dimensions");
72     const auto& awhDimParams = awhBiasParams.dimParams();
73     for (int d = 0; d < gmx::ssize(awhDimParams); d++)
74     {
75         coordValue_[d] = dimParams[d].scaleUserInputToInternal(awhDimParams[d].initialCoordinate());
76     }
77
78     /* The grid-point index is always the nearest point to the coordinate.
79      * We initialize the umbrella location to the same nearest point.
80      * More correctly one would sample from the biased distribution,
81      * but it doesn't really matter, as this will happen after a few steps.
82      */
83     gridpointIndex_    = grid.nearestIndex(coordValue_);
84     umbrellaGridpoint_ = gridpointIndex_;
85 }
86
87 namespace
88 {
89
90 /*! \brief Generate a sample from a discrete probability distribution defined on [0, distr.size() - 1].
91  *
92  * The pair (indexSeed0,indexSeed1) should be different for every invocation.
93  *
94  * \param[in] distr       Normalized probability distribution to generate a sample from.
95  * \param[in] seed        Random seed for initializing the random number generator.
96  * \param[in] indexSeed0  Random seed needed by the random number generator.
97  * \param[in] indexSeed1  Random seed needed by the random number generator.
98  * \returns a sample index in [0, distr.size() - 1]
99  */
100 int getSampleFromDistribution(ArrayRef<const double> distr, int64_t seed, int64_t indexSeed0, int64_t indexSeed1)
101 {
102     gmx::ThreeFry2x64<0>               rng(seed, gmx::RandomDomain::AwhBiasing);
103     gmx::UniformRealDistribution<real> uniformRealDistr;
104
105     GMX_RELEASE_ASSERT(!distr.empty(), "We need a non-zero length distribution to sample from");
106
107     /* Generate the cumulative probability distribution function */
108     std::vector<double> cumulativeDistribution(distr.size());
109
110     cumulativeDistribution[0] = distr[0];
111
112     for (gmx::index i = 1; i < distr.ssize(); i++)
113     {
114         cumulativeDistribution[i] = cumulativeDistribution[i - 1] + distr[i];
115     }
116
117     GMX_RELEASE_ASSERT(gmx_within_tol(cumulativeDistribution.back(), 1.0, 0.01),
118                        "Attempt to get sample from non-normalized/zero distribution");
119
120     /* Use binary search to convert the real value to an integer in [0, ndistr - 1] distributed according to distr. */
121     rng.restart(indexSeed0, indexSeed1);
122
123     double value = uniformRealDistr(rng);
124     int sample = std::upper_bound(cumulativeDistribution.begin(), cumulativeDistribution.end() - 1, value)
125                  - cumulativeDistribution.begin();
126
127     return sample;
128 }
129
130 } // namespace
131
132 void CoordState::sampleUmbrellaGridpoint(const BiasGrid&             grid,
133                                          int                         gridpointIndex,
134                                          gmx::ArrayRef<const double> probWeightNeighbor,
135                                          int64_t                     step,
136                                          int64_t                     seed,
137                                          int                         indexSeed)
138 {
139     /* Sample new umbrella reference value from the probability distribution
140      * which is defined for the neighboring points of the current coordinate.
141      */
142     const std::vector<int>& neighbor = grid.point(gridpointIndex).neighbor;
143
144     /* In order to use the same seed for all AWH biases and get independent
145        samples we use the index of the bias. */
146     int localIndex = getSampleFromDistribution(probWeightNeighbor, seed, step, indexSeed);
147
148     umbrellaGridpoint_ = neighbor[localIndex];
149 }
150
151 void CoordState::setCoordValue(const BiasGrid& grid, const awh_dvec coordValue)
152 {
153     /* We need to check for valid (probable) coordinate values, to give
154      * a clear error message instead of a low-level assertion failure.
155      * We allow values up to 10*sigma beyond the bounds. For points at
156      * the bounds this means a chance of less than 2*10-45 of a false positive
157      * in the usual case that the PMF is flat or goes up beyond the bounds.
158      * In the worst reasonable case of the PMF curving down with a curvature
159      * of half the harmonic force constant, the chance is 1.5*10-12.
160      */
161     constexpr int c_marginInSigma = 10;
162
163     for (int dim = 0; dim < grid.numDimensions(); dim++)
164     {
165         const GridAxis& axis = grid.axis(dim);
166         /* We do not check periodic coordinates, since that is more complicated
167          * and those cases are less likely to cause problems.
168          */
169         if (!axis.isPeriodic())
170         {
171             const double margin = axis.spacing() * c_marginInSigma / BiasGrid::c_numPointsPerSigma;
172             if (coordValue[dim] < axis.origin() - margin
173                 || coordValue[dim] > axis.origin() + axis.length() + margin)
174             {
175                 std::string mesg = gmx::formatString(
176                         "Coordinate %d of an AWH bias has a value %f which is more than %d sigma "
177                         "out of the AWH range of [%f, %f]. You seem to have an unstable reaction "
178                         "coordinate setup or an unequilibrated system.",
179                         dim + 1,
180                         coordValue[dim],
181                         c_marginInSigma,
182                         axis.origin(),
183                         axis.origin() + axis.length());
184                 GMX_THROW(SimulationInstabilityError(mesg));
185             }
186         }
187
188         coordValue_[dim] = coordValue[dim];
189     }
190
191     /* The grid point closest to the coordinate value defines the current
192      * neighborhood of points. Besides at steps when global updates and/or
193      * checks are performed, only the neighborhood will be touched.
194      */
195     gridpointIndex_ = grid.nearestIndex(coordValue_);
196 }
197
198 void CoordState::restoreFromHistory(const AwhBiasStateHistory& stateHistory)
199 {
200     umbrellaGridpoint_ = stateHistory.umbrellaGridpoint;
201 }
202
203 void CoordState::setUmbrellaGridpointToGridpoint()
204 {
205     umbrellaGridpoint_ = gridpointIndex_;
206 }
207
208 } // namespace gmx