From b61516a874f3c067c0a2fce81d76bd50d6bdd72d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 18 Dec 2023 11:51:30 +0100 Subject: [PATCH] Disable NpT propagation of angular velocities --- doc/sphinx/integration.rst | 1 + src/core/integrators/velocity_verlet_npt.cpp | 9 ++++---- testsuite/python/integrator_npt.py | 23 +++++++++----------- 3 files changed, 15 insertions(+), 18 deletions(-) diff --git a/doc/sphinx/integration.rst b/doc/sphinx/integration.rst index cf521a64d6..9e2a4e3967 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 21a3c9fab0..f73592cb93 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/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 diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index 5c86768694..919c926d36 100644 --- a/testsuite/python/integrator_npt.py +++ b/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