2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
38 * Implements the CoordState class.
40 * \author Viveca Lindahl
41 * \author Berk Hess <hess@kth.se>
47 #include "coordstate.h"
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/gmxassert.h"
64 CoordState::CoordState(const AwhBiasParams
&awhBiasParams
,
65 const std::vector
<DimParams
> &dimParams
,
68 for (size_t d
= 0; d
< dimParams
.size(); d
++)
70 coordValue_
[d
] = dimParams
[d
].scaleUserInputToInternal(awhBiasParams
.dimParams
[d
].coordValueInit
);
73 /* The grid-point index is always the nearest point to the coordinate.
74 * We initialize the umbrella location to the same nearest point.
75 * More correctly one would sample from the biased distribution,
76 * but it doesn't really matter, as this will happen after a few steps.
78 gridpointIndex_
= grid
.nearestIndex(coordValue_
);
79 umbrellaGridpoint_
= gridpointIndex_
;
85 /*! \brief Generate a sample from a discrete probability distribution defined on [0, distr.size() - 1].
87 * The pair (indexSeed0,indexSeed1) should be different for every invocation.
89 * \param[in] distr Normalized probability distribution to generate a sample from.
90 * \param[in] seed Random seed for initializing the random number generator.
91 * \param[in] indexSeed0 Random seed needed by the random number generator.
92 * \param[in] indexSeed1 Random seed needed by the random number generator.
93 * \returns a sample index in [0, distr.size() - 1]
95 int getSampleFromDistribution(ArrayRef
<const double> distr
,
97 gmx_int64_t indexSeed0
,
98 gmx_int64_t indexSeed1
)
100 gmx::ThreeFry2x64
<0> rng(seed
, gmx::RandomDomain::AwhBiasing
);
101 gmx::UniformRealDistribution
<real
> uniformRealDistr
;
103 GMX_RELEASE_ASSERT(distr
.size() > 0, "We need a non-zero length distribution to sample from");
105 /* Generate the cumulative probability distribution function */
106 std::vector
<double> cumulativeDistribution(distr
.size());
108 cumulativeDistribution
[0] = distr
[0];
110 for (size_t i
= 1; i
< distr
.size(); i
++)
112 cumulativeDistribution
[i
] = cumulativeDistribution
[i
- 1] + distr
[i
];
115 GMX_RELEASE_ASSERT(gmx_within_tol(cumulativeDistribution
.back(), 1.0, 0.01), "Attempt to get sample from non-normalized/zero distribution");
117 /* Use binary search to convert the real value to an integer in [0, ndistr - 1] distributed according to distr. */
118 rng
.restart(indexSeed0
, indexSeed1
);
120 double value
= uniformRealDistr(rng
);
121 int sample
= std::upper_bound(cumulativeDistribution
.begin(), cumulativeDistribution
.end() - 1, value
) - cumulativeDistribution
.begin();
129 CoordState::sampleUmbrellaGridpoint(const Grid
&grid
,
131 gmx::ArrayRef
<const double> probWeightNeighbor
,
136 /* Sample new umbrella reference value from the probability distribution
137 * which is defined for the neighboring points of the current coordinate.
139 const std::vector
<int> &neighbor
= grid
.point(gridpointIndex
).neighbor
;
141 /* In order to use the same seed for all AWH biases and get independent
142 samples we use the index of the bias. */
143 int localIndex
= getSampleFromDistribution(probWeightNeighbor
,
144 seed
, step
, indexSeed
);
146 umbrellaGridpoint_
= neighbor
[localIndex
];
149 void CoordState::setCoordValue(const Grid
&grid
,
150 const awh_dvec coordValue
)
152 for (int dim
= 0; dim
< grid
.numDimensions(); dim
++)
154 coordValue_
[dim
] = coordValue
[dim
];
157 /* The grid point closest to the coordinate value defines the current
158 * neighborhood of points. Besides at steps when global updates and/or
159 * checks are performed, only the neighborhood will be touched.
161 gridpointIndex_
= grid
.nearestIndex(coordValue_
);
165 CoordState::restoreFromHistory(const AwhBiasStateHistory
&stateHistory
)
167 umbrellaGridpoint_
= stateHistory
.umbrellaGridpoint
;