Add lambda state to dhdl.xvg with AWH fep

[alexxy/gromacs.git] / src / gromacs / applied_forces / awh / read_params.cpp

diff --git a/src/gromacs/applied_forces/awh/read_params.cpp b/src/gromacs/applied_forces/awh/read_params.cpp
index 926f989286d202292b47343356fb9de7eeb57469..f2b1805209eb113ecb0642b9f1eb82b543d6668e 100644 (file)
--- a/src/gromacs/applied_forces/awh/read_params.cpp
+++ b/src/gromacs/applied_forces/awh/read_params.cpp
@@ -35,10 +35,13 @@
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+
 #include "gmxpre.h"
 
 #include "read_params.h"
 
+#include <algorithm>
+
 #include "gromacs/applied_forces/awh/awh.h"
 #include "gromacs/fileio/readinp.h"
 #include "gromacs/fileio/warninp.h"
@@ -1220,7 +1223,7 @@ static void setStateDependentAwhPullDimParams(AwhDimParams*        dimParams,
     }
 
     /* The initial coordinate value, converted to external user units. */
-    double initialCoordinate = get_pull_coord_value(pull_work, dimParams->coordinateIndex(), &pbc);
+    double initialCoordinate = get_pull_coord_value(pull_work, dimParams->coordinateIndex(), pbc);
     initialCoordinate *= pull_conversion_factor_internal2userinput(pullCoordParams);
     dimParams->setInitialCoordinate(initialCoordinate);
 }
@@ -1310,10 +1313,42 @@ void checkAwhParams(const AwhParams& awhParams, const t_inputrec& ir, warninp_t
         warning_error(wi, opt + " needs to be an integer > 0");
     }
 
-    for (int k = 0; k < awhParams.numBias(); k++)
+    bool       haveFepLambdaDim = false;
+    const auto awhBiasParams    = awhParams.awhBiasParams();
+    for (int k = 0; k < awhParams.numBias() && !haveFepLambdaDim; k++)
     {
         std::string prefixawh = formatString("awh%d", k + 1);
-        checkBiasParams(awhParams.awhBiasParams()[k], prefixawh, ir, wi);
+        checkBiasParams(awhBiasParams[k], prefixawh, ir, wi);
+        /* Check if there is a FEP lambda dimension. */
+        const auto dimParams = awhBiasParams[k].dimParams();
+        haveFepLambdaDim = std::any_of(dimParams.begin(), dimParams.end(), [](const auto& dimParam) {
+            return dimParam.coordinateProvider() == AwhCoordinateProviderType::FreeEnergyLambda;
+        });
+    }
+
+    if (haveFepLambdaDim)
+    {
+        if (awhParams.nstSampleCoord() % ir.nstcalcenergy != 0)
+        {
+            opt          = "awh-nstsample";
+            auto message = formatString(
+                    "%s (%d) should be a multiple of nstcalcenergy (%d) when using AWH for "
+                    "sampling an FEP lambda dimension",
+                    opt.c_str(),
+                    awhParams.nstSampleCoord(),
+                    ir.nstcalcenergy);
+            warning_error(wi, message);
+        }
+        if (awhParams.potential() != AwhPotentialType::Umbrella)
+        {
+            opt          = "awh-potential";
+            auto message = formatString(
+                    "%s (%s) must be set to %s when using AWH for sampling an FEP lambda dimension",
+                    opt.c_str(),
+                    enumValueToString(awhParams.potential()),
+                    enumValueToString(AwhPotentialType::Umbrella));
+            warning_error(wi, message);
+        }
     }
 
     if (ir.init_step != 0)
@@ -1322,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