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