2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2015,2016,2017,2018,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.
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.
39 * Declares and defines the PointState class.
41 * Since nearly all operations on PointState objects occur in loops over
42 * (parts of) the grid of an AWH bias, all these methods should be inlined.
43 * Only samplePmf() is called only once per step and is thus not inlined.
45 * \author Viveca Lindahl
46 * \author Berk Hess <hess@kth.se>
50 #ifndef GMX_AWH_POINTSTATE_H
51 #define GMX_AWH_POINTSTATE_H
55 #include "gromacs/mdtypes/awh_history.h"
56 #include "gromacs/mdtypes/awh_params.h"
57 #include "gromacs/utility/gmxassert.h"
59 #include "biasparams.h"
67 //! A value that can be passed to exp() with result 0, also with SIMD
68 constexpr double c_largeNegativeExponent = -10000.0;
70 //! The largest acceptable positive exponent for variables that are passed to exp().
71 constexpr double c_largePositiveExponent = 700.0;
76 * \brief The state of a coordinate point.
78 * This class contains all the state variables of a coordinate point
79 * (on the bias grid) and methods to update the state of a point.
84 /*! \brief Constructs a point state with default values. */
89 targetConstantWeight_(1),
90 weightSumIteration_(0),
95 numVisitsIteration_(0),
101 * Set all values in the state to those from a history.
103 * \param[in] psh Coordinate point history to copy from.
105 void setFromHistory(const AwhPointStateHistory& psh)
107 target_ = psh.target;
108 freeEnergy_ = psh.free_energy;
110 weightSumIteration_ = psh.weightsum_iteration;
111 weightSumTot_ = psh.weightsum_tot;
112 weightSumRef_ = psh.weightsum_ref;
113 lastUpdateIndex_ = psh.last_update_index;
114 logPmfSum_ = psh.log_pmfsum;
115 numVisitsIteration_ = psh.visits_iteration;
116 numVisitsTot_ = psh.visits_tot;
120 * Store the state of a point in a history struct.
122 * \param[in,out] psh Coordinate point history to copy to.
124 void storeState(AwhPointStateHistory* psh) const
126 psh->target = target_;
127 psh->free_energy = freeEnergy_;
129 psh->weightsum_iteration = weightSumIteration_;
130 psh->weightsum_tot = weightSumTot_;
131 psh->weightsum_ref = weightSumRef_;
132 psh->last_update_index = lastUpdateIndex_;
133 psh->log_pmfsum = logPmfSum_;
134 psh->visits_iteration = numVisitsIteration_;
135 psh->visits_tot = numVisitsTot_;
139 * Query if the point is in the target region.
141 * \returns true if the point is in the target region.
143 bool inTargetRegion() const { return target_ > 0; }
145 /*! \brief Return the bias function estimate. */
146 double bias() const { return bias_; }
148 /*! \brief Set the target to zero and the bias to minus infinity. */
149 void setTargetToZero()
152 /* the bias = log(target) + const = -infty */
153 bias_ = detail::c_largeNegativeExponent;
156 /*! \brief Return the free energy. */
157 double freeEnergy() const { return freeEnergy_; }
159 /*! \brief Set the free energy, only to be used at initialization.
161 * \param[in] freeEnergy The free energy.
163 void setFreeEnergy(double freeEnergy) { freeEnergy_ = freeEnergy; }
165 /*! \brief Return the target distribution value. */
166 double target() const { return target_; }
168 /*! \brief Return the weight accumulated since the last update. */
169 double weightSumIteration() const { return weightSumIteration_; }
171 /*! \brief Increases the weight accumulated since the last update.
173 * \param[in] weight The amount to add to the weight
175 void increaseWeightSumIteration(double weight) { weightSumIteration_ += weight; }
177 /*! \brief Returns the accumulated weight */
178 double weightSumTot() const { return weightSumTot_; }
180 /*! \brief Return the reference weight histogram. */
181 double weightSumRef() const { return weightSumRef_; }
183 /*! \brief Return log(PmfSum). */
184 double logPmfSum() const { return logPmfSum_; }
186 /*! \brief Set log(PmfSum).
188 * TODO: Replace this setter function with a more elegant solution.
190 * \param[in] logPmfSum The log(PmfSum).
192 void setLogPmfSum(double logPmfSum) { logPmfSum_ = logPmfSum; }
194 /*! \brief Return the number of visits since the last update */
195 double numVisitsIteration() const { return numVisitsIteration_; }
197 /*! \brief Return the total number of visits */
198 double numVisitsTot() const { return numVisitsTot_; }
200 /*! \brief Set the constant target weight factor.
202 * \param[in] targetConstantWeight The target weight factor.
204 void setTargetConstantWeight(double targetConstantWeight)
206 targetConstantWeight_ = targetConstantWeight;
209 /*! \brief Updates the bias of a point. */
212 GMX_ASSERT(target_ > 0, "AWH target distribution must be > 0 to calculate the point bias.");
214 bias_ = freeEnergy() + std::log(target_);
217 /*! \brief Set the initial reference weighthistogram.
219 * \param[in] histogramSize The weight histogram size.
221 void setInitialReferenceWeightHistogram(double histogramSize)
223 weightSumRef_ = histogramSize * target_;
226 /*! \brief Correct free energy and PMF sum for the change in minimum.
228 * \param[in] minimumFreeEnergy The free energy at the minimum;
230 void normalizeFreeEnergyAndPmfSum(double minimumFreeEnergy)
232 if (inTargetRegion())
234 /* The sign of the free energy and PMF constants are opposite
235 * because the PMF samples are reweighted with the negative
236 * bias e^(-bias) ~ e^(-free energy).
238 freeEnergy_ -= minimumFreeEnergy;
239 logPmfSum_ += minimumFreeEnergy;
243 /*! \brief Apply previous updates that were skipped.
245 * An update can only be skipped if the parameters needed for the update are constant or
246 * deterministic so that the same update can be performed at a later time.
247 * Here, the necessary parameters are the sampled weight and scaling factors for the
248 * histograms. The scaling factors are provided as arguments only to avoid recalculating
249 * them for each point
251 * The last update index is also updated here.
253 * \param[in] params The AWH bias parameters.
254 * \param[in] numUpdates The global number of updates.
255 * \param[in] weighthistScaling Scale factor for the reference weight histogram.
256 * \param[in] logPmfSumScaling Scale factor for the reference PMF histogram.
257 * \returns true if at least one update was applied.
259 bool performPreviouslySkippedUpdates(const BiasParams& params,
261 double weighthistScaling,
262 double logPmfSumScaling)
264 GMX_ASSERT(params.skipUpdates(),
265 "Calling function for skipped updates when skipping updates is not allowed");
267 if (!inTargetRegion())
272 /* The most current past update */
273 int64_t lastUpdateIndex = numUpdates;
274 int64_t numUpdatesSkipped = lastUpdateIndex - lastUpdateIndex_;
276 if (numUpdatesSkipped == 0)
278 /* Was not updated */
282 for (int64_t i = 0; i < numUpdatesSkipped; i++)
284 /* This point was non-local at the time of the update meaning no weight */
285 updateFreeEnergyAndWeight(params, 0, weighthistScaling, logPmfSumScaling);
288 /* Only past updates are applied here. */
289 lastUpdateIndex_ = lastUpdateIndex;
294 /*! \brief Apply a point update with new sampling.
296 * \note The last update index is also updated here.
297 * \note The new sampling containers are cleared here.
299 * \param[in] params The AWH bias parameters.
300 * \param[in] numUpdates The global number of updates.
301 * \param[in] weighthistScaling Scaling factor for the reference weight histogram.
302 * \param[in] logPmfSumScaling Log of the scaling factor for the PMF histogram.
304 void updateWithNewSampling(const BiasParams& params, int64_t numUpdates, double weighthistScaling, double logPmfSumScaling)
306 GMX_RELEASE_ASSERT(lastUpdateIndex_ == numUpdates,
307 "When doing a normal update, the point update index should match the "
308 "global index, otherwise we lost (skipped?) updates.");
310 updateFreeEnergyAndWeight(params, weightSumIteration_, weighthistScaling, logPmfSumScaling);
311 lastUpdateIndex_ += 1;
313 /* Clear the iteration collection data */
314 weightSumIteration_ = 0;
315 numVisitsIteration_ = 0;
319 /*! \brief Update the PMF histogram with the current coordinate value.
321 * \param[in] convolvedBias The convolved bias.
323 void samplePmf(double convolvedBias);
325 /*! \brief Update the PMF histogram of unvisited coordinate values
326 * (along a lambda axis)
328 * \param[in] bias The bias to update with.
330 void updatePmfUnvisited(double bias);
333 /*! \brief Update the free energy estimate of a point.
335 * The free energy update here is inherently local, i.e. it just depends on local sampling and
336 * on constant AWH parameters. This assumes that the variables used here are kept constant, at
337 * least in between global updates.
339 * \param[in] params The AWH bias parameters.
340 * \param[in] weightAtPoint Sampled probability weight at this point.
342 void updateFreeEnergy(const BiasParams& params, double weightAtPoint)
344 double weighthistSampled = weightSumRef() + weightAtPoint;
345 double weighthistTarget = weightSumRef() + params.updateWeight * target_;
347 double df = -std::log(weighthistSampled / weighthistTarget);
350 GMX_RELEASE_ASSERT(std::abs(freeEnergy_) < detail::c_largePositiveExponent,
351 "Very large free energy differences or badly normalized free energy in "
355 /*! \brief Update the reference weight histogram of a point.
357 * \param[in] params The AWH bias parameters.
358 * \param[in] weightAtPoint Sampled probability weight at this point.
359 * \param[in] scaleFactor Factor to rescale the histogram with.
361 void updateWeightHistogram(const BiasParams& params, double weightAtPoint, double scaleFactor)
363 if (params.idealWeighthistUpdate)
365 /* Grow histogram using the target distribution. */
366 weightSumRef_ += target_ * params.updateWeight * params.localWeightScaling;
370 /* Grow using the actual samples (which are distributed ~ as target). */
371 weightSumRef_ += weightAtPoint * params.localWeightScaling;
374 weightSumRef_ *= scaleFactor;
377 /*! \brief Apply a point update.
379 * This updates local properties that can be updated without
380 * accessing or affecting all points.
381 * This excludes updating the size of reference weight histogram and
382 * the target distribution. The bias update is excluded only because
383 * if updates have been skipped this function will be called multiple
384 * times, while the bias only needs to be updated once (last).
386 * Since this function only performs the update with the given
387 * arguments and does not know anything about the time of the update,
388 * the last update index is not updated here. The caller should take
389 * care of updating the update index.
391 * \param[in] params The AWH bias parameters.
392 * \param[in] weightAtPoint Sampled probability weight at this point.
393 * \param[in] weighthistScaling Scaling factor for the reference weight histogram.
394 * \param[in] logPmfSumScaling Log of the scaling factor for the PMF histogram.
396 void updateFreeEnergyAndWeight(const BiasParams& params,
397 double weightAtPoint,
398 double weighthistScaling,
399 double logPmfSumScaling)
401 updateFreeEnergy(params, weightAtPoint);
402 updateWeightHistogram(params, weightAtPoint, weighthistScaling);
403 logPmfSum_ += logPmfSumScaling;
408 /*! \brief Update the target weight of a point.
410 * Note that renormalization over all points is needed after the update.
412 * \param[in] params The AWH bias parameters.
413 * \param[in] freeEnergyCutoff The cut-off for the free energy for target type "cutoff".
414 * \returns the updated value of the target.
416 double updateTargetWeight(const BiasParams& params, double freeEnergyCutoff)
418 switch (params.eTarget)
420 case eawhtargetCONSTANT: target_ = 1; break;
421 case eawhtargetCUTOFF:
423 double df = freeEnergy_ - freeEnergyCutoff;
424 target_ = 1 / (1 + std::exp(df));
427 case eawhtargetBOLTZMANN:
428 target_ = std::exp(-params.temperatureScaleFactor * freeEnergy_);
430 case eawhtargetLOCALBOLTZMANN: target_ = weightSumRef_; break;
433 /* All target types can be modulated by a constant factor. */
434 target_ *= targetConstantWeight_;
439 /*! \brief Set the weight and count accumulated since the last update.
441 * \param[in] weightSum The weight-sum value
442 * \param[in] numVisits The number of visits
444 void setPartialWeightAndCount(double weightSum, double numVisits)
446 weightSumIteration_ = weightSum;
447 numVisitsIteration_ = numVisits;
450 /*! \brief Add the weights and counts accumulated between updates. */
451 void addPartialWeightAndCount()
453 weightSumTot_ += weightSumIteration_;
454 numVisitsTot_ += numVisitsIteration_;
457 /*! \brief Scale the target weight of the point.
459 * \param[in] scaleFactor Factor to scale with.
461 void scaleTarget(double scaleFactor) { target_ *= scaleFactor; }
464 double bias_; /**< Current biasing function estimate */
465 double freeEnergy_; /**< Current estimate of the convolved free energy/PMF. */
466 double target_; /**< Current target distribution, normalized to 1 */
467 double targetConstantWeight_; /**< Constant target weight, from user data. */
468 double weightSumIteration_; /**< Accumulated weight this iteration; note: only contains data for this Bias, even when sharing biases. */
469 double weightSumTot_; /**< Accumulated weights, never reset */
470 double weightSumRef_; /**< The reference weight histogram determining the free energy updates */
471 int64_t lastUpdateIndex_; /**< The last update that was performed at this point, in units of number of updates. */
472 double logPmfSum_; /**< Logarithm of the PMF histogram */
473 double numVisitsIteration_; /**< Visits to this bin this iteration; note: only contains data for this Bias, even when sharing biases. */
474 double numVisitsTot_; /**< Accumulated visits to this bin */
479 #endif /* GMX_AWH_POINTSTATE_H */