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

@@ -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:
 

src/core/integrators/velocity_verlet_npt.cpp

@@ -30,7 +30,6 @@
 #include "errorhandling.hpp"
 #include "integrate.hpp"
 #include "npt.hpp"
-#include "rotation.hpp"
 #include "system/System.hpp"
 #include "thermostat.hpp"
 #include "thermostats/npt_inline.hpp"
@@ -170,7 +169,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
@@ -201,9 +203,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 // NPT

testsuite/python/integrator_npt.py

@@ -18,7 +18,6 @@
 #
 import unittest as ut
 import unittest_decorators as utx
-import tests_common
 import numpy as np
 import espressomd
 import espressomd.integrate
@@ -123,6 +122,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, container, p3m, method):
         system = self.system
         npt_kwargs = {"ext_pressure": 0.001, "piston": 0.001}
@@ -135,8 +142,6 @@ def run_with_p3m(self, container, p3m, method):
             pos=np.random.uniform(0, system.box_l[0], (11, 3)))
         if espressomd.has_features("P3M"):
             partcl.q = np.sign(np.arange(-5, 6))
-        if espressomd.has_features("DP3M"):
-            partcl.dip = tests_common.random_dipoles(11)
         system.non_bonded_inter[0, 0].lennard_jones.set_params(
             epsilon=1, sigma=1, cutoff=2**(1 / 6), shift=0.25)
         system.integrator.set_steepest_descent(
@@ -150,8 +155,8 @@ def run_with_p3m(self, container, 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)
@@ -164,14 +169,6 @@ def run_with_p3m(self, container, p3m, method):
         with self.assertRaisesRegex(Exception, err_msg):
             system.integrator.run(0, recalc_forces=True)
 
-    @utx.skipIfMissingFeatures(["DP3M"])
-    def test_npt_dp3m_cpu(self):
-        import espressomd.magnetostatics
-        dp3m = espressomd.magnetostatics.DipolarP3M(
-            prefactor=1.0, accuracy=1e-2, mesh=3 * [36], cao=7, r_cut=1.0,
-            alpha=2.995, tune=False)
-        self.run_with_p3m(self.system.magnetostatics, dp3m, "magnetostatics")
-
     @utx.skipIfMissingFeatures(["P3M"])
     def test_npt_p3m_cpu(self):
         import espressomd.electrostatics