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>
57 #include "gromacs/fileio/enxio.h"
58 #include "gromacs/gmxlib/network.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/mdrunutility/multisim.h"
61 #include "gromacs/mdtypes/awh_history.h"
62 #include "gromacs/mdtypes/awh_params.h"
63 #include "gromacs/mdtypes/commrec.h"
64 #include "gromacs/mdtypes/forceoutput.h"
65 #include "gromacs/mdtypes/inputrec.h"
66 #include "gromacs/mdtypes/pull_params.h"
67 #include "gromacs/mdtypes/state.h"
68 #include "gromacs/pbcutil/pbc.h"
69 #include "gromacs/pulling/pull.h"
70 #include "gromacs/timing/wallcycle.h"
71 #include "gromacs/trajectory/energyframe.h"
72 #include "gromacs/utility/exceptions.h"
73 #include "gromacs/utility/gmxassert.h"
74 #include "gromacs/utility/pleasecite.h"
77 #include "biassharing.h"
78 #include "correlationgrid.h"
79 #include "pointstate.h"
85 * \brief A bias and its coupling to the system.
87 * This struct is used to separate the bias machinery in the Bias class,
88 * which should be independent from the reaction coordinate, from the
89 * obtaining of the reaction coordinate values and passing the computed forces.
90 * Currently the AWH method couples to the system by mapping each
91 * AWH bias to a pull coordinate. This can easily be generalized here.
93 struct BiasCoupledToSystem
95 /*! \brief Constructor, couple a bias to a set of pull coordinates.
97 * \param[in] bias The bias.
98 * \param[in] pullCoordIndex The pull coordinate indices.
100 BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex);
102 Bias bias_; /**< The bias. */
103 const std::vector<int> pullCoordIndex_; /**< The pull coordinates this bias acts on. */
105 /* Here AWH can be extended to work on other coordinates than pull. */
108 BiasCoupledToSystem::BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex) :
109 bias_(std::move(bias)),
110 pullCoordIndex_(pullCoordIndex)
112 /* We already checked for this in grompp, but check again here. */
113 GMX_RELEASE_ASSERT(static_cast<size_t>(bias_.ndim()) == pullCoordIndex_.size(),
114 "The bias dimensionality should match the number of pull coordinates.");
117 Awh::Awh(FILE* fplog,
118 const t_inputrec& inputRecord,
119 const t_commrec* commRecord,
120 const gmx_multisim_t* multiSimRecord,
121 const AwhParams& awhParams,
122 const std::string& biasInitFilename,
124 seed_(awhParams.seed),
125 nstout_(awhParams.nstOut),
126 commRecord_(commRecord),
127 multiSimRecord_(multiSimRecord),
131 /* We already checked for this in grompp, but check again here. */
132 GMX_RELEASE_ASSERT(inputRecord.pull != nullptr, "With AWH we should have pull parameters");
133 GMX_RELEASE_ASSERT(pull_work != nullptr,
134 "With AWH pull should be initialized before initializing AWH");
136 if (fplog != nullptr)
138 please_cite(fplog, "Lindahl2014");
141 if (haveBiasSharingWithinSimulation(awhParams))
143 /* This has likely been checked by grompp, but throw anyhow. */
145 InvalidInputError("Biases within a simulation are shared, currently sharing of "
146 "biases is only supported between simulations"));
149 int numSharingSimulations = 1;
150 if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
152 numSharingSimulations = multiSimRecord_->numSimulations_;
155 /* Initialize all the biases */
156 const double beta = 1 / (BOLTZ * inputRecord.opts.ref_t[0]);
157 for (int k = 0; k < awhParams.numBias; k++)
159 const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
161 std::vector<int> pullCoordIndex;
162 std::vector<DimParams> dimParams;
163 for (int d = 0; d < awhBiasParams.ndim; d++)
165 const AwhDimParams& awhDimParams = awhBiasParams.dimParams[d];
166 GMX_RELEASE_ASSERT(awhDimParams.eCoordProvider == eawhcoordproviderPULL,
167 "Currently only the pull code is supported as coordinate provider");
168 const t_pull_coord& pullCoord = inputRecord.pull->coord[awhDimParams.coordIndex];
169 GMX_RELEASE_ASSERT(pullCoord.eGeom != epullgDIRPBC,
170 "Pull geometry 'direction-periodic' is not supported by AWH");
171 double conversionFactor = pull_coordinate_is_angletype(&pullCoord) ? DEG2RAD : 1;
172 dimParams.emplace_back(conversionFactor, awhDimParams.forceConstant, beta);
174 pullCoordIndex.push_back(awhDimParams.coordIndex);
177 /* Construct the bias and couple it to the system. */
178 Bias::ThisRankWillDoIO thisRankWillDoIO =
179 (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
180 biasCoupledToSystem_.emplace_back(
181 Bias(k, awhParams, awhParams.awhBiasParams[k], dimParams, beta, inputRecord.delta_t,
182 numSharingSimulations, biasInitFilename, thisRankWillDoIO),
185 biasCoupledToSystem_.back().bias_.printInitializationToLog(fplog);
188 /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
189 registerAwhWithPull(awhParams, pull_);
191 if (numSharingSimulations > 1 && MASTER(commRecord_))
193 std::vector<size_t> pointSize;
194 for (auto const& biasCts : biasCoupledToSystem_)
196 pointSize.push_back(biasCts.bias_.state().points().size());
198 /* Ensure that the shared biased are compatible between simulations */
199 biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, multiSimRecord_);
203 Awh::~Awh() = default;
205 bool Awh::isOutputStep(int64_t step) const
207 return (nstout_ > 0 && step % nstout_ == 0);
210 real Awh::applyBiasForcesAndUpdateBias(PbcType pbcType,
213 gmx::ForceWithVirial* forceWithVirial,
216 gmx_wallcycle* wallcycle,
219 GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
221 wallcycle_start(wallcycle, ewcAWH);
224 set_pbc(&pbc, pbcType, box);
226 /* During the AWH update the potential can instantaneously jump due to either
227 an bias update or moving the umbrella. The jumps are kept track of and
228 subtracted from the potential in order to get a useful conserved energy quantity. */
229 double awhPotential = potentialOffset_;
231 for (auto& biasCts : biasCoupledToSystem_)
233 /* Update the AWH coordinate values with those of the corresponding
236 awh_dvec coordValue = { 0, 0, 0, 0 };
237 for (int d = 0; d < biasCts.bias_.ndim(); d++)
239 coordValue[d] = get_pull_coord_value(pull_, biasCts.pullCoordIndex_[d], &pbc);
242 /* Perform an AWH biasing step: this means, at regular intervals,
243 * sampling observables based on the input pull coordinate value,
244 * setting the bias force and/or updating the AWH bias state.
246 double biasPotential;
247 double biasPotentialJump;
248 /* Note: In the near future this call will be split in calls
249 * to supports bias sharing within a single simulation.
251 gmx::ArrayRef<const double> biasForce = biasCts.bias_.calcForceAndUpdateBias(
252 coordValue, &biasPotential, &biasPotentialJump, commRecord_, multiSimRecord_, t,
255 awhPotential += biasPotential;
257 /* Keep track of the total potential shift needed to remove the potential jumps. */
258 potentialOffset_ -= biasPotentialJump;
260 /* Communicate the bias force to the pull struct.
261 * The bias potential is returned at the end of this function,
262 * so that it can be added externally to the correct energy data block.
264 for (int d = 0; d < biasCts.bias_.ndim(); d++)
266 apply_external_pull_coord_force(pull_, biasCts.pullCoordIndex_[d], biasForce[d], masses,
270 if (isOutputStep(step))
272 /* We might have skipped updates for part of the grid points.
273 * Ensure all points are updated before writing out their data.
275 biasCts.bias_.doSkippedUpdatesForAllPoints();
279 wallcycle_stop(wallcycle, ewcAWH);
281 return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
284 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
286 if (MASTER(commRecord_))
288 std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
289 awhHistory->bias.clear();
290 awhHistory->bias.resize(biasCoupledToSystem_.size());
292 for (size_t k = 0; k < awhHistory->bias.size(); k++)
294 biasCoupledToSystem_[k].bias_.initHistoryFromState(&awhHistory->bias[k]);
301 /* Return an empty pointer */
302 return std::shared_ptr<AwhHistory>();
306 void Awh::restoreStateFromHistory(const AwhHistory* awhHistory)
308 /* Restore the history to the current state */
309 if (MASTER(commRecord_))
311 GMX_RELEASE_ASSERT(awhHistory != nullptr,
312 "The master rank should have a valid awhHistory when restoring the "
313 "state from history.");
315 if (awhHistory->bias.size() != biasCoupledToSystem_.size())
317 GMX_THROW(InvalidInputError(
318 "AWH state and history contain different numbers of biases. Likely you "
319 "provided a checkpoint from a different simulation."));
322 potentialOffset_ = awhHistory->potentialOffset;
324 if (PAR(commRecord_))
326 gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_->mpi_comm_mygroup);
329 for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
331 biasCoupledToSystem_[k].bias_.restoreStateFromHistory(
332 awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
336 void Awh::updateHistory(AwhHistory* awhHistory) const
338 if (!MASTER(commRecord_))
343 /* This assert will also catch a non-master rank calling this function. */
344 GMX_RELEASE_ASSERT(awhHistory->bias.size() == biasCoupledToSystem_.size(),
345 "AWH state and history bias count should match");
347 awhHistory->potentialOffset = potentialOffset_;
349 for (size_t k = 0; k < awhHistory->bias.size(); k++)
351 biasCoupledToSystem_[k].bias_.updateHistory(&awhHistory->bias[k]);
355 const char* Awh::externalPotentialString()
360 void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work)
362 GMX_RELEASE_ASSERT(pull_work, "Need a valid pull object");
364 for (int k = 0; k < awhParams.numBias; k++)
366 const AwhBiasParams& biasParams = awhParams.awhBiasParams[k];
368 for (int d = 0; d < biasParams.ndim; d++)
370 register_external_pull_potential(pull_work, biasParams.dimParams[d].coordIndex,
371 Awh::externalPotentialString());
376 /* Fill the AWH data block of an energy frame with data (if there is any). */
377 void Awh::writeToEnergyFrame(int64_t step, t_enxframe* frame) const
379 GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
380 GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
382 if (!isOutputStep(step))
384 /* This is not an AWH output step, don't write any AWH data */
388 /* Get the total number of energy subblocks that AWH needs */
389 int numSubblocks = 0;
390 for (auto& biasCoupledToSystem : biasCoupledToSystem_)
392 numSubblocks += biasCoupledToSystem.bias_.numEnergySubblocksToWrite();
394 GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
396 /* Add 1 energy block */
397 add_blocks_enxframe(frame, frame->nblock + 1);
399 /* Take the block that was just added and set the number of subblocks. */
400 t_enxblock* awhEnergyBlock = &(frame->block[frame->nblock - 1]);
401 add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
403 /* Claim it as an AWH block. */
404 awhEnergyBlock->id = enxAWH;
406 /* Transfer AWH data blocks to energy sub blocks */
407 int energySubblockCount = 0;
408 for (auto& biasCoupledToSystem : biasCoupledToSystem_)
410 energySubblockCount += biasCoupledToSystem.bias_.writeToEnergySubblocks(
411 &(awhEnergyBlock->sub[energySubblockCount]));
415 std::unique_ptr<Awh> prepareAwhModule(FILE* fplog,
416 const t_inputrec& inputRecord,
417 t_state* stateGlobal,
418 const t_commrec* commRecord,
419 const gmx_multisim_t* multiSimRecord,
420 const bool startingFromCheckpoint,
421 const bool usingShellParticles,
422 const std::string& biasInitFilename,
425 if (!inputRecord.bDoAwh)
429 if (usingShellParticles)
431 GMX_THROW(InvalidInputError("AWH biasing does not support shell particles."));
434 auto awh = std::make_unique<Awh>(fplog, inputRecord, commRecord, multiSimRecord,
435 *inputRecord.awhParams, biasInitFilename, pull_work);
437 if (startingFromCheckpoint)
439 // Restore the AWH history read from checkpoint
440 awh->restoreStateFromHistory(MASTER(commRecord) ? stateGlobal->awhHistory.get() : nullptr);
442 else if (MASTER(commRecord))
444 // Initialize the AWH history here
445 stateGlobal->awhHistory = awh->initHistoryFromState();