Add lambda state to dhdl.xvg with AWH fep
[alexxy/gromacs.git] / src / gromacs / mdlib / energyoutput.cpp
index 571edc54c8feb7bec602550012ebf1ec82312406..2a6c91e9b6e5f1eedbb9065689b836a337f26325 100644 (file)
@@ -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"
@@ -82,7 +83,7 @@
 #include "gromacs/utility/arraysize.h"
 #include "gromacs/utility/enumerationhelpers.h"
 #include "gromacs/utility/fatalerror.h"
-#include "gromacs/utility/mdmodulenotification.h"
+#include "gromacs/utility/mdmodulesnotifiers.h"
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/stringutil.h"
 
 
 //! Labels for energy file quantities
 //! \{
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 static const char* conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
 
-static std::array<const char*, 3> boxs_nm = { "Box-X", "Box-Y", "Box-Z" };
+static constexpr std::array<const char*, 3> boxs_nm = { "Box-X", "Box-Y", "Box-Z" };
 
-static std::array<const char*, 6> tricl_boxs_nm = { "Box-XX", "Box-YY", "Box-ZZ",
-                                                    "Box-YX", "Box-ZX", "Box-ZY" };
+static constexpr std::array<const char*, 6> tricl_boxs_nm = { "Box-XX", "Box-YY", "Box-ZZ",
+                                                              "Box-YX", "Box-ZX", "Box-ZY" };
 
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 static const char* vol_nm[] = { "Volume" };
 
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 static const char* dens_nm[] = { "Density" };
 
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 static const char* pv_nm[] = { "pV" };
 
+// NOLINTNEXTLINE(cppcoreguidelines-avoid-non-const-global-variables)
 static const char* enthalpy_nm[] = { "Enthalpy" };
 
-static std::array<const char*, 6> boxvel_nm = { "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
-                                                "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY" };
+static constexpr std::array<const char*, 6> boxvel_nm = { "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
+                                                          "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY" };
+
+const char* enumValueToString(NonBondedEnergyTerms enumValue)
+{
+    static constexpr gmx::EnumerationArray<NonBondedEnergyTerms, const char*> nonBondedEnergyTermTypeNames = {
+        "Coul-SR", "LJ-SR", "Buck-SR", "Coul-14", "LJ-14"
+    };
+    return nonBondedEnergyTermTypeNames[enumValue];
+}
 
-const char* egrp_nm[egNR + 1] = { "Coul-SR", "LJ-SR", "Buck-SR", "Coul-14", "LJ-14", nullptr };
 //! \}
 
