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