Merge branch release-2021 into merge-2021-into-master
[alexxy/gromacs.git] / src / gromacs / mdlib / enerdata_utils.cpp
index 0ebef33d1e74741c0de981adab7041fe84cc41aa..f76dca56a1e3fafd5398f69f6b992b750f26d014 100644 (file)
 #include "gromacs/mdtypes/commrec.h"
 #include "gromacs/mdtypes/enerdata.h"
 #include "gromacs/mdtypes/inputrec.h"
+#include "gromacs/utility/enumerationhelpers.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
 ForeignLambdaTerms::ForeignLambdaTerms(int numLambdas) :
-    numLambdas_(numLambdas),
-    energies_(1 + numLambdas),
-    dhdl_(1 + numLambdas)
+    numLambdas_(numLambdas), energies_(1 + numLambdas), dhdl_(1 + numLambdas)
 {
 }
 
@@ -81,9 +80,7 @@ void ForeignLambdaTerms::zeroAllTerms()
 }
 
 gmx_enerdata_t::gmx_enerdata_t(int numEnergyGroups, int numFepLambdas) :
-    grpp(numEnergyGroups),
-    foreignLambdaTerms(numFepLambdas),
-    foreign_grpp(numEnergyGroups)
+    grpp(numEnergyGroups), foreignLambdaTerms(numFepLambdas), foreign_grpp(numEnergyGroups)
 {
 }
 
@@ -106,15 +103,15 @@ void sum_epot(const gmx_grppairener_t& grpp, real* epot)
     int i;
 
     /* Accumulate energies */
-    epot[F_COUL_SR] = sum_v(grpp.nener, grpp.ener[egCOULSR]);
-    epot[F_LJ]      = sum_v(grpp.nener, grpp.ener[egLJSR]);
-    epot[F_LJ14]    = sum_v(grpp.nener, grpp.ener[egLJ14]);
-    epot[F_COUL14]  = sum_v(grpp.nener, grpp.ener[egCOUL14]);
+    epot[F_COUL_SR] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::CoulombSR]);
+    epot[F_LJ]      = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJSR]);
+    epot[F_LJ14]    = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::LJ14]);
+    epot[F_COUL14]  = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::Coulomb14]);
 
     /* lattice part of LR doesnt belong to any group
      * and has been added earlier
      */
-    epot[F_BHAM] = sum_v(grpp.nener, grpp.ener[egBHAMSR]);
+    epot[F_BHAM] = sum_v(grpp.nener, grpp.energyGroupPairTerms[NonBondedEnergyTerms::BuckinghamSR]);
 
     epot[F_EPOT] = 0;
     for (i = 0; (i < F_EPOT); i++)
@@ -127,7 +124,7 @@ void sum_epot(const gmx_grppairener_t& grpp, real* epot)
 }
 
 // Adds computed dH/dlambda contribution i to the requested output
