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