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 3ddb5a2cb5aeb1d9e25d34427a6b15761908fe41..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"
@@ -626,11 +629,7 @@ void checkInputConsistencyAwh(const AwhParams& awhParams, warninp_t wi)
 
 } // namespace
 
-AwhDimParams::AwhDimParams(std::vector<t_inpfile>* inp,
-                           const std::string&      prefix,
-                           const t_inputrec&       ir,
-                           warninp_t               wi,
-                           bool                    bComment)
+AwhDimParams::AwhDimParams(std::vector<t_inpfile>* inp, const std::string& prefix, warninp_t wi, bool bComment)
 {
     std::string opt;
     if (bComment)
@@ -710,8 +709,6 @@ AwhDimParams::AwhDimParams(std::vector<t_inpfile>* inp,
     }
     opt            = prefix + "-cover-diameter";
     coverDiameter_ = get_ereal(inp, opt, 0, wi);
-
-    checkDimParams(prefix, *this, ir, wi);
 }
 
 AwhDimParams::AwhDimParams(ISerializer* serializer)
@@ -744,11 +741,7 @@ void AwhDimParams::serialize(ISerializer* serializer)
     serializer->doDouble(&coverDiameter_);
 }
 
-AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp,
-                             const std::string&      prefix,
-                             const t_inputrec&       ir,
-                             warninp_t               wi,
-                             bool                    bComment)
+AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp, const std::string& prefix, warninp_t wi, bool bComment)
 {
     if (bComment)
     {
@@ -830,9 +823,8 @@ AwhBiasParams::AwhBiasParams(std::vector<t_inpfile>* inp,
     {
         bComment              = bComment && d == 0;
         std::string prefixdim = prefix + formatString("-dim%d", d + 1);
-        dimParams_.emplace_back(inp, prefixdim, ir, wi, bComment);
+        dimParams_.emplace_back(inp, prefixdim, wi, bComment);
     }
-    checkBiasParams(*this, prefix, ir, wi);
 }
 
 AwhBiasParams::AwhBiasParams(ISerializer* serializer)
@@ -882,7 +874,7 @@ void AwhBiasParams::serialize(ISerializer* serializer)
     }
 }
 
-AwhParams::AwhParams(std::vector<t_inpfile>* inp, const t_inputrec& ir, warninp_t wi)
+AwhParams::AwhParams(std::vector<t_inpfile>* inp, warninp_t wi)
 {
     std::string opt;
 
@@ -934,9 +926,8 @@ AwhParams::AwhParams(std::vector<t_inpfile>* inp, const t_inputrec& ir, warninp_
     {
         bool        bComment  = (k == 0);
         std::string prefixawh = formatString("awh%d", k + 1);
-        awhBiasParams_.emplace_back(inp, prefixawh, ir, wi, bComment);
+        awhBiasParams_.emplace_back(inp, prefixawh, wi, bComment);
     }
-    checkAwhParams(*this, ir, wi);
     checkInputConsistencyAwh(*this, wi);
 }
 
@@ -1232,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);
 }
@@ -1322,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)
@@ -1334,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