src/gromacs/modularsimulator/energyelement.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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.
   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 /*! \internal \file
  36  * \brief Defines the microstate for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "energyelement.h"
  45
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/mdlib/compute_io.h"
  48 #include "gromacs/mdlib/enerdata_utils.h"
  49 #include "gromacs/mdlib/energyoutput.h"
  50 #include "gromacs/mdlib/mdatoms.h"
  51 #include "gromacs/mdlib/mdoutf.h"
  52 #include "gromacs/mdlib/stat.h"
  53 #include "gromacs/mdlib/update.h"
  54 #include "gromacs/mdrunutility/handlerestart.h"
  55 #include "gromacs/mdtypes/enerdata.h"
  56 #include "gromacs/mdtypes/energyhistory.h"
  57 #include "gromacs/mdtypes/inputrec.h"
  58 #include "gromacs/mdtypes/mdatom.h"
  59 #include "gromacs/mdtypes/observableshistory.h"
  60 #include "gromacs/mdtypes/pullhistory.h"
  61 #include "gromacs/mdtypes/state.h"
  62 #include "gromacs/topology/topology.h"
  63
  64 #include "freeenergyperturbationelement.h"
  65 #include "parrinellorahmanbarostat.h"
  66 #include "statepropagatordata.h"
  67 #include "vrescalethermostat.h"
  68
  69 struct pull_t;
  70 class t_state;
  71
  72 namespace gmx
  73 {
  74 class Awh;
  75
  76 EnergyElement::EnergyElement(StatePropagatorData*           statePropagatorData,
  77                              FreeEnergyPerturbationElement* freeEnergyPerturbationElement,
  78                              const gmx_mtop_t*              globalTopology,
  79                              const t_inputrec*              inputrec,
  80                              const MDAtoms*                 mdAtoms,
  81                              gmx_enerdata_t*                enerd,
  82                              gmx_ekindata_t*                ekind,
  83                              const Constraints*             constr,
  84                              FILE*                          fplog,
  85                              t_fcdata*                      fcd,
  86                              const MdModulesNotifier&       mdModulesNotifier,
  87                              bool                           isMasterRank,
  88                              ObservablesHistory*            observablesHistory,
  89                              StartingBehavior               startingBehavior) :
  90     isMasterRank_(isMasterRank),
  91     energyWritingStep_(-1),
  92     energyCalculationStep_(-1),
  93     freeEnergyCalculationStep_(-1),
  94     forceVirialStep_(-1),
  95     shakeVirialStep_(-1),
  96     totalVirialStep_(-1),
  97     pressureStep_(-1),
  98     needToSumEkinhOld_(false),
  99     startingBehavior_(startingBehavior),
 100     statePropagatorData_(statePropagatorData),
 101     freeEnergyPerturbationElement_(freeEnergyPerturbationElement),
 102     vRescaleThermostat_(nullptr),
 103     parrinelloRahmanBarostat_(nullptr),
 104     inputrec_(inputrec),
 105     top_global_(globalTopology),
 106     mdAtoms_(mdAtoms),
 107     enerd_(enerd),
 108     ekind_(ekind),
 109     constr_(constr),
 110     fplog_(fplog),
 111     fcd_(fcd),
 112     mdModulesNotifier_(mdModulesNotifier),
 113     groups_(&globalTopology->groups),
 114     observablesHistory_(observablesHistory)
 115 {
 116     clear_mat(forceVirial_);
 117     clear_mat(shakeVirial_);
 118     clear_mat(totalVirial_);
 119     clear_mat(pressure_);
 120     clear_rvec(muTot_);
 121
 122     if (freeEnergyPerturbationElement_)
 123     {
 124         dummyLegacyState_.flags = (1U << estFEPSTATE);
 125     }
 126 }
 127
 128 void EnergyElement::scheduleTask(Step step, Time time, const RegisterRunFunctionPtr& registerRunFunction)
 129 {
 130     if (!isMasterRank_)
 131     {
 132         return;
 133     }
 134     auto writeEnergy                 = energyWritingStep_ == step;
 135     auto isEnergyCalculationStep     = energyCalculationStep_ == step;
 136     auto isFreeEnergyCalculationStep = freeEnergyCalculationStep_ == step;
 137     if (isEnergyCalculationStep || writeEnergy)
 138     {
 139         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 140                 [this, time, isEnergyCalculationStep, isFreeEnergyCalculationStep]() {
 141                     doStep(time, isEnergyCalculationStep, isFreeEnergyCalculationStep);
 142                 }));
 143     }
 144     else
 145     {
 146         (*registerRunFunction)(std::make_unique<SimulatorRunFunction>(
 147                 [this]() { energyOutput_->recordNonEnergyStep(); }));
 148     }
 149 }
 150
 151 void EnergyElement::elementTeardown()
 152 {
 153     if (inputrec_->nstcalcenergy > 0 && isMasterRank_)
 154     {
 155         energyOutput_->printAverages(fplog_, groups_);
 156     }
 157 }
 158
 159 void EnergyElement::trajectoryWriterSetup(gmx_mdoutf* outf)
 160 {
 161     pull_t* pull_work = nullptr;
 162     energyOutput_ = std::make_unique<EnergyOutput>(mdoutf_get_fp_ene(outf), top_global_, inputrec_,
 163                                                    pull_work, mdoutf_get_fp_dhdl(outf), false,
 164                                                    startingBehavior_, mdModulesNotifier_);
 165
 166     if (!isMasterRank_)
 167     {
 168         return;
 169     }
 170
 171     initializeEnergyHistory(startingBehavior_, observablesHistory_, energyOutput_.get());
 172
 173     // TODO: This probably doesn't really belong here...
 174     //       but we have all we need in this element,
 175     //       so we'll leave it here for now!
 176     double io = compute_io(inputrec_, top_global_->natoms, *groups_, energyOutput_->numEnergyTerms(), 1);
 177     if ((io > 2000) && isMasterRank_)
 178     {
 179         fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io);
 180     }
 181     if (!inputrec_->bContinuation)
 182     {
 183         real temp = enerd_->term[F_TEMP];
 184         if (inputrec_->eI != eiVV)
 185         {
 186             /* Result of Ekin averaged over velocities of -half
 187              * and +half step, while we only have -half step here.
 188              */
 189             temp *= 2;
 190         }
 191         fprintf(fplog_, "Initial temperature: %g K\n", temp);
 192     }
 193 }
 194
 195 ITrajectoryWriterCallbackPtr EnergyElement::registerTrajectoryWriterCallback(TrajectoryEvent event)
 196 {
 197     if (event == TrajectoryEvent::EnergyWritingStep && isMasterRank_)
 198     {
 199         return std::make_unique<ITrajectoryWriterCallback>(
 200                 [this](gmx_mdoutf* mdoutf, Step step, Time time, bool writeTrajectory,
 201                        bool writeLog) { write(mdoutf, step, time, writeTrajectory, writeLog); });
 202     }
 203     return nullptr;
 204 }
 205
 206 SignallerCallbackPtr EnergyElement::registerTrajectorySignallerCallback(gmx::TrajectoryEvent event)
 207 {
 208     if (event == TrajectoryEvent::EnergyWritingStep && isMasterRank_)
 209     {
 210         return std::make_unique<SignallerCallback>(
 211                 [this](Step step, Time /*unused*/) { energyWritingStep_ = step; });
 212     }
 213     return nullptr;
 214 }
 215
 216 SignallerCallbackPtr EnergyElement::registerEnergyCallback(EnergySignallerEvent event)
 217 {
 218     if (event == EnergySignallerEvent::EnergyCalculationStep && isMasterRank_)
 219     {
 220         return std::make_unique<SignallerCallback>(
 221                 [this](Step step, Time /*unused*/) { energyCalculationStep_ = step; });
 222     }
 223     if (event == EnergySignallerEvent::FreeEnergyCalculationStep && isMasterRank_)
 224     {
 225         return std::make_unique<SignallerCallback>(
 226                 [this](Step step, Time /*unused*/) { freeEnergyCalculationStep_ = step; });
 227     }
 228     return nullptr;
 229 }
 230
 231 void EnergyElement::doStep(Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep)
 232 {
 233     enerd_->term[F_ETOT] = enerd_->term[F_EPOT] + enerd_->term[F_EKIN];
 234     if (vRescaleThermostat_)
 235     {
 236         dummyLegacyState_.therm_integral = vRescaleThermostat_->thermostatIntegral();
 237     }
 238     if (freeEnergyPerturbationElement_)
 239     {
 240         sum_dhdl(enerd_, freeEnergyPerturbationElement_->constLambdaView(), *inputrec_->fepvals);
 241         dummyLegacyState_.fep_state = freeEnergyPerturbationElement_->currentFEPState();
 242     }
 243     if (parrinelloRahmanBarostat_)
 244     {
 245         copy_mat(parrinelloRahmanBarostat_->boxVelocities(), dummyLegacyState_.boxv);
 246         copy_mat(statePropagatorData_->constBox(), dummyLegacyState_.box);
 247     }
 248     if (integratorHasConservedEnergyQuantity(inputrec_))
 249     {
 250         enerd_->term[F_ECONSERVED] =
 251                 enerd_->term[F_ETOT] + NPT_energy(inputrec_, &dummyLegacyState_, nullptr);
 252     }
 253     energyOutput_->addDataAtEnergyStep(isFreeEnergyCalculationStep, isEnergyCalculationStep, time,
 254                                        mdAtoms_->mdatoms()->tmass, enerd_, &dummyLegacyState_,
 255                                        inputrec_->fepvals, inputrec_->expandedvals,
 256                                        statePropagatorData_->constPreviousBox(), shakeVirial_,
 257                                        forceVirial_, totalVirial_, pressure_, ekind_, muTot_, constr_);
 258 }
 259
 260 void EnergyElement::write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog)
 261 {
 262     if (writeLog)
 263     {
 264         energyOutput_->printHeader(fplog_, step, time);
 265     }
 266
 267     bool do_dr = do_per_step(step, inputrec_->nstdisreout);
 268     bool do_or = do_per_step(step, inputrec_->nstorireout);
 269
 270     // energyOutput_->printAnnealingTemperatures(writeLog ? fplog_ : nullptr, groups_, &(inputrec_->opts));
 271     Awh* awh = nullptr;
 272     energyOutput_->printStepToEnergyFile(mdoutf_get_fp_ene(outf), writeTrajectory, do_dr, do_or,
 273                                          writeLog ? fplog_ : nullptr, step, time, fcd_, awh);
 274 }
 275
 276 void EnergyElement::addToForceVirial(const tensor virial, Step step)
 277 {
 278     if (step > forceVirialStep_)
 279     {
 280         forceVirialStep_ = step;
 281         clear_mat(forceVirial_);
 282     }
 283     m_add(forceVirial_, virial, forceVirial_);
 284 }
 285
 286 void EnergyElement::addToConstraintVirial(const tensor virial, Step step)
 287 {
 288     if (step > shakeVirialStep_)
 289     {
 290         shakeVirialStep_ = step;
 291         clear_mat(shakeVirial_);
 292     }
 293     m_add(shakeVirial_, virial, shakeVirial_);
 294 }
 295
 296 rvec* EnergyElement::forceVirial(Step gmx_unused step)
 297 {
 298     if (step > forceVirialStep_)
 299     {
 300         forceVirialStep_ = step;
 301         clear_mat(forceVirial_);
 302     }
 303     GMX_ASSERT(step >= forceVirialStep_ || forceVirialStep_ == -1,
 304                "Asked for force virial of previous step.");
 305     return forceVirial_;
 306 }
 307
 308 rvec* EnergyElement::constraintVirial(Step gmx_unused step)
 309 {
 310     if (step > shakeVirialStep_)
 311     {
 312         shakeVirialStep_ = step;
 313         clear_mat(shakeVirial_);
 314     }
 315     GMX_ASSERT(step >= shakeVirialStep_ || shakeVirialStep_ == -1,
 316                "Asked for constraint virial of previous step.");
 317     return shakeVirial_;
 318 }
 319
 320 rvec* EnergyElement::totalVirial(Step gmx_unused step)
 321 {
 322     if (step > totalVirialStep_)
 323     {
 324         totalVirialStep_ = step;
 325         clear_mat(totalVirial_);
 326     }
 327     GMX_ASSERT(step >= totalVirialStep_ || totalVirialStep_ == -1,
 328                "Asked for total virial of previous step.");
 329     return totalVirial_;
 330 }
 331
 332 rvec* EnergyElement::pressure(Step gmx_unused step)
 333 {
 334     if (step > pressureStep_)
 335     {
 336         pressureStep_ = step;
 337         clear_mat(pressure_);
 338     }
 339     GMX_ASSERT(step >= pressureStep_ || pressureStep_ == -1,
 340                "Asked for pressure of previous step.");
 341     return pressure_;
 342 }
 343
 344 real* EnergyElement::muTot()
 345 {
 346     return muTot_;
 347 }
 348
 349 gmx_enerdata_t* EnergyElement::enerdata()
 350 {
 351     return enerd_;
 352 }
 353
 354 gmx_ekindata_t* EnergyElement::ekindata()
 355 {
 356     return ekind_;
 357 }
 358
 359 bool* EnergyElement::needToSumEkinhOld()
 360 {
 361     return &needToSumEkinhOld_;
 362 }
 363
 364 void EnergyElement::writeCheckpoint(t_state gmx_unused* localState, t_state* globalState)
 365 {
 366     if (isMasterRank_)
 367     {
 368         if (needToSumEkinhOld_)
 369         {
 370             globalState->ekinstate.bUpToDate = false;
 371         }
 372         else
 373         {
 374             update_ekinstate(&globalState->ekinstate, ekind_);
 375             globalState->ekinstate.bUpToDate = true;
 376         }
 377         energyOutput_->fillEnergyHistory(observablesHistory_->energyHistory.get());
 378     }
 379 }
 380
 381 void EnergyElement::initializeEnergyHistory(StartingBehavior    startingBehavior,
 382                                             ObservablesHistory* observablesHistory,
 383                                             EnergyOutput*       energyOutput)
 384 {
 385     if (startingBehavior != StartingBehavior::NewSimulation)
 386     {
 387         /* Restore from energy history if appending to output files */
 388         if (startingBehavior == StartingBehavior::RestartWithAppending)
 389         {
 390             /* If no history is available (because a checkpoint is from before
 391              * it was written) make a new one later, otherwise restore it.
 392              */
 393             if (observablesHistory->energyHistory)
 394             {
 395                 energyOutput->restoreFromEnergyHistory(*observablesHistory->energyHistory);
 396             }
 397         }
 398         else if (observablesHistory->energyHistory)
 399         {
 400             /* We might have read an energy history from checkpoint.
 401              * As we are not appending, we want to restart the statistics.
 402              * Free the allocated memory and reset the counts.
 403              */
 404             observablesHistory->energyHistory = {};
 405             /* We might have read a pull history from checkpoint.
 406              * We will still want to keep the statistics, so that the files
 407              * can be joined and still be meaningful.
 408              * This means that observablesHistory_->pullHistory
 409              * should not be reset.
 410              */
 411         }
 412     }
 413     if (!observablesHistory->energyHistory)
 414     {
 415         observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
 416     }
 417     if (!observablesHistory->pullHistory)
 418     {
 419         observablesHistory->pullHistory = std::make_unique<PullHistory>();
 420     }
 421     /* Set the initial energy history */
 422     energyOutput->fillEnergyHistory(observablesHistory->energyHistory.get());
 423 }
 424
 425 void EnergyElement::setVRescaleThermostat(const gmx::VRescaleThermostat* vRescaleThermostat)
 426 {
 427     vRescaleThermostat_ = vRescaleThermostat;
 428     if (vRescaleThermostat_)
 429     {
 430         dummyLegacyState_.flags |= (1U << estTHERM_INT);
 431     }
 432 }
 433
 434 void EnergyElement::setParrinelloRahamnBarostat(const gmx::ParrinelloRahmanBarostat* parrinelloRahmanBarostat)
 435 {
 436     parrinelloRahmanBarostat_ = parrinelloRahmanBarostat;
 437     if (parrinelloRahmanBarostat_)
 438     {
 439         dummyLegacyState_.flags |= (1U << estBOX) | (1U << estBOXV);
 440     }
 441 }
 442
 443 } // namespace gmx