From f25c404f340028dafc8446cf07d556d02b476b60 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Thu, 21 Oct 2021 11:50:37 +0000 Subject: [PATCH] Add lambda state to dhdl.xvg with AWH fep --- .../applied_forces/awh/read_params.cpp | 16 ++++++++++ src/gromacs/applied_forces/awh/read_params.h | 3 ++ src/gromacs/mdlib/energyoutput.cpp | 32 ++++++++++++------- src/gromacs/mdlib/energyoutput.h | 5 ++- src/gromacs/mdlib/tests/energyoutput.cpp | 1 - src/gromacs/mdrun/md.cpp | 1 - src/gromacs/mdrun/mimic.cpp | 1 - src/gromacs/mdrun/minimize.cpp | 5 --- src/gromacs/mdrun/rerun.cpp | 1 - src/gromacs/modularsimulator/energydata.cpp | 1 - 10 files changed, 41 insertions(+), 25 deletions(-) diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp index 6fa98dc4c3..f2b1805209 100644 --- a/src/gromacs/applied_forces/awh/read_params.cpp +++ b/src/gromacs/applied_forces/awh/read_params.cpp @@ -1357,4 +1357,20 @@ void checkAwhParams(const AwhParams& awhParams, const t_inputrec& ir, warninp_t } } +bool awhHasFepLambdaDimension(const AwhParams& awhParams) +{ + for (const auto& biasParams : awhParams.awhBiasParams()) + { + for (const auto& dimParams : biasParams.dimParams()) + { + if (dimParams.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda) + { + return true; + } + } + } + + return false; +} + } // namespace gmx diff --git a/src/gromacs/applied_forces/awh/read_params.h b/src/gromacs/applied_forces/awh/read_params.h index 22247ae480..9e33649a60 100644 --- a/src/gromacs/applied_forces/awh/read_params.h +++ b/src/gromacs/applied_forces/awh/read_params.h @@ -100,6 +100,9 @@ void setStateDependentAwhParams(AwhParams* awhParams, const gmx_mtop_t& mtop, warninp_t wi); +//! Returns true when AWH has a bias with a free energy lambda state dimension +bool awhHasFepLambdaDimension(const AwhParams& awhParams); + } // namespace gmx #endif /* GMX_AWH_READPARAMS_H */ diff --git a/src/gromacs/mdlib/energyoutput.cpp b/src/gromacs/mdlib/energyoutput.cpp index 469fac5336..2a6c91e9b6 100644 --- a/src/gromacs/mdlib/energyoutput.cpp +++ b/src/gromacs/mdlib/energyoutput.cpp @@ -56,6 +56,7 @@ #include #include "gromacs/applied_forces/awh/awh.h" +#include "gromacs/applied_forces/awh/read_params.h" #include "gromacs/fileio/enxio.h" #include "gromacs/fileio/gmxfio.h" #include "gromacs/fileio/xvgr.h" @@ -123,6 +124,13 @@ const char* enumValueToString(NonBondedEnergyTerms enumValue) //! \} +static bool haveFepLambdaMoves(const t_inputrec& inputrec) +{ + return (inputrec.bExpanded && inputrec.expandedvals->elmcmove > LambdaMoveCalculation::No) + || (inputrec.efep != FreeEnergyPerturbationType::No && inputrec.bDoAwh + && awhHasFepLambdaDimension(*inputrec.awhParams)); +} + namespace gmx { @@ -142,7 +150,8 @@ EnergyOutput::EnergyOutput(ener_file* fp_ene, bool isRerun, const StartingBehavior startingBehavior, const bool simulationsShareState, - const MDModulesNotifiers& mdModulesNotifiers) + const MDModulesNotifiers& mdModulesNotifiers) : + haveFepLambdaMoves_(haveFepLambdaMoves(inputrec)) { const char* ener_nm[F_NRE]; static const char* vir_nm[] = { "Vir-XX", "Vir-XY", "Vir-XZ", "Vir-YX", "Vir-YY", @@ -635,11 +644,10 @@ FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env FILE* fp; const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda", *lambdastate = "\\lambda state"; - int i, nsets, nsets_de, nsetsbegin; - int n_lambda_terms = 0; - t_lambda* fep = ir->fepvals.get(); /* for simplicity */ - t_expanded* expand = ir->expandedvals.get(); - char lambda_vec_str[STRLEN], lambda_name_str[STRLEN]; + int i, nsets, nsets_de, nsetsbegin; + int n_lambda_terms = 0; + t_lambda* fep = ir->fepvals.get(); /* for simplicity */ + char lambda_vec_str[STRLEN], lambda_name_str[STRLEN]; int nsets_dhdl = 0; int s = 0; @@ -678,7 +686,8 @@ FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env buf = gmx::formatString("T = %g (K) ", ir->opts.ref_t[0]); } if ((ir->efep != FreeEnergyPerturbationType::SlowGrowth) - && (ir->efep != FreeEnergyPerturbationType::Expanded)) + && (ir->efep != FreeEnergyPerturbationType::Expanded) + && !(ir->bDoAwh && awhHasFepLambdaDimension(*ir->awhParams))) { if ((fep->init_lambda >= 0) && (n_lambda_terms == 1)) { @@ -706,7 +715,7 @@ FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env nsets = nsets_dhdl + nsets_de; /* dhdl + fep differences */ - if (fep->n_lambda > 0 && (expand->elmcmove > LambdaMoveCalculation::No)) + if (haveFepLambdaMoves(*ir)) { nsets += 1; /*add fep state for expanded ensemble */ } @@ -728,7 +737,7 @@ FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env } std::vector setname(nsetsextend); - if (expand->elmcmove > LambdaMoveCalculation::No) + if (haveFepLambdaMoves(*ir)) { /* state for the fep_vals, if we have alchemical sampling */ setname[s++] = "Thermodynamic state"; @@ -781,7 +790,7 @@ FILE* open_dhdl(const char* filename, const t_inputrec* ir, const gmx_output_env * from this xvg legend. */ - if (expand->elmcmove > LambdaMoveCalculation::No) + if (haveFepLambdaMoves(*ir)) { nsetsbegin = 1; /* for including the expanded ensemble */ } @@ -838,7 +847,6 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, real tmass, const gmx_enerdata_t* enerd, const t_lambda* fep, - const t_expanded* expand, const matrix box, PTCouplingArrays ptCouplingArrays, int fep_state, @@ -1048,7 +1056,7 @@ void EnergyOutput::addDataAtEnergyStep(bool bDoDHDL, /* the current free energy state */ /* print the current state if we are doing expanded ensemble */ - if (expand->elmcmove > LambdaMoveCalculation::No) + if (haveFepLambdaMoves_) { fprintf(fp_dhdl_, " %4d", fep_state); } diff --git a/src/gromacs/mdlib/energyoutput.h b/src/gromacs/mdlib/energyoutput.h index 6c5805514e..54ae8b01b2 100644 --- a/src/gromacs/mdlib/energyoutput.h +++ b/src/gromacs/mdlib/energyoutput.h @@ -65,7 +65,6 @@ struct gmx_mtop_t; struct gmx_output_env_t; struct pull_t; struct t_ebin; -struct t_expanded; struct t_fcdata; struct t_grpopts; struct t_inputrec; @@ -168,7 +167,6 @@ public: * \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. @@ -184,7 +182,6 @@ public: real tmass, const gmx_enerdata_t* enerd, const t_lambda* fep, - const t_expanded* expand, const matrix lastbox, PTCouplingArrays ptCouplingArrays, int fep_state, @@ -402,6 +399,8 @@ private: //! The dhdl.xvg output file FILE* fp_dhdl_ = nullptr; + //! Whether the free-energy lambda moves dynamically between lambda states + bool haveFepLambdaMoves_; //! Energy components for dhdl.xvg output std::vector dE_; //! The delta U components (raw data + histogram) diff --git a/src/gromacs/mdlib/tests/energyoutput.cpp b/src/gromacs/mdlib/tests/energyoutput.cpp index c5b714a908..51aed47c64 100644 --- a/src/gromacs/mdlib/tests/energyoutput.cpp +++ b/src/gromacs/mdlib/tests/energyoutput.cpp @@ -615,7 +615,6 @@ TEST_P(EnergyOutputTest, CheckOutput) tmass_, enerdata_.get(), nullptr, - nullptr, box_, PTCouplingArrays({ state_.boxv, state_.nosehoover_xi, diff --git a/src/gromacs/mdrun/md.cpp b/src/gromacs/mdrun/md.cpp index 4f6f0263f2..0a59a4f253 100644 --- a/src/gromacs/mdrun/md.cpp +++ b/src/gromacs/mdrun/md.cpp @@ -1831,7 +1831,6 @@ void gmx::LegacySimulator::do_md() md->tmass, enerd, ir->fepvals.get(), - ir->expandedvals.get(), lastbox, PTCouplingArrays{ state->boxv, state->nosehoover_xi, diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp index 96b3666369..92d183f083 100644 --- a/src/gromacs/mdrun/mimic.cpp +++ b/src/gromacs/mdrun/mimic.cpp @@ -691,7 +691,6 @@ void gmx::LegacySimulator::do_mimic() mdatoms->tmass, enerd, ir->fepvals.get(), - ir->expandedvals.get(), state->box, PTCouplingArrays({ state->boxv, state->nosehoover_xi, diff --git a/src/gromacs/mdrun/minimize.cpp b/src/gromacs/mdrun/minimize.cpp index 42b34a45eb..eec144faf4 100644 --- a/src/gromacs/mdrun/minimize.cpp +++ b/src/gromacs/mdrun/minimize.cpp @@ -1391,7 +1391,6 @@ void LegacySimulator::do_cg() mdatoms->tmass, enerd, nullptr, - nullptr, nullBox, PTCouplingArrays(), 0, @@ -1842,7 +1841,6 @@ void LegacySimulator::do_cg() mdatoms->tmass, enerd, nullptr, - nullptr, nullBox, PTCouplingArrays(), 0, @@ -2149,7 +2147,6 @@ void LegacySimulator::do_lbfgs() mdatoms->tmass, enerd, nullptr, - nullptr, nullBox, PTCouplingArrays(), 0, @@ -2671,7 +2668,6 @@ void LegacySimulator::do_lbfgs() mdatoms->tmass, enerd, nullptr, - nullptr, nullBox, PTCouplingArrays(), 0, @@ -2964,7 +2960,6 @@ void LegacySimulator::do_steep() mdatoms->tmass, enerd, nullptr, - nullptr, nullBox, PTCouplingArrays(), 0, diff --git a/src/gromacs/mdrun/rerun.cpp b/src/gromacs/mdrun/rerun.cpp index 3ee37edb73..3b39e7a6bb 100644 --- a/src/gromacs/mdrun/rerun.cpp +++ b/src/gromacs/mdrun/rerun.cpp @@ -779,7 +779,6 @@ void gmx::LegacySimulator::do_rerun() mdatoms->tmass, enerd, ir->fepvals.get(), - ir->expandedvals.get(), state->box, PTCouplingArrays({ state->boxv, state->nosehoover_xi, diff --git a/src/gromacs/modularsimulator/energydata.cpp b/src/gromacs/modularsimulator/energydata.cpp index fd24ac18ac..8999abfc39 100644 --- a/src/gromacs/modularsimulator/energydata.cpp +++ b/src/gromacs/modularsimulator/energydata.cpp @@ -259,7 +259,6 @@ void EnergyData::doStep(Step step, Time time, bool isEnergyCalculationStep, bool mdAtoms_->mdatoms()->tmass, enerd_, inputrec_->fepvals.get(), - inputrec_->expandedvals.get(), statePropagatorData_->constPreviousBox(), PTCouplingArrays({ parrinelloRahmanBoxVelocities_ ? parrinelloRahmanBoxVelocities_() : nullMatrix, {}, -- 2.22.0