src/gromacs/applied_forces/awh/awh.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019,2020,2021, 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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35
  36 /*! \internal \file
  37  * \brief
  38  * Implements the Awh class.
  39  *
  40  * \author Viveca Lindahl
  41  * \author Berk Hess <hess@kth.se>
  42  * \author Magnus Lundborg
  43  * \ingroup module_awh
  44  */
  45
  46 #include "gmxpre.h"
  47
  48 #include "awh.h"
  49
  50 #include <cassert>
  51 #include <cmath>
  52 #include <cstdio>
  53 #include <cstdlib>
  54 #include <cstring>
  55
  56 #include <algorithm>
  57 #include <vector>
  58
  59 #include "gromacs/fileio/enxio.h"
  60 #include "gromacs/gmxlib/network.h"
  61 #include "gromacs/math/units.h"
  62 #include "gromacs/math/utilities.h"
  63 #include "gromacs/mdrunutility/multisim.h"
  64 #include "gromacs/mdtypes/awh_history.h"
  65 #include "gromacs/mdtypes/awh_params.h"
  66 #include "gromacs/mdtypes/commrec.h"
  67 #include "gromacs/mdtypes/forceoutput.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/pull_params.h"
  70 #include "gromacs/mdtypes/state.h"
  71 #include "gromacs/pbcutil/pbc.h"
  72 #include "gromacs/pulling/pull.h"
  73 #include "gromacs/timing/wallcycle.h"
  74 #include "gromacs/trajectory/energyframe.h"
  75 #include "gromacs/utility/arrayref.h"
  76 #include "gromacs/utility/exceptions.h"
  77 #include "gromacs/utility/gmxassert.h"
  78 #include "gromacs/utility/pleasecite.h"
  79
  80 #include "bias.h"
  81 #include "biassharing.h"
  82 #include "correlationgrid.h"
  83 #include "pointstate.h"
  84
  85 namespace gmx
  86 {
  87
  88 /*! \internal
  89  * \brief A bias and its coupling to the system.
  90  *
  91  * This struct is used to separate the bias machinery in the Bias class,
  92  * which should be independent from the reaction coordinate, from the
  93  * obtaining of the reaction coordinate values and passing the computed forces.
  94  * Currently the AWH method couples to the system by mapping each
  95  * AWH bias to a pull coordinate. This can easily be generalized here.
  96  */
  97 struct BiasCoupledToSystem
  98 {
  99     /*! \brief Constructor, couple a bias to a set of pull coordinates.
 100      *
 101      * \param[in] bias            The bias.
 102      * \param[in] pullCoordIndex  The pull coordinate indices.
 103      */
 104     BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex);
 105
 106     Bias                   bias_;           /**< The bias. */
 107     const std::vector<int> pullCoordIndex_; /**< The pull coordinates this bias acts on. */
 108
 109     /* Here AWH can be extended to work on other coordinates than pull. */
 110 };
 111
 112 /*! \brief Checks whether any dimension uses the given coordinate provider type.
 113  *
 114  * \param[in] awhBiasParams The bias params to check.
 115  * \param[in] awhCoordProvider The type of coordinate provider
 116  * \returns true if any dimension of the bias is linked to the given provider
 117  */
 118 static bool anyDimUsesProvider(const AwhBiasParams&            awhBiasParams,
 119                                const AwhCoordinateProviderType awhCoordProvider)
 120 {
 121     return std::any_of(awhBiasParams.dimParams().begin(),
 122                        awhBiasParams.dimParams().end(),
 123                        [&awhCoordProvider](const auto& awhDimParam) {
 124                            return awhDimParam.coordinateProvider() == awhCoordProvider;
 125                        });
 126 }
 127
 128 /*! \brief Checks whether any dimension uses the given coordinate provider type.
 129  *
 130  * \param[in] awhParams The AWH params to check.
 131  * \param[in] awhCoordProvider The type of coordinate provider
 132  * \returns true if any dimension of awh is linked to the given provider type.
 133  */
 134 static bool anyDimUsesProvider(const AwhParams& awhParams, const AwhCoordinateProviderType awhCoordProvider)
 135 {
 136     return std::any_of(awhParams.awhBiasParams().begin(),
 137                        awhParams.awhBiasParams().end(),
 138                        [&awhCoordProvider](const auto& awhBiasParam) {
 139                            return anyDimUsesProvider(awhBiasParam, awhCoordProvider);
 140                        });
 141 }
 142
 143 /*! \brief Checks whether any dimension uses pulling as a coordinate provider.
 144  *
 145  * \param[in] biasCoupledToSystem The AWH biases to check.
 146  * \returns true if any dimension of the provided biases is linked to pulling.
 147  */
 148 static bool anyDimUsesPull(const ArrayRef<BiasCoupledToSystem> biasCoupledToSystem)
 149 {
 150     return std::any_of(biasCoupledToSystem.begin(), biasCoupledToSystem.end(), [](const auto& biasCts) {
 151         return !biasCts.pullCoordIndex_.empty();
 152     });
 153 }
 154
 155 BiasCoupledToSystem::BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex) :
 156     bias_(std::move(bias)), pullCoordIndex_(pullCoordIndex)
 157 {
 158     /* We already checked for this in grompp, but check again here. */
 159     GMX_RELEASE_ASSERT(
 160             static_cast<size_t>(bias_.ndim()) == pullCoordIndex_.size() + bias_.hasFepLambdaDimension() ? 1 : 0,
 161             "The bias dimensionality should match the number of pull and lambda coordinates.");
 162 }
 163
 164 Awh::Awh(FILE*                 fplog,
 165          const t_inputrec&     inputRecord,
 166          const t_commrec*      commRecord,
 167          const gmx_multisim_t* multiSimRecord,
 168          const AwhParams&      awhParams,
 169          const std::string&    biasInitFilename,
 170          pull_t*               pull_work,
 171          int                   numFepLambdaStates,
 172          int                   fepLambdaState) :
 173     seed_(awhParams.seed()),
 174     nstout_(awhParams.nstout()),
 175     commRecord_(commRecord),
 176     pull_(pull_work),
 177     potentialOffset_(0),
 178     numFepLambdaStates_(numFepLambdaStates),
 179     fepLambdaState_(fepLambdaState)
 180 {
 181     if (anyDimUsesProvider(awhParams, AwhCoordinateProviderType::Pull))
 182     {
 183         GMX_RELEASE_ASSERT(inputRecord.pull != nullptr, "With AWH we should have pull parameters");
 184         GMX_RELEASE_ASSERT(pull_work != nullptr,
 185                            "With AWH pull should be initialized before initializing AWH");
 186     }
 187
 188     if (fplog != nullptr)
 189     {
 190         please_cite(fplog, "Lindahl2014");
 191
 192         if (anyDimUsesProvider(awhParams, AwhCoordinateProviderType::FreeEnergyLambda))
 193         {
 194             please_cite(fplog, "Lundborg2021");
 195         }
 196     }
 197
 198     if (haveBiasSharingWithinSimulation(awhParams))
 199     {
 200         /* This has likely been checked by grompp, but throw anyhow. */
 201         GMX_THROW(
 202                 InvalidInputError("Biases within a simulation are shared, currently sharing of "
 203                                   "biases is only supported between simulations"));
 204     }
 205
 206     if (awhParams.shareBiasMultisim() && multiSimRecord != nullptr)
 207     {
 208         GMX_RELEASE_ASSERT(commRecord, "Need a valid commRecord");
 209         biasSharing_ = std::make_unique<BiasSharing>(awhParams, *commRecord, multiSimRecord->mastersComm_);
 210         if (fplog)
 211         {
 212             for (int k = 0; k < awhParams.numBias(); k++)
 213             {
 214                 const int shareGroup = awhParams.awhBiasParams()[k].shareGroup();
 215                 if (shareGroup > 0)
 216                 {
 217                     fprintf(fplog,
 218                             "awh%d: bias with share group %d is shared between %d simulations\n",
 219                             1 + k,
 220                             shareGroup,
 221                             biasSharing_->numSharingSimulations(k));
 222                 }
 223                 else
 224                 {
 225                     fprintf(fplog, "awh%d: bias is not shared between simulations\n", 1 + k);
 226                 }
 227             }
 228         }
 229     }
 230
 231     /* Initialize all the biases */
 232     const double beta          = 1 / (gmx::c_boltz * inputRecord.opts.ref_t[0]);
 233     const auto&  awhBiasParams = awhParams.awhBiasParams();
 234     for (int k = 0; k < gmx::ssize(awhBiasParams); k++)
 235     {
 236         std::vector<int>       pullCoordIndex;
 237         std::vector<DimParams> dimParams;
 238         for (const auto& awhDimParam : awhBiasParams[k].dimParams())
 239         {
 240             if (awhDimParam.coordinateProvider() != AwhCoordinateProviderType::Pull
 241                 && awhDimParam.coordinateProvider() != AwhCoordinateProviderType::FreeEnergyLambda)
 242             {
 243                 GMX_THROW(
 244                         InvalidInputError("Currently only the pull code and lambda are supported "
 245                                           "as coordinate providers"));
 246             }
 247             if (awhDimParam.coordinateProvider() == AwhCoordinateProviderType::Pull)
 248             {
 249                 const t_pull_coord& pullCoord = inputRecord.pull->coord[awhDimParam.coordinateIndex()];
 250                 if (pullCoord.eGeom == PullGroupGeometry::DirectionPBC)
 251                 {
 252                     GMX_THROW(InvalidInputError(
 253                             "Pull geometry 'direction-periodic' is not supported by AWH"));
 254                 }
 255                 double conversionFactor = pull_conversion_factor_userinput2internal(pullCoord);
 256                 pullCoordIndex.push_back(awhDimParam.coordinateIndex());
 257                 dimParams.emplace_back(DimParams::pullDimParams(
 258                         conversionFactor, awhDimParam.forceConstant(), beta));
 259             }
 260             else
 261             {
 262                 dimParams.push_back(DimParams::fepLambdaDimParams(numFepLambdaStates_, beta));
 263             }
 264         }
 265
 266         /* Construct the bias and couple it to the system. */
 267         Bias::ThisRankWillDoIO thisRankWillDoIO =
 268                 (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
 269         biasCoupledToSystem_.emplace_back(Bias(k,
 270                                                awhParams,
 271                                                awhBiasParams[k],
 272                                                dimParams,
 273                                                beta,
 274                                                inputRecord.delta_t,
 275                                                biasSharing_.get(),
 276                                                biasInitFilename,
 277                                                thisRankWillDoIO),
 278                                           pullCoordIndex);
 279
 280         biasCoupledToSystem_.back().bias_.printInitializationToLog(fplog);
 281     }
 282
 283     /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
 284     registerAwhWithPull(awhParams, pull_);
 285
 286     if (biasSharing_ && MASTER(commRecord_))
 287     {
 288         std::vector<size_t> pointSize;
 289         for (auto const& biasCts : biasCoupledToSystem_)
 290         {
 291             pointSize.push_back(biasCts.bias_.state().points().size());
 292         }
 293         /* Ensure that the shared biased are compatible between simulations */
 294         biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, *biasSharing_);
 295     }
 296 }
 297
 298 Awh::~Awh() = default;
 299
 300 bool Awh::isOutputStep(int64_t step) const
 301 {
 302     return (nstout_ > 0 && step % nstout_ == 0);
 303 }
 304
 305 real Awh::applyBiasForcesAndUpdateBias(PbcType                pbcType,
 306                                        ArrayRef<const real>   masses,
 307                                        ArrayRef<const double> neighborLambdaEnergies,
 308                                        ArrayRef<const double> neighborLambdaDhdl,
 309                                        const matrix           box,
 310                                        gmx::ForceWithVirial*  forceWithVirial,
 311                                        double                 t,
 312                                        int64_t                step,
 313                                        gmx_wallcycle*         wallcycle,
 314                                        FILE*                  fplog)
 315 {
 316     if (anyDimUsesPull(biasCoupledToSystem_))
 317     {
 318         GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
 319     }
 320
 321     wallcycle_start(wallcycle, WallCycleCounter::Awh);
 322
 323     t_pbc pbc;
 324     set_pbc(&pbc, pbcType, box);
 325
 326     /* During the AWH update the potential can instantaneously jump due to either
 327        an bias update or moving the umbrella. The jumps are kept track of and
 328        subtracted from the potential in order to get a useful conserved energy quantity. */
 329     double awhPotential = potentialOffset_;
 330
 331     for (auto& biasCts : biasCoupledToSystem_)
 332     {
 333         /* Update the AWH coordinate values with those of the corresponding
 334          * pull coordinates.
 335          */
 336         awh_dvec coordValue           = { 0, 0, 0, 0 };
 337         int      numLambdaDimsCounted = 0;
 338         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 339         {
 340             if (biasCts.bias_.dimParams()[d].isPullDimension())
 341             {
 342                 coordValue[d] = get_pull_coord_value(
 343                         pull_, biasCts.pullCoordIndex_[d - numLambdaDimsCounted], pbc);
 344             }
 345             else
 346             {
 347                 coordValue[d] = fepLambdaState_;
 348                 numLambdaDimsCounted += 1;
 349             }
 350         }
 351
 352         /* Perform an AWH biasing step: this means, at regular intervals,
 353          * sampling observables based on the input pull coordinate value,
 354          * setting the bias force and/or updating the AWH bias state.
 355          */
 356         double biasPotential;
 357         double biasPotentialJump;
 358         /* Note: In the near future this call will be split in calls
 359          *       to supports bias sharing within a single simulation.
 360          */
 361         gmx::ArrayRef<const double> biasForce =
 362                 biasCts.bias_.calcForceAndUpdateBias(coordValue,
 363                                                      neighborLambdaEnergies,
 364                                                      neighborLambdaDhdl,
 365                                                      &biasPotential,
 366                                                      &biasPotentialJump,
 367                                                      t,
 368                                                      step,
 369                                                      seed_,
 370                                                      fplog);
 371
 372         awhPotential += biasPotential;
 373
 374         /* Keep track of the total potential shift needed to remove the potential jumps. */
 375         potentialOffset_ -= biasPotentialJump;
 376
 377         /* Communicate the bias force to the pull struct.
 378          * The bias potential is returned at the end of this function,
 379          * so that it can be added externally to the correct energy data block.
 380          */
 381         numLambdaDimsCounted = 0;
 382         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 383         {
 384             if (biasCts.bias_.dimParams()[d].isPullDimension())
 385             {
 386                 apply_external_pull_coord_force(pull_,
 387                                                 biasCts.pullCoordIndex_[d - numLambdaDimsCounted],
 388                                                 biasForce[d],
 389                                                 masses,
 390                                                 forceWithVirial);
 391             }
 392             else
 393             {
 394                 int umbrellaGridpointIndex = biasCts.bias_.state().coordState().umbrellaGridpoint();
 395                 fepLambdaState_ = biasCts.bias_.getGridCoordValue(umbrellaGridpointIndex)[d];
 396                 numLambdaDimsCounted += 1;
 397             }
 398         }
 399
 400         if (isOutputStep(step))
 401         {
 402             /* We might have skipped updates for part of the grid points.
 403              * Ensure all points are updated before writing out their data.
 404              */
 405             biasCts.bias_.doSkippedUpdatesForAllPoints();
 406         }
 407     }
 408
 409     wallcycle_stop(wallcycle, WallCycleCounter::Awh);
 410
 411     return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
 412 }
 413
 414 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
 415 {
 416     if (MASTER(commRecord_))
 417     {
 418         std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
 419         awhHistory->bias.clear();
 420         awhHistory->bias.resize(biasCoupledToSystem_.size());
 421
 422         for (size_t k = 0; k < awhHistory->bias.size(); k++)
 423         {
 424             biasCoupledToSystem_[k].bias_.initHistoryFromState(&awhHistory->bias[k]);
 425         }
 426
 427         return awhHistory;
 428     }
 429     else
 430     {
 431         /* Return an empty pointer */
 432         return std::shared_ptr<AwhHistory>();
 433     }
 434 }
 435
 436 void Awh::restoreStateFromHistory(const AwhHistory* awhHistory)
 437 {
 438     /* Restore the history to the current state */
 439     if (MASTER(commRecord_))
 440     {
 441         GMX_RELEASE_ASSERT(awhHistory != nullptr,
 442                            "The master rank should have a valid awhHistory when restoring the "
 443                            "state from history.");
 444
 445         if (awhHistory->bias.size() != biasCoupledToSystem_.size())
 446         {
 447             GMX_THROW(InvalidInputError(
 448                     "AWH state and history contain different numbers of biases. Likely you "
 449                     "provided a checkpoint from a different simulation."));
 450         }
 451
 452         potentialOffset_ = awhHistory->potentialOffset;
 453     }
 454     if (PAR(commRecord_))
 455     {
 456         gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_->mpi_comm_mygroup);
 457     }
 458
 459     for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
 460     {
 461         biasCoupledToSystem_[k].bias_.restoreStateFromHistory(
 462                 awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
 463     }
 464 }
 465
 466 void Awh::updateHistory(AwhHistory* awhHistory) const
 467 {
 468     if (!MASTER(commRecord_))
 469     {
 470         return;
 471     }
 472
 473     /* This assert will also catch a non-master rank calling this function. */
 474     GMX_RELEASE_ASSERT(awhHistory->bias.size() == biasCoupledToSystem_.size(),
 475                        "AWH state and history bias count should match");
 476
 477     awhHistory->potentialOffset = potentialOffset_;
 478
 479     for (size_t k = 0; k < awhHistory->bias.size(); k++)
 480     {
 481         biasCoupledToSystem_[k].bias_.updateHistory(&awhHistory->bias[k]);
 482     }
 483 }
 484
 485 const char* Awh::externalPotentialString()
 486 {
 487     return "AWH";
 488 }
 489
 490 void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work)
 491 {
 492     GMX_RELEASE_ASSERT(!anyDimUsesProvider(awhParams, AwhCoordinateProviderType::Pull) || pull_work,
 493                        "Need a valid pull object");
 494
 495     for (const auto& biasParam : awhParams.awhBiasParams())
 496     {
 497         for (const auto& dimParam : biasParam.dimParams())
 498         {
 499             if (dimParam.coordinateProvider() == AwhCoordinateProviderType::Pull)
 500             {
 501                 register_external_pull_potential(
 502                         pull_work, dimParam.coordinateIndex(), Awh::externalPotentialString());
 503             }
 504         }
 505     }
 506 }
 507
 508 /* Fill the AWH data block of an energy frame with data (if there is any). */
 509 void Awh::writeToEnergyFrame(int64_t step, t_enxframe* frame) const
 510 {
 511     GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
 512     GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
 513
 514     if (!isOutputStep(step))
 515     {
 516         /* This is not an AWH output step, don't write any AWH data */
 517         return;
 518     }
 519
 520     /* Get the total number of energy subblocks that AWH needs */
 521     int numSubblocks = 0;
 522     for (const auto& biasCoupledToSystem : biasCoupledToSystem_)
 523     {
 524         numSubblocks += biasCoupledToSystem.bias_.numEnergySubblocksToWrite();
 525     }
 526     GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
 527
 528     /* Add 1 energy block */
 529     add_blocks_enxframe(frame, frame->nblock + 1);
 530
 531     /* Take the block that was just added and set the number of subblocks. */
 532     t_enxblock* awhEnergyBlock = &(frame->block[frame->nblock - 1]);
 533     add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
 534
 535     /* Claim it as an AWH block. */
 536     awhEnergyBlock->id = enxAWH;
 537
 538     /* Transfer AWH data blocks to energy sub blocks */
 539     int energySubblockCount = 0;
 540     for (const auto& biasCoupledToSystem : biasCoupledToSystem_)
 541     {
 542         energySubblockCount += biasCoupledToSystem.bias_.writeToEnergySubblocks(
 543                 &(awhEnergyBlock->sub[energySubblockCount]));
 544     }
 545 }
 546
 547 bool Awh::hasFepLambdaDimension() const
 548 {
 549     return std::any_of(
 550             std::begin(biasCoupledToSystem_),
 551             std::end(biasCoupledToSystem_),
 552             [](const auto& coupledBias) { return coupledBias.bias_.hasFepLambdaDimension(); });
 553 }
 554
 555 bool Awh::needForeignEnergyDifferences(const int64_t step) const
 556 {
 557     /* If there is no FEP lambda dimension at all in any bias there will be no need for foreign
 558      * energy differences */
 559     if (!hasFepLambdaDimension())
 560     {
 561         return false;
 562     }
 563     if (step == 0)
 564     {
 565         return true;
 566     }
 567     /* Check whether the bias(es) that has/have a FEP lambda dimension should sample coordinates
 568      * this step. Since the biases may have different sampleCoordStep it is necessary to check
 569      * this combination. */
 570     return std::any_of(biasCoupledToSystem_.begin(), biasCoupledToSystem_.end(), [step](const auto& biasCts) {
 571         return biasCts.bias_.hasFepLambdaDimension() && biasCts.bias_.isSampleCoordStep(step);
 572     });
 573 }
 574
 575 std::unique_ptr<Awh> prepareAwhModule(FILE*                 fplog,
 576                                       const t_inputrec&     inputRecord,
 577                                       t_state*              stateGlobal,
 578                                       const t_commrec*      commRecord,
 579                                       const gmx_multisim_t* multiSimRecord,
 580                                       const bool            startingFromCheckpoint,
 581                                       const bool            usingShellParticles,
 582                                       const std::string&    biasInitFilename,
 583                                       pull_t*               pull_work)
 584 {
 585     if (!inputRecord.bDoAwh)
 586     {
 587         return nullptr;
 588     }
 589     if (usingShellParticles)
 590     {
 591         GMX_THROW(InvalidInputError("AWH biasing does not support shell particles."));
 592     }
 593
 594     auto awh = std::make_unique<Awh>(fplog,
 595                                      inputRecord,
 596                                      commRecord,
 597                                      multiSimRecord,
 598                                      *inputRecord.awhParams,
 599                                      biasInitFilename,
 600                                      pull_work,
 601                                      inputRecord.fepvals->n_lambda,
 602                                      inputRecord.fepvals->init_fep_state);
 603
 604     if (startingFromCheckpoint)
 605     {
 606         // Restore the AWH history read from checkpoint
 607         awh->restoreStateFromHistory(MASTER(commRecord) ? stateGlobal->awhHistory.get() : nullptr);
 608     }
 609     else if (MASTER(commRecord))
 610     {
 611         // Initialize the AWH history here
 612         stateGlobal->awhHistory = awh->initHistoryFromState();
 613     }
 614     return awh;
 615 }
 616
 617 } // namespace gmx