Bugfix/mshake

CHANGELOG.md

@@ -12,6 +12,7 @@ All notable changes to this project will be documented in this file.
 ### Bug Fixes
 
 - Eigen version finally fixed to 5.0.0 (latest aka master broken on 28.09.25)
+- Fix problems with convergence and indices in MShake
 
 <!-- insertion marker -->
 ## [v0.6.2](https://github.com/MolarVerse/PQ/releases/tag/v0.6.2) - 2025-08-22

include/constraints/mShake.hpp

@@ -56,7 +56,7 @@ namespace constraints
 
         void initMShake();
         void initMShakeReferences();
-        void applyMShake(const double, pq::SimBox &);
+        void applyMShake(const double, const size_t, pq::SimBox &);
         void applyMRattle(pq::SimBox &);
 
         [[nodiscard]] size_t calcNumberOfMShakeMolecules(pq::SimBox &) const;

src/constraints/constraints.cpp

@@ -147,7 +147,7 @@ void Constraints::_applyShake(SimulationBox &simBox)
 void Constraints::_applyMShake(SimulationBox &simulationBox)
 {
     startTimingsSection("MShake - Shake");
-    _mShake.applyMShake(_shakeTolerance, simulationBox);
+    _mShake.applyMShake(_shakeTolerance, _shakeMaxIter, simulationBox);
     stopTimingsSection("MShake - Shake");
 }
 

src/constraints/mShake.cpp

@@ -129,7 +129,7 @@ void MShake::initMShakeReferences()
  * @param simBox
  *
  */
-void MShake::applyMShake(const double shakeTolerance, SimulationBox &simBox)
+void MShake::applyMShake(const double shakeTolerance, const size_t shakeIterations, SimulationBox &simBox)
 {
     auto &molecules = simBox.getMolecules();
 
@@ -212,11 +212,15 @@ void MShake::applyMShake(const double shakeTolerance, SimulationBox &simBox)
                 ++index_ij;
             }
 
+        size_t iteration = 0;
+
         while (true)
         {
             auto converged = true;
             index_ij       = 0;
 
+            iteration++;
+
             /*****************************************
              * fill (nBonds x nBonds) m-Shake matrix *
              *****************************************/
@@ -232,7 +236,7 @@ void MShake::applyMShake(const double shakeTolerance, SimulationBox &simBox)
 
                     for (size_t k = 0; k < nAtoms - 1; ++k)
                     {
-                        for (size_t l = 0; l < nAtoms - 1; ++l)
+                        for (size_t l = k + 1; l < nAtoms; ++l)
                         {
                             const auto mShakeElement = calcMatrixElement(
                                 {i, j, k, l},
@@ -331,8 +335,19 @@ void MShake::applyMShake(const double shakeTolerance, SimulationBox &simBox)
                 }
             }
 
-            if (!converged)
+            if (converged)
                 break;
+
+            if (iteration >= shakeIterations)
+            {
+                throw customException::MShakeException(
+                    std::format(
+                        "M-Shake did not converge within {} iterations for molecule type {}",
+                        shakeIterations,
+                        moltype
+                    )
+                );
+            }
         }
 
         for (size_t i = 0; i < nAtoms; ++i)