+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
 {
 
@@ -120,28 +140,22 @@ namespace gmx
  * be written out to the .edr file.
  *
  * \todo Use more std containers.
- * \todo Remove GMX_CONSTRAINTVIR
  * \todo Write free-energy output also to energy file (after adding more tests)
  */
-EnergyOutput::EnergyOutput(ener_file*               fp_ene,
-                           const gmx_mtop_t&        mtop,
-                           const t_inputrec&        inputrec,
-                           const pull_t*            pull_work,
-                           FILE*                    fp_dhdl,
-                           bool                     isRerun,
-                           const StartingBehavior   startingBehavior,
-                           const bool               simulationsShareState,
-                           const MdModulesNotifier& mdModulesNotifier)
+EnergyOutput::EnergyOutput(ener_file*                fp_ene,
+                           const gmx_mtop_t&         mtop,
+                           const t_inputrec&         inputrec,
+                           const pull_t*             pull_work,
+                           FILE*                     fp_dhdl,
+                           bool                      isRerun,
+                           const StartingBehavior    startingBehavior,
+                           const bool                simulationsShareState,
+                           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",
                                     "Vir-YZ", "Vir-ZX", "Vir-ZY", "Vir-ZZ" };
-    static const char* sv_nm[]    = { "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
-                                   "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
-                                   "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ" };
-    static const char* fv_nm[]    = { "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
-                                   "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
-                                   "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ" };
     static const char* pres_nm[]  = { "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
                                      "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ" };
     static const char* surft_nm[] = { "#Surf*SurfTen" };
@@ -154,7 +168,7 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     char**                  gnm;
     char                    buf[256];
     const char*             bufi;
-    int                     i, j, ni, nj, n, k, kk, ncon, nset;
+    int                     i, j, ni, nj, n, ncon, nset;
     bool                    bBHAM, b14;
 
     if (EI_DYNAMICS(inputrec.eI))
@@ -174,7 +188,6 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     ncon         = gmx_mtop_ftype_count(mtop, F_CONSTR);
     nset         = gmx_mtop_ftype_count(mtop, F_SETTLE);
     bool bConstr = (ncon > 0 || nset > 0) && !isRerun;
-    bConstrVir_  = false;
     nCrmsd_      = 0;
     if (bConstr)
     {
@@ -182,7 +195,6 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
         {
             nCrmsd_ = 1;
         }
-        bConstrVir_ = (getenv("GMX_CONSTRAINTVIR") != nullptr);
     }
     else
     {
@@ -190,9 +202,9 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     }
 
     /* Energy monitoring */
-    for (i = 0; i < egNR; i++)
+    for (auto& term : bEInd_)
     {
-        bEInd_[i] = false;
+        term = false;
     }
 
     // Setting true only to those energy terms, that have active interactions and
@@ -216,7 +228,6 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
 
     bEner_[F_LJ]   = !bBHAM;
     bEner_[F_BHAM] = bBHAM;
-    bEner_[F_EQM]  = inputrec.bQMMM;
     bEner_[F_RF_EXCL] = (EEL_RF(inputrec.coulombtype) && inputrec.cutoff_scheme == CutoffScheme::Group);
     bEner_[F_COUL_RECIP]   = EEL_FULL(inputrec.coulombtype);
     bEner_[F_LJ_RECIP]     = EVDW_PME(inputrec.vdwtype);
@@ -252,11 +263,16 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     bEner_[F_ORIRESDEV]  = (gmx_mtop_ftype_count(mtop, F_ORIRES) > 0);
     bEner_[F_COM_PULL]   = ((inputrec.bPull && pull_have_potential(*pull_work)) || inputrec.bRot);
 
-    MdModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest;
-    mdModulesNotifier.simulationSetupNotifications_.notify(&mdModulesAddOutputToDensityFittingFieldRequest);
+    // Check MDModules for any energy output
+    MDModulesEnergyOutputToDensityFittingRequestChecker mdModulesAddOutputToDensityFittingFieldRequest;
+    mdModulesNotifiers.simulationSetupNotifier_.notify(&mdModulesAddOutputToDensityFittingFieldRequest);
 
     bEner_[F_DENSITYFITTING] = mdModulesAddOutputToDensityFittingFieldRequest.energyOutputToDensityFitting_;
 
+    MDModulesEnergyOutputToQMMMRequestChecker mdModulesAddOutputToQMMMFieldRequest;
+    mdModulesNotifiers.simulationSetupNotifier_.notify(&mdModulesAddOutputToQMMMFieldRequest);
+
+    bEner_[F_EQM] = mdModulesAddOutputToQMMMFieldRequest.energyOutputToQMMM_;
 
     // Counting the energy terms that will be printed and saving their names
     f_nre_ = 0;
@@ -307,11 +323,6 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
             ienthalpy_ = get_ebin_space(ebin_, 1, enthalpy_nm, unit_energy);
         }
     }
-    if (bConstrVir_)
-    {
-        isvir_ = get_ebin_space(ebin_, asize(sv_nm), sv_nm, unit_energy);
-        ifvir_ = get_ebin_space(ebin_, asize(fv_nm), fv_nm, unit_energy);
-    }
     if (bPres_)
     {
         ivir_   = get_ebin_space(ebin_, asize(vir_nm), vir_nm, unit_energy);
@@ -333,27 +344,27 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     }
 
     /* Energy monitoring */
-    for (i = 0; i < egNR; i++)
+    for (auto& term : bEInd_)
     {
-        bEInd_[i] = false;
+        term = false;
     }
-    bEInd_[egCOULSR] = true;
-    bEInd_[egLJSR]   = true;
+    bEInd_[NonBondedEnergyTerms::CoulombSR] = true;
+    bEInd_[NonBondedEnergyTerms::LJSR]      = true;
 
     if (bBHAM)
     {
-        bEInd_[egLJSR]   = false;
-        bEInd_[egBHAMSR] = true;
+        bEInd_[NonBondedEnergyTerms::LJSR]         = false;
+        bEInd_[NonBondedEnergyTerms::BuckinghamSR] = true;
     }
     if (b14)
     {
-        bEInd_[egLJ14]   = true;
-        bEInd_[egCOUL14] = true;
+        bEInd_[NonBondedEnergyTerms::LJ14]      = true;
+        bEInd_[NonBondedEnergyTerms::Coulomb14] = true;
     }
     nEc_ = 0;
