src/gromacs/modularsimulator/energyelement.cpp

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