/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/compute_io.h"
-#include "gromacs/mdlib/coupling.h"
#include "gromacs/mdlib/enerdata_utils.h"
#include "gromacs/mdlib/energyoutput.h"
#include "gromacs/mdlib/mdatoms.h"
#include "gromacs/mdlib/mdoutf.h"
#include "gromacs/mdlib/stat.h"
+#include "gromacs/mdlib/tgroup.h"
#include "gromacs/mdlib/update.h"
-#include "gromacs/mdrunutility/handlerestart.h"
-#include "gromacs/mdtypes/checkpointdata.h"
#include "gromacs/mdtypes/commrec.h"
#include "gromacs/mdtypes/enerdata.h"
#include "gromacs/mdtypes/energyhistory.h"
+#include "gromacs/mdtypes/group.h"
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/observableshistory.h"
#include "freeenergyperturbationdata.h"
#include "modularsimulator.h"
-#include "parrinellorahmanbarostat.h"
#include "simulatoralgorithm.h"
#include "statepropagatordata.h"
-#include "velocityscalingtemperaturecoupling.h"
struct pull_t;
class t_state;
EnergyData::EnergyData(StatePropagatorData* statePropagatorData,
FreeEnergyPerturbationData* freeEnergyPerturbationData,
- const gmx_mtop_t* globalTopology,
+ const gmx_mtop_t& globalTopology,
const t_inputrec* inputrec,
const MDAtoms* mdAtoms,
gmx_enerdata_t* enerd,
const Constraints* constr,
FILE* fplog,
t_fcdata* fcd,
- const MdModulesNotifier& mdModulesNotifier,
+ const MDModulesNotifiers& mdModulesNotifiers,
bool isMasterRank,
ObservablesHistory* observablesHistory,
StartingBehavior startingBehavior,
- bool simulationsShareState) :
+ bool simulationsShareState,
+ pull_t* pullWork) :
element_(std::make_unique<Element>(this, isMasterRank)),
isMasterRank_(isMasterRank),
forceVirialStep_(-1),
startingBehavior_(startingBehavior),
statePropagatorData_(statePropagatorData),
freeEnergyPerturbationData_(freeEnergyPerturbationData),
- velocityScalingTemperatureCoupling_(nullptr),
- parrinelloRahmanBarostat_(nullptr),
inputrec_(inputrec),
top_global_(globalTopology),
mdAtoms_(mdAtoms),
constr_(constr),
fplog_(fplog),
fcd_(fcd),
- mdModulesNotifier_(mdModulesNotifier),
- groups_(&globalTopology->groups),
+ mdModulesNotifiers_(mdModulesNotifiers),
+ groups_(&globalTopology.groups),
observablesHistory_(observablesHistory),
- simulationsShareState_(simulationsShareState)
+ simulationsShareState_(simulationsShareState),
+ pullWork_(pullWork)
{
clear_mat(forceVirial_);
clear_mat(shakeVirial_);
auto isFreeEnergyCalculationStep = freeEnergyCalculationStep_ == step;
if (isEnergyCalculationStep || writeEnergy)
{
- registerRunFunction([this, time, isEnergyCalculationStep, isFreeEnergyCalculationStep]() {
- energyData_->doStep(time, isEnergyCalculationStep, isFreeEnergyCalculationStep);
+ registerRunFunction([this, step, time, isEnergyCalculationStep, isFreeEnergyCalculationStep]() {
+ energyData_->doStep(step, time, isEnergyCalculationStep, isFreeEnergyCalculationStep);
});
}
else
{
if (inputrec_->nstcalcenergy > 0 && isMasterRank_)
{
+ energyOutput_->printEnergyConservation(fplog_, inputrec_->simulation_part, EI_MD(inputrec_->eI));
energyOutput_->printAverages(fplog_, groups_);
}
}
void EnergyData::setup(gmx_mdoutf* outf)
{
- pull_t* pull_work = nullptr;
- energyOutput_ = std::make_unique<EnergyOutput>(
- mdoutf_get_fp_ene(outf), top_global_, inputrec_, pull_work, mdoutf_get_fp_dhdl(outf),
- false, startingBehavior_, simulationsShareState_, mdModulesNotifier_);
+ energyOutput_ = std::make_unique<EnergyOutput>(mdoutf_get_fp_ene(outf),
+ top_global_,
+ *inputrec_,
+ pullWork_,
+ mdoutf_get_fp_dhdl(outf),
+ false,
+ startingBehavior_,
+ simulationsShareState_,
+ mdModulesNotifiers_);
if (!isMasterRank_)
{
// TODO: This probably doesn't really belong here...
// but we have all we need in this element,
// so we'll leave it here for now!
- double io = compute_io(inputrec_, top_global_->natoms, *groups_, energyOutput_->numEnergyTerms(), 1);
+ double io = compute_io(inputrec_, top_global_.natoms, *groups_, energyOutput_->numEnergyTerms(), 1);
if ((io > 2000) && isMasterRank_)
{
fprintf(stderr, "\nWARNING: This run will generate roughly %.0f Mb of data\n\n", io);
if (!inputrec_->bContinuation)
{
real temp = enerd_->term[F_TEMP];
- if (inputrec_->eI != eiVV)
+ if (inputrec_->eI != IntegrationAlgorithm::VV)
{
/* Result of Ekin averaged over velocities of -half
* and +half step, while we only have -half step here.
return std::nullopt;
}
-void EnergyData::doStep(Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep)
+void EnergyData::doStep(Step step, Time time, bool isEnergyCalculationStep, bool isFreeEnergyCalculationStep)
{
enerd_->term[F_ETOT] = enerd_->term[F_EPOT] + enerd_->term[F_EKIN];
if (freeEnergyPerturbationData_)
{
- accumulateKineticLambdaComponents(enerd_, freeEnergyPerturbationData_->constLambdaView(),
- *inputrec_->fepvals);
+ accumulateKineticLambdaComponents(
+ enerd_, freeEnergyPerturbationData_->constLambdaView(), *inputrec_->fepvals);
}
if (integratorHasConservedEnergyQuantity(inputrec_))
{
- enerd_->term[F_ECONSERVED] =
- enerd_->term[F_ETOT]
- + (velocityScalingTemperatureCoupling_
- ? velocityScalingTemperatureCoupling_->conservedEnergyContribution()
- : 0)
- + (parrinelloRahmanBarostat_ ? parrinelloRahmanBarostat_->conservedEnergyContribution() : 0);
+ enerd_->term[F_ECONSERVED] = enerd_->term[F_ETOT];
+ for (const auto& energyContibution : conservedEnergyContributions_)
+ {
+ enerd_->term[F_ECONSERVED] += energyContibution(step, time);
+ }
}
matrix nullMatrix = {};
energyOutput_->addDataAtEnergyStep(
- isFreeEnergyCalculationStep, isEnergyCalculationStep, time, mdAtoms_->mdatoms()->tmass, enerd_,
- inputrec_->fepvals, inputrec_->expandedvals, statePropagatorData_->constPreviousBox(),
- PTCouplingArrays({ parrinelloRahmanBarostat_ ? parrinelloRahmanBarostat_->boxVelocities() : nullMatrix,
+ isFreeEnergyCalculationStep,
+ isEnergyCalculationStep,
+ time,
+ mdAtoms_->mdatoms()->tmass,
+ enerd_,
+ inputrec_->fepvals.get(),
+ statePropagatorData_->constPreviousBox(),
+ PTCouplingArrays({ parrinelloRahmanBoxVelocities_ ? parrinelloRahmanBoxVelocities_() : nullMatrix,
{},
{},
{},
{} }),
freeEnergyPerturbationData_ ? freeEnergyPerturbationData_->currentFEPState() : 0,
- shakeVirial_, forceVirial_, totalVirial_, pressure_, ekind_, muTot_, constr_);
+ totalVirial_,
+ pressure_,
+ ekind_,
+ muTot_,
+ constr_);
}
void EnergyData::write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog)
// energyOutput_->printAnnealingTemperatures(writeLog ? fplog_ : nullptr, groups_, &(inputrec_->opts));
Awh* awh = nullptr;
- energyOutput_->printStepToEnergyFile(mdoutf_get_fp_ene(outf), writeTrajectory, do_dr, do_or,
- writeLog ? fplog_ : nullptr, step, time, fcd_, awh);
+ energyOutput_->printStepToEnergyFile(
+ mdoutf_get_fp_ene(outf), writeTrajectory, do_dr, do_or, writeLog ? fplog_ : nullptr, step, time, fcd_, awh);
}
void EnergyData::addToForceVirial(const tensor virial, Step step)
pressureStep_ = step;
clear_mat(pressure_);
}
- GMX_ASSERT(step >= pressureStep_ || pressureStep_ == -1,
- "Asked for pressure of previous step.");
+ GMX_ASSERT(step >= pressureStep_ || pressureStep_ == -1, "Asked for pressure of previous step.");
return pressure_;
}
return enerd_;
}
+const gmx_enerdata_t* EnergyData::enerdata() const
+{
+ return enerd_;
+}
+
gmx_ekindata_t* EnergyData::ekindata()
{
return ekind_;
energyData_->hasReadEkinFromCheckpoint_ = MASTER(cr) ? energyData_->ekinstate_.bUpToDate : false;
if (PAR(cr))
{
- gmx_bcast(sizeof(hasReadEkinFromCheckpoint_), &energyData_->hasReadEkinFromCheckpoint_,
+ gmx_bcast(sizeof(hasReadEkinFromCheckpoint_),
+ &energyData_->hasReadEkinFromCheckpoint_,
cr->mpi_comm_mygroup);
}
if (energyData_->hasReadEkinFromCheckpoint_)
energyOutput->fillEnergyHistory(observablesHistory->energyHistory.get());
}
-void EnergyData::setVelocityScalingTemperatureCoupling(const VelocityScalingTemperatureCoupling* velocityScalingTemperatureCoupling)
+void EnergyData::addConservedEnergyContribution(EnergyContribution&& energyContribution)
{
- velocityScalingTemperatureCoupling_ = velocityScalingTemperatureCoupling;
+ conservedEnergyContributions_.emplace_back(std::move(energyContribution));
}
-void EnergyData::setParrinelloRahamnBarostat(const gmx::ParrinelloRahmanBarostat* parrinelloRahmanBarostat)
+void EnergyData::setParrinelloRahmanBoxVelocities(std::function<const rvec*()>&& parrinelloRahmanBoxVelocities)
{
- parrinelloRahmanBarostat_ = parrinelloRahmanBarostat;
+ GMX_RELEASE_ASSERT(!parrinelloRahmanBoxVelocities_,
+ "Received a second callback to the Parrinello-Rahman velocities");
+ parrinelloRahmanBoxVelocities_ = parrinelloRahmanBoxVelocities;
+}
+
+void EnergyData::updateKineticEnergy()
+{
+ // The legacy sum_ekin function does not offer named types, so define variables for readability
+ // dEkin/dlambda is not handled here
+ real* dEkinDLambda = nullptr;
+ // Whether we use the full step kinetic energy (vs the average of half step KEs)
+ const bool useFullStepKineticEnergy = (inputrec_->eI == IntegrationAlgorithm::VV);
+ /* Whether we're ignoring the NHC scaling factor, only used if useFullStepKineticEnergy
+ * is true. (This parameter is confusing, as it is named `bScaleEkin`, but prompts the
+ * function to ignore scaling. There is no use case within modular simulator to ignore
+ * these, so we set this to false.) */
+ const bool ignoreScalingFactor = false;
+
+ enerd_->term[F_TEMP] = sum_ekin(
+ &(inputrec_->opts), ekind_, dEkinDLambda, useFullStepKineticEnergy, ignoreScalingFactor);
+ enerd_->term[F_EKIN] = trace(ekind_->ekin);
}
EnergyData::Element* EnergyData::element()
StatePropagatorData gmx_unused* statePropagatorData,
EnergyData* energyData,
FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
- GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
+ GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
+ ObservablesReducer* /*observablesReducer*/)
{
return energyData->element();
}