src/gromacs/modularsimulator/energydata.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 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 /*! \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 "energydata.h"
  45
  46 #include "gromacs/gmxlib/network.h"
  47 #include "gromacs/math/vec.h"
  48 #include "gromacs/mdlib/compute_io.h"
  49 #include "gromacs/mdlib/enerdata_utils.h"
  50 #include "gromacs/mdlib/energyoutput.h"
  51 #include "gromacs/mdlib/mdatoms.h"
  52 #include "gromacs/mdlib/mdoutf.h"
  53 #include "gromacs/mdlib/stat.h"
  54 #include "gromacs/mdlib/update.h"
  55 #include "gromacs/mdtypes/commrec.h"
  56 #include "gromacs/mdtypes/enerdata.h"
  57 #include "gromacs/mdtypes/energyhistory.h"
  58 #include "gromacs/mdtypes/inputrec.h"
  59 #include "gromacs/mdtypes/mdatom.h"
  60 #include "gromacs/mdtypes/observableshistory.h"
  61 #include "gromacs/mdtypes/pullhistory.h"
  62 #include "gromacs/topology/topology.h"
  63
  64 #include "freeenergyperturbationdata.h"
  65 #include "modularsimulator.h"
  66 #include "simulatoralgorithm.h"
  67 #include "statepropagatordata.h"
  68
  69 struct pull_t;
  70 class t_state;
  71
  72 namespace gmx
  73 {
  74 class Awh;
  75
  76 EnergyData::EnergyData(StatePropagatorData*        statePropagatorData,
  77                        FreeEnergyPerturbationData* freeEnergyPerturbationData,
  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 MDModulesNotifiers&   mdModulesNotifiers,
  87                        bool                        isMasterRank,
  88                        ObservablesHistory*         observablesHistory,
  89                        StartingBehavior            startingBehavior,
  90                        bool                        simulationsShareState) :
  91     element_(std::make_unique<Element>(this, isMasterRank)),
  92     isMasterRank_(isMasterRank),
  93     forceVirialStep_(-1),
  94     shakeVirialStep_(-1),
  95     totalVirialStep_(-1),
  96     pressureStep_(-1),
  97     needToSumEkinhOld_(false),
  98     hasReadEkinFromCheckpoint_(false),
  99     startingBehavior_(startingBehavior),
 100     statePropagatorData_(statePropagatorData),
 101     freeEnergyPerturbationData_(freeEnergyPerturbationData),
 102     inputrec_(inputrec),
 103     top_global_(globalTopology),
 104     mdAtoms_(mdAtoms),
 105     enerd_(enerd),
 106     ekind_(ekind),
 107     constr_(constr),
 108     fplog_(fplog),
 109     fcd_(fcd),
 110     mdModulesNotifiers_(mdModulesNotifiers),
 111     groups_(&globalTopology.groups),
 112     observablesHistory_(observablesHistory),
 113     simulationsShareState_(simulationsShareState)
 114 {
 115     clear_mat(forceVirial_);
 116     clear_mat(shakeVirial_);
 117     clear_mat(totalVirial_);
 118     clear_mat(pressure_);
 119     clear_rvec(muTot_);
 120
 121     init_ekinstate(&ekinstate_, inputrec_);
 122     observablesHistory_->energyHistory = std::make_unique<energyhistory_t>();
 123 }
 124
 125 void EnergyData::Element::scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction)
 126 {
 127     if (!isMasterRank_)
 128     {
 129         return;
 130     }
 131     auto writeEnergy                 = energyWritingStep_ == step;
 132     auto isEnergyCalculationStep     = energyCalculationStep_ == step;
 133     auto isFreeEnergyCalculationStep = freeEnergyCalculationStep_ == step;
 134     if (isEnergyCalculationStep || writeEnergy)
 135     {
 136         registerRunFunction([this, step, time, isEnergyCalculationStep, isFreeEnergyCalculationStep]() {
 137             energyData_->doStep(step, time, isEnergyCalculationStep, isFreeEnergyCalculationStep);
 138         });
 139     }
 140     else
 141     {
 142         registerRunFunction([this]() { energyData_->energyOutput_->recordNonEnergyStep(); });
 143     }
 144 }
 145
 146 void EnergyData::teardown()
 147 {
 148     if (inputrec_->nstcalcenergy > 0 && isMasterRank_)
 149     {
 150         energyOutput_->printEnergyConservation(fplog_, inputrec_->simulation_part, EI_MD(inputrec_->eI));
 151         energyOutput_->printAverages(fplog_, groups_);
 152     }
 153 }
 154
 155 void EnergyData::Element::trajectoryWriterSetup(gmx_mdoutf* outf)
 156 {
 157     energyData_->setup(outf);
 158 }
 159
 160 void EnergyData::setup(gmx_mdoutf* outf)
 161 {
 162     pull_t* pull_work = nullptr;
 163     energyOutput_     = std::make_unique<EnergyOutput>(mdoutf_get_fp_ene(outf),
 164                                                    top_global_,
 165                                                    *inputrec_,
 166                                                    pull_work,
 167                                                    mdoutf_get_fp_dhdl(outf),
 168                                                    false,
 169                                                    startingBehavior_,
 170                                                    simulationsShareState_,
 171                                                    mdModulesNotifiers_);
 172
 173     if (!isMasterRank_)
 174     {
 175         return;
 176     }
 177
 178     initializeEnergyHistory(startingBehavior_, observablesHistory_, energyOutput_.get());
 179
 180     // TODO: This probably doesn't really belong here...
 181     //       but we have all we need in this element,
 182     //       so we'll leave it here for now!
 183     double io = compute_io(inputrec_, top_global_.natoms, *groups_, energyOutput_->numEnergyTerms(), 1);
 184     if ((io > 2000) && isMasterRank_)
 185     {
 186         fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io);
 187     }
 188     if (!inputrec_->bContinuation)
 189     {
 190         real temp = enerd_->term[F_TEMP];
 191         if (inputrec_->eI != IntegrationAlgorithm::VV)
 192         {
 193             /* Result of Ekin averaged over velocities of -half
 194              * and +half step, while we only have -half step here.
 195              */
 196             temp *= 2;
 197         }
 198         fprintf(fplog_, "Initial temperature: %g K\n", temp);
 199     }
 200 }
 201
 202 std::optional<ITrajectoryWriterCallback> EnergyData::Element::registerTrajectoryWriterCallback(TrajectoryEvent event)
 203 {
 204     if (event == TrajectoryEvent::EnergyWritingStep && isMasterRank_)
 205     {
 206         return [this](gmx_mdoutf* mdoutf, Step step, Time time, bool writeTrajectory, bool writeLog) {
 207             energyData_->write(mdoutf, step, time, writeTrajectory, writeLog);
 208         };
 209     }
 210     return std::nullopt;
 211 }
 212
 213 std::optional<SignallerCallback> EnergyData::Element::registerTrajectorySignallerCallback(gmx::TrajectoryEvent event)
 214 {
 215     if (event == TrajectoryEvent::EnergyWritingStep && isMasterRank_)
 216     {
 217         return [this](Step step, Time /*unused*/) { energyWritingStep_ = step; };
 218     }
 219     return std::nullopt;
 220 }
 221
 222 std::optional<SignallerCallback> EnergyData::Element::registerEnergyCallback(EnergySignallerEvent event)
 223 {
 224     if (event == EnergySignallerEvent::EnergyCalculationStep && isMasterRank_)
 225     {
 226         return [this](Step step, Time /*unused*/) { energyCalculationStep_ = step; };
 227     }
 228     if (event == EnergySignallerEvent::FreeEnergyCalculationStep && isMasterRank_)
 229     {
 230         return [this](Step step, Time /*unused*/) { freeEnergyCalculationStep_ = step; };
 231     }
 232     return std::nullopt;
 233 }
 234
 235 void EnergyData::doStep(Step step, Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep)
 236 {
 237     enerd_->term[F_ETOT] = enerd_->term[F_EPOT] + enerd_->term[F_EKIN];
 238     if (freeEnergyPerturbationData_)
 239     {
 240         accumulateKineticLambdaComponents(
 241                 enerd_, freeEnergyPerturbationData_->constLambdaView(), *inputrec_->fepvals);
 242     }
 243     if (integratorHasConservedEnergyQuantity(inputrec_))
 244     {
 245         enerd_->term[F_ECONSERVED] = enerd_->term[F_ETOT];
 246         for (const auto& energyContibution : conservedEnergyContributions_)
 247         {
 248             enerd_->term[F_ECONSERVED] += energyContibution(step, time);
 249         }
 250     }
 251     matrix nullMatrix = {};
 252     energyOutput_->addDataAtEnergyStep(
 253             isFreeEnergyCalculationStep,
 254             isEnergyCalculationStep,
 255             time,
 256             mdAtoms_->mdatoms()->tmass,
 257             enerd_,
 258             inputrec_->fepvals.get(),
 259             inputrec_->expandedvals.get(),
 260             statePropagatorData_->constPreviousBox(),
 261             PTCouplingArrays({ parrinelloRahmanBoxVelocities_ ? parrinelloRahmanBoxVelocities_() : nullMatrix,
 262                                {},
 263                                {},
 264                                {},
 265                                {} }),
 266             freeEnergyPerturbationData_ ? freeEnergyPerturbationData_->currentFEPState() : 0,
 267             shakeVirial_,
 268             forceVirial_,
 269             totalVirial_,
 270             pressure_,
 271             ekind_,
 272             muTot_,
 273             constr_);
 274 }
 275
 276 void EnergyData::write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog)
 277 {
 278     if (writeLog)
 279     {
 280         energyOutput_->printHeader(fplog_, step, time);
 281     }
 282
 283     bool do_dr = do_per_step(step, inputrec_->nstdisreout);
 284     bool do_or = do_per_step(step, inputrec_->nstorireout);
 285
 286     // energyOutput_->printAnnealingTemperatures(writeLog ? fplog_ : nullptr, groups_, &(inputrec_->opts));
 287     Awh* awh = nullptr;
 288     energyOutput_->printStepToEnergyFile(
 289             mdoutf_get_fp_ene(outf), writeTrajectory, do_dr, do_or, writeLog ? fplog_ : nullptr, step, time, fcd_, awh);
 290 }
 291
 292 void EnergyData::addToForceVirial(const tensor virial, Step step)
 293 {
 294     if (step > forceVirialStep_)
 295     {
 296         forceVirialStep_ = step;
 297         clear_mat(forceVirial_);
 298     }
 299     m_add(forceVirial_, virial, forceVirial_);
 300 }
 301
 302 void EnergyData::addToConstraintVirial(const tensor virial, Step step)
 303 {
 304     if (step > shakeVirialStep_)
 305     {
 306         shakeVirialStep_ = step;
 307         clear_mat(shakeVirial_);
 308     }
 309     m_add(shakeVirial_, virial, shakeVirial_);
 310 }
 311
 312 rvec* EnergyData::forceVirial(Step gmx_unused step)
 313 {
 314     if (step > forceVirialStep_)
 315     {
 316         forceVirialStep_ = step;
 317         clear_mat(forceVirial_);
 318     }
 319     GMX_ASSERT(step >= forceVirialStep_ || forceVirialStep_ == -1,
 320                "Asked for force virial of previous step.");
 321     return forceVirial_;
 322 }
 323
 324 rvec* EnergyData::constraintVirial(Step gmx_unused step)
 325 {
 326     if (step > shakeVirialStep_)
 327     {
 328         shakeVirialStep_ = step;
 329         clear_mat(shakeVirial_);
 330     }
 331     GMX_ASSERT(step >= shakeVirialStep_ || shakeVirialStep_ == -1,
 332                "Asked for constraint virial of previous step.");
 333     return shakeVirial_;
 334 }
 335
 336 rvec* EnergyData::totalVirial(Step gmx_unused step)
 337 {
 338     if (step > totalVirialStep_)
 339     {
 340         totalVirialStep_ = step;
 341         clear_mat(totalVirial_);
 342     }
 343     GMX_ASSERT(step >= totalVirialStep_ || totalVirialStep_ == -1,
 344                "Asked for total virial of previous step.");
 345     return totalVirial_;
 346 }
 347
 348 rvec* EnergyData::pressure(Step gmx_unused step)
 349 {
 350     if (step > pressureStep_)
 351     {
 352         pressureStep_ = step;
 353         clear_mat(pressure_);
 354     }
 355     GMX_ASSERT(step >= pressureStep_ || pressureStep_ == -1, "Asked for pressure of previous step.");
 356     return pressure_;
 357 }
 358
 359 real* EnergyData::muTot()
 360 {
 361     return muTot_;
 362 }
 363
 364 gmx_enerdata_t* EnergyData::enerdata()
 365 {
 366     return enerd_;
 367 }
 368
 369 gmx_ekindata_t* EnergyData::ekindata()
 370 {
 371     return ekind_;
 372 }
 373
 374 bool* EnergyData::needToSumEkinhOld()
 375 {
 376     return &needToSumEkinhOld_;
 377 }
 378
 379 bool EnergyData::hasReadEkinFromCheckpoint() const
 380 {
 381     return hasReadEkinFromCheckpoint_;
 382 }
 383
 384 namespace
 385 {
 386 /*!
 387  * \brief Enum describing the contents EnergyData::Element writes to modular checkpoint
 388  *
 389  * When changing the checkpoint content, add a new element just above Count, and adjust the
 390  * checkpoint functionality.
 391  */
 392 enum class CheckpointVersion
 393 {
 394     Base, //!< First version of modular checkpointing
 395     Count //!< Number of entries. Add new versions right above this!
 396 };
 397 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 398 } // namespace
 399
 400 template<CheckpointDataOperation operation>
 401 void EnergyData::Element::doCheckpointData(CheckpointData<operation>* checkpointData)
 402 {
 403     checkpointVersion(checkpointData, "EnergyData version", c_currentVersion);
 404
 405     energyData_->observablesHistory_->energyHistory->doCheckpoint<operation>(
 406             checkpointData->subCheckpointData("energy history"));
 407     energyData_->ekinstate_.doCheckpoint<operation>(checkpointData->subCheckpointData("ekinstate"));
 408 }
 409
 410 void EnergyData::Element::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 411                                               const t_commrec*                   cr)
 412 {
 413     if (MASTER(cr))
 414     {
 415         if (energyData_->needToSumEkinhOld_)
 416         {
 417             energyData_->ekinstate_.bUpToDate = false;
 418         }
 419         else
 420         {
 421             update_ekinstate(&energyData_->ekinstate_, energyData_->ekind_);
 422             energyData_->ekinstate_.bUpToDate = true;
 423         }
 424         energyData_->energyOutput_->fillEnergyHistory(
 425                 energyData_->observablesHistory_->energyHistory.get());
 426         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 427     }
 428 }
 429
 430 void EnergyData::Element::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 431                                                  const t_commrec*                  cr)
 432 {
 433     if (MASTER(cr))
 434     {
 435         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 436     }
 437     energyData_->hasReadEkinFromCheckpoint_ = MASTER(cr) ? energyData_->ekinstate_.bUpToDate : false;
 438     if (PAR(cr))
 439     {
 440         gmx_bcast(sizeof(hasReadEkinFromCheckpoint_),
 441                   &energyData_->hasReadEkinFromCheckpoint_,
 442                   cr->mpi_comm_mygroup);
 443     }
 444     if (energyData_->hasReadEkinFromCheckpoint_)
 445     {
 446         // this takes care of broadcasting from master to agents
 447         restore_ekinstate_from_state(cr, energyData_->ekind_, &energyData_->ekinstate_);
 448     }
 449 }
 450
 451 const std::string& EnergyData::Element::clientID()
 452 {
 453     return identifier_;
 454 }
 455
 456 void EnergyData::initializeEnergyHistory(StartingBehavior    startingBehavior,
 457                                          ObservablesHistory* observablesHistory,
 458                                          EnergyOutput*       energyOutput)
 459 {
 460     if (startingBehavior != StartingBehavior::NewSimulation)
 461     {
 462         /* Restore from energy history if appending to output files */
 463         if (startingBehavior == StartingBehavior::RestartWithAppending)
 464         {
 465             /* If no history is available (because a checkpoint is from before
 466              * it was written) make a new one later, otherwise restore it.
 467              */
 468             if (observablesHistory->energyHistory)
 469             {
 470                 energyOutput->restoreFromEnergyHistory(*observablesHistory->energyHistory);
 471             }
 472         }
 473         else if (observablesHistory->energyHistory)
 474         {
 475             /* We might have read an energy history from checkpoint.
 476              * As we are not appending, we want to restart the statistics.
 477              * Free the allocated memory and reset the counts.
 478              */
 479             observablesHistory->energyHistory = {};
 480             /* We might have read a pull history from checkpoint.
 481              * We will still want to keep the statistics, so that the files
 482              * can be joined and still be meaningful.
 483              * This means that observablesHistory_->pullHistory
 484              * should not be reset.
 485              */
 486         }
 487     }
 488     if (!observablesHistory->energyHistory)
 489     {
 490         observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
 491     }
 492     if (!observablesHistory->pullHistory)
 493     {
 494         observablesHistory->pullHistory = std::make_unique<PullHistory>();
 495     }
 496     /* Set the initial energy history */
 497     energyOutput->fillEnergyHistory(observablesHistory->energyHistory.get());
 498 }
 499
 500 void EnergyData::addConservedEnergyContribution(EnergyContribution&& energyContribution)
 501 {
 502     conservedEnergyContributions_.emplace_back(std::move(energyContribution));
 503 }
 504
 505 void EnergyData::setParrinelloRahmanBoxVelocities(std::function<const rvec*()>&& parrinelloRahmanBoxVelocities)
 506 {
 507     GMX_RELEASE_ASSERT(!parrinelloRahmanBoxVelocities_,
 508                        "Received a second callback to the Parrinello-Rahman velocities");
 509     parrinelloRahmanBoxVelocities_ = parrinelloRahmanBoxVelocities;
 510 }
 511
 512 EnergyData::Element* EnergyData::element()
 513 {
 514     return element_.get();
 515 }
 516
 517 EnergyData::Element::Element(EnergyData* energyData, bool isMasterRank) :
 518     energyData_(energyData),
 519     isMasterRank_(isMasterRank),
 520     energyWritingStep_(-1),
 521     energyCalculationStep_(-1),
 522     freeEnergyCalculationStep_(-1)
 523 {
 524 }
 525
 526 ISimulatorElement* EnergyData::Element::getElementPointerImpl(
 527         LegacySimulatorData gmx_unused*        legacySimulatorData,
 528         ModularSimulatorAlgorithmBuilderHelper gmx_unused* builderHelper,
 529         StatePropagatorData gmx_unused* statePropagatorData,
 530         EnergyData*                     energyData,
 531         FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
 532         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
 533 {
 534     return energyData->element();
 535 }
 536
 537 } // namespace gmx