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,
 101                         const std::vector<int> &pullCoordIndex);
 102
 103     Bias                   bias_;           /**< The bias. */
 104     const std::vector<int> pullCoordIndex_; /**< The pull coordinates this bias acts on. */
 105
 106     /* Here AWH can be extended to work on other coordinates than pull. */
 107 };
 108
 109 BiasCoupledToSystem::BiasCoupledToSystem(Bias                    bias,
 110                                          const std::vector<int> &pullCoordIndex) :
 111     bias_(std::move(bias)),
 112     pullCoordIndex_(pullCoordIndex)
 113 {
 114     /* We already checked for this in grompp, but check again here. */
 115     GMX_RELEASE_ASSERT(static_cast<size_t>(bias_.ndim()) == pullCoordIndex_.size(), "The bias dimensionality should match the number of pull coordinates.");
 116 }
 117
 118 Awh::Awh(FILE                 *fplog,
 119          const t_inputrec     &inputRecord,
 120          const t_commrec      *commRecord,
 121          const gmx_multisim_t *multiSimRecord,
 122          const AwhParams      &awhParams,
 123          const std::string    &biasInitFilename,
 124          pull_t               *pull_work) :
 125     seed_(awhParams.seed),
 126     nstout_(awhParams.nstOut),
 127     commRecord_(commRecord),
 128     multiSimRecord_(multiSimRecord),
 129     pull_(pull_work),
 130     potentialOffset_(0)
 131 {
 132     /* We already checked for this in grompp, but check again here. */
 133     GMX_RELEASE_ASSERT(inputRecord.pull != nullptr, "With AWH we should have pull parameters");
 134     GMX_RELEASE_ASSERT(pull_work != nullptr, "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(InvalidInputError("Biases within a simulation are shared, currently sharing of biases is only supported between simulations"));
 145     }
 146
 147     int numSharingSimulations = 1;
 148     if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
 149     {
 150         numSharingSimulations = multiSimRecord_->nsim;
 151     }
 152
 153     /* Initialize all the biases */
 154     const double beta = 1/(BOLTZ*inputRecord.opts.ref_t[0]);
 155     for (int k = 0; k < awhParams.numBias; k++)
 156     {
 157         const AwhBiasParams    &awhBiasParams = awhParams.awhBiasParams[k];
 158
 159         std::vector<int>        pullCoordIndex;
 160         std::vector<DimParams>  dimParams;
 161         for (int d = 0; d < awhBiasParams.ndim; d++)
 162         {
 163             const AwhDimParams &awhDimParams      = awhBiasParams.dimParams[d];
 164             GMX_RELEASE_ASSERT(awhDimParams.eCoordProvider == eawhcoordproviderPULL, "Currently only the pull code is supported as coordinate provider");
 165             const t_pull_coord &pullCoord         = inputRecord.pull->coord[awhDimParams.coordIndex];
 166             GMX_RELEASE_ASSERT(pullCoord.eGeom != epullgDIRPBC, "Pull geometry 'direction-periodic' is not supported by AWH");
 167             double              conversionFactor  = pull_coordinate_is_angletype(&pullCoord) ? DEG2RAD : 1;
 168             dimParams.emplace_back(conversionFactor, awhDimParams.forceConstant, beta);
 169
 170             pullCoordIndex.push_back(awhDimParams.coordIndex);
 171         }
 172
 173         /* Construct the bias and couple it to the system. */
 174         Bias::ThisRankWillDoIO thisRankWillDoIO = (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
 175         biasCoupledToSystem_.emplace_back(Bias(k, awhParams, awhParams.awhBiasParams[k], dimParams, beta, inputRecord.delta_t, numSharingSimulations, biasInitFilename, thisRankWillDoIO),
 176                                           pullCoordIndex);
 177
 178         biasCoupledToSystem_.back().bias_.printInitializationToLog(fplog);
 179     }
 180
 181     /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
 182     registerAwhWithPull(awhParams, pull_);
 183
 184     if (numSharingSimulations > 1 && MASTER(commRecord_))
 185     {
 186         std::vector<size_t> pointSize;
 187         for (auto const &biasCts : biasCoupledToSystem_)
 188         {
 189             pointSize.push_back(biasCts.bias_.state().points().size());
 190         }
 191         /* Ensure that the shared biased are compatible between simulations */
 192         biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, multiSimRecord_);
 193     }
 194 }
 195
 196 Awh::~Awh() = default;
 197
 198 bool Awh::isOutputStep(int64_t step) const
 199 {
 200     return (nstout_ > 0 && step % nstout_ == 0);
 201 }
 202
 203 real Awh::applyBiasForcesAndUpdateBias(int                     ePBC,
 204                                        const t_mdatoms        &mdatoms,
 205                                        const matrix            box,
 206                                        gmx::ForceWithVirial   *forceWithVirial,
 207                                        double                  t,
 208                                        int64_t                 step,
 209                                        gmx_wallcycle          *wallcycle,
 210                                        FILE                   *fplog)
 211 {
 212     GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
 213
 214     wallcycle_start(wallcycle, ewcAWH);
 215
 216     t_pbc  pbc;
 217     set_pbc(&pbc, ePBC, box);
 218
 219     /* During the AWH update the potential can instantaneously jump due to either
 220        an bias update or moving the umbrella. The jumps are kept track of and
 221        subtracted from the potential in order to get a useful conserved energy quantity. */
 222     double awhPotential = potentialOffset_;
 223
 224     for (auto &biasCts : biasCoupledToSystem_)
 225     {
 226         /* Update the AWH coordinate values with those of the corresponding
 227          * pull coordinates.
 228          */
 229         awh_dvec coordValue = { 0, 0, 0, 0 };
 230         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 231         {
 232             coordValue[d] = get_pull_coord_value(pull_, biasCts.pullCoordIndex_[d], &pbc);
 233         }
 234
 235         /* Perform an AWH biasing step: this means, at regular intervals,
 236          * sampling observables based on the input pull coordinate value,
 237          * setting the bias force and/or updating the AWH bias state.
 238          */
 239         double              biasPotential;
 240         double              biasPotentialJump;
 241         /* Note: In the near future this call will be split in calls
 242          *       to supports bias sharing within a single simulation.
 243          */
 244         gmx::ArrayRef<const double> biasForce =
 245             biasCts.bias_.calcForceAndUpdateBias(coordValue,
 246                                                  &biasPotential, &biasPotentialJump,
 247                                                  commRecord_,
 248                                                  multiSimRecord_,
 249                                                  t, step, seed_, fplog);
 250
 251         awhPotential += biasPotential;
 252
 253         /* Keep track of the total potential shift needed to remove the potential jumps. */
 254         potentialOffset_ -= biasPotentialJump;
 255
 256         /* Communicate the bias force to the pull struct.
 257          * The bias potential is returned at the end of this function,
 258          * so that it can be added externally to the correct energy data block.
 259          */
 260         for (int d = 0; d < biasCts.bias_.ndim(); d++)
 261         {
 262             apply_external_pull_coord_force(pull_, biasCts.pullCoordIndex_[d],
 263                                             biasForce[d], &mdatoms,
 264                                             forceWithVirial);
 265         }
 266
 267         if (isOutputStep(step))
 268         {
 269             /* We might have skipped updates for part of the grid points.
 270              * Ensure all points are updated before writing out their data.
 271              */
 272             biasCts.bias_.doSkippedUpdatesForAllPoints();
 273         }
 274     }
 275
 276     wallcycle_stop(wallcycle, ewcAWH);
 277
 278     return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
 279 }
 280
 281 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
 282 {
 283     if (MASTER(commRecord_))
 284     {
 285         std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
 286         awhHistory->bias.clear();
 287         awhHistory->bias.resize(biasCoupledToSystem_.size());
 288
 289         for (size_t k = 0; k < awhHistory->bias.size(); k++)
 290         {
 291             biasCoupledToSystem_[k].bias_.initHistoryFromState(&awhHistory->bias[k]);
 292         }
 293
 294         return awhHistory;
 295     }
 296     else
 297     {
 298         /* Return an empty pointer */
 299         return std::shared_ptr<AwhHistory>();
 300     }
 301 }
 302
 303 void Awh::restoreStateFromHistory(const AwhHistory *awhHistory)
 304 {
 305     /* Restore the history to the current state */
 306     if (MASTER(commRecord_))
 307     {
 308         GMX_RELEASE_ASSERT(awhHistory != nullptr, "The master rank should have a valid awhHistory when restoring the state from history.");
 309
 310         if (awhHistory->bias.size() != biasCoupledToSystem_.size())
 311         {
 312             GMX_THROW(InvalidInputError("AWH state and history contain different numbers of biases. Likely you provided a checkpoint from a different simulation."));
 313         }
 314
 315         potentialOffset_ = awhHistory->potentialOffset;
 316     }
 317     if (PAR(commRecord_))
 318     {
 319         gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_);
 320     }
 321
 322     for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
 323     {
 324         biasCoupledToSystem_[k].bias_.restoreStateFromHistory(awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
 325     }
 326 }
 327
 328 void Awh::updateHistory(AwhHistory *awhHistory) const
 329 {
 330     if (!MASTER(commRecord_))
 331     {
 332         return;
 333     }
 334
 335     /* This assert will also catch a non-master rank calling this function. */
 336     GMX_RELEASE_ASSERT(awhHistory->bias.size() == biasCoupledToSystem_.size(), "AWH state and history bias count should match");
 337
 338     awhHistory->potentialOffset = potentialOffset_;
 339
 340     for (size_t k = 0; k < awhHistory->bias.size(); k++)
 341     {
 342         biasCoupledToSystem_[k].bias_.updateHistory(&awhHistory->bias[k]);
 343     }
 344 }
 345
 346 const char * Awh::externalPotentialString()
 347 {
 348     return "AWH";
 349 }
 350
 351 void Awh::registerAwhWithPull(const AwhParams &awhParams,
 352                               pull_t          *pull_work)
 353 {
 354     GMX_RELEASE_ASSERT(pull_work, "Need a valid pull object");
 355
 356     for (int k = 0; k < awhParams.numBias; k++)
 357     {
 358         const AwhBiasParams &biasParams = awhParams.awhBiasParams[k];
 359
 360         for (int d = 0; d < biasParams.ndim; d++)
 361         {
 362             register_external_pull_potential(pull_work, biasParams.dimParams[d].coordIndex, Awh::externalPotentialString());
 363         }
 364     }
 365 }
 366
 367 /* Fill the AWH data block of an energy frame with data (if there is any). */
 368 void Awh::writeToEnergyFrame(int64_t      step,
 369                              t_enxframe  *frame) const
 370 {
 371     GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
 372     GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
 373
 374     if (!isOutputStep(step))
 375     {
 376         /* This is not an AWH output step, don't write any AWH data */
 377         return;
 378     }
 379
 380     /* Get the total number of energy subblocks that AWH needs */
 381     int numSubblocks  = 0;
 382     for (auto &biasCoupledToSystem : biasCoupledToSystem_)
 383     {
 384         numSubblocks += biasCoupledToSystem.bias_.numEnergySubblocksToWrite();
 385     }
 386     GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
 387
 388     /* Add 1 energy block */
 389     add_blocks_enxframe(frame, frame->nblock + 1);
 390
 391     /* Take the block that was just added and set the number of subblocks. */
 392     t_enxblock *awhEnergyBlock = &(frame->block[frame->nblock - 1]);
 393     add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
 394
 395     /* Claim it as an AWH block. */
 396     awhEnergyBlock->id = enxAWH;
 397
 398     /* Transfer AWH data blocks to energy sub blocks */
 399     int energySubblockCount = 0;
 400     for (auto &biasCoupledToSystem : biasCoupledToSystem_)
 401     {
 402         energySubblockCount += biasCoupledToSystem.bias_.writeToEnergySubblocks(&(awhEnergyBlock->sub[energySubblockCount]));
 403     }
 404 }
 405
 406 std::unique_ptr<Awh>
 407 prepareAwhModule(FILE                 *fplog,
 408                  const t_inputrec     &inputRecord,
 409                  t_state              *stateGlobal,
 410                  const t_commrec      *commRecord,
 411                  const gmx_multisim_t *multiSimRecord,
 412                  const bool            startingFromCheckpoint,
 413                  const bool            usingShellParticles,
 414                  const std::string    &biasInitFilename,
 415                  pull_t               *pull_work)
 416 {
 417     if (!inputRecord.bDoAwh)
 418     {
 419         return nullptr;
 420     }
 421     if (usingShellParticles)
 422     {
 423         GMX_THROW(InvalidInputError("AWH biasing does not support shell particles."));
 424     }
 425
 426     auto awh = std::make_unique<Awh>(fplog, inputRecord, commRecord, multiSimRecord, *inputRecord.awhParams,
 427                                      biasInitFilename, pull_work);
 428
 429     if (startingFromCheckpoint)
 430     {
 431         // Restore the AWH history read from checkpoint
 432         awh->restoreStateFromHistory(MASTER(commRecord) ? stateGlobal->awhHistory.get() : nullptr);
 433     }
 434     else if (MASTER(commRecord))
 435     {
 436         // Initialize the AWH history here
 437         stateGlobal->awhHistory = awh->initHistoryFromState();
 438     }
 439     return awh;
 440 }
 441
 442 } // namespace gmx