src/gromacs/awh/awh.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2015,2016,2017,2018,2019, 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  * \ingroup module_awh
  43  */
  44
  45 #include "gmxpre.h"
  46
  47 #include "awh.h"
  48
  49 #include <cassert>
  50 #include <cmath>
  51 #include <cstdio>
  52 #include <cstdlib>
  53 #include <cstring>
  54
  55 #include <algorithm>
  56
  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"
  75
  76 #include "bias.h"
  77 #include "biassharing.h"
  78 #include "correlationgrid.h"
  79 #include "pointstate.h"
  80
  81 namespace gmx
  82 {
  83
  84 /*! \internal
  85  * \brief A bias and its coupling to the system.
  86  *
  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.
  92  */
  93 struct BiasCoupledToSystem
  94 {
  95     /*! \brief Constructor, couple a bias to a set of pull coordinates.
  96      *
  97      * \param[in] bias            The bias.
  98      * \param[in] pullCoordIndex  The pull coordinate indices.
  99      */
 100     BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex);
 101
 102     Bias                   bias_;           /**< The bias. */
 103     const std::vector<int> pullCoordIndex_; /**< The pull coordinates this bias acts on. */
 104
 105     /* Here AWH can be extended to work on other coordinates than pull. */
 106 };
 107
 108 BiasCoupledToSystem::BiasCoupledToSystem(Bias bias, const std::vector<int>& pullCoordIndex) :
 109     bias_(std::move(bias)),
 110     pullCoordIndex_(pullCoordIndex)
 111 {
 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.");
 115 }
 116
 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,
 123          pull_t*               pull_work) :
 124     seed_(awhParams.seed),
 125     nstout_(awhParams.nstOut),
 126     commRecord_(commRecord),
 127     multiSimRecord_(multiSimRecord),
 128     pull_(pull_work),
 129     potentialOffset_(0)
 130 {
 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");
 135
 136     if (fplog != nullptr)
 137     {
 138         please_cite(fplog, "Lindahl2014");
 139     }
 140
 141     if (haveBiasSharingWithinSimulation(awhParams))
 142     {
 143         /* This has likely been checked by grompp, but throw anyhow. */
 144         GMX_THROW(
 145                 InvalidInputError("Biases within a simulation are shared, currently sharing of "
 146                                   "biases is only supported between simulations"));
 147     }
 148
 149     int numSharingSimulations = 1;
 150     if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
 151     {
 152         numSharingSimulations = multiSimRecord_->nsim;
 153     }
 154
 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++)
 158     {
 159         const AwhBiasParams& awhBiasParams = awhParams.awhBiasParams[k];
 160
 161         std::vector<int>       pullCoordIndex;
 162         std::vector<DimParams> dimParams;
 163         for (int d = 0; d < awhBiasParams.ndim; d++)
 164         {
 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);
 173
 174             pullCoordIndex.push_back(awhDimParams.coordIndex);
 175         }
 176
 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),
 183                 pullCoordIndex);
 184
 185         biasCoupledToSystem_.back().bias_.printInitializationToLog(fplog);
 186     }
 187
 188     /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
 189     registerAwhWithPull(awhParams, pull_);
 190
 191     if (numSharingSimulations > 1 && MASTER(commRecord_))
 192     {
 193         std::vector<size_t> pointSize;
 194         for (auto const& biasCts : biasCoupledToSystem_)
 195         {
 196             pointSize.push_back(biasCts.bias_.state().points().size());
 197         }
 198         /* Ensure that the shared biased are compatible between simulations */
 199         biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, multiSimRecord_);
 200     }
 201 }
 202
 203 Awh::~Awh() = default;
 204
 205 bool Awh::isOutputStep(int64_t step) const
 206 {
 207     return (nstout_ > 0 && step % nstout_ == 0);
 208 }
 209
 210 real Awh::applyBiasForcesAndUpdateBias(int                   ePBC,
 211                                        const t_mdatoms&      mdatoms,
 212                                        const matrix          box,
 213                                        gmx::ForceWithVirial* forceWithVirial,
 214                                        double                t,
 215                                        int64_t               step,
 216                                        gmx_wallcycle*        wallcycle,
 217                                        FILE*                 fplog)
 218 {
 219     GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
 220
 221     wallcycle_start(wallcycle, ewcAWH);
 222
 223     t_pbc pbc;
 224     set_pbc(&pbc, ePBC, box);
 225
 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_;
 230
 231     for (auto& biasCts : biasCoupledToSystem_)
 232     {
 233         /* Update the AWH coordinate values with those of the corresponding
 234          * pull coordinates.
 235          */
 236         awh_dvec coordValue = { 0, 0, 0, 0 };
 237         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 238         {
 239             coordValue[d] = get_pull_coord_value(pull_, biasCts.pullCoordIndex_[d], &pbc);
 240         }
 241
 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.
 245          */
 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.
 250          */
 251         gmx::ArrayRef<const double> biasForce = biasCts.bias_.calcForceAndUpdateBias(
 252                 coordValue, &biasPotential, &biasPotentialJump, commRecord_, multiSimRecord_, t,
 253                 step, seed_, fplog);
 254
 255         awhPotential += biasPotential;
 256
 257         /* Keep track of the total potential shift needed to remove the potential jumps. */
 258         potentialOffset_ -= biasPotentialJump;
 259
 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.
 263          */
 264         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 265         {
 266             apply_external_pull_coord_force(pull_, biasCts.pullCoordIndex_[d], biasForce[d],
 267                                             &mdatoms, forceWithVirial);
 268         }
 269
 270         if (isOutputStep(step))
 271         {
 272             /* We might have skipped updates for part of the grid points.
 273              * Ensure all points are updated before writing out their data.
 274              */
 275             biasCts.bias_.doSkippedUpdatesForAllPoints();
 276         }
 277     }
 278
 279     wallcycle_stop(wallcycle, ewcAWH);
 280
 281     return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
 282 }
 283
 284 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
 285 {
 286     if (MASTER(commRecord_))
 287     {
 288         std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
 289         awhHistory->bias.clear();
 290         awhHistory->bias.resize(biasCoupledToSystem_.size());
 291
 292         for (size_t k = 0; k < awhHistory->bias.size(); k++)
 293         {
 294             biasCoupledToSystem_[k].bias_.initHistoryFromState(&awhHistory->bias[k]);
 295         }
 296
 297         return awhHistory;
 298     }
 299     else
 300     {
 301         /* Return an empty pointer */
 302         return std::shared_ptr<AwhHistory>();
 303     }
 304 }
 305
 306 void Awh::restoreStateFromHistory(const AwhHistory* awhHistory)
 307 {
 308     /* Restore the history to the current state */
 309     if (MASTER(commRecord_))
 310     {
 311         GMX_RELEASE_ASSERT(awhHistory != nullptr,
 312                            "The master rank should have a valid awhHistory when restoring the "
 313                            "state from history.");
 314
 315         if (awhHistory->bias.size() != biasCoupledToSystem_.size())
 316         {
 317             GMX_THROW(InvalidInputError(
 318                     "AWH state and history contain different numbers of biases. Likely you "
 319                     "provided a checkpoint from a different simulation."));
 320         }
 321
 322         potentialOffset_ = awhHistory->potentialOffset;
 323     }
 324     if (PAR(commRecord_))
 325     {
 326         gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_);
 327     }
 328
 329     for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
 330     {
 331         biasCoupledToSystem_[k].bias_.restoreStateFromHistory(
 332                 awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
 333     }
 334 }
 335
 336 void Awh::updateHistory(AwhHistory* awhHistory) const
 337 {
 338     if (!MASTER(commRecord_))
 339     {
 340         return;
 341     }
 342
 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");
 346
 347     awhHistory->potentialOffset = potentialOffset_;
 348
 349     for (size_t k = 0; k < awhHistory->bias.size(); k++)
 350     {
 351         biasCoupledToSystem_[k].bias_.updateHistory(&awhHistory->bias[k]);
 352     }
 353 }
 354
 355 const char* Awh::externalPotentialString()
 356 {
 357     return "AWH";
 358 }
 359
 360 void Awh::registerAwhWithPull(const AwhParams& awhParams, pull_t* pull_work)
 361 {
 362     GMX_RELEASE_ASSERT(pull_work, "Need a valid pull object");
 363
 364     for (int k = 0; k < awhParams.numBias; k++)
 365     {
 366         const AwhBiasParams& biasParams = awhParams.awhBiasParams[k];
 367
 368         for (int d = 0; d < biasParams.ndim; d++)
 369         {
 370             register_external_pull_potential(pull_work, biasParams.dimParams[d].coordIndex,
 371                                              Awh::externalPotentialString());
 372         }
 373     }
 374 }
 375
 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
 378 {
 379     GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
 380     GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
 381
 382     if (!isOutputStep(step))
 383     {
 384         /* This is not an AWH output step, don't write any AWH data */
 385         return;
 386     }
 387
 388     /* Get the total number of energy subblocks that AWH needs */
 389     int numSubblocks = 0;
 390     for (auto& biasCoupledToSystem : biasCoupledToSystem_)
 391     {
 392         numSubblocks += biasCoupledToSystem.bias_.numEnergySubblocksToWrite();
 393     }
 394     GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
 395
 396     /* Add 1 energy block */
 397     add_blocks_enxframe(frame, frame->nblock + 1);
 398
 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);
 402
 403     /* Claim it as an AWH block. */
 404     awhEnergyBlock->id = enxAWH;
 405
 406     /* Transfer AWH data blocks to energy sub blocks */
 407     int energySubblockCount = 0;
 408     for (auto& biasCoupledToSystem : biasCoupledToSystem_)
 409     {
 410         energySubblockCount += biasCoupledToSystem.bias_.writeToEnergySubblocks(
 411                 &(awhEnergyBlock->sub[energySubblockCount]));
 412     }
 413 }
 414
 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,
 423                                       pull_t*               pull_work)
 424 {
 425     if (!inputRecord.bDoAwh)
 426     {
 427         return nullptr;
 428     }
 429     if (usingShellParticles)
 430     {
 431         GMX_THROW(InvalidInputError("AWH biasing does not support shell particles."));
 432     }
 433
 434     auto awh = std::make_unique<Awh>(fplog, inputRecord, commRecord, multiSimRecord,
 435                                      *inputRecord.awhParams, biasInitFilename, pull_work);
 436
 437     if (startingFromCheckpoint)
 438     {
 439         // Restore the AWH history read from checkpoint
 440         awh->restoreStateFromHistory(MASTER(commRecord) ? stateGlobal->awhHistory.get() : nullptr);
 441     }
 442     else if (MASTER(commRecord))
 443     {
 444         // Initialize the AWH history here
 445         stateGlobal->awhHistory = awh->initHistoryFromState();
 446     }
 447     return awh;
 448 }
 449
 450 } // namespace gmx