Add lambda state to dhdl.xvg with AWH fep

diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp
index 6fa98dc4c38b96b53bcf8475569e4809c9f2b0c2..f2b1805209eb113ecb0642b9f1eb82b543d6668e 100644 (file)
--- 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 22247ae48065f2032f2b99b6a286719026d6e606..9e33649a60b7025990317aaf5c9c3c2bcccb5ba0 100644 (file)
--- 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 469fac5336bf840d1959f598996ad726c57473d2..2a6c91e9b6e5f1eedbb9065689b836a337f26325 100644 (file)
--- a/src/gromacs/mdlib/energyoutput.cpp
+++ b/src/gromacs/mdlib/energyoutput.cpp
@@ -56,6 +56,7 @@
 #include <string>
 
 #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<std::string> 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 6c5805514e5a7ed919e8f8e2a90ede3dbf033209..54ae8b01b278139a37aa3a74bf52cae7ac09b217 100644 (file)
--- 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<double> 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 c5b714a90829adea034729e25270516748254b15..51aed47c64def29dd2bdbc518983c654fada1b64 100644 (file)
--- 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 4f6f0263f27a5faab674c33718711b7148cbfeda..0a59a4f2535bca973a34366ef336b2dce519ee23 100644 (file)
--- 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 96b3666369ff0775d219c8b00b1150d6aeabd5a8..92d183f08328fb497629d4f8d73c2118a0432a33 100644 (file)
--- 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 42b34a45ebb65dfad3a75720547581f460664f7d..eec144faf42604fa6ae6a41a0607abc4bf580253 100644 (file)
--- 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 3ee37edb73684a536091234b68cba7c93d402561..3b39e7a6bbb7bafebfc5acb08892165944feb079 100644 (file)
--- 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 fd24ac18ac6662177a3efc57913870e670e1be7e..8999abfc39920ba6f6549697776a673732d31157 100644 (file)
--- 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,
                                {},