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.
38 * Implements the Awh class.
40 * \author Viveca Lindahl
41 * \author Berk Hess <hess@kth.se>
42 * \author Magnus Lundborg
58 #include "gromacs/fileio/enxio.h"
59 #include "gromacs/gmxlib/network.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/mdrunutility/multisim.h"
62 #include "gromacs/mdtypes/awh_history.h"
63 #include "gromacs/mdtypes/awh_params.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/forceoutput.h"
66 #include "gromacs/mdtypes/inputrec.h"
67 #include "gromacs/mdtypes/pull_params.h"
68 #include "gromacs/mdtypes/state.h"
69 #include "gromacs/pbcutil/pbc.h"
70 #include "gromacs/pulling/pull.h"
71 #include "gromacs/timing/wallcycle.h"
72 #include "gromacs/trajectory/energyframe.h"
73 #include "gromacs/utility/arrayref.h"
74 #include "gromacs/utility/exceptions.h"
75 #include "gromacs/utility/gmxassert.h"
76 #include "gromacs/utility/pleasecite.h"
79 #include "biassharing.h"
80 #include "correlationgrid.h"
81 #include "pointstate.h"
87 * \brief A bias and its coupling to the system.
89 * This struct is used to separate the bias machinery in the Bias class,
90 * which should be independent from the reaction coordinate, from the
91 * obtaining of the reaction coordinate values and passing the computed forces.
92 * Currently the AWH method couples to the system by mapping each
93 * AWH bias to a pull coordinate. This can easily be generalized here.
95 struct BiasCoupledToSystem
97 /*! \brief Constructor, couple a bias to a set of pull coordinates.
99 * \param[in] bias The bias.
100 * \param[in] pullCoordIndex The pull coordinate indices.
102 BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex);
104 Bias bias_; /**< The bias. */
105 const std::vector<int> pullCoordIndex_; /**< The pull coordinates this bias acts on. */
107 /* Here AWH can be extended to work on other coordinates than pull. */
110 /*! \brief Checks whether any dimension uses pulling as a coordinate provider.
112 * \param[in] awhBiasParams The bias params to check.
113 * \returns true if any dimension of the bias is linked to pulling.
115 static bool anyDimUsesPull(const AwhBiasParams& awhBiasParams)
117 for (int d = 0; d < awhBiasParams.ndim; d++)
119 const AwhDimParams& awhDimParams = awhBiasParams.dimParams[d];
120 if (awhDimParams.eCoordProvider == eawhcoordproviderPULL)
128 /*! \brief Checks whether any dimension uses pulling as a coordinate provider.
130 * \param[in] awhParams The AWH params to check.
131 * \returns true if any dimension of awh is linked to pulling.
133 static bool anyDimUsesPull(const AwhParams& awhParams)
135 for (int k = 0; k < awhParams.numBias; k++)
137 const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
138 if (anyDimUsesPull(awhBiasParams))
146 /*! \brief Checks whether any dimension uses pulling as a coordinate provider.
148 * \param[in] biasCoupledToSystem The AWH biases to check.
149 * \returns true if any dimension of the provided biases is linked to pulling.
151 static bool anyDimUsesPull(const ArrayRef<BiasCoupledToSystem> biasCoupledToSystem)
153 for (const auto& biasCts : biasCoupledToSystem)
155 if (!biasCts.pullCoordIndex_.empty())
163 BiasCoupledToSystem::BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex) :
164 bias_(std::move(bias)),
165 pullCoordIndex_(pullCoordIndex)
167 /* We already checked for this in grompp, but check again here. */
169 static_cast<size_t>(bias_.ndim()) == pullCoordIndex_.size() + bias_.hasFepLambdaDimension() ? 1 : 0,
170 "The bias dimensionality should match the number of pull and lambda coordinates.");
173 Awh::Awh(FILE* fplog,
174 const t_inputrec& inputRecord,
175 const t_commrec* commRecord,
176 const gmx_multisim_t* multiSimRecord,
177 const AwhParams& awhParams,
178 const std::string& biasInitFilename,
180 int numFepLambdaStates,
181 int fepLambdaState) :
182 seed_(awhParams.seed),
183 nstout_(awhParams.nstOut),
184 commRecord_(commRecord),
185 multiSimRecord_(multiSimRecord),
188 numFepLambdaStates_(numFepLambdaStates),
189 fepLambdaState_(fepLambdaState)
191 if (anyDimUsesPull(awhParams))
193 GMX_RELEASE_ASSERT(inputRecord.pull != nullptr, "With AWH we should have pull parameters");
194 GMX_RELEASE_ASSERT(pull_work != nullptr,
195 "With AWH pull should be initialized before initializing AWH");
198 if (fplog != nullptr)
200 please_cite(fplog, "Lindahl2014");
203 if (haveBiasSharingWithinSimulation(awhParams))
205 /* This has likely been checked by grompp, but throw anyhow. */
207 InvalidInputError("Biases within a simulation are shared, currently sharing of "
208 "biases is only supported between simulations"));
211 int numSharingSimulations = 1;
212 if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
214 numSharingSimulations = multiSimRecord_->numSimulations_;
217 /* Initialize all the biases */
218 const double beta = 1 / (BOLTZ * inputRecord.opts.ref_t[0]);
219 for (int k = 0; k < awhParams.numBias; k++)
221 const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
223 std::vector<int> pullCoordIndex;
224 std::vector<DimParams> dimParams;
225 for (int d = 0; d < awhBiasParams.ndim; d++)
227 const AwhDimParams& awhDimParams = awhBiasParams.dimParams[d];
228 if (awhDimParams.eCoordProvider != eawhcoordproviderPULL
229 && awhDimParams.eCoordProvider != eawhcoordproviderFREE_ENERGY_LAMBDA)
232 InvalidInputError("Currently only the pull code and lambda are supported "
233 "as coordinate providers"));
235 if (awhDimParams.eCoordProvider == eawhcoordproviderPULL)
237 const t_pull_coord& pullCoord = inputRecord.pull->coord[awhDimParams.coordIndex];
238 if (pullCoord.eGeom == epullgDIRPBC)
240 GMX_THROW(InvalidInputError(
241 "Pull geometry 'direction-periodic' is not supported by AWH"));
243 double conversionFactor = pull_coordinate_is_angletype(&pullCoord) ? DEG2RAD : 1;
244 pullCoordIndex.push_back(awhDimParams.coordIndex);
245 dimParams.emplace_back(conversionFactor, awhDimParams.forceConstant, beta);
249 dimParams.emplace_back(awhDimParams.forceConstant, beta, numFepLambdaStates_);
253 /* Construct the bias and couple it to the system. */
254 Bias::ThisRankWillDoIO thisRankWillDoIO =
255 (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
256 biasCoupledToSystem_.emplace_back(
257 Bias(k, awhParams, awhParams.awhBiasParams[k], dimParams, beta, inputRecord.delta_t,
258 numSharingSimulations, biasInitFilename, thisRankWillDoIO),
261 biasCoupledToSystem_.back().bias_.printInitializationToLog(fplog);
264 /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
265 registerAwhWithPull(awhParams, pull_);
267 if (numSharingSimulations > 1 && MASTER(commRecord_))
269 std::vector<size_t> pointSize;
270 for (auto const& biasCts : biasCoupledToSystem_)
272 pointSize.push_back(biasCts.bias_.state().points().size());
274 /* Ensure that the shared biased are compatible between simulations */
275 biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, multiSimRecord_);
279 Awh::~Awh() = default;
281 bool Awh::isOutputStep(int64_t step) const
283 return (nstout_ > 0 && step % nstout_ == 0);
286 real Awh::applyBiasForcesAndUpdateBias(PbcType pbcType,
288 ArrayRef<const double> neighborLambdaEnergies,
289 ArrayRef<const double> neighborLambdaDhdl,
291 gmx::ForceWithVirial* forceWithVirial,
294 gmx_wallcycle* wallcycle,
297 if (anyDimUsesPull(biasCoupledToSystem_))
299 GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
302 wallcycle_start(wallcycle, ewcAWH);
305 set_pbc(&pbc, pbcType, box);
307 /* During the AWH update the potential can instantaneously jump due to either
308 an bias update or moving the umbrella. The jumps are kept track of and
309 subtracted from the potential in order to get a useful conserved energy quantity. */
310 double awhPotential = potentialOffset_;
312 for (auto& biasCts : biasCoupledToSystem_)
314 /* Update the AWH coordinate values with those of the corresponding
317 awh_dvec coordValue = { 0, 0, 0, 0 };
318 int numLambdaDimsCounted = 0;
319 for (int d = 0; d < biasCts.bias_.ndim(); d++)
321 if (!biasCts.bias_.isFepLambdaDimension(d))
323 coordValue[d] = get_pull_coord_value(
324 pull_, biasCts.pullCoordIndex_[d - numLambdaDimsCounted], &pbc);
328 coordValue[d] = fepLambdaState_;
329 numLambdaDimsCounted += 1;
333 /* Perform an AWH biasing step: this means, at regular intervals,
334 * sampling observables based on the input pull coordinate value,
335 * setting the bias force and/or updating the AWH bias state.
337 double biasPotential;
338 double biasPotentialJump;
339 /* Note: In the near future this call will be split in calls
340 * to supports bias sharing within a single simulation.
342 gmx::ArrayRef<const double> biasForce = biasCts.bias_.calcForceAndUpdateBias(
343 coordValue, neighborLambdaEnergies, neighborLambdaDhdl, &biasPotential,
344 &biasPotentialJump, commRecord_, multiSimRecord_, t, step, seed_, fplog);
346 awhPotential += biasPotential;
348 /* Keep track of the total potential shift needed to remove the potential jumps. */
349 potentialOffset_ -= biasPotentialJump;
351 /* Communicate the bias force to the pull struct.
352 * The bias potential is returned at the end of this function,
353 * so that it can be added externally to the correct energy data block.
355 numLambdaDimsCounted = 0;
356 for (int d = 0; d < biasCts.bias_.ndim(); d++)
358 if (!biasCts.bias_.dimParams()[d].isFepLambdaDimension())
360 apply_external_pull_coord_force(pull_, biasCts.pullCoordIndex_[d - numLambdaDimsCounted],
361 biasForce[d], masses, forceWithVirial);
365 int umbrellaGridpointIndex = biasCts.bias_.state().coordState().umbrellaGridpoint();
366 fepLambdaState_ = biasCts.bias_.getGridCoordValue(umbrellaGridpointIndex)[d];
367 numLambdaDimsCounted += 1;
371 if (isOutputStep(step))
373 /* We might have skipped updates for part of the grid points.
374 * Ensure all points are updated before writing out their data.
376 biasCts.bias_.doSkippedUpdatesForAllPoints();
380 wallcycle_stop(wallcycle, ewcAWH);
382 return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
385 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
387 if (MASTER(commRecord_))
389 std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
390 awhHistory->bias.clear();
391 awhHistory->bias.resize(biasCoupledToSystem_.size());
393 for (size_t k = 0; k < awhHistory->bias.size(); k++)
395 biasCoupledToSystem_[k].bias_.initHistoryFromState(&awhHistory->bias[k]);
402 /* Return an empty pointer */
403 return std::shared_ptr<AwhHistory>();
407 void Awh::restoreStateFromHistory(const AwhHistory* awhHistory)
409 /* Restore the history to the current state */
410 if (MASTER(commRecord_))
412 GMX_RELEASE_ASSERT(awhHistory != nullptr,
413 "The master rank should have a valid awhHistory when restoring the "
414 "state from history.");
416 if (awhHistory->bias.size() != biasCoupledToSystem_.size())
418 GMX_THROW(InvalidInputError(
419 "AWH state and history contain different numbers of biases. Likely you "
420 "provided a checkpoint from a different simulation."));
423 potentialOffset_ = awhHistory->potentialOffset;
425 if (PAR(commRecord_))
427 gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_->mpi_comm_mygroup);
430 for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
432 biasCoupledToSystem_[k].bias_.restoreStateFromHistory(
433 awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
437 void Awh::updateHistory(AwhHistory* awhHistory) const
439 if (!MASTER(commRecord_))
444 /* This assert will also catch a non-master rank calling this function. */
445 GMX_RELEASE_ASSERT(awhHistory->bias.size() == biasCoupledToSystem_.size(),
446 "AWH state and history bias count should match");
448 awhHistory->potentialOffset = potentialOffset_;
450 for (size_t k = 0; k < awhHistory->bias.size(); k++)
452 biasCoupledToSystem_[k].bias_.updateHistory(&awhHistory->bias[k]);
456 const char* Awh::externalPotentialString()
461 void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work)
463 GMX_RELEASE_ASSERT(!anyDimUsesPull(awhParams) || pull_work, "Need a valid pull object");
465 for (int k = 0; k < awhParams.numBias; k++)
467 const AwhBiasParams& biasParams = awhParams.awhBiasParams[k];
469 for (int d = 0; d < biasParams.ndim; d++)
471 if (biasParams.dimParams[d].eCoordProvider == eawhcoordproviderPULL)
473 register_external_pull_potential(pull_work, biasParams.dimParams[d].coordIndex,
474 Awh::externalPotentialString());
480 /* Fill the AWH data block of an energy frame with data (if there is any). */
481 void Awh::writeToEnergyFrame(int64_t step, t_enxframe* frame) const
483 GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
484 GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
486 if (!isOutputStep(step))
488 /* This is not an AWH output step, don't write any AWH data */
492 /* Get the total number of energy subblocks that AWH needs */
493 int numSubblocks = 0;
494 for (auto& biasCoupledToSystem : biasCoupledToSystem_)
496 numSubblocks += biasCoupledToSystem.bias_.numEnergySubblocksToWrite();
498 GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
500 /* Add 1 energy block */
501 add_blocks_enxframe(frame, frame->nblock + 1);
503 /* Take the block that was just added and set the number of subblocks. */
504 t_enxblock* awhEnergyBlock = &(frame->block[frame->nblock - 1]);
505 add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
507 /* Claim it as an AWH block. */
508 awhEnergyBlock->id = enxAWH;
510 /* Transfer AWH data blocks to energy sub blocks */
511 int energySubblockCount = 0;
512 for (auto& biasCoupledToSystem : biasCoupledToSystem_)
514 energySubblockCount += biasCoupledToSystem.bias_.writeToEnergySubblocks(
515 &(awhEnergyBlock->sub[energySubblockCount]));
519 bool Awh::hasFepLambdaDimension() const
522 std::begin(biasCoupledToSystem_), std::end(biasCoupledToSystem_),
523 [](const auto& coupledBias) { return coupledBias.bias_.hasFepLambdaDimension(); });
526 bool Awh::needForeignEnergyDifferences(const int64_t step) const
528 /* If there is no FEP lambda dimension at all in any bias there will be no need for foreign
529 * energy differences */
530 if (!hasFepLambdaDimension())
538 /* Check whether the bias(es) that has/have a FEP lambda dimension should sample coordinates
539 * this step. Since the biases may have different sampleCoordStep it is necessary to check
540 * this combination. */
541 for (const auto& biasCts : biasCoupledToSystem_)
543 if (biasCts.bias_.hasFepLambdaDimension() && biasCts.bias_.isSampleCoordStep(step))
551 std::unique_ptr<Awh> prepareAwhModule(FILE* fplog,
552 const t_inputrec& inputRecord,
553 t_state* stateGlobal,
554 const t_commrec* commRecord,
555 const gmx_multisim_t* multiSimRecord,
556 const bool startingFromCheckpoint,
557 const bool usingShellParticles,
558 const std::string& biasInitFilename,
561 if (!inputRecord.bDoAwh)
565 if (usingShellParticles)
567 GMX_THROW(InvalidInputError("AWH biasing does not support shell particles."));
570 auto awh = std::make_unique<Awh>(
571 fplog, inputRecord, commRecord, multiSimRecord, *inputRecord.awhParams, biasInitFilename,
572 pull_work, inputRecord.fepvals->n_lambda, inputRecord.fepvals->init_fep_state);
574 if (startingFromCheckpoint)
576 // Restore the AWH history read from checkpoint
577 awh->restoreStateFromHistory(MASTER(commRecord) ? stateGlobal->awhHistory.get() : nullptr);
579 else if (MASTER(commRecord))
581 // Initialize the AWH history here
582 stateGlobal->awhHistory = awh->initHistoryFromState();