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