-    for (i = 0; (i < egNR); i++)
+    for (auto term : bEInd_)
     {
-        if (bEInd_[i])
+        if (term)
         {
             nEc_++;
         }
@@ -362,12 +373,12 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
     nEg_ = n;
     nE_  = (n * (n + 1)) / 2;
 
-    snew(igrp_, nE_);
+    igrp_.resize(nE_);
     if (nE_ > 1)
     {
         n = 0;
         snew(gnm, nEc_);
-        for (k = 0; (k < nEc_); k++)
+        for (int k = 0; (k < nEc_); k++)
         {
             snew(gnm[k], STRLEN);
         }
@@ -376,24 +387,25 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
             ni = groups->groups[SimulationAtomGroupType::EnergyOutput][i];
             for (j = i; (j < gmx::ssize(groups->groups[SimulationAtomGroupType::EnergyOutput])); j++)
             {
-                nj = groups->groups[SimulationAtomGroupType::EnergyOutput][j];
-                for (k = kk = 0; (k < egNR); k++)
+                nj    = groups->groups[SimulationAtomGroupType::EnergyOutput][j];
+                int k = 0;
+                for (auto key : keysOf(bEInd_))
                 {
-                    if (bEInd_[k])
+                    if (bEInd_[key])
                     {
-                        sprintf(gnm[kk],
+                        sprintf(gnm[k],
                                 "%s:%s-%s",
-                                egrp_nm[k],
+                                enumValueToString(key),
                                 *(groups->groupNames[ni]),
                                 *(groups->groupNames[nj]));
-                        kk++;
+                        k++;
                     }
                 }
                 igrp_[n] = get_ebin_space(ebin_, nEc_, gnm, unit_energy);
                 n++;
             }
         }
-        for (k = 0; (k < nEc_); k++)
+        for (int k = 0; (k < nEc_); k++)
         {
             sfree(gnm[k]);
         }
@@ -437,7 +449,7 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
         mdeb_n_ = 0;
     }
 
-    snew(tmp_r_, mde_n_);
+    tmp_r_.resize(mde_n_);
     // TODO redo the group name memory management to make it more clear
     char** grpnms;
     snew(grpnms, std::max(mde_n_, mdeb_n_)); // Just in case mdeb_n_ > mde_n_
@@ -540,31 +552,19 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
         /* Currently dh histograms are only written with dynamics */
         if (EI_DYNAMICS(inputrec.eI))
         {
-            snew(dhc_, 1);
-
-            mde_delta_h_coll_init(dhc_, inputrec);
+            dhc_ = std::make_unique<t_mde_delta_h_coll>(inputrec);
         }
         fp_dhdl_ = nullptr;
-        snew(dE_, inputrec.fepvals->n_lambda);
+        dE_.resize(inputrec.fepvals->n_lambda);
     }
     else
     {
         fp_dhdl_ = fp_dhdl;
-        snew(dE_, inputrec.fepvals->n_lambda);
+        dE_.resize(inputrec.fepvals->n_lambda);
     }
     if (inputrec.bSimTemp)
     {
-        int i;
-        snew(temperatures_, inputrec.fepvals->n_lambda);
-        numTemperatures_ = inputrec.fepvals->n_lambda;
-        for (i = 0; i < inputrec.fepvals->n_lambda; i++)
-        {
-            temperatures_[i] = inputrec.simtempvals->temperatures[i];
-        }
-    }
-    else
-    {
-        numTemperatures_ = 0;
+        temperatures_ = inputrec.simtempvals->temperatures;
     }
 
     if (EI_MD(inputrec.eI) && !simulationsShareState)
@@ -575,16 +575,7 @@ EnergyOutput::EnergyOutput(ener_file*               fp_ene,
 
 EnergyOutput::~EnergyOutput()
 {
-    sfree(igrp_);
-    sfree(tmp_r_);
-    sfree(tmp_v_);
     done_ebin(ebin_);
-    done_mde_delta_h_coll(dhc_);
-    sfree(dE_);
-    if (numTemperatures_ > 0)
-    {
-        sfree(temperatures_);
-    }
 }
 
 } // namespace gmx
