From ad3ac5d65d474147384675610c1c901fc18258e4 Mon Sep 17 00:00:00 2001
From: "kodiakhq[bot]" <49736102+kodiakhq[bot]@users.noreply.github.com>
Date: Mon, 18 Dec 2023 16:10:18 +0000
Subject: [PATCH] Prevent propagation of angular velocities with NpT (#4843)

Fixes #4842

Description of changes:
- The current implementation of the velocity Verlet NpT propagator doesn't apply friction and noise on angular velocities. ESPResSo now throws an error when NpT encounters a rotating particle.
---
 doc/sphinx/integration.rst                   |  1 +
 src/core/integrators/velocity_verlet_npt.cpp |  9 ++++-----
 testsuite/python/integrator_npt.py           | 12 ++++++++++--
 3 files changed, 15 insertions(+), 7 deletions(-)

diff --git a/doc/sphinx/integration.rst b/doc/sphinx/integration.rst
index dfe66078b02..5517bb01abd 100644
--- a/doc/sphinx/integration.rst
+++ b/doc/sphinx/integration.rst
@@ -199,6 +199,7 @@ Notes:
 * The particle forces :math:`F` include interactions as well as a friction (:math:`\gamma^0`) and noise term (:math:`\sqrt{k_B T \gamma^0 dt} \overline{\eta}`) analogous to the terms in the :ref:`Langevin thermostat`.
 * The particle forces are only calculated in step 5 and then reused in step 1 of the next iteration. See :ref:`Velocity Verlet Algorithm` for the implications of that.
 * The NpT algorithm doesn't support :ref:`Lees-Edwards boundary conditions`.
+* The NpT algorithm doesn't support propagation of angular velocities.
 
 .. _Steepest descent:
 
diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp
index 8758b147f36..5d482b87209 100644
--- a/src/core/integrators/velocity_verlet_npt.cpp
+++ b/src/core/integrators/velocity_verlet_npt.cpp
@@ -31,7 +31,6 @@
 #include "grid.hpp"
 #include "integrate.hpp"
 #include "npt.hpp"
-#include "rotation.hpp"
 #include "thermostat.hpp"
 #include "thermostats/npt_inline.hpp"
 
@@ -164,7 +163,10 @@ void velocity_verlet_npt_propagate_vel(const ParticleRange &particles,
 
   for (auto &p : particles) {
 #ifdef ROTATION
-    propagate_omega_quat_particle(p, time_step);
+    if (p.can_rotate()) {
+      runtimeErrorMsg() << "The isotropic NpT integrator doesn't propagate "
+                           "angular velocities";
+    }
 #endif
 
     // Don't propagate translational degrees of freedom of vs
@@ -196,9 +198,6 @@ void velocity_verlet_npt_step_1(const ParticleRange &particles,
 void velocity_verlet_npt_step_2(const ParticleRange &particles,
                                 double time_step) {
   velocity_verlet_npt_propagate_vel_final(particles, time_step);
-#ifdef ROTATION
-  convert_torques_propagate_omega(particles, time_step);
-#endif
   velocity_verlet_npt_finalize_p_inst(time_step);
 }
 #endif
diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py
index 4371f9b7703..e8be111a3f1 100644
--- a/testsuite/python/integrator_npt.py
+++ b/testsuite/python/integrator_npt.py
@@ -120,6 +120,14 @@ def test_00_integrator_recovery(self):
             np.copy(system.part.all().pos),
             positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))
 
+        if espressomd.has_features("ROTATION"):
+            # propagator doesn't handle angular velocities
+            system.part.all().rotation = [3 * [False], [True, False, False]]
+            system.thermostat.set_npt(kT=1., gamma0=2., gammav=0.04, seed=42)
+            system.integrator.set_isotropic_npt(**npt_params)
+            with self.assertRaisesRegex(Exception, "The isotropic NpT integrator doesn't propagate angular velocities"):
+                system.integrator.run(1)
+
     def run_with_p3m(self, p3m, method):
         system = self.system
         npt_kwargs = {"ext_pressure": 0.001, "piston": 0.001}
@@ -147,8 +155,8 @@ def run_with_p3m(self, p3m, method):
         system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
         system.integrator.run(10)
         # check runtime warnings
-        self.system.thermostat.turn_off()
-        self.system.integrator.set_vv()
+        system.thermostat.turn_off()
+        system.integrator.set_vv()
         err_msg = f"If {method} is being used you must use the cubic box NpT"
         with self.assertRaisesRegex(RuntimeError, err_msg):
             system.integrator.set_isotropic_npt(**npt_kwargs_rectangular)