2 * This file is part of the GROMACS molecular simulation package.
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.
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 Bias class.
40 * \author Viveca Lindahl
41 * \author Berk Hess <hess@kth.se>
58 #include "gromacs/compat/make_unique.h"
59 #include "gromacs/fileio/gmxfio.h"
60 #include "gromacs/gmxlib/network.h"
61 #include "gromacs/math/utilities.h"
62 #include "gromacs/mdtypes/awh-history.h"
63 #include "gromacs/mdtypes/awh-params.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/utility/exceptions.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/stringutil.h"
69 #include "correlationgrid.h"
70 #include "correlationhistory.h"
71 #include "pointstate.h"
76 void Bias::warnForHistogramAnomalies(double t, int64_t step, FILE *fplog)
78 const int maxNumWarningsInCheck = 1; /* The maximum number of warnings to print per check */
79 const int maxNumWarningsInRun = 10; /* The maximum number of warnings to print in a run */
81 if (fplog == nullptr || numWarningsIssued_ >= maxNumWarningsInRun || state_.inInitialStage() ||
82 !params_.isCheckHistogramForAnomaliesStep(step))
88 state_.warnForHistogramAnomalies(grid_, biasIndex(), t, fplog,
89 maxNumWarningsInCheck);
91 if (numWarningsIssued_ >= maxNumWarningsInRun)
93 fprintf(fplog, "\nawh%d: suppressing future AWH warnings.\n", biasIndex() + 1);
97 void Bias::doSkippedUpdatesForAllPoints()
99 if (params_.skipUpdates())
101 state_.doSkippedUpdatesForAllPoints(params_);
105 gmx::ArrayRef<const double>
106 Bias::calcForceAndUpdateBias(const awh_dvec coordValue,
107 double *awhPotential,
108 double *potentialJump,
109 const t_commrec *commRecord,
110 const gmx_multisim_t *ms,
118 GMX_THROW(InvalidInputError("The step number is negative which is not supported by the AWH code."));
121 state_.setCoordValue(grid_, coordValue);
123 std::vector < double, AlignedAllocator < double>> &probWeightNeighbor = alignedTempWorkSpace_;
125 /* If the convolved force is needed or this is a sampling step,
126 * the bias in the current neighborhood needs to be up-to-date
127 * and the probablity weights need to be calculated.
129 const bool sampleCoord = params_.isSampleCoordStep(step);
130 const bool moveUmbrella = (sampleCoord || step == 0);
131 double convolvedBias = 0;
132 if (params_.convolveForce || moveUmbrella || sampleCoord)
134 if (params_.skipUpdates())
136 state_.doSkippedUpdatesInNeighborhood(params_, grid_);
139 convolvedBias = state_.updateProbabilityWeightsAndConvolvedBias(dimParams_, grid_, &probWeightNeighbor);
143 updateForceCorrelationGrid(probWeightNeighbor, t);
145 state_.sampleCoordAndPmf(grid_, probWeightNeighbor, convolvedBias);
149 const CoordState &coordState = state_.coordState();
151 /* Set the bias force and get the potential contribution from this bias.
152 * The potential jump occurs at different times depending on how
153 * the force is applied (and how the potential is normalized).
154 * For the convolved force it happens when the bias is updated,
155 * for the umbrella when the umbrella is moved.
159 if (params_.convolveForce)
161 state_.calcConvolvedForce(dimParams_, grid_, probWeightNeighbor,
162 tempForce_, biasForce_);
164 potential = -convolvedBias*params_.invBeta;
169 GMX_RELEASE_ASSERT(state_.points()[coordState.umbrellaGridpoint()].inTargetRegion(),
170 "AWH bias grid point for the umbrella reference value is outside of the target region.");
172 state_.calcUmbrellaForceAndPotential(dimParams_, grid_, coordState.umbrellaGridpoint(), biasForce_);
174 /* Moving the umbrella results in a force correction and
175 * a new potential. The umbrella center is sampled as often as
176 * the coordinate so we know the probability weights needed
177 * for moving the umbrella are up-to-date.
181 double newPotential = state_.moveUmbrella(dimParams_, grid_, probWeightNeighbor, biasForce_, step, seed, params_.biasIndex);
182 *potentialJump = newPotential - potential;
186 /* Update the free energy estimates and bias and other history dependent method parameters */
187 if (params_.isUpdateFreeEnergyStep(step))
189 state_.updateFreeEnergyAndAddSamplesToHistogram(dimParams_, grid_,
195 if (params_.convolveForce)
197 /* The update results in a potential jump, so we need the new convolved potential. */
198 double newPotential = -calcConvolvedBias(coordState.coordValue())*params_.invBeta;
199 *potentialJump = newPotential - potential;
203 /* Return the potential. */
204 *awhPotential = potential;
206 /* Check the sampled histograms and potentially warn user if something is suspicious */
207 warnForHistogramAnomalies(t, step, fplog);
213 * Count the total number of samples / sample weight over all grid points.
215 * \param[in] pointState The state of the points in a bias.
216 * \returns the total sample count.
218 static int64_t countSamples(const std::vector<PointState> &pointState)
220 double numSamples = 0;
221 for (const PointState &point : pointState)
223 numSamples += point.weightSumTot();
226 return static_cast<int64_t>(numSamples + 0.5);
230 * Check if the state (loaded from checkpoint) and the run are consistent.
232 * When the state and the run setup are inconsistent, an exception is thrown.
234 * \param[in] params The parameters of the bias.
235 * \param[in] state The state of the bias.
237 static void ensureStateAndRunConsistency(const BiasParams ¶ms,
238 const BiasState &state)
240 int64_t numSamples = countSamples(state.points());
241 int64_t numUpdatesFromSamples = numSamples/(params.numSamplesUpdateFreeEnergy_*params.numSharedUpdate);
242 int64_t numUpdatesExpected = state.histogramSize().numUpdates();
243 if (numUpdatesFromSamples != numUpdatesExpected)
245 std::string mesg = gmx::formatString("The number of AWH updates in the checkpoint file (%ld) does not match the total number of AWH samples divided by the number of samples per update for %d sharing AWH bias(es) (%ld/%d=%ld)",
247 params.numSharedUpdate,
249 params.numSamplesUpdateFreeEnergy_*params.numSharedUpdate,
250 numUpdatesFromSamples);
251 mesg += " Maybe you changed AWH parameters.";
252 /* Unfortunately we currently do not store the number of simulations
253 * sharing the bias or the state to checkpoint. But we can hint at
254 * a mismatch in the number of sharing simulations.
256 if (numUpdatesFromSamples % state.histogramSize().numUpdates() == 0)
258 mesg += gmx::formatString(" Or the run you continued from used %ld sharing simulations, whereas you now specified %d sharing simulations.",
259 numUpdatesFromSamples/state.histogramSize().numUpdates(),
260 params.numSharedUpdate);
262 GMX_THROW(InvalidInputError(mesg));
266 void Bias::restoreStateFromHistory(const AwhBiasHistory *biasHistory,
269 GMX_RELEASE_ASSERT(thisRankDoesIO_ == MASTER(cr), "The master rank should do I/O, the other ranks should not");
273 GMX_RELEASE_ASSERT(biasHistory != nullptr, "On the master rank we need a valid history object to restore from");
274 state_.restoreFromHistory(*biasHistory, grid_);
276 /* Ensure that the state is consistent with our current run setup,
277 * since the user can have changed AWH parameters or the number
278 * of simulations sharing the bias.
280 ensureStateAndRunConsistency(params_, state_);
282 if (forceCorrelationGrid_ != nullptr)
284 forceCorrelationGrid_->restoreStateFromHistory(biasHistory->forceCorrelationGrid);
290 state_.broadcast(cr);
294 void Bias::initHistoryFromState(AwhBiasHistory *biasHistory) const
296 GMX_RELEASE_ASSERT(biasHistory != nullptr, "Need a valid biasHistory");
298 state_.initHistoryFromState(biasHistory);
300 if (forceCorrelationGrid_ != nullptr)
302 biasHistory->forceCorrelationGrid = initCorrelationGridHistoryFromState(forceCorrelationGrid());
306 void Bias::updateHistory(AwhBiasHistory *biasHistory) const
308 GMX_RELEASE_ASSERT(biasHistory != nullptr, "Need a valid biasHistory");
310 state_.updateHistory(biasHistory, grid_);
312 if (forceCorrelationGrid_ != nullptr)
314 updateCorrelationGridHistory(&biasHistory->forceCorrelationGrid, forceCorrelationGrid());
318 Bias::Bias(int biasIndexInCollection,
319 const AwhParams &awhParams,
320 const AwhBiasParams &awhBiasParams,
321 const std::vector<DimParams> &dimParamsInit,
324 int numSharingSimulations,
325 const std::string &biasInitFilename,
326 ThisRankWillDoIO thisRankWillDoIO,
327 BiasParams::DisableUpdateSkips disableUpdateSkips) :
328 dimParams_(dimParamsInit),
329 grid_(dimParamsInit, awhBiasParams.dimParams),
330 params_(awhParams, awhBiasParams, dimParams_, beta, mdTimeStep, disableUpdateSkips, numSharingSimulations, grid_.axis(), biasIndexInCollection),
331 state_(awhBiasParams, params_.initialHistogramSize, dimParams_, grid_),
332 thisRankDoesIO_(thisRankWillDoIO == ThisRankWillDoIO::Yes),
336 numWarningsIssued_(0)
338 /* For a global update updateList covers all points, so reserve that */
339 updateList_.reserve(grid_.numPoints());
341 state_.initGridPointState(awhBiasParams, dimParams_, grid_, params_, biasInitFilename, awhParams.numBias);
345 /* Set up the force correlation object. */
347 /* We let the correlation init function set its parameters
348 * to something useful for now.
350 double blockLength = 0;
351 /* Construct the force correlation object. */
352 forceCorrelationGrid_ =
353 compat::make_unique<CorrelationGrid>(state_.points().size(), ndim(),
354 blockLength, CorrelationGrid::BlockLengthMeasure::Time,
355 awhParams.nstSampleCoord*mdTimeStep);
357 writer_ = compat::make_unique<BiasWriter>(*this);
361 void Bias::printInitializationToLog(FILE *fplog) const
363 if (fplog != nullptr && forceCorrelationGrid_ != nullptr)
366 gmx::formatString("\nawh%d:", params_.biasIndex + 1);
369 "%s initial force correlation block length = %g %s"
370 "%s force correlation number of blocks = %d",
371 prefix.c_str(), forceCorrelationGrid().getBlockLength(),
372 forceCorrelationGrid().blockLengthMeasure == CorrelationGrid::BlockLengthMeasure::Weight ? "" : "ps",
373 prefix.c_str(), forceCorrelationGrid().getNumBlocks());
377 void Bias::updateForceCorrelationGrid(gmx::ArrayRef<const double> probWeightNeighbor,
380 if (forceCorrelationGrid_ == nullptr)
385 const std::vector<int> &neighbor = grid_.point(state_.coordState().gridpointIndex()).neighbor;
387 gmx::ArrayRef<double> forceFromNeighbor = tempForce_;
388 for (size_t n = 0; n < neighbor.size(); n++)
390 double weightNeighbor = probWeightNeighbor[n];
391 int indexNeighbor = neighbor[n];
393 /* Add the force data of this neighbor point. Note: the sum of these forces is the convolved force.
395 We actually add the force normalized by beta which has the units of 1/length. This means that the
396 resulting correlation time integral is directly in units of friction time/length^2 which is really what
397 we're interested in. */
398 state_.calcUmbrellaForceAndPotential(dimParams_, grid_, indexNeighbor, forceFromNeighbor);
400 /* Note: we might want to give a whole list of data to add instead and have this loop in the data adding function */
401 forceCorrelationGrid_->addData(indexNeighbor, weightNeighbor, forceFromNeighbor, t);
405 /* Return the number of data blocks that have been prepared for writing. */
406 int Bias::numEnergySubblocksToWrite() const
408 GMX_RELEASE_ASSERT(writer_ != nullptr, "Should only request data from an initialized writer");
410 return writer_->numBlocks();
413 /* Write bias data blocks to energy subblocks. */
414 int Bias::writeToEnergySubblocks(t_enxsubblock *subblock) const
416 GMX_RELEASE_ASSERT(writer_ != nullptr, "Should only request data from an initialized writer");
418 return writer_->writeToEnergySubblocks(*this, subblock);