@@ -653,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;
@@ -696,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))
         {
@@ -724,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 */
     }
@@ -746,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";
@@ -799,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 */
         }
@@ -856,12 +847,9 @@ 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,
-                                       const tensor            svir,
-                                       const tensor            fvir,
                                        const tensor            vir,
                                        const tensor            pres,
                                        const gmx_ekindata_t*   ekind,
@@ -871,7 +859,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
     int  j, k, kk, n, gid;
     real crmsd[2], tmp6[6];
     real bs[tricl_boxs_nm.size()], vol, dens, enthalpy;
-    real eee[egNR];
+    real eee[static_cast<int>(NonBondedEnergyTerms::Count)];
     gmx::EnumerationArray<FreeEnergyPerturbationCouplingType, double> store_dhdl;
     real                                                              store_energy = 0;
     real                                                              tmp;
@@ -908,7 +896,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
             nboxs = boxs_nm.size();
         }
         vol  = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
-        dens = (tmass * AMU) / (vol * NANO * NANO * NANO);
+        dens = (tmass * gmx::c_amu) / (vol * gmx::c_nano * gmx::c_nano * gmx::c_nano);
         add_ebin(ebin_, ib_, nboxs, bs, bSum);
         add_ebin(ebin_, ivol_, 1, &vol, bSum);
         add_ebin(ebin_, idens_, 1, &dens, bSum);
@@ -917,18 +905,13 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
         {
             /* This is pV (in kJ/mol).  The pressure is the reference pressure,
                not the instantaneous pressure */
-            pv = vol * ref_p_ / PRESFAC;
+            pv = vol * ref_p_ / gmx::c_presfac;
 
             add_ebin(ebin_, ipv_, 1, &pv, bSum);
             enthalpy = pv + enerd->term[F_ETOT];
             add_ebin(ebin_, ienthalpy_, 1, &enthalpy, bSum);
         }
     }
-    if (bConstrVir_)
-    {
-        add_ebin(ebin_, isvir_, 9, svir[0], bSum);
-        add_ebin(ebin_, ifvir_, 9, fvir[0], bSum);
-    }
     if (bPres_)
     {
         add_ebin(ebin_, ivir_, 9, vir[0], bSum);
@@ -953,12 +936,12 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
     if (ekind && ekind->cosacc.cos_accel != 0)
     {
         vol  = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
-        dens = (tmass * AMU) / (vol * NANO * NANO * NANO);
+        dens = (tmass * gmx::c_amu) / (vol * gmx::c_nano * gmx::c_nano * gmx::c_nano);
         add_ebin(ebin_, ivcos_, 1, &(ekind->cosacc.vcos), bSum);
         /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
         tmp = 1
-              / (ekind->cosacc.cos_accel / (ekind->cosacc.vcos * PICO) * dens
-                 * gmx::square(box[ZZ][ZZ] * NANO / (2 * M_PI)));
+              / (ekind->cosacc.cos_accel / (ekind->cosacc.vcos * gmx::c_pico) * dens
+                 * gmx::square(box[ZZ][ZZ] * gmx::c_nano / (2 * M_PI)));
         add_ebin(ebin_, ivisc_, 1, &tmp, bSum);
     }
     if (nE_ > 1)
@@ -969,11 +952,11 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
             for (j = i; (j < nEg_); j++)
             {
                 gid = GID(i, j, nEg_);
-                for (k = kk = 0; (k < egNR); k++)
+                for (k = kk = 0; (k < static_cast<int>(NonBondedEnergyTerms::Count)); k++)
                 {
                     if (bEInd_[k])
                     {
-                        eee[kk++] = enerd->grpp.ener[k][gid];
+                        eee[kk++] = enerd->grpp.energyGroupPairTerms[k][gid];
                     }
                 }
                 add_ebin(ebin_, igrp_[n], nEc_, eee, bSum);
@@ -988,7 +971,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
         {
             tmp_r_[i] = ekind->tcstat[i].T;
         }
-        add_ebin(ebin_, itemp_, nTC_, tmp_r_, bSum);
+        add_ebin(ebin_, itemp_, nTC_, tmp_r_.data(), bSum);
 
         if (etc_ == TemperatureCoupling::NoseHoover)
         {
@@ -1006,7 +989,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
                             tmp_r_[2 * k + 1] = ptCouplingArrays.nosehoover_vxi[k];
                         }
                     }
-                    add_ebin(ebin_, itc_, mde_n_, tmp_r_, bSum);
+                    add_ebin(ebin_, itc_, mde_n_, tmp_r_.data(), bSum);
 
                     if (bMTTK_)
                     {
@@ -1019,7 +1002,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
                                 tmp_r_[2 * k + 1] = ptCouplingArrays.nhpres_vxi[k];
                             }
                         }
