/*
* 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 "energydata.h"
+#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/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 "gromacs/mdtypes/pullhistory.h"
-#include "gromacs/mdtypes/state.h"
#include "gromacs/topology/topology.h"
#include "freeenergyperturbationdata.h"
#include "modularsimulator.h"
-#include "parrinellorahmanbarostat.h"
#include "simulatoralgorithm.h"
#include "statepropagatordata.h"
-#include "vrescalethermostat.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) :
+ StartingBehavior startingBehavior,
+ bool simulationsShareState,
+ pull_t* pullWork) :
element_(std::make_unique<Element>(this, isMasterRank)),
isMasterRank_(isMasterRank),
forceVirialStep_(-1),
totalVirialStep_(-1),
pressureStep_(-1),
needToSumEkinhOld_(false),
+ hasReadEkinFromCheckpoint_(false),
startingBehavior_(startingBehavior),
statePropagatorData_(statePropagatorData),
freeEnergyPerturbationData_(freeEnergyPerturbationData),
- vRescaleThermostat_(nullptr),
- parrinelloRahmanBarostat_(nullptr),
inputrec_(inputrec),
top_global_(globalTopology),
mdAtoms_(mdAtoms),
constr_(constr),
fplog_(fplog),
fcd_(fcd),
- mdModulesNotifier_(mdModulesNotifier),
- groups_(&globalTopology->groups),
- observablesHistory_(observablesHistory)
+ mdModulesNotifiers_(mdModulesNotifiers),
+ groups_(&globalTopology.groups),
+ observablesHistory_(observablesHistory),
+ simulationsShareState_(simulationsShareState),
+ pullWork_(pullWork)
{
clear_mat(forceVirial_);
clear_mat(shakeVirial_);
clear_mat(pressure_);
clear_rvec(muTot_);
- if (freeEnergyPerturbationData_)
- {
- dummyLegacyState_.flags = (1U << estFEPSTATE);
- }
+ init_ekinstate(&ekinstate_, inputrec_);
+ observablesHistory_->energyHistory = std::make_unique<energyhistory_t>();
}
void EnergyData::Element::scheduleTask(Step step, Time time, const RegisterRunFunction& registerRunFunction)
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_, 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 (vRescaleThermostat_)
- {
- dummyLegacyState_.therm_integral = vRescaleThermostat_->thermostatIntegral();
- }
if (freeEnergyPerturbationData_)
{
- accumulateKineticLambdaComponents(enerd_, freeEnergyPerturbationData_->constLambdaView(),
- *inputrec_->fepvals);
- dummyLegacyState_.fep_state = freeEnergyPerturbationData_->currentFEPState();
- }
- if (parrinelloRahmanBarostat_)
- {
- copy_mat(parrinelloRahmanBarostat_->boxVelocities(), dummyLegacyState_.boxv);
- copy_mat(statePropagatorData_->constBox(), dummyLegacyState_.box);
+ accumulateKineticLambdaComponents(
+ enerd_, freeEnergyPerturbationData_->constLambdaView(), *inputrec_->fepvals);
}
if (integratorHasConservedEnergyQuantity(inputrec_))
{
- enerd_->term[F_ECONSERVED] =
- enerd_->term[F_ETOT] + NPT_energy(inputrec_, &dummyLegacyState_, nullptr);
+ enerd_->term[F_ECONSERVED] = enerd_->term[F_ETOT];
+ for (const auto& energyContibution : conservedEnergyContributions_)
+ {
+ enerd_->term[F_ECONSERVED] += energyContibution(step, time);
+ }
}
- energyOutput_->addDataAtEnergyStep(isFreeEnergyCalculationStep, isEnergyCalculationStep, time,
- mdAtoms_->mdatoms()->tmass, enerd_, &dummyLegacyState_,
- inputrec_->fepvals, inputrec_->expandedvals,
- statePropagatorData_->constPreviousBox(), shakeVirial_,
- forceVirial_, totalVirial_, pressure_, ekind_, muTot_, constr_);
+ matrix nullMatrix = {};
+ energyOutput_->addDataAtEnergyStep(
+ isFreeEnergyCalculationStep,
+ isEnergyCalculationStep,
+ time,
+ mdAtoms_->mdatoms()->tmass,
+ enerd_,
+ inputrec_->fepvals.get(),
+ statePropagatorData_->constPreviousBox(),
+ PTCouplingArrays({ parrinelloRahmanBoxVelocities_ ? parrinelloRahmanBoxVelocities_() : nullMatrix,
+ {},
+ {},
+ {},
+ {} }),
+ freeEnergyPerturbationData_ ? freeEnergyPerturbationData_->currentFEPState() : 0,
+ 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_;
return &needToSumEkinhOld_;
}
-void EnergyData::Element::writeCheckpoint(t_state gmx_unused* localState, t_state* globalState)
+bool EnergyData::hasReadEkinFromCheckpoint() const
+{
+ return hasReadEkinFromCheckpoint_;
+}
+
+namespace
{
- if (isMasterRank_)
+/*!
+ * \brief Enum describing the contents EnergyData::Element writes to modular checkpoint
+ *
+ * When changing the checkpoint content, add a new element just above Count, and adjust the
+ * checkpoint functionality.
+ */
+enum class CheckpointVersion
+{
+ Base, //!< First version of modular checkpointing
+ Count //!< Number of entries. Add new versions right above this!
+};
+constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
+} // namespace
+
+template<CheckpointDataOperation operation>
+void EnergyData::Element::doCheckpointData(CheckpointData<operation>* checkpointData)
+{
+ checkpointVersion(checkpointData, "EnergyData version", c_currentVersion);
+
+ energyData_->observablesHistory_->energyHistory->doCheckpoint<operation>(
+ checkpointData->subCheckpointData("energy history"));
+ energyData_->ekinstate_.doCheckpoint<operation>(checkpointData->subCheckpointData("ekinstate"));
+}
+
+void EnergyData::Element::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
+ const t_commrec* cr)
+{
+ if (MASTER(cr))
{
if (energyData_->needToSumEkinhOld_)
{
- globalState->ekinstate.bUpToDate = false;
+ energyData_->ekinstate_.bUpToDate = false;
}
else
{
- update_ekinstate(&globalState->ekinstate, energyData_->ekind_);
- globalState->ekinstate.bUpToDate = true;
+ update_ekinstate(&energyData_->ekinstate_, energyData_->ekind_);
+ energyData_->ekinstate_.bUpToDate = true;
}
energyData_->energyOutput_->fillEnergyHistory(
energyData_->observablesHistory_->energyHistory.get());
+ doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
+ }
+}
+
+void EnergyData::Element::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
+ const t_commrec* cr)
+{
+ if (MASTER(cr))
+ {
+ doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
+ }
+ energyData_->hasReadEkinFromCheckpoint_ = MASTER(cr) ? energyData_->ekinstate_.bUpToDate : false;
+ if (PAR(cr))
+ {
+ gmx_bcast(sizeof(hasReadEkinFromCheckpoint_),
+ &energyData_->hasReadEkinFromCheckpoint_,
+ cr->mpi_comm_mygroup);
+ }
+ if (energyData_->hasReadEkinFromCheckpoint_)
+ {
+ // this takes care of broadcasting from master to agents
+ restore_ekinstate_from_state(cr, energyData_->ekind_, &energyData_->ekinstate_);
}
}
+const std::string& EnergyData::Element::clientID()
+{
+ return identifier_;
+}
+
void EnergyData::initializeEnergyHistory(StartingBehavior startingBehavior,
ObservablesHistory* observablesHistory,
EnergyOutput* energyOutput)
energyOutput->fillEnergyHistory(observablesHistory->energyHistory.get());
}
-void EnergyData::setVRescaleThermostat(const gmx::VRescaleThermostat* vRescaleThermostat)
+void EnergyData::addConservedEnergyContribution(EnergyContribution&& energyContribution)
{
- vRescaleThermostat_ = vRescaleThermostat;
- if (vRescaleThermostat_)
- {
- dummyLegacyState_.flags |= (1U << estTHERM_INT);
- }
+ conservedEnergyContributions_.emplace_back(std::move(energyContribution));
}
-void EnergyData::setParrinelloRahamnBarostat(const gmx::ParrinelloRahmanBarostat* parrinelloRahmanBarostat)
+void EnergyData::setParrinelloRahmanBoxVelocities(std::function<const rvec*()>&& parrinelloRahmanBoxVelocities)
{
- parrinelloRahmanBarostat_ = parrinelloRahmanBarostat;
- if (parrinelloRahmanBarostat_)
- {
- dummyLegacyState_.flags |= (1U << estBOX) | (1U << estBOXV);
- }
+ 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();
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff