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