-                        add_ebin(ebin_, itcb_, mdeb_n_, tmp_r_, bSum);
+                        add_ebin(ebin_, itcb_, mdeb_n_, tmp_r_.data(), bSum);
                     }
                 }
                 else
@@ -1029,7 +1012,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
                         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);
+                    add_ebin(ebin_, itc_, mde_n_, tmp_r_.data(), bSum);
                 }
             }
         }
@@ -1040,7 +1023,7 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
             {
                 tmp_r_[i] = ekind->tcstat[i].lambda;
             }
-            add_ebin(ebin_, itc_, nTC_, tmp_r_, bSum);
+            add_ebin(ebin_, itc_, nTC_, tmp_r_.data(), bSum);
         }
     }
 
@@ -1054,12 +1037,12 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
         {
             /* zero for simulated tempering */
             dE_[i] = foreignTerms.deltaH(i);
-            if (numTemperatures_ > 0)
+            if (!temperatures_.empty())
             {
-                GMX_RELEASE_ASSERT(numTemperatures_ > fep_state,
+                GMX_RELEASE_ASSERT(gmx::ssize(temperatures_) > fep_state,
                                    "Number of lambdas in state is bigger then in input record");
                 GMX_RELEASE_ASSERT(
-                        numTemperatures_ >= foreignTerms.numLambdas(),
+                        gmx::ssize(temperatures_) >= 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? */
@@ -1073,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);
             }
@@ -1133,8 +1116,13 @@ void EnergyOutput::addDataAtEnergyStep(bool                    bDoDHDL,
             }
             store_energy = enerd->term[F_ETOT];
             /* store_dh is dE */
-            mde_delta_h_coll_add_dh(
-                    dhc_, static_cast<double>(fep_state), store_energy, pv, store_dhdl, dE_ + fep->lambda_start_n, time);
+            mde_delta_h_coll_add_dh(dhc_.get(),
+                                    static_cast<double>(fep_state),
+                                    store_energy,
+                                    pv,
+                                    store_dhdl,
+                                    dE_.data() + fep->lambda_start_n,
+                                    time);
         }
     }
 
@@ -1194,18 +1182,20 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene,
         nr[i] = 0;
     }
 
-    if (bOR && fcd->orires->nr > 0)
+    if (bOR && fcd->orires)
     {
         t_oriresdata& orires = *fcd->orires;
         diagonalize_orires_tensors(&orires);
-        nr[enxOR]     = orires.nr;
-        block[enxOR]  = orires.otav;
+        nr[enxOR]     = orires.numRestraints;
+        block[enxOR]  = orires.orientationsTimeAndEnsembleAv.data();
         id[enxOR]     = enxOR;
-        nr[enxORI]    = (orires.oinsl != orires.otav) ? orires.nr : 0;
-        block[enxORI] = orires.oinsl;
+        nr[enxORI]    = (orires.orientations.data() != orires.orientationsTimeAndEnsembleAv.data())
+                                ? orires.numRestraints
+                                : 0;
+        block[enxORI] = orires.orientations.data();
         id[enxORI]    = enxORI;
-        nr[enxORT]    = orires.nex * 12;
-        block[enxORT] = orires.eig;
+        nr[enxORT]    = ssize(orires.eigenOutput);
+        block[enxORT] = orires.eigenOutput.data();
         id[enxORT]    = enxORT;
     }
 
@@ -1228,10 +1218,10 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene,
             fr.block[b].id        = id[b];
             fr.block[b].sub[0].nr = nr[b];
 #if !GMX_DOUBLE
-            fr.block[b].sub[0].type = xdr_datatype_float;
+            fr.block[b].sub[0].type = XdrDataType::Float;
             fr.block[b].sub[0].fval = block[b];
 #else
-            fr.block[b].sub[0].type  = xdr_datatype_double;
+            fr.block[b].sub[0].type  = XdrDataType::Double;
             fr.block[b].sub[0].dval  = block[b];
 #endif
         }
