From 00e8ec253bcf966ae356a0a8f9f38e9a32a206e9 Mon Sep 17 00:00:00 2001 From: Pascal Merz Date: Mon, 7 Sep 2020 20:23:36 +0000 Subject: [PATCH] Make EnergyOutput::addDataAtEnergyStep independent of t_state EnergyOutput::addDataAtEnergyStep took a pointer to a full t_state, although it only uses a small part of the object. This makes the dependencies explicit. This is a necessary step to make EnergyData independent of t_state. Refs #3419 --- src/gromacs/mdlib/energyoutput.cpp | 35 +++++++------- src/gromacs/mdlib/energyoutput.h | 53 ++++++++++++++------- src/gromacs/mdlib/tests/energyoutput.cpp | 9 ++-- src/gromacs/mdrun/md.cpp | 9 ++-- src/gromacs/mdrun/mimic.cpp | 10 ++-- src/gromacs/mdrun/minimize.cpp | 23 ++++----- src/gromacs/mdrun/rerun.cpp | 10 ++-- src/gromacs/modularsimulator/energydata.cpp | 17 ++++--- 8 files changed, 101 insertions(+), 65 deletions(-) diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index 120591f0a7..f658ea377f 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -855,10 +855,11 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, double time, real tmass, const gmx_enerdata_t* enerd, - const t_state* state, const t_lambda* fep, const t_expanded* expand, const matrix box, + PTCouplingArrays ptCouplingArrays, + int fep_state, const tensor svir, const tensor fvir, const tensor vir, @@ -936,12 +937,12 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, } if (epc_ == epcPARRINELLORAHMAN || epc_ == epcMTTK) { - tmp6[0] = state->boxv[XX][XX]; - tmp6[1] = state->boxv[YY][YY]; - tmp6[2] = state->boxv[ZZ][ZZ]; - tmp6[3] = state->boxv[YY][XX]; - tmp6[4] = state->boxv[ZZ][XX]; - tmp6[5] = state->boxv[ZZ][YY]; + tmp6[0] = ptCouplingArrays.boxv[XX][XX]; + tmp6[1] = ptCouplingArrays.boxv[YY][YY]; + tmp6[2] = ptCouplingArrays.boxv[ZZ][ZZ]; + tmp6[3] = ptCouplingArrays.boxv[YY][XX]; + tmp6[4] = ptCouplingArrays.boxv[ZZ][XX]; + tmp6[5] = ptCouplingArrays.boxv[ZZ][YY]; add_ebin(ebin_, ipc_, bTricl_ ? 6 : 3, tmp6, bSum); } if (bMu_) @@ -1000,8 +1001,8 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, for (j = 0; j < nNHC_; j++) { k = i * nNHC_ + j; - tmp_r_[2 * k] = state->nosehoover_xi[k]; - tmp_r_[2 * k + 1] = state->nosehoover_vxi[k]; + tmp_r_[2 * k] = ptCouplingArrays.nosehoover_xi[k]; + tmp_r_[2 * k + 1] = ptCouplingArrays.nosehoover_vxi[k]; } } add_ebin(ebin_, itc_, mde_n_, tmp_r_, bSum); @@ -1013,8 +1014,8 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, for (j = 0; j < nNHC_; j++) { k = i * nNHC_ + j; - tmp_r_[2 * k] = state->nhpres_xi[k]; - tmp_r_[2 * k + 1] = state->nhpres_vxi[k]; + tmp_r_[2 * k] = ptCouplingArrays.nhpres_xi[k]; + tmp_r_[2 * k + 1] = ptCouplingArrays.nhpres_vxi[k]; } } add_ebin(ebin_, itcb_, mdeb_n_, tmp_r_, bSum); @@ -1024,8 +1025,8 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, { for (int i = 0; (i < nTC_); i++) { - tmp_r_[2 * i] = state->nosehoover_xi[i]; - tmp_r_[2 * i + 1] = state->nosehoover_vxi[i]; + tmp_r_[2 * i] = ptCouplingArrays.nosehoover_xi[i]; + tmp_r_[2 * i + 1] = ptCouplingArrays.nosehoover_vxi[i]; } add_ebin(ebin_, itc_, mde_n_, tmp_r_, bSum); } @@ -1062,14 +1063,14 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, dE_[i] = foreignTerms.deltaH(i); if (numTemperatures_ > 0) { - GMX_RELEASE_ASSERT(numTemperatures_ > state->fep_state, + GMX_RELEASE_ASSERT(numTemperatures_ > fep_state, "Number of lambdas in state is bigger then in input record"); GMX_RELEASE_ASSERT( numTemperatures_ >= foreignTerms.numLambdas(), "Number of lambdas in energy data is bigger then in input record"); /* MRS: is this right, given the way we have defined the exchange probabilities? */ /* is this even useful to have at all? */ - dE_[i] += (temperatures_[i] / temperatures_[state->fep_state] - 1.0) * enerd->term[F_EKIN]; + dE_[i] += (temperatures_[i] / temperatures_[fep_state] - 1.0) * enerd->term[F_EKIN]; } } @@ -1081,7 +1082,7 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, /* print the current state if we are doing expanded ensemble */ if (expand->elmcmove > elmcmoveNO) { - fprintf(fp_dhdl_, " %4d", state->fep_state); + fprintf(fp_dhdl_, " %4d", fep_state); } /* total energy (for if the temperature changes */ @@ -1137,7 +1138,7 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, } store_energy = enerd->term[F_ETOT]; /* store_dh is dE */ - mde_delta_h_coll_add_dh(dhc_, static_cast(state->fep_state), store_energy, pv, + mde_delta_h_coll_add_dh(dhc_, static_cast(fep_state), store_energy, pv, store_dhdl, dE_ + fep->lambda_start_n, time); } } diff --git a/src/gromacs/mdlib/energyoutput.h b/src/gromacs/mdlib/energyoutput.h index 6bfff750de..18bfa98de9 100644 --- a/src/gromacs/mdlib/energyoutput.h +++ b/src/gromacs/mdlib/energyoutput.h @@ -98,6 +98,23 @@ enum namespace gmx { +/*! \internal + * \brief Arrays connected to Pressure and Temperature coupling + */ +struct PTCouplingArrays +{ + //! Box velocities for Parrinello-Rahman P-coupling. + const rvec* boxv; + //! Nose-Hoover coordinates. + ArrayRef nosehoover_xi; + //! Nose-Hoover velocities. + ArrayRef nosehoover_vxi; + //! Pressure Nose-Hoover coordinates. + ArrayRef nhpres_xi; + //! Pressure Nose-Hoover velocities. + ArrayRef nhpres_vxi; +}; + /* Energy output class * * This is the collection of energy averages collected during mdrun, and to @@ -136,32 +153,34 @@ public: * * Called every step on which the thermodynamic values are evaluated. * - * \param[in] bDoDHDL Whether the FEP is enabled. - * \param[in] bSum If this stepshould be recorded to compute sums and averaes. - * \param[in] time Current simulation time. - * \param[in] tmass Total mass - * \param[in] enerd Energy data object. - * \param[in] state System state. - * \param[in] fep FEP data. - * \param[in] expand Expanded ensemble (for FEP). - * \param[in] lastbox PBC data. - * \param[in] svir Constraint virial. - * \param[in] fvir Force virial. - * \param[in] vir Total virial. - * \param[in] pres Pressure. - * \param[in] ekind Kinetic energy data. - * \param[in] mu_tot Total dipole. - * \param[in] constr Constraints object to get RMSD from (for LINCS only). + * \param[in] bDoDHDL Whether the FEP is enabled. + * \param[in] bSum If this stepshould be recorded to compute sums and averages. + * \param[in] time Current simulation time. + * \param[in] tmass Total mass + * \param[in] enerd Energy data object. + * \param[in] fep FEP data. + * \param[in] expand Expanded ensemble (for FEP). + * \param[in] lastbox PBC data. + * \param[in] ptCouplingArrays Arrays connected to pressure and temperature coupling. + * \param[in] fep_state The current alchemical state we are in. + * \param[in] svir Constraint virial. + * \param[in] fvir Force virial. + * \param[in] vir Total virial. + * \param[in] pres Pressure. + * \param[in] ekind Kinetic energy data. + * \param[in] mu_tot Total dipole. + * \param[in] constr Constraints object to get RMSD from (for LINCS only). */ void addDataAtEnergyStep(bool bDoDHDL, bool bSum, double time, real tmass, const gmx_enerdata_t* enerd, - const t_state* state, const t_lambda* fep, const t_expanded* expand, const matrix lastbox, + PTCouplingArrays ptCouplingArrays, + int fep_state, const tensor svir, const tensor fvir, const tensor vir, diff --git a/src/gromacs/mdlib/tests/energyoutput.cpp b/src/gromacs/mdlib/tests/energyoutput.cpp index 422213bb85..bb4ce4ce00 100644 --- a/src/gromacs/mdlib/tests/energyoutput.cpp +++ b/src/gromacs/mdlib/tests/energyoutput.cpp @@ -634,9 +634,12 @@ TEST_P(EnergyOutputTest, CheckOutput) for (int frame = 0; frame < parameters.numFrames; frame++) { setStepData(&testValue); - energyOutput->addDataAtEnergyStep(false, true, time_, tmass_, enerdata_.get(), &state_, nullptr, - nullptr, box_, constraintsVirial_, forceVirial_, totalVirial_, - pressure_, &ekindata_, muTotal_, constraints_.get()); + energyOutput->addDataAtEnergyStep( + false, true, time_, tmass_, enerdata_.get(), nullptr, nullptr, box_, + PTCouplingArrays({ state_.boxv, state_.nosehoover_xi, state_.nosehoover_vxi, + state_.nhpres_xi, state_.nhpres_vxi }), + state_.fep_state, constraintsVirial_, forceVirial_, totalVirial_, pressure_, + &ekindata_, muTotal_, constraints_.get()); energyOutput->printAnnealingTemperatures(log_, &mtop_.groups, &inputrec_.opts); energyOutput->printStepToEnergyFile(energyFile_, true, false, false, log_, 100 * frame, diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 3d1c95a411..a97e0eeb75 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -1509,9 +1509,12 @@ void gmx::LegacySimulator::do_md() } if (bCalcEner) { - energyOutput.addDataAtEnergyStep(bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, state, - ir->fepvals, ir->expandedvals, lastbox, shake_vir, - force_vir, total_vir, pres, ekind, mu_tot, constr); + energyOutput.addDataAtEnergyStep( + bDoDHDL, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals, + ir->expandedvals, lastbox, + PTCouplingArrays{ state->boxv, state->nosehoover_xi, state->nosehoover_vxi, + state->nhpres_xi, state->nhpres_vxi }, + state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr); } else { diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index e6c924c2eb..edd2814e58 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -515,10 +515,12 @@ void gmx::LegacySimulator::do_mimic() if (MASTER(cr)) { const bool bCalcEnerStep = true; - energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, bCalcEnerStep, t, - mdatoms->tmass, enerd, state, ir->fepvals, - ir->expandedvals, state->box, shake_vir, force_vir, - total_vir, pres, ekind, mu_tot, constr); + energyOutput.addDataAtEnergyStep( + doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals, + ir->expandedvals, state->box, + PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi, + state->nhpres_xi, state->nhpres_vxi }), + state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr); const bool do_ene = true; const bool do_log = true; diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 95815d24b7..9c2b75760f 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -1126,8 +1126,8 @@ void LegacySimulator::do_cg() /* Copy stuff to the energy bin for easy printing etc. */ matrix nullBox = {}; energyOutput.addDataAtEnergyStep(false, false, static_cast(step), mdatoms->tmass, - enerd, nullptr, nullptr, nullptr, nullBox, nullptr, - nullptr, vir, pres, nullptr, mu_tot, constr); + enerd, nullptr, nullptr, nullBox, PTCouplingArrays(), 0, + nullptr, nullptr, vir, pres, nullptr, mu_tot, constr); EnergyOutput::printHeader(fplog, step, step); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, @@ -1535,8 +1535,8 @@ void LegacySimulator::do_cg() /* Store the new (lower) energies */ matrix nullBox = {}; energyOutput.addDataAtEnergyStep(false, false, static_cast(step), mdatoms->tmass, - enerd, nullptr, nullptr, nullptr, nullBox, nullptr, - nullptr, vir, pres, nullptr, mu_tot, constr); + enerd, nullptr, nullptr, nullBox, PTCouplingArrays(), 0, + nullptr, nullptr, vir, pres, nullptr, mu_tot, constr); do_log = do_per_step(step, inputrec->nstlog); do_ene = do_per_step(step, inputrec->nstenergy); @@ -1776,8 +1776,8 @@ void LegacySimulator::do_lbfgs() /* Copy stuff to the energy bin for easy printing etc. */ matrix nullBox = {}; energyOutput.addDataAtEnergyStep(false, false, static_cast(step), mdatoms->tmass, - enerd, nullptr, nullptr, nullptr, nullBox, nullptr, - nullptr, vir, pres, nullptr, mu_tot, constr); + enerd, nullptr, nullptr, nullBox, PTCouplingArrays(), 0, + nullptr, nullptr, vir, pres, nullptr, mu_tot, constr); EnergyOutput::printHeader(fplog, step, step); energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), TRUE, FALSE, FALSE, fplog, step, @@ -2261,8 +2261,8 @@ void LegacySimulator::do_lbfgs() /* Store the new (lower) energies */ matrix nullBox = {}; energyOutput.addDataAtEnergyStep(false, false, static_cast(step), mdatoms->tmass, - enerd, nullptr, nullptr, nullptr, nullBox, nullptr, - nullptr, vir, pres, nullptr, mu_tot, constr); + enerd, nullptr, nullptr, nullBox, PTCouplingArrays(), 0, + nullptr, nullptr, vir, pres, nullptr, mu_tot, constr); do_log = do_per_step(step, inputrec->nstlog); do_ene = do_per_step(step, inputrec->nstenergy); @@ -2473,9 +2473,10 @@ void LegacySimulator::do_steep() { /* Store the new (lower) energies */ matrix nullBox = {}; - energyOutput.addDataAtEnergyStep(false, false, static_cast(count), mdatoms->tmass, - enerd, nullptr, nullptr, nullptr, nullBox, nullptr, - nullptr, vir, pres, nullptr, mu_tot, constr); + energyOutput.addDataAtEnergyStep(false, false, static_cast(count), + mdatoms->tmass, enerd, nullptr, nullptr, nullBox, + PTCouplingArrays(), 0, nullptr, nullptr, vir, pres, + nullptr, mu_tot, constr); imdSession->fillEnergyRecord(count, TRUE); diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 12dcb98b7b..8230157ded 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -594,10 +594,12 @@ void gmx::LegacySimulator::do_rerun() if (MASTER(cr)) { const bool bCalcEnerStep = true; - energyOutput.addDataAtEnergyStep(doFreeEnergyPerturbation, bCalcEnerStep, t, - mdatoms->tmass, enerd, state, ir->fepvals, - ir->expandedvals, state->box, shake_vir, force_vir, - total_vir, pres, ekind, mu_tot, constr); + energyOutput.addDataAtEnergyStep( + doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals, + ir->expandedvals, state->box, + PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi, + state->nhpres_xi, state->nhpres_vxi }), + state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr); const bool do_ene = true; const bool do_log = true; diff --git a/src/gromacs/modularsimulator/energydata.cpp b/src/gromacs/modularsimulator/energydata.cpp index 88383cc0f3..d8ec87cb3d 100644 --- a/src/gromacs/modularsimulator/energydata.cpp +++ b/src/gromacs/modularsimulator/energydata.cpp @@ -241,7 +241,6 @@ void EnergyData::doStep(Time time, bool isEnergyCalculationStep, bool isFreeEner { accumulateKineticLambdaComponents(enerd_, freeEnergyPerturbationData_->constLambdaView(), *inputrec_->fepvals); - dummyLegacyState_.fep_state = freeEnergyPerturbationData_->currentFEPState(); } if (parrinelloRahmanBarostat_) { @@ -253,11 +252,17 @@ void EnergyData::doStep(Time time, bool isEnergyCalculationStep, bool isFreeEner enerd_->term[F_ECONSERVED] = enerd_->term[F_ETOT] + NPT_energy(inputrec_, &dummyLegacyState_, nullptr); } - 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, inputrec_->expandedvals, statePropagatorData_->constPreviousBox(), + PTCouplingArrays({ parrinelloRahmanBarostat_ ? parrinelloRahmanBarostat_->boxVelocities() : nullMatrix, + {}, + {}, + {}, + {} }), + freeEnergyPerturbationData_ ? freeEnergyPerturbationData_->currentFEPState() : 0, + shakeVirial_, forceVirial_, totalVirial_, pressure_, ekind_, muTot_, constr_); } void EnergyData::write(gmx_mdoutf* outf, Step step, Time time, bool writeTrajectory, bool writeLog) -- 2.22.0