Bugfix: Stochastic rescaling with Pade approximation

include/linearAlgebra/staticMatrix/staticMatrix3x3.hpp

@@ -128,6 +128,9 @@ namespace linearAlgebra
     template <typename T>
     [[nodiscard]] StaticMatrix3x3<T> exp(const StaticMatrix3x3<T> &mat);
 
+    template <typename T>
+    [[nodiscard]] StaticMatrix3x3<T> expPade(const StaticMatrix3x3<T> &mat);
+
 }   // namespace linearAlgebra
 
 #include "staticMatrix3x3.tpp.hpp"   // DO NOT MOVE THIS LINE

include/linearAlgebra/staticMatrix/staticMatrix3x3.tpp.hpp

@@ -26,6 +26,9 @@
 
 #include "staticMatrix3x3.hpp"
 
+#include "staticMatrix3x3Class.hpp"   // for StaticMatrix3x3
+#include "vector3d.hpp"               // for Vector3D
+
 namespace linearAlgebra
 {
     /**
@@ -458,7 +461,32 @@ namespace linearAlgebra
         auto result = StaticMatrix3x3<T>(0.0);
 
         for (size_t i = 0; i < 3; ++i)
-            for (size_t j = i + 1; j < 3; ++j) result[i][j] = ::exp(mat[i][j]);
+            for (size_t j = 0; j < 3; ++j) result[i][j] = ::exp(mat[i][j]);
+
+        return result;
+    }
+
+    /**
+     * @brief Pade approximation of the exponential of a StaticMatrix3x3
+     *
+     * @link https://en.wikipedia.org/wiki/Matrix_exponential
+     * @link https://en.wikipedia.org/wiki/Pad%C3%A9_table
+     * @link https://en.wikipedia.org/wiki/Pad%C3%A9_approximant
+     * @link
+     * https://github.com/bussilab/crescale/blob/master/simplemd_anisotropic/simplemd.cpp#L351
+     */
+    template <typename T>
+    [[nodiscard]] StaticMatrix3x3<T> expPade(const StaticMatrix3x3<T> &mat)
+    {
+        auto result = StaticMatrix3x3<T>(0.0);
+
+        auto mat2 = mat * mat;
+        auto mat3 = mat2 * mat;
+
+        auto den = diagonalMatrix(1.0) - 0.5 * mat + 0.1 * mat2 - mat3 / 120.0;
+        auto num = diagonalMatrix(1.0) + 0.5 * mat + 0.1 * mat2 + mat3 / 120.0;
+
+        result = inverse(den) * num;
 
         return result;
     }

include/manostat/manostat.hpp

@@ -50,6 +50,7 @@ namespace manostat
         virtual ~Manostat() = default;
 
         void         calculatePressure(const pq::SimBox &, pq::PhysicalData &);
+        void         rotateMu(linearAlgebra::tensor3D &mu);
         virtual void applyManostat(pq::SimBox &, pq::PhysicalData &);
 
         virtual pq::ManostatType getManostatType() const;

src/manostat/berendsenManostat.cpp

@@ -73,7 +73,7 @@ SemiIsotropicBerendsenManostat::SemiIsotropicBerendsenManostat(
 )
     : BerendsenManostat(targetPressure, tau, compressibility),
       _2DAnisotropicAxis(anisotropicAxis),
-      _2DIsotropicAxes(isotropicAxes){};
+      _2DIsotropicAxes(isotropicAxes) {};
 
 /**
  * @brief apply Berendsen manostat for NPT ensemble
@@ -178,13 +178,7 @@ tensor3D FullAnisotropicBerendsenManostat::calculateMu() const
 
     auto mu = kronecker - preFactor * (pTarget - _pressureTensor);
 
-    mu[0][1] += mu[1][0];
-    mu[0][2] += mu[2][0];
-    mu[1][2] += mu[2][1];
-
-    mu[1][0] = 0.0;
-    mu[2][0] = 0.0;
-    mu[2][1] = 0.0;
+    rotateMu(mu);
 
     return mu;
 }

src/manostat/manostat.cpp

@@ -33,6 +33,7 @@ using namespace simulationBox;
 using namespace physicalData;
 using namespace constants;
 using namespace settings;
+using namespace linearAlgebra;
 
 /**
  * @brief Construct a new Manostat:: Manostat object
@@ -66,6 +67,25 @@ void Manostat::calculatePressure(const SimulationBox &box, PhysicalData &data)
     data.setPressure(_pressure);
 }
 
+/**
+ * @brief rotate mu back into upper diagonal space
+ *
+ * @details first order approximation of mu rotation according to gromacs
+ * @link https://manual.gromacs.org/current/reference-manual/algorithms/molecular-dynamics.html
+ *
+ * @param mu
+ */
+void Manostat::rotateMu(tensor3D &mu)
+{
+    mu[0][1] += mu[1][0];
+    mu[0][2] += mu[2][0];
+    mu[1][2] += mu[2][1];
+
+    mu[1][0] = 0.0;
+    mu[2][0] = 0.0;
+    mu[2][1] = 0.0;
+}
+
 /**
  * @brief apply dummy manostat for NVT ensemble
  *

src/manostat/stochasticRescalingManostat.cpp

@@ -56,7 +56,7 @@ StochasticRescalingManostat::StochasticRescalingManostat(
     : Manostat(other),
       _tau(other._tau),
       _compressibility(other._compressibility),
-      _dt(other._dt){};
+      _dt(other._dt) {};
 
 /**
  * @brief Construct a new Stochastic Rescaling Manostat:: Stochastic Rescaling
@@ -78,7 +78,7 @@ SemiIsotropicStochasticRescalingManostat::
     )
     : StochasticRescalingManostat(targetPressure, tau, compressibility),
       _2DAnisotropicAxis(anisotropicAxis),
-      _2DIsotropicAxes(isotropicAxes){};
+      _2DIsotropicAxes(isotropicAxes) {};
 
 /**
  * @brief Construct a new Stochastic Rescaling Manostat:: Stochastic Rescaling
@@ -251,8 +251,11 @@ tensor3D FullAnisotropicStochasticRescalingManostat::calculateMu(
     stochasticFactor       = ::sqrt(stochasticFactor) * random;
 
     const auto deltaP = diagonalMatrix(_targetPressure) - _pressureTensor;
+    auto       mu     = expPade(-compress * deltaP / 3.0 + stochasticFactor);
 
-    return exp(-compress * deltaP / 3.0 + stochasticFactor);
+    rotateMu(mu);
+
+    return mu;
 }
 
 /***************************