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, 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 <assert.h>
  50
  51 #include <cmath>
  52 #include <cstdio>
  53 #include <cstdlib>
  54 #include <cstring>
  55
  56 #include <algorithm>
  57
  58 #include "gromacs/fileio/enxio.h"
  59 #include "gromacs/gmxlib/network.h"
  60 #include "gromacs/math/units.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/utility/exceptions.h"
  72 #include "gromacs/utility/gmxassert.h"
  73 #include "gromacs/utility/pleasecite.h"
  74
  75 #include "bias.h"
  76 #include "biassharing.h"
  77 #include "correlationgrid.h"
  78 #include "pointstate.h"
  79
  80 namespace gmx
  81 {
  82
  83 /*! \internal
  84  * \brief A bias and its coupling to the system.
  85  *
  86  * This struct is used to separate the bias machinery in the Bias class,
  87  * which should be independent from the reaction coordinate, from the
  88  * obtaining of the reaction coordinate values and passing the computed forces.
  89  * Currently the AWH method couples to the system by mapping each
  90  * AWH bias to a pull coordinate. This can easily be generalized here.
  91  */
  92 struct BiasCoupledToSystem
  93 {
  94     /*! \brief Constructor, couple a bias to a set of pull coordinates.
  95      *
  96      * \param[in] bias            The bias.
  97      * \param[in] pullCoordIndex  The pull coordinate indices.
  98      */
  99     BiasCoupledToSystem(Bias                    bias,
 100                         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,
 109                                          const std::vector<int> &pullCoordIndex) :
 110     bias(std::move(bias)),
 111     pullCoordIndex(pullCoordIndex)
 112 {
 113     /* We already checked for this in grompp, but check again here. */
 114     GMX_RELEASE_ASSERT(static_cast<size_t>(bias.ndim()) == pullCoordIndex.size(), "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, "With AWH pull should be initialized before initializing AWH");
 134
 135     if (fplog != nullptr)
 136     {
 137         please_cite(fplog, "Lindahl2014");
 138     }
 139
 140     if (haveBiasSharingWithinSimulation(awhParams))
 141     {
 142         /* This has likely been checked by grompp, but throw anyhow. */
 143         GMX_THROW(InvalidInputError("Biases within a simulation are shared, currently sharing of biases is only supported between simulations"));
 144     }
 145
 146     int numSharingSimulations = 1;
 147     if (awhParams.shareBiasMultisim && isMultiSim(multiSimRecord_))
 148     {
 149         numSharingSimulations = multiSimRecord_->nsim;
 150     }
 151
 152     /* Initialize all the biases */
 153     const double beta = 1/(BOLTZ*inputRecord.opts.ref_t[0]);
 154     for (int k = 0; k < awhParams.numBias; k++)
 155     {
 156         const AwhBiasParams    &awhBiasParams = awhParams.awhBiasParams[k];
 157
 158         std::vector<int>        pullCoordIndex;
 159         std::vector<DimParams>  dimParams;
 160         for (int d = 0; d < awhBiasParams.ndim; d++)
 161         {
 162             const AwhDimParams &awhDimParams      = awhBiasParams.dimParams[d];
 163             GMX_RELEASE_ASSERT(awhDimParams.eCoordProvider == eawhcoordproviderPULL, "Currently only the pull code is supported as coordinate provider");
 164             const t_pull_coord &pullCoord         = inputRecord.pull->coord[awhDimParams.coordIndex];
 165             double              conversionFactor  = pull_coordinate_is_angletype(&pullCoord) ? DEG2RAD : 1;
 166             dimParams.push_back(DimParams(conversionFactor, awhDimParams.forceConstant, beta));
 167
 168             pullCoordIndex.push_back(awhDimParams.coordIndex);
 169         }
 170
 171         /* Construct the bias and couple it to the system. */
 172         Bias::ThisRankWillDoIO thisRankWillDoIO = (MASTER(commRecord_) ? Bias::ThisRankWillDoIO::Yes : Bias::ThisRankWillDoIO::No);
 173         biasCoupledToSystem_.emplace_back(Bias(k, awhParams, awhParams.awhBiasParams[k], dimParams, beta, inputRecord.delta_t, numSharingSimulations, biasInitFilename, thisRankWillDoIO),
 174                                           pullCoordIndex);
 175
 176         biasCoupledToSystem_.back().bias.printInitializationToLog(fplog);
 177     }
 178
 179     /* Need to register the AWH coordinates to be allowed to apply forces to the pull coordinates. */
 180     registerAwhWithPull(awhParams, pull_);
 181
 182     if (numSharingSimulations > 1 && MASTER(commRecord_))
 183     {
 184         std::vector<size_t> pointSize;
 185         for (auto const &biasCts : biasCoupledToSystem_)
 186         {
 187             pointSize.push_back(biasCts.bias.state().points().size());
 188         }
 189         /* Ensure that the shared biased are compatible between simulations */
 190         biasesAreCompatibleForSharingBetweenSimulations(awhParams, pointSize, multiSimRecord_);
 191     }
 192 }
 193
 194 Awh::~Awh() = default;
 195
 196 bool Awh::isOutputStep(gmx_int64_t step) const
 197 {
 198     return (nstout_ > 0 && step % nstout_ == 0);
 199 }
 200
 201 real Awh::applyBiasForcesAndUpdateBias(int                     ePBC,
 202                                        const t_mdatoms        &mdatoms,
 203                                        const matrix            box,
 204                                        gmx::ForceWithVirial   *forceWithVirial,
 205                                        double                  t,
 206                                        gmx_int64_t             step,
 207                                        gmx_wallcycle          *wallcycle,
 208                                        FILE                   *fplog)
 209 {
 210     GMX_ASSERT(forceWithVirial, "Need a valid ForceWithVirial object");
 211
 212     wallcycle_start(wallcycle, ewcAWH);
 213
 214     t_pbc  pbc;
 215     set_pbc(&pbc, ePBC, box);
 216
 217     /* During the AWH update the potential can instantaneously jump due to either
 218        an bias update or moving the umbrella. The jumps are kept track of and
 219        subtracted from the potential in order to get a useful conserved energy quantity. */
 220     double awhPotential = potentialOffset_;
 221
 222     for (auto &biasCts : biasCoupledToSystem_)
 223     {
 224         /* Update the AWH coordinate values with those of the corresponding
 225          * pull coordinates.
 226          */
 227         awh_dvec coordValue = { 0, 0, 0, 0 };
 228         for (int d = 0; d < biasCts.bias.ndim(); d++)
 229         {
 230             coordValue[d] = get_pull_coord_value(pull_, biasCts.pullCoordIndex[d], &pbc);
 231         }
 232
 233         /* Perform an AWH biasing step: this means, at regular intervals,
 234          * sampling observables based on the input pull coordinate value,
 235          * setting the bias force and/or updating the AWH bias state.
 236          */
 237         double              biasPotential;
 238         double              biasPotentialJump;
 239         /* Note: In the near future this call will be split in calls
 240          *       to supports bias sharing within a single simulation.
 241          */
 242         gmx::ArrayRef<const double> biasForce =
 243             biasCts.bias.calcForceAndUpdateBias(coordValue,
 244                                                 &biasPotential, &biasPotentialJump,
 245                                                 commRecord_,
 246                                                 multiSimRecord_,
 247                                                 t, step, seed_, fplog);
 248
 249         awhPotential += biasPotential;
 250
 251         /* Keep track of the total potential shift needed to remove the potential jumps. */
 252         potentialOffset_ -= biasPotentialJump;
 253
 254         /* Communicate the bias force to the pull struct.
 255          * The bias potential is returned at the end of this function,
 256          * so that it can be added externally to the correct energy data block.
 257          */
 258         for (int d = 0; d < biasCts.bias.ndim(); d++)
 259         {
 260             apply_external_pull_coord_force(pull_, biasCts.pullCoordIndex[d],
 261                                             biasForce[d], &mdatoms,
 262                                             forceWithVirial);
 263         }
 264
 265         if (isOutputStep(step))
 266         {
 267             /* We might have skipped updates for part of the grid points.
 268              * Ensure all points are updated before writing out their data.
 269              */
 270             biasCts.bias.doSkippedUpdatesForAllPoints();
 271         }
 272     }
 273
 274     wallcycle_stop(wallcycle, ewcAWH);
 275
 276     return MASTER(commRecord_) ? static_cast<real>(awhPotential) : 0;
 277 }
 278
 279 std::shared_ptr<AwhHistory> Awh::initHistoryFromState() const
 280 {
 281     if (MASTER(commRecord_))
 282     {
 283         std::shared_ptr<AwhHistory> awhHistory(new AwhHistory);
 284         awhHistory->bias.clear();
 285         awhHistory->bias.resize(biasCoupledToSystem_.size());
 286
 287         for (size_t k = 0; k < awhHistory->bias.size(); k++)
 288         {
 289             biasCoupledToSystem_[k].bias.initHistoryFromState(&awhHistory->bias[k]);
 290         }
 291
 292         return awhHistory;
 293     }
 294     else
 295     {
 296         /* Return an empty pointer */
 297         return std::shared_ptr<AwhHistory>();
 298     }
 299 }
 300
 301 void Awh::restoreStateFromHistory(const AwhHistory *awhHistory)
 302 {
 303     /* Restore the history to the current state */
 304     if (MASTER(commRecord_))
 305     {
 306         GMX_RELEASE_ASSERT(awhHistory != nullptr, "The master rank should have a valid awhHistory when restoring the state from history.");
 307
 308         if (awhHistory->bias.size() != biasCoupledToSystem_.size())
 309         {
 310             GMX_THROW(InvalidInputError("AWH state and history contain different numbers of biases. Likely you provided a checkpoint from a different simulation."));
 311         }
 312
 313         potentialOffset_ = awhHistory->potentialOffset;
 314     }
 315     if (PAR(commRecord_))
 316     {
 317         gmx_bcast(sizeof(potentialOffset_), &potentialOffset_, commRecord_);
 318     }
 319
 320     for (size_t k = 0; k < biasCoupledToSystem_.size(); k++)
 321     {
 322         biasCoupledToSystem_[k].bias.restoreStateFromHistory(awhHistory ? &awhHistory->bias[k] : nullptr, commRecord_);
 323     }
 324 }
 325
 326 void Awh::updateHistory(AwhHistory *awhHistory) const
 327 {
 328     if (!MASTER(commRecord_))
 329     {
 330         return;
 331     }
 332
 333     /* This assert will also catch a non-master rank calling this function. */
 334     GMX_RELEASE_ASSERT(awhHistory->bias.size() == biasCoupledToSystem_.size(), "AWH state and history bias count should match");
 335
 336     awhHistory->potentialOffset = potentialOffset_;
 337
 338     for (size_t k = 0; k < awhHistory->bias.size(); k++)
 339     {
 340         biasCoupledToSystem_[k].bias.updateHistory(&awhHistory->bias[k]);
 341     }
 342 }
 343
 344 const char * Awh::externalPotentialString()
 345 {
 346     return "AWH";
 347 }
 348
 349 void Awh::registerAwhWithPull(const AwhParams &awhParams,
 350                               pull_t          *pull_work)
 351 {
 352     GMX_RELEASE_ASSERT(pull_work, "Need a valid pull object");
 353
 354     for (int k = 0; k < awhParams.numBias; k++)
 355     {
 356         const AwhBiasParams &biasParams = awhParams.awhBiasParams[k];
 357
 358         for (int d = 0; d < biasParams.ndim; d++)
 359         {
 360             register_external_pull_potential(pull_work, biasParams.dimParams[d].coordIndex, Awh::externalPotentialString());
 361         }
 362     }
 363 }
 364
 365 /* Fill the AWH data block of an energy frame with data (if there is any). */
 366 void Awh::writeToEnergyFrame(gmx_int64_t  step,
 367                              t_enxframe  *frame) const
 368 {
 369     GMX_ASSERT(MASTER(commRecord_), "writeToEnergyFrame should only be called on the master rank");
 370     GMX_ASSERT(frame != nullptr, "Need a valid energy frame");
 371
 372     if (!isOutputStep(step))
 373     {
 374         /* This is not an AWH output step, don't write any AWH data */
 375         return;
 376     }
 377
 378     /* Get the total number of energy subblocks that AWH needs */
 379     int numSubblocks  = 0;
 380     for (auto &biasCoupledToSystem : biasCoupledToSystem_)
 381     {
 382         numSubblocks += biasCoupledToSystem.bias.numEnergySubblocksToWrite();
 383     }
 384     GMX_ASSERT(numSubblocks > 0, "We should always have data to write");
 385
 386     /* Add 1 energy block */
 387     add_blocks_enxframe(frame, frame->nblock + 1);
 388
 389     /* Take the block that was just added and set the number of subblocks. */
 390     t_enxblock *awhEnergyBlock = &(frame->block[frame->nblock - 1]);
 391     add_subblocks_enxblock(awhEnergyBlock, numSubblocks);
 392
 393     /* Claim it as an AWH block. */
 394     awhEnergyBlock->id = enxAWH;
 395
 396     /* Transfer AWH data blocks to energy sub blocks */
 397     int energySubblockCount = 0;
 398     for (auto &biasCoupledToSystem : biasCoupledToSystem_)
 399     {
 400         energySubblockCount += biasCoupledToSystem.bias.writeToEnergySubblocks(&(awhEnergyBlock->sub[energySubblockCount]));
 401     }
 402 }
 403
 404 } // namespace gmx