Make it possible to use FEP lambda states as a reaction coordinate in AWH. Atom masse...

[alexxy/gromacs.git] / src / gromacs / awh / biasparams.cpp

diff --git a/src/gromacs/awh/biasparams.cpp b/src/gromacs/awh/biasparams.cpp
index d5afec9f5fae4a6098232c3cedfdfdfa974f1a7b..ab3b6474962ff4c977a2d1d8b321d74ee32aa6b9 100644 (file)
--- a/src/gromacs/awh/biasparams.cpp
+++ b/src/gromacs/awh/biasparams.cpp
@@ -124,20 +124,28 @@ int64_t calcCheckCoveringInterval(const AwhParams&              awhParams,
     int minNumSamplesCover = 0;
     for (size_t d = 0; d < gridAxis.size(); d++)
     {
-        GMX_RELEASE_ASSERT(dimParams[d].betak > 0,
-                           "Inverse temperature (beta) and force constant (k) should be positive.");
-        double sigma = 1.0 / std::sqrt(dimParams[d].betak);
-
-        /* The additional sample is here because to cover a discretized
-           axis of length sigma one needs two samples, one for each
-           end point. */
-        GMX_RELEASE_ASSERT(gridAxis[d].length() / sigma < std::numeric_limits<int>::max(),
-                           "The axis length in units of sigma should fit in an int");
-        int numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length() / sigma)) + 1;
-
+        int numSamplesCover;
+        if (!dimParams[d].isFepLambdaDimension())
+        {
+            GMX_RELEASE_ASSERT(
+                    dimParams[d].betak > 0,
+                    "Inverse temperature (beta) and force constant (k) should be positive.");
+            double sigma = 1.0 / std::sqrt(dimParams[d].betak);
+
+            /* The additional sample is here because to cover a discretized
+            axis of length sigma one needs two samples, one for each
+            end point. */
+            GMX_RELEASE_ASSERT(gridAxis[d].length() / sigma < std::numeric_limits<int>::max(),
+                               "The axis length in units of sigma should fit in an int");
+            numSamplesCover = static_cast<int>(std::ceil(gridAxis[d].length() / sigma)) + 1;
+        }
+        else
+        {
+            numSamplesCover = dimParams[d].numFepLambdaStates;
+        }
         /* The minimum number of samples needed for simultaneously
-           covering all axes is limited by the axis requiring most
-           samples. */
+        covering all axes is limited by the axis requiring most
+        samples. */
         minNumSamplesCover = std::max(minNumSamplesCover, numSamplesCover);
     }
 
@@ -255,12 +263,16 @@ double getInitialHistogramSizeEstimate(const std::vector<DimParams>& dimParams,
     std::vector<double> x;
     for (size_t d = 0; d < gridAxis.size(); d++)
     {
-        double axisLength = gridAxis[d].length();
+        double axisLength = gridAxis[d].isFepLambdaAxis() ? 1.0 : gridAxis[d].length();
         if (axisLength > 0)
         {
             crossingTime += awhBiasParams.dimParams[d].diffusion / (axisLength * axisLength);
             /* The sigma of the Gaussian distribution in the umbrella */
-            double sigma = 1. / std::sqrt(dimParams[d].betak);
+            double sigma = 1.;
+            if (dimParams[d].betak != 0)
+            {
+                sigma /= std::sqrt(dimParams[d].betak);
+            }
             x.push_back(sigma / axisLength);
         }
     }
@@ -334,6 +346,8 @@ BiasParams::BiasParams(const AwhParams&              awhParams,
 
     for (int d = 0; d < awhBiasParams.ndim; d++)
     {
+        /* The spacing in FEP dimensions is 1. The calculated coverRadius will be in lambda states
+         * (cf points in other dimensions). */
         double coverRadiusInNm = 0.5 * awhBiasParams.dimParams[d].coverDiameter;
         double spacing         = gridAxis[d].spacing();
         coverRadius_[d] = spacing > 0 ? static_cast<int>(std::round(coverRadiusInNm / spacing)) : 0;