@@ -1249,13 +1239,13 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene,
             fr.block[db].sub[0].nr     = ndisre;
             fr.block[db].sub[1].nr     = ndisre;
 #if !GMX_DOUBLE
-            fr.block[db].sub[0].type = xdr_datatype_float;
-            fr.block[db].sub[1].type = xdr_datatype_float;
+            fr.block[db].sub[0].type = XdrDataType::Float;
+            fr.block[db].sub[1].type = XdrDataType::Float;
             fr.block[db].sub[0].fval = disres.rt;
             fr.block[db].sub[1].fval = disres.rm3tav;
 #else
-            fr.block[db].sub[0].type = xdr_datatype_double;
-            fr.block[db].sub[1].type = xdr_datatype_double;
+            fr.block[db].sub[0].type = XdrDataType::Double;
+            fr.block[db].sub[1].type = XdrDataType::Double;
             fr.block[db].sub[0].dval = disres.rt;
             fr.block[db].sub[1].dval = disres.rm3tav;
 #endif
@@ -1265,13 +1255,13 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene,
         /* Free energy perturbation blocks */
         if (dhc_)
         {
-            mde_delta_h_coll_handle_block(dhc_, &fr, fr.nblock);
+            mde_delta_h_coll_handle_block(dhc_.get(), &fr, fr.nblock);
         }
 
         /* we can now free & reset the data in the blocks */
         if (dhc_)
         {
-            mde_delta_h_coll_reset(dhc_);
+            mde_delta_h_coll_reset(dhc_.get());
         }
 
         /* AWH bias blocks. */
@@ -1291,9 +1281,9 @@ void EnergyOutput::printStepToEnergyFile(ener_file* fp_ene,
     free_enxframe(&fr);
     if (log)
     {
-        if (bOR && fcd->orires->nr > 0)
+        if (bOR && fcd->orires)
         {
-            print_orires_log(log, fcd->orires);
+            print_orires_log(log, fcd->orires.get());
         }
 
         fprintf(log, "   Energies (%s)\n", unit_energy);
@@ -1357,15 +1347,6 @@ void EnergyOutput::printAverages(FILE* log, const SimulationGroups* groups)
             pr_ebin(log, ebin_, ib_, bTricl_ ? tricl_boxs_nm.size() : boxs_nm.size(), 5, eprAVER, true);
             fprintf(log, "\n");
         }
-        if (bConstrVir_)
-        {
-            fprintf(log, "   Constraint Virial (%s)\n", unit_energy);
-            pr_ebin(log, ebin_, isvir_, 9, 3, eprAVER, false);
-            fprintf(log, "\n");
-            fprintf(log, "   Force Virial (%s)\n", unit_energy);
-            pr_ebin(log, ebin_, ifvir_, 9, 3, eprAVER, false);
-            fprintf(log, "\n");
-        }
         if (bPres_)
         {
             fprintf(log, "   Total Virial (%s)\n", unit_energy);
@@ -1386,11 +1367,11 @@ void EnergyOutput::printAverages(FILE* log, const SimulationGroups* groups)
         {
             int padding = 8 - strlen(unit_energy);
             fprintf(log, "%*sEpot (%s)   ", padding, "", unit_energy);
-            for (int i = 0; (i < egNR); i++)
+            for (auto key : keysOf(bEInd_))
             {
-                if (bEInd_[i])
+                if (bEInd_[key])
                 {
-                    fprintf(log, "%12s   ", egrp_nm[i]);
+                    fprintf(log, "%12s   ", enumValueToString(key));
                 }
             }
             fprintf(log, "\n");
@@ -1453,7 +1434,7 @@ void EnergyOutput::fillEnergyHistory(energyhistory_t* enerhist) const
     }
     if (dhc_)
     {
-        mde_delta_h_coll_update_energyhistory(dhc_, enerhist);
+        mde_delta_h_coll_update_energyhistory(dhc_.get(), enerhist);
     }
 }
 
@@ -1485,7 +1466,7 @@ void EnergyOutput::restoreFromEnergyHistory(const energyhistory_t& enerhist)
     }
     if (dhc_)
     {
-        mde_delta_h_coll_restore_energyhistory(dhc_, enerhist.deltaHForeignLambdas.get());
+        mde_delta_h_coll_restore_energyhistory(dhc_.get(), enerhist.deltaHForeignLambdas.get());
     }
 }