-static void set_dhdl_output(gmx_enerdata_t* enerd, int i, const t_lambda& fepvals)
+static void set_dhdl_output(gmx_enerdata_t* enerd, FreeEnergyPerturbationCouplingType i, const t_lambda& fepvals)
 {
     if (fepvals.separate_dvdl[i])
     {
@@ -137,18 +134,23 @@ static void set_dhdl_output(gmx_enerdata_t* enerd, int i, const t_lambda& fepval
         int index;
         switch (i)
         {
-            case (efptMASS): index = F_DKDL; break;
-            case (efptCOUL): index = F_DVDL_COUL; break;
-            case (efptVDW): index = F_DVDL_VDW; break;
-            case (efptBONDED): index = F_DVDL_BONDED; break;
-            case (efptRESTRAINT): index = F_DVDL_RESTRAINT; break;
+            case (FreeEnergyPerturbationCouplingType::Mass): index = F_DKDL; break;
+            case (FreeEnergyPerturbationCouplingType::Coul): index = F_DVDL_COUL; break;
+            case (FreeEnergyPerturbationCouplingType::Vdw): index = F_DVDL_VDW; break;
+            case (FreeEnergyPerturbationCouplingType::Bonded): index = F_DVDL_BONDED; break;
+            case (FreeEnergyPerturbationCouplingType::Restraint): index = F_DVDL_RESTRAINT; break;
             default: index = F_DVDL; break;
         }
         enerd->term[index] = enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
         if (debug)
         {
-            fprintf(debug, "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n", efpt_names[i], i,
-                    enerd->term[index], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
+            fprintf(debug,
+                    "dvdl-%s[%2d]: %f: non-linear %f + linear %f\n",
+                    enumValueToString(i),
+                    static_cast<int>(i),
+                    enerd->term[index],
+                    enerd->dvdl_nonlin[i],
+                    enerd->dvdl_lin[i]);
         }
     }
     else
@@ -156,8 +158,13 @@ static void set_dhdl_output(gmx_enerdata_t* enerd, int i, const t_lambda& fepval
         enerd->term[F_DVDL] += enerd->dvdl_lin[i] + enerd->dvdl_nonlin[i];
         if (debug)
         {
-            fprintf(debug, "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n", efpt_names[0], i,
-                    enerd->term[F_DVDL], enerd->dvdl_nonlin[i], enerd->dvdl_lin[i]);
+            fprintf(debug,
+                    "dvd-%sl[%2d]: %f: non-linear %f + linear %f\n",
+                    enumValueToString(FreeEnergyPerturbationCouplingType::Fep),
+                    static_cast<int>(i),
+                    enerd->term[F_DVDL],
+                    enerd->dvdl_nonlin[i],
+                    enerd->dvdl_lin[i]);
         }
     }
 }
@@ -180,7 +187,7 @@ void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const doub
     }
 
     double dvdl_lin = 0;
-    for (int i = 0; i < efptNR; i++)
+    for (int i = 0; i < static_cast<int>(FreeEnergyPerturbationCouplingType::Count); i++)
     {
         dvdl_lin += dvdlLinear[i];
     }
@@ -218,15 +225,15 @@ void ForeignLambdaTerms::finalizePotentialContributions(gmx::ArrayRef<const doub
 
 void accumulatePotentialEnergies(gmx_enerdata_t* enerd, gmx::ArrayRef<const real> lambda, const t_lambda* fepvals)
 {
-    sum_epot(enerd->grpp, enerd->term);
+    sum_epot(enerd->grpp, enerd->term.data());
 
     if (fepvals)
     {
         enerd->term[F_DVDL] = 0.0;
-        for (int i = 0; i < efptNR; i++)
+        for (auto i : gmx::EnumerationWrapper<FreeEnergyPerturbationCouplingType>{})
         {
             // Skip kinetic terms here, as those are not available here yet
-            if (i != efptMASS)
+            if (i != FreeEnergyPerturbationCouplingType::Mass)
             {
                 set_dhdl_output(enerd, i, *fepvals);
             }
@@ -252,7 +259,7 @@ void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real>
 
     // Treat current lambda, the deltaH contribution is 0 as delta-lambda=0 for the current lambda
     accumulateKinetic(0, 0.0, energyTerms[F_DVDL_CONSTR]);
-    if (!fepvals.separate_dvdl[efptMASS])
+    if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
     {
         accumulateKinetic(0, 0.0, energyTerms[F_DKDL]);
     }
@@ -265,13 +272,17 @@ void ForeignLambdaTerms::finalizeKineticContributions(gmx::ArrayRef<const real>
          * a linear extrapolation. This is an approximation, but usually
          * quite accurate since constraints change little between lambdas.
          */
-        const int    lambdaIndex = (fepvals.separate_dvdl[efptBONDED] ? efptBONDED : efptFEP);
-        const double dlam        = fepvals.all_lambda[lambdaIndex][i] - lambda[lambdaIndex];
+        const FreeEnergyPerturbationCouplingType lambdaIndex =
+                (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded]
+                         ? FreeEnergyPerturbationCouplingType::Bonded
+                         : FreeEnergyPerturbationCouplingType::Fep);
+        const double dlam = fepvals.all_lambda[lambdaIndex][i] - lambda[static_cast<int>(lambdaIndex)];
         accumulateKinetic(1 + i, dlam * energyTerms[F_DVDL_CONSTR], energyTerms[F_DVDL_CONSTR]);
 
-        if (!fepvals.separate_dvdl[efptMASS])
+        if (!fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Mass])
         {
-            const double dlam = fepvals.all_lambda[efptMASS][i] - lambda[efptMASS];
+            const double dlam = fepvals.all_lambda[FreeEnergyPerturbationCouplingType::Mass][i]
+                                - lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Mass)];
             accumulateKinetic(1 + i, dlam * energyTerms[F_DKDL], energyTerms[F_DKDL]);
         }
     }
@@ -281,7 +292,7 @@ void accumulateKineticLambdaComponents(gmx_enerdata_t*           enerd,
                                        gmx::ArrayRef<const real> lambda,
                                        const t_lambda&           fepvals)
 {
-    if (fepvals.separate_dvdl[efptBONDED])
+    if (fepvals.separate_dvdl[FreeEnergyPerturbationCouplingType::Bonded])
     {
         enerd->term[F_DVDL_BONDED] += enerd->term[F_DVDL_CONSTR];
     }
@@ -291,10 +302,10 @@ void accumulateKineticLambdaComponents(gmx_enerdata_t*           enerd,
     }
 
     // Add computed mass dH/dlambda contribution to the requested output
-    set_dhdl_output(enerd, efptMASS, fepvals);
+    set_dhdl_output(enerd, FreeEnergyPerturbationCouplingType::Mass, fepvals);
 
-    enerd->foreignLambdaTerms.finalizeKineticContributions(enerd->term, enerd->dvdl_lin[efptMASS],
-                                                           lambda, fepvals);
+    enerd->foreignLambdaTerms.finalizeKineticContributions(
+            enerd->term, enerd->dvdl_lin[FreeEnergyPerturbationCouplingType::Mass], lambda, fepvals);
 
     /* The constrain contribution is now included in other terms, so clear it */
     enerd->term[F_DVDL_CONSTR] = 0;
@@ -306,11 +317,11 @@ void reset_foreign_enerdata(gmx_enerdata_t* enerd)
 
     /* First reset all foreign energy components.  Foreign energies always called on
        neighbor search steps */
-    for (i = 0; (i < egNR); i++)
+    for (i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
     {
         for (j = 0; (j < enerd->grpp.nener); j++)
         {
-            enerd->foreign_grpp.ener[i][j] = 0.0;
+            enerd->foreign_grpp.energyGroupPairTerms[i][j] = 0.0;
         }
     }
 
@@ -323,7 +334,7 @@ void reset_foreign_enerdata(gmx_enerdata_t* enerd)
 
 void reset_dvdl_enerdata(gmx_enerdata_t* enerd)
 {
-    for (int i = 0; i < efptNR; i++)
+    for (auto i : keysOf(enerd->dvdl_lin))
     {
         enerd->dvdl_lin[i]    = 0.0;
         enerd->dvdl_nonlin[i] = 0.0;
@@ -335,11 +346,11 @@ void reset_enerdata(gmx_enerdata_t* enerd)
     int i, j;
 
     /* First reset all energy components. */
-    for (i = 0; (i < egNR); i++)
+    for (i = 0; (i < static_cast<int>(NonBondedEnergyTerms::Count)); i++)
     {
         for (j = 0; (j < enerd->grpp.nener); j++)
         {
-            enerd->grpp.ener[i][j] = 0.0_real;
+            enerd->grpp.energyGroupPairTerms[i][j] = 0.0_real;
         }
     }