From ca4b75d5fa3fcccbd3c66b3be0d2c26e7c5e66b5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Thu, 17 Mar 2022 17:34:52 +0100
Subject: [PATCH 01/10] testsuite: Fix regressions in test cases

Check Particle setters/getters on empty myconfig files. Fix
incorrect position of the LB boundary in the VTK test. Make
VirtualSitesRelative tests independent from each other and
faster using less particles (runtime scales with n_part^2).
Fix MAX_NUM_PROC regression in the checkpoint test CMake logic.
---
 samples/grand_canonical.py                    |  2 +-
 src/core/unit_tests/Particle_test.cpp         | 20 +++++---
 testsuite/python/CMakeLists.txt               | 14 ++---
 testsuite/python/field_test.py                | 12 ++---
 .../interactions_non-bonded_interface.py      |  7 ++-
 testsuite/python/lb_vtk.py                    |  2 +-
 testsuite/python/particle.py                  | 46 ++++++++---------
 testsuite/python/random_pairs.py              |  2 +-
 testsuite/python/reaction_methods.py          | 51 ++++++-------------
 testsuite/python/virtual_sites_relative.py    | 27 +++++-----
 10 files changed, 83 insertions(+), 100 deletions(-)

diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py
index 8326e9c8d76..7fe1c2b4ad1 100644
--- a/samples/grand_canonical.py
+++ b/samples/grand_canonical.py
@@ -38,7 +38,7 @@
 import espressomd.reaction_ensemble
 import espressomd.electrostatics
 
-required_features = ["P3M", "EXTERNAL_FORCES", "WCA"]
+required_features = ["P3M", "WCA"]
 espressomd.assert_features(required_features)
 
 parser = argparse.ArgumentParser(epilog=__doc__ + epilog)
diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp
index af08ba21579..2726a7b7621 100644
--- a/src/core/unit_tests/Particle_test.cpp
+++ b/src/core/unit_tests/Particle_test.cpp
@@ -243,8 +243,6 @@ BOOST_AUTO_TEST_CASE(rattle_constructors) {
 }
 #endif // BOND_CONSTRAINT
 
-#ifdef EXTERNAL_FORCES
-#ifdef ROTATION
 BOOST_AUTO_TEST_CASE(particle_bitfields) {
   auto p = Particle();
 
@@ -255,24 +253,32 @@ BOOST_AUTO_TEST_CASE(particle_bitfields) {
   BOOST_CHECK(not p.can_rotate_around(1));
 
   // check setting of one axis
+#ifdef EXTERNAL_FORCES
   p.set_fixed_along(1, true);
-  p.set_can_rotate_around(1, true);
   BOOST_CHECK(p.is_fixed_along(1));
-  BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.has_fixed_coordinates());
+#endif
+#ifdef ROTATION
+  p.set_can_rotate_around(1, true);
+  BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.can_rotate());
+#endif
 
   // check that unsetting is properly registered
+#ifdef EXTERNAL_FORCES
   p.set_fixed_along(1, false);
-  p.set_can_rotate_around(1, false);
   BOOST_CHECK(not p.has_fixed_coordinates());
+#endif
+#ifdef ROTATION
+  p.set_can_rotate_around(1, false);
   BOOST_CHECK(not p.can_rotate());
+#endif
 
   // check setting of all flags at once
+#ifdef ROTATION
   p.set_can_rotate_all_axes();
   BOOST_CHECK(p.can_rotate_around(0));
   BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.can_rotate_around(2));
+#endif
 }
-#endif // ROTATION
-#endif // EXTERNAL_FORCES
\ No newline at end of file
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index 77cd2112bec..38a9b15527f 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -14,7 +14,7 @@ function(PYTHON_TEST)
   set(TEST_FILE ${TEST_FILE_CONFIGURED})
 
   if(NOT DEFINED TEST_MAX_NUM_PROC)
-    set(TEST_MAX_NUM_PROC 1)
+    set(TEST_MAX_NUM_PROC 4)
   endif()
 
   if(${TEST_MAX_NUM_PROC} GREATER ${TEST_NP})
@@ -61,20 +61,20 @@ function(CHECKPOINT_TEST)
   endif()
   if(TEST_SUFFIX)
     set(TEST_ARGUMENTS "Test_suffix_${TEST_SUFFIX}__${TEST_MODES}")
-    set(TEST_SUFFIX "${TEST_MODES}_${TEST_SUFFIX}")
+    set(TEST_SUFFIX "_${TEST_MODES}_${TEST_SUFFIX}")
   else()
     set(TEST_ARGUMENTS "Test__${TEST_MODES}")
-    set(TEST_SUFFIX "${TEST_MODES}")
+    set(TEST_SUFFIX "_${TEST_MODES}")
   endif()
   python_test(
-    FILE save_checkpoint.py MAX_NUM_PROC ${TEST_NPROCS} LABELS ${TEST_LABELS}
-    SUFFIX ${TEST_SUFFIX} ARGUMENTS ${TEST_ARGUMENTS} DEPENDENCIES
-    unittest_generator.py)
+    FILE save_checkpoint.py MAX_NUM_PROC ${TEST_MAX_NUM_PROC} LABELS
+    ${TEST_LABELS} SUFFIX ${TEST_SUFFIX} ARGUMENTS ${TEST_ARGUMENTS}
+    DEPENDENCIES unittest_generator.py)
   python_test(
     FILE
     test_checkpoint.py
     MAX_NUM_PROC
-    ${TEST_NPROCS}
+    ${TEST_MAX_NUM_PROC}
     LABELS
     ${TEST_LABELS}
     SUFFIX
diff --git a/testsuite/python/field_test.py b/testsuite/python/field_test.py
index 327b8ebaaf9..27b5aa13e11 100644
--- a/testsuite/python/field_test.py
+++ b/testsuite/python/field_test.py
@@ -149,11 +149,10 @@ def test_homogeneous_flow_field(self):
 
     def test_potential_field(self):
         h = np.array([.2, .2, .2])
-        box = np.array([10., 10., 10.])
         scaling = 2.6
 
         field_data = espressomd.constraints.PotentialField.field_from_fn(
-            box, h, self.potential)
+            self.system.box_l, h, self.potential)
 
         F = espressomd.constraints.PotentialField(
             field=field_data,
@@ -186,10 +185,9 @@ def test_potential_field(self):
     @utx.skipIfMissingFeatures("ELECTROSTATICS")
     def test_electric_potential_field(self):
         h = np.array([.2, .2, .2])
-        box = np.array([10., 10., 10.])
 
         field_data = espressomd.constraints.ElectricPotential.field_from_fn(
-            box, h, self.potential)
+            self.system.box_l, h, self.potential)
 
         F = espressomd.constraints.ElectricPotential(
             field=field_data, grid_spacing=h)
@@ -213,11 +211,10 @@ def test_electric_potential_field(self):
 
     def test_force_field(self):
         h = np.array([.8, .8, .8])
-        box = np.array([10., 10., 10.])
         scaling = 2.6
 
         field_data = espressomd.constraints.ForceField.field_from_fn(
-            box, h, self.force)
+            self.system.box_l, h, self.force)
 
         F = espressomd.constraints.ForceField(
             field=field_data,
@@ -245,11 +242,10 @@ def test_force_field(self):
 
     def test_flow_field(self):
         h = np.array([.8, .8, .8])
-        box = np.array([10., 10., 10.])
         gamma = 2.6
 
         field_data = espressomd.constraints.FlowField.field_from_fn(
-            box, h, self.force)
+            self.system.box_l, h, self.force)
 
         F = espressomd.constraints.FlowField(
             field=field_data,
diff --git a/testsuite/python/interactions_non-bonded_interface.py b/testsuite/python/interactions_non-bonded_interface.py
index d893f28ef79..b842fb8394b 100644
--- a/testsuite/python/interactions_non-bonded_interface.py
+++ b/testsuite/python/interactions_non-bonded_interface.py
@@ -25,7 +25,12 @@
 
 
 class Non_bonded_interactionsTests(ut.TestCase):
-    system = espressomd.System(box_l=[20.0, 20.0, 20.0])
+    system = espressomd.System(box_l=[30.0, 30.0, 30.0])
+
+    def tearDown(self):
+        if espressomd.has_features(["LENNARD_JONES"]):
+            self.system.non_bonded_inter[0, 0].lennard_jones.set_params(
+                epsilon=0., sigma=0., cutoff=0., shift=0.)
 
     def intersMatch(self, inType, outInter, inParams, outParams, msg_long):
         """Check, if the interaction type set and gotten back as well as the
diff --git a/testsuite/python/lb_vtk.py b/testsuite/python/lb_vtk.py
index 7ed32ebef82..8cd8e065c02 100644
--- a/testsuite/python/lb_vtk.py
+++ b/testsuite/python/lb_vtk.py
@@ -56,7 +56,7 @@ def set_lbf(self):
             self.system.lbboundaries.add(espressomd.lbboundaries.LBBoundary(
                 shape=espressomd.shapes.Wall(normal=[1, 0, 0], dist=1.5)))
             self.system.lbboundaries.add(espressomd.lbboundaries.LBBoundary(
-                shape=espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-10.5)))
+                shape=espressomd.shapes.Wall(normal=[-1, 0, 0], dist=-8.5)))
         return lbf
 
     def parse_vtk(self, filepath, name, shape):
diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py
index 812b60e0d0c..43f704c105e 100644
--- a/testsuite/python/particle.py
+++ b/testsuite/python/particle.py
@@ -182,29 +182,29 @@ def test_gamma_rot_single(self):
 
     if espressomd.has_features(["VIRTUAL_SITES"]):
         test_virtual = generateTestForScalarProperty("virtual", 1)
-    if espressomd.has_features(["VIRTUAL_SITES_RELATIVE"]):
-        def test_yy_vs_relative(self):
-            self.system.part.add(id=0, pos=(0, 0, 0))
-            p1 = self.system.part.add(id=1, pos=(0, 0, 0))
-            p1.vs_relative = (0, 5.0, (0.5, -0.5, -0.5, -0.5))
-            p1.vs_quat = [1, 2, 3, 4]
-            np.testing.assert_array_equal(
-                p1.vs_quat, [1, 2, 3, 4])
-            res = p1.vs_relative
-            self.assertEqual(res[0], 0, "vs_relative: " + res.__str__())
-            self.assertEqual(res[1], 5.0, "vs_relative: " + res.__str__())
-            np.testing.assert_allclose(
-                res[2], np.array((0.5, -0.5, -0.5, -0.5)),
-                err_msg="vs_relative: " + res.__str__(), atol=self.tol)
-            # check exceptions
-            with self.assertRaisesRegex(ValueError, "needs input in the form"):
-                p1.vs_relative = (0, 5.0)
-            with self.assertRaisesRegex(ValueError, "particle id has to be given as an int"):
-                p1.vs_relative = ('0', 5.0, (1, 0, 0, 0))
-            with self.assertRaisesRegex(ValueError, "distance has to be given as a float"):
-                p1.vs_relative = (0, '5', (1, 0, 0, 0))
-            with self.assertRaisesRegex(ValueError, "quaternion has to be given as a tuple of 4 floats"):
-                p1.vs_relative = (0, 5.0, (1, 0, 0))
+
+    @utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE"])
+    def test_vs_relative(self):
+        self.system.part.add(id=0, pos=(0, 0, 0))
+        p1 = self.system.part.add(id=1, pos=(0, 0, 0))
+        p1.vs_relative = (0, 5.0, (0.5, -0.5, -0.5, -0.5))
+        p1.vs_quat = [1, 2, 3, 4]
+        np.testing.assert_array_equal(p1.vs_quat, [1, 2, 3, 4])
+        res = p1.vs_relative
+        self.assertEqual(res[0], 0, f"vs_relative: {res}")
+        self.assertEqual(res[1], 5.0, f"vs_relative: {res}")
+        np.testing.assert_allclose(
+            res[2], np.array((0.5, -0.5, -0.5, -0.5)),
+            err_msg=f"vs_relative: {res}", atol=self.tol)
+        # check exceptions
+        with self.assertRaisesRegex(ValueError, "needs input in the form"):
+            p1.vs_relative = (0, 5.0)
+        with self.assertRaisesRegex(ValueError, "particle id has to be given as an int"):
+            p1.vs_relative = ('0', 5.0, (1, 0, 0, 0))
+        with self.assertRaisesRegex(ValueError, "distance has to be given as a float"):
+            p1.vs_relative = (0, '5', (1, 0, 0, 0))
+        with self.assertRaisesRegex(ValueError, "quaternion has to be given as a tuple of 4 floats"):
+            p1.vs_relative = (0, 5.0, (1, 0, 0))
 
     @utx.skipIfMissingFeatures("DIPOLES")
     def test_contradicting_properties_dip_dipm(self):
diff --git a/testsuite/python/random_pairs.py b/testsuite/python/random_pairs.py
index 5c4b19d2e7c..2d9a31354ef 100644
--- a/testsuite/python/random_pairs.py
+++ b/testsuite/python/random_pairs.py
@@ -91,7 +91,7 @@ def check_n_squared(self, n2_pairs):
 
     def test(self):
         periods = [0, 1]
-        self.system.periodicity = True, True, True
+        self.system.periodicity = [True, True, True]
         check_non_bonded_loop_trace(self.system)
 
         for periodicity in itertools.product(periods, periods, periods):
diff --git a/testsuite/python/reaction_methods.py b/testsuite/python/reaction_methods.py
index 4a568aea974..d82a5119698 100644
--- a/testsuite/python/reaction_methods.py
+++ b/testsuite/python/reaction_methods.py
@@ -72,28 +72,20 @@ def check_reaction_parameters(reactions, parameters):
             method.exclusion_range,
             exclusion_range,
             delta=1e-10)
-        if exclusion_radius_per_type:  # Avoid this test for widom method
+        if not isinstance(method, espressomd.reaction_ensemble.WidomInsertion):
             self.assertEqual(
-                list(
-                    method.exclusion_radius_per_type.keys()),
-                [1])
+                list(method.exclusion_radius_per_type.keys()), [1])
             self.assertAlmostEqual(
                 method.exclusion_radius_per_type[1],
                 exclusion_radius_per_type[1],
                 delta=1e-10)
             method.exclusion_radius_per_type = {2: 0.2}
             self.assertEqual(
-                list(
-                    method.exclusion_radius_per_type.keys()),
-                [2])
+                list(method.exclusion_radius_per_type.keys()), [2])
             self.assertAlmostEqual(
-                method.exclusion_radius_per_type[2],
-                0.2,
-                delta=1e-10)
+                method.exclusion_radius_per_type[2], 0.2, delta=1e-10)
         self.assertAlmostEqual(
-            method.get_volume(),
-            self.system.volume(),
-            delta=1e-10)
+            method.get_volume(), self.system.volume(), delta=1e-10)
         method.set_volume(volume=1.)
         self.assertAlmostEqual(method.get_volume(), 1., delta=1e-10)
         self.assertEqual(method.get_non_interacting_type(), 100)
@@ -165,36 +157,23 @@ def check_reaction_parameters(reactions, parameters):
         self.assertEqual(len(method.reactions), 0)
 
     def test_interface(self):
+        params = {'exclusion_range': 0.8,
+                  'exclusion_radius_per_type': {1: 0.1}}
+
         # reaction ensemble
         method = espressomd.reaction_ensemble.ReactionEnsemble(
-            kT=1.5, exclusion_range=0.8, seed=12, exclusion_radius_per_type={1: 0.1})
-        self.check_interface(
-            method,
-            kT=1.5,
-            exclusion_range=0.8,
-            gamma=1.2,
-            exclusion_radius_per_type={
-                1: 0.1})
+            kT=1.4, seed=12, **params)
+        self.check_interface(method, kT=1.4, gamma=1.2, **params)
 
         # constant pH ensemble
         method = espressomd.reaction_ensemble.ConstantpHEnsemble(
-            kT=1.5, exclusion_range=0.8, seed=12, constant_pH=10, exclusion_radius_per_type={1: 0.1})
-        self.check_interface(
-            method,
-            kT=1.5,
-            exclusion_range=0.8,
-            gamma=1.2,
-            exclusion_radius_per_type={
-                1: 0.1})
+            kT=1.5, seed=14, constant_pH=10, **params)
+        self.check_interface(method, kT=1.5, gamma=1.2, **params)
 
         # Widom insertion
-        method = espressomd.reaction_ensemble.WidomInsertion(kT=1.6, seed=12)
-        self.check_interface(
-            method,
-            kT=1.6,
-            exclusion_range=0.,
-            gamma=1.,
-            exclusion_radius_per_type={})
+        method = espressomd.reaction_ensemble.WidomInsertion(kT=1.6, seed=16)
+        self.check_interface(method, kT=1.6, gamma=1., exclusion_range=0.,
+                             exclusion_radius_per_type={})
 
     def test_exceptions(self):
         single_reaction_params = {
diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py
index cefc2d2b86c..d48fb5f16eb 100644
--- a/testsuite/python/virtual_sites_relative.py
+++ b/testsuite/python/virtual_sites_relative.py
@@ -24,14 +24,21 @@
 import tests_common
 
 
-@utx.skipIfMissingFeatures("VIRTUAL_SITES_RELATIVE")
+@utx.skipIfMissingFeatures(["VIRTUAL_SITES_RELATIVE", "LENNARD_JONES"])
 class VirtualSites(ut.TestCase):
     system = espressomd.System(box_l=[1.0, 1.0, 1.0])
 
     np.random.seed(42)
 
+    def setUp(self):
+        self.system.box_l = [10.0, 10.0, 10.0]
+
     def tearDown(self):
         self.system.part.clear()
+        self.system.thermostat.turn_off()
+        self.system.integrator.set_vv()
+        self.system.non_bonded_inter[0, 0].lennard_jones.set_params(
+            epsilon=0., sigma=0., cutoff=0., shift=0.)
 
     def multiply_quaternions(self, a, b):
         return np.array(
@@ -85,8 +92,10 @@ def test_aa_method_switching(self):
             espressomd.virtual_sites.VirtualSitesRelative)
 
     def test_vs_quat(self):
+        self.system.time_step = 0.01
+        self.system.min_global_cut = 0.23
         # First check that quaternion of virtual particle is unchanged if
-        # have_quaterion is false.
+        # have_quaternion is false.
         self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(
             have_quaternion=False)
         self.assertFalse(self.system.virtual_sites.have_quaternion)
@@ -129,20 +138,12 @@ def test_vs_quat(self):
         # Check exceptions.
         with self.assertRaisesRegex(ValueError, "Argument of vs_auto_relate_to has to be of type ParticleHandle or int"):
             p2.vs_auto_relate_to('0')
-        try:
-            p2.vs_auto_relate_to(p1)
-        except BaseException:
-            self.fail('Failed to set a vs from a particle handle')
 
     def test_pos_vel_forces(self):
         system = self.system
         system.cell_system.skin = 0.3
         system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
-        system.box_l = [10, 10, 10]
         system.time_step = 0.004
-        system.thermostat.turn_off()
-        system.non_bonded_inter[0, 0].lennard_jones.set_params(
-            epsilon=0, sigma=0, cutoff=0, shift=0)
 
         # Check setting of min_global_cut
         system.min_global_cut = 0.23
@@ -215,7 +216,7 @@ def run_test_lj(self):
         system = self.system
         system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
         # Parameters
-        n = 90
+        n = 40
         phi = 0.6
         sigma = 1.
         eps = .025
@@ -296,10 +297,6 @@ def run_test_lj(self):
         self.assertNotAlmostEqual(enegry_pre_change, enegry_post_change)
         self.assertNotAlmostEqual(pressure_pre_change, pressure_post_change)
 
-        # Turn off lj interaction
-        system.non_bonded_inter[0, 0].lennard_jones.set_params(
-            epsilon=0, sigma=0, cutoff=0, shift=0)
-
     def test_lj(self):
         """Run LJ fluid test for different cell systems."""
         system = self.system

From 0e09fbf733d737fed45dc31d0d7dcc2bebd26f26 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Thu, 17 Mar 2022 17:46:51 +0100
Subject: [PATCH 02/10] testsuite: Simplify python tests

Use a temporary directory in I/O tests to automate file cleanup.
Improve testing of importlib_wrapper and reduce code duplication.
---
 testsuite/python/h5md.py                    | 64 +++++++++------------
 testsuite/scripts/test_importlib_wrapper.py | 15 +++++
 testsuite/scripts/tutorials/test_convert.py | 49 ++++++----------
 3 files changed, 60 insertions(+), 68 deletions(-)

diff --git a/testsuite/python/h5md.py b/testsuite/python/h5md.py
index de8a7427e89..57ede3d2ff2 100644
--- a/testsuite/python/h5md.py
+++ b/testsuite/python/h5md.py
@@ -28,6 +28,7 @@
 import espressomd
 import espressomd.interactions
 import espressomd.io.writer
+import tempfile
 try:
     import h5py  # h5py has to be imported *after* espressomd (MPI)
     skipIfMissingPythonPackage = utx.no_skip
@@ -46,18 +47,15 @@ class H5mdTests(ut.TestCase):
     Test the core implementation of writing hdf5 files.
 
     """
-    system = espressomd.System(box_l=[1.0, 1.0, 1.0])
-    # avoid particles to be set outside of the main box, otherwise particle
-    # positions are folded in the core when writing out and we cannot directly
-    # compare positions in the dataset and where particles were set. One would
-    # need to unfold the positions of the hdf5 file.
-    box_l = N_PART / 2.0
-    system.box_l = [box_l, box_l, box_l]
+    box_l = N_PART // 2
+    box_l = [box_l, box_l + 1, box_l + 2]
+    system = espressomd.System(box_l=box_l)
     system.cell_system.skin = 0.4
     system.time_step = 0.01
 
+    # set particles outside the main box to verify position folding
     for i in range(N_PART):
-        p = system.part.add(id=i, pos=np.array(3 * [i], dtype=float),
+        p = system.part.add(id=i, pos=np.array(3 * [i - 4], dtype=float),
                             v=np.array([1.0, 2.0, 3.0]), type=23)
         if espressomd.has_features(['MASS']):
             p.mass = 2.3
@@ -77,17 +75,17 @@ class H5mdTests(ut.TestCase):
 
     @classmethod
     def setUpClass(cls):
-        if os.path.isfile('test.h5'):
-            os.remove('test.h5')
+        cls.temp_dir = tempfile.TemporaryDirectory()
+        cls.temp_file = os.path.join(cls.temp_dir.name, 'test.h5')
         h5_units = espressomd.io.writer.h5md.UnitSystem(
             time='ps', mass='u', length='m', charge='e')
         h5 = espressomd.io.writer.h5md.H5md(
-            file_path="test.h5", unit_system=h5_units)
+            file_path=cls.temp_file, unit_system=h5_units)
         h5.write()
         h5.write()
         h5.flush()
         h5.close()
-        cls.py_file = h5py.File("test.h5", 'r')
+        cls.py_file = h5py.File(cls.temp_file, 'r')
         cls.py_pos = cls.py_file['particles/atoms/position/value'][1]
         cls.py_img = cls.py_file['particles/atoms/image/value'][1]
         cls.py_mass = cls.py_file['particles/atoms/mass/value'][1]
@@ -102,10 +100,10 @@ def setUpClass(cls):
 
     @classmethod
     def tearDownClass(cls):
-        os.remove("test.h5")
+        cls.temp_dir.cleanup()
 
     def test_opening(self):
-        h5 = espressomd.io.writer.h5md.H5md(file_path="test.h5")
+        h5 = espressomd.io.writer.h5md.H5md(file_path=self.temp_file)
         h5.close()
 
     def test_box(self):
@@ -126,9 +124,10 @@ def test_metadata(self):
 
     def test_pos(self):
         """Test if positions have been written properly."""
-        np.testing.assert_allclose(
-            np.array([3 * [float(i) % self.box_l] for i in range(N_PART)]),
-            np.array([x for (_, x) in sorted(zip(self.py_id, self.py_pos))]))
+        pos_ref = np.outer(np.arange(N_PART) - 4, [1, 1, 1])
+        pos_ref = np.mod(pos_ref, self.box_l)
+        pos_read = [x for (_, x) in sorted(zip(self.py_id, self.py_pos))]
+        np.testing.assert_allclose(pos_read, pos_ref)
 
     def test_time(self):
         """Test for time dataset."""
@@ -136,11 +135,10 @@ def test_time(self):
 
     def test_img(self):
         """Test if images have been written properly."""
-        images = np.append(np.zeros((int(N_PART / 2), 3)),
-                           np.ones((int(N_PART / 2), 3)))
-        images = images.reshape(N_PART, 3)
-        np.testing.assert_allclose(
-            [x for (_, x) in sorted(zip(self.py_id, self.py_img))], images)
+        pos_ref = np.outer(np.arange(N_PART) - 4, [1, 1, 1])
+        images_ref = np.floor_divide(pos_ref, self.box_l)
+        images_read = [x for (_, x) in sorted(zip(self.py_id, self.py_img))]
+        np.testing.assert_allclose(images_read, images_ref)
 
     @utx.skipIfMissingFeatures("MASS")
     def test_mass(self):
@@ -185,22 +183,16 @@ def test_script(self):
         self.assertEqual(data, ref)
 
     def test_units(self):
-        self.assertEqual(
-            self.py_file['particles/atoms/id/time'].attrs['unit'], b'ps')
-        self.assertEqual(
-            self.py_file['particles/atoms/position/value'].attrs['unit'], b'm')
+        def get_unit(path):
+            return self.py_file[path].attrs['unit']
+        self.assertEqual(get_unit('particles/atoms/id/time'), b'ps')
+        self.assertEqual(get_unit('particles/atoms/position/value'), b'm')
         if espressomd.has_features(['ELECTROSTATICS']):
-            self.assertEqual(
-                self.py_file['particles/atoms/charge/value'].attrs['unit'], b'e')
+            self.assertEqual(get_unit('particles/atoms/charge/value'), b'e')
         if espressomd.has_features(['MASS']):
-            self.assertEqual(
-                self.py_file['particles/atoms/mass/value'].attrs['unit'], b'u')
-        self.assertEqual(
-            self.py_file['particles/atoms/force/value'].attrs['unit'],
-            b'm u ps-2')
-        self.assertEqual(
-            self.py_file['particles/atoms/velocity/value'].attrs['unit'],
-            b'm ps-1')
+            self.assertEqual(get_unit('particles/atoms/mass/value'), b'u')
+        self.assertEqual(get_unit('particles/atoms/force/value'), b'm u ps-2')
+        self.assertEqual(get_unit('particles/atoms/velocity/value'), b'm ps-1')
 
     def test_links(self):
         time_ref = self.py_id_time
diff --git a/testsuite/scripts/test_importlib_wrapper.py b/testsuite/scripts/test_importlib_wrapper.py
index 9d959edae26..e21904ecb6f 100644
--- a/testsuite/scripts/test_importlib_wrapper.py
+++ b/testsuite/scripts/test_importlib_wrapper.py
@@ -380,6 +380,21 @@ def test_delimit_statements(self):
         linenos_out = iw.delimit_statements(source_code)
         self.assertEqual(linenos_out, linenos_exp)
 
+    def test_ipython_magics(self):
+        lines = [
+            '%matplotlib inline',
+            '%%matplotlib notebook',
+            'import matplotlib.pyplot as plt',
+            '__doc__="%matplotlib inline"',
+        ]
+        code = '\n'.join(lines)
+        code_protected = iw.protect_ipython_magics(code)
+        code_protected_ref = f'\n{code}'.replace(
+            '\n%', '\n#_IPYTHON_MAGIC_%')[1:]
+        self.assertEqual(code_protected, code_protected_ref)
+        code_deprotected = iw.deprotect_ipython_magics(code_protected)
+        self.assertEqual(code_deprotected, code)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/scripts/tutorials/test_convert.py b/testsuite/scripts/tutorials/test_convert.py
index 1386ab93826..0f24519cdc2 100644
--- a/testsuite/scripts/tutorials/test_convert.py
+++ b/testsuite/scripts/tutorials/test_convert.py
@@ -85,11 +85,18 @@ class HtmlRunner(ut.TestCase):
             "version": ".".join(map(str, sys.version_info[:3]))}
     }
 
-    def failed_to_run(self, cmd):
-        traceback.print_exc()
-        self.fail('Could not run @CMAKE_BINARY_DIR@/pypresso '
-                  '@CMAKE_BINARY_DIR@/doc/tutorials/convert.py ' +
-                  ' '.join(cmd))
+    def run_command(self, cmd, output=None):
+        error_msg = ("Could not run @CMAKE_BINARY_DIR@/pypresso "
+                     "@CMAKE_BINARY_DIR@/doc/tutorials/convert.py "
+                     f"{' '.join(cmd)}")
+        try:
+            args = convert.parser.parse_args(cmd)
+            args.callback(args)
+        except BaseException:
+            traceback.print_exc()
+            self.fail(error_msg)
+        if output is not None:
+            self.assertTrue(os.path.isfile(output), f"{output} not created")
 
     def test_html_wrapper(self):
         f_input = '@CMAKE_CURRENT_BINARY_DIR@/test_convert_notebook.ipynb'
@@ -114,12 +121,7 @@ def test_html_wrapper(self):
                '--scripts', f_script,
                '--substitutions', 'global_var=20',
                '--execute']
-        try:
-            args = convert.parser.parse_args(cmd)
-            args.callback(args)
-        except BaseException:
-            self.failed_to_run(cmd)
-        self.assertTrue(os.path.isfile(f_output), f_output + ' not created')
+        self.run_command(cmd, f_output)
         # read processed notebook
         with open(f_output, encoding='utf-8') as f:
             nb_output = nbformat.read(f, as_version=4)
@@ -193,12 +195,7 @@ def test_exercise2_plugin(self):
                '--output', f_output,
                '--substitutions', 'global_var=20',
                '--exercise2', '--remove-empty-cells']
-        try:
-            args = convert.parser.parse_args(cmd)
-            args.callback(args)
-        except BaseException:
-            self.failed_to_run(cmd)
-        self.assertTrue(os.path.isfile(f_output), f_output + ' not created')
+        self.run_command(cmd, f_output)
         # read processed notebook
         with open(f_output, encoding='utf-8') as f:
             nb_output = nbformat.read(f, as_version=4)
@@ -246,11 +243,7 @@ def test_exercise2_conversion(self):
             nbformat.write(nb, f)
         # run command and check for errors
         cmd = ['exercise2', '--to-py', f_input]
-        try:
-            args = convert.parser.parse_args(cmd)
-            args.callback(args)
-        except BaseException:
-            self.failed_to_run(cmd)
+        self.run_command(cmd, f_input)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f:
             nb_output = nbformat.read(f, as_version=4)
@@ -267,11 +260,7 @@ def test_exercise2_conversion(self):
         self.assertEqual(next(cells, 'EOF'), 'EOF')
         # run command and check for errors
         cmd = ['exercise2', '--to-md', f_input]
-        try:
-            args = convert.parser.parse_args(cmd)
-            args.callback(args)
-        except BaseException:
-            self.failed_to_run(cmd)
+        self.run_command(cmd, f_input)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f:
             nb_output = nbformat.read(f, as_version=4)
@@ -305,11 +294,7 @@ def test_exercise2_autopep8(self):
             nbformat.write(nb, f)
         # run command and check for errors
         cmd = ['exercise2', '--pep8', f_input]
-        try:
-            args = convert.parser.parse_args(cmd)
-            args.callback(args)
-        except BaseException:
-            self.failed_to_run(cmd)
+        self.run_command(cmd)
         # read processed notebook
         with open(f_input, encoding='utf-8') as f:
             nb_output = nbformat.read(f, as_version=4)

From aa3fc5d04cec7c2d07cbfbbd35421b36f1f82931 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Thu, 17 Mar 2022 18:23:04 +0100
Subject: [PATCH 03/10] core: Unit test the rotation code

---
 src/core/rotation.cpp                 |   3 +-
 src/core/unit_tests/CMakeLists.txt    |   2 +
 src/core/unit_tests/Particle_test.cpp |   2 +
 src/core/unit_tests/rotation_test.cpp | 219 ++++++++++++++++++++++++++
 4 files changed, 225 insertions(+), 1 deletion(-)
 create mode 100644 src/core/unit_tests/rotation_test.cpp

diff --git a/src/core/rotation.cpp b/src/core/rotation.cpp
index d402130f5a1..f9adc15c32c 100644
--- a/src/core/rotation.cpp
+++ b/src/core/rotation.cpp
@@ -117,7 +117,8 @@ static void define_Qdd(Particle const &p, Utils::Quaternion<double> &Qd,
  *
  *  For very high angular velocities (e.g. if the product of @p time_step
  *  with the largest component of @ref ParticleMomentum::omega "p.omega()"
- *  is superior to ~2.0), the calculation might fail.
+ *  is superior to ~2.0) and for @p time_step superior or equal to unity,
+ *  the calculation might fail.
  *
  *  \todo implement for fixed_coord_flag
  */
diff --git a/src/core/unit_tests/CMakeLists.txt b/src/core/unit_tests/CMakeLists.txt
index 95560b0931a..fdc28c5c5c7 100644
--- a/src/core/unit_tests/CMakeLists.txt
+++ b/src/core/unit_tests/CMakeLists.txt
@@ -37,6 +37,8 @@ unit_test(NAME p3m_test SRC p3m_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME link_cell_test SRC link_cell_test.cpp DEPENDS EspressoUtils)
 unit_test(NAME Particle_test SRC Particle_test.cpp DEPENDS EspressoUtils
           Boost::serialization)
+unit_test(NAME rotation_test SRC rotation_test.cpp DEPENDS EspressoUtils
+          EspressoCore)
 unit_test(NAME field_coupling_couplings SRC field_coupling_couplings_test.cpp
           DEPENDS EspressoUtils)
 unit_test(NAME field_coupling_fields SRC field_coupling_fields_test.cpp DEPENDS
diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp
index 2726a7b7621..86a0af10410 100644
--- a/src/core/unit_tests/Particle_test.cpp
+++ b/src/core/unit_tests/Particle_test.cpp
@@ -280,5 +280,7 @@ BOOST_AUTO_TEST_CASE(particle_bitfields) {
   BOOST_CHECK(p.can_rotate_around(0));
   BOOST_CHECK(p.can_rotate_around(1));
   BOOST_CHECK(p.can_rotate_around(2));
+  p.set_cannot_rotate_all_axes();
+  BOOST_CHECK(not p.can_rotate());
 #endif
 }
diff --git a/src/core/unit_tests/rotation_test.cpp b/src/core/unit_tests/rotation_test.cpp
new file mode 100644
index 00000000000..5ff5d98cb6c
--- /dev/null
+++ b/src/core/unit_tests/rotation_test.cpp
@@ -0,0 +1,219 @@
+/*
+ * Copyright (C) 2022 The ESPResSo project
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#define BOOST_TEST_MODULE rotation test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include "Particle.hpp"
+#include "config.hpp"
+#include "rotation.hpp"
+
+#include <utils/Vector.hpp>
+#include <utils/constants.hpp>
+#include <utils/quaternion.hpp>
+
+#include <limits>
+#include <stdexcept>
+#include <tuple>
+
+#ifdef ROTATION
+
+auto constexpr tol = 5. * 100. * std::numeric_limits<double>::epsilon();
+
+namespace Testing {
+std::tuple<Utils::Quaternion<double>, Utils::Vector3d>
+setup_trivial_quat(int i, Utils::Vector3d const &v_in) {
+  auto quat = Utils::Quaternion<double>{{0., 0., 0., 0.}};
+  quat[i] = 1.;
+  auto v_ref = v_in;
+  if (i) {
+    v_ref *= -1.;
+    v_ref[i - 1] *= -1.;
+  }
+  return std::make_tuple(quat, v_ref);
+}
+} // namespace Testing
+
+BOOST_AUTO_TEST_CASE(convert_vector_space_to_body_test) {
+  auto const t_in = Utils::Vector3d{{1., 2., 3.}};
+  for (int i : {0, 1, 2, 3}) {
+    auto p = Particle();
+    Utils::Vector3d t_ref;
+    std::tie(p.quat(), t_ref) = Testing::setup_trivial_quat(i, t_in);
+    auto const t_out = convert_vector_space_to_body(p, t_in);
+    for (int j : {0, 1, 2}) {
+      BOOST_CHECK_CLOSE(t_out[j], t_ref[j], tol);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(convert_torque_to_body_frame_apply_fix_test) {
+  auto const t_in = Utils::Vector3d{{1., 2., 3.}};
+  {
+    // test particle torque conversion
+    for (int i : {0, 1, 2, 3}) {
+      auto p = Particle();
+      p.set_can_rotate_all_axes();
+      Utils::Vector3d t_ref;
+      std::tie(p.quat(), t_ref) = Testing::setup_trivial_quat(i, t_in);
+      p.torque() = t_in;
+      convert_torque_to_body_frame_apply_fix(p);
+      auto const t_out = p.torque();
+      for (int j : {0, 1, 2}) {
+        BOOST_CHECK_CLOSE(t_out[j], t_ref[j], tol);
+      }
+    }
+  }
+  {
+    // torque is set to zero for axes without rotation
+    for (int j : {0, 1, 2}) {
+      auto p = Particle();
+      p.set_can_rotate_all_axes();
+      p.set_can_rotate_around(j, false);
+      p.quat() = Utils::Quaternion<double>::identity();
+      p.torque() = t_in;
+      auto t_ref = t_in;
+      t_ref[j] = 0.;
+      convert_torque_to_body_frame_apply_fix(p);
+      BOOST_TEST(p.torque() == t_ref, boost::test_tools::per_element());
+    }
+  }
+  {
+    // torque is always zero for non-rotatable particles
+    auto p = Particle();
+    p.set_cannot_rotate_all_axes();
+    p.quat() = Utils::Quaternion<double>::identity();
+    p.torque() = t_in;
+    convert_torque_to_body_frame_apply_fix(p);
+    auto const t_ref = Utils::Vector3d{};
+    BOOST_TEST(p.torque() == t_ref, boost::test_tools::per_element());
+  }
+}
+
+BOOST_AUTO_TEST_CASE(rotate_particle_body_test) {
+  auto p = Particle();
+  p.quat() = {1., 2., 3., 4.};
+  {
+    // fixed particles are unaffected, quaternion is identical to original
+    p.set_cannot_rotate_all_axes();
+    auto const phi = 2.;
+    auto const quat = local_rotate_particle_body(p, {0., 0., 1.}, phi);
+    BOOST_TEST((quat == p.quat()));
+  }
+  {
+    // edge case: null rotation throws an exception
+    p.set_can_rotate_around(2, true);
+    auto const phi = 2.;
+    BOOST_CHECK_THROW(local_rotate_particle_body(p, {1., 1., 0.}, phi),
+                      std::exception);
+  }
+  {
+    // an angle of zero has no effect, quaternion is identical to original
+    p.set_can_rotate_all_axes();
+    auto const phi = 0.;
+    auto const quat = local_rotate_particle_body(p, {1., 2., 3.}, phi);
+    BOOST_TEST((quat == p.quat()));
+  }
+  {
+    // an angle of pi around the z-axis flips the quaternion sequence
+    p.set_can_rotate_all_axes();
+    auto const phi = Utils::pi<double>();
+    auto const quat = local_rotate_particle_body(p, {0., 0., 1.}, phi);
+    auto const quat_ref = Utils::Vector4d{{-4., 3., -2., 1.}};
+    for (int i : {0, 1, 2, 3}) {
+      BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+    }
+  }
+  {
+    // an angle of 2 pi around the z-axis flips the quaternion sign
+    p.set_can_rotate_all_axes();
+    auto const phi = 2. * Utils::pi<double>();
+    auto const quat = local_rotate_particle_body(p, {0., 0., 1.}, phi);
+    auto const quat_ref = Utils::Vector4d{{-1., -2., -3., -4.}};
+    for (int i : {0, 1, 2, 3}) {
+      BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+    }
+  }
+}
+
+BOOST_AUTO_TEST_CASE(propagate_omega_quat_particle_test) {
+  auto p = Particle();
+  p.set_can_rotate_all_axes();
+  {
+    // test edge case: null quaternion and no rotation
+    p.quat() = {0., 0., 0., 0.};
+    p.omega() = {0., 0., 0.};
+    propagate_omega_quat_particle(p, 0.01);
+    auto const quat = p.quat();
+    auto const quat_ref = Utils::Quaternion<double>::identity();
+    for (int i : {0, 1, 2, 3}) {
+      BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+    }
+  }
+  {
+    // test trivial cases with parameters extremely close to the limit:
+    // time step almost 1.0 and product of time step with omega almost 2.0
+    auto const time_step = 0.99;
+    for (int j : {0, 1, 2}) {
+      p.quat() = Utils::Quaternion<double>::identity();
+      p.omega() = {0., 0., 0.};
+      p.omega()[j] = 2.;
+      propagate_omega_quat_particle(p, time_step);
+      auto const quat = p.quat();
+      auto quat_ref = Utils::Quaternion<double>::identity();
+      quat_ref[1 + j] = time_step;
+      quat_ref[0] = std::sqrt(1. - time_step * time_step);
+      for (int i : {0, 1, 2, 3}) {
+        BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+      }
+    }
+  }
+}
+
+#ifdef DIPOLES
+BOOST_AUTO_TEST_CASE(convert_dip_to_quat_test) {
+  auto const quat_to_vector4d = [](Utils::Quaternion<double> const &quat) {
+    return Utils::Vector4d{quat.data(), quat.data() + 4};
+  };
+  auto p = Particle();
+  p.quat() = {1., 2., 3., 4.};
+  {
+    auto const dipm = 0.8;
+    auto const pair = convert_dip_to_quat({0., 0., dipm});
+    auto const quat = quat_to_vector4d(pair.first);
+    auto const quat_ref = Utils::Vector4d{{1., 0., 0., 0.}};
+    for (int i : {0, 1, 2, 3}) {
+      BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+    }
+    BOOST_CHECK_CLOSE(pair.second, dipm, tol);
+  }
+  {
+    auto const dipm = 1.6;
+    auto const pair = convert_dip_to_quat({dipm, 0., 0.});
+    auto const quat = quat_to_vector4d(pair.first);
+    auto const quat_ref = Utils::Vector4d{{0.5, -0.5, 0.5, -0.5}};
+    for (int i : {0, 1, 2, 3}) {
+      BOOST_CHECK_CLOSE(quat[i], quat_ref[i], tol);
+    }
+    BOOST_CHECK_CLOSE(pair.second, dipm, tol);
+  }
+}
+#endif // DIPOLES
+#endif // ROTATION

From 11939cce7a2e31d30f9a57857d7124d06178c317 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Thu, 17 Mar 2022 18:55:59 +0100
Subject: [PATCH 04/10] core: MPI-safe virtual sites relative exceptions

When a real particle is not found, queue a runtime error instead
of throwing a fatal error.
---
 .../virtual_sites/VirtualSitesRelative.cpp    | 84 +++++++++----------
 testsuite/python/virtual_sites_relative.py    | 14 +++-
 2 files changed, 51 insertions(+), 47 deletions(-)

diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp
index c88fb2a33fd..b80996bbef9 100644
--- a/src/core/virtual_sites/VirtualSitesRelative.cpp
+++ b/src/core/virtual_sites/VirtualSitesRelative.cpp
@@ -23,6 +23,7 @@
 
 #include "Particle.hpp"
 #include "cells.hpp"
+#include "errorhandling.hpp"
 #include "forces.hpp"
 #include "grid.hpp"
 #include "integrate.hpp"
@@ -33,80 +34,74 @@
 #include <utils/math/tensor_product.hpp>
 #include <utils/quaternion.hpp>
 
-#include <stdexcept>
-
-namespace {
 /**
  * @brief Vector pointing from the real particle to the virtual site.
  *
  * @return Relative distance.
  */
-auto connection_vector(
-    Particle const &p_ref,
-    ParticleProperties::VirtualSitesRelativeParameters const &vs_rel) {
+static auto connection_vector(Particle const &p_ref, Particle const &p_vs) {
   // Calculate the quaternion defining the orientation of the vector connecting
   // the virtual site and the real particle
   // This is obtained, by multiplying the quaternion representing the director
   // of the real particle with the quaternion of the virtual particle, which
   // specifies the relative orientation.
   auto const director = Utils::convert_quaternion_to_director(
-                            p_ref.quat() * vs_rel.rel_orientation)
+                            p_ref.quat() * p_vs.vs_relative().rel_orientation)
                             .normalize();
 
-  return vs_rel.distance * director;
+  return p_vs.vs_relative().distance * director;
 }
 
 /**
  * @brief Velocity of the virtual site
  * @param p_ref Reference particle for the virtual site.
- * @param vs_rel Parameters for the virtual site.
+ * @param p_vs Virtual site.
  * @return Velocity of the virtual site.
  */
-Utils::Vector3d
-velocity(Particle const &p_ref,
-         ParticleProperties::VirtualSitesRelativeParameters const &vs_rel) {
-  auto const d = connection_vector(p_ref, vs_rel);
+static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs) {
+  auto const d = connection_vector(p_ref, p_vs);
 
   // Get omega of real particle in space-fixed frame
   auto const omega_space_frame =
       convert_vector_body_to_space(p_ref, p_ref.omega());
-  // Obtain velocity from v=v_real particle + omega_real_particle \times
-  // director
+  // Obtain velocity from v = v_real particle + omega_real_particle * director
   return vector_product(omega_space_frame, d) + p_ref.v();
 }
 
 /**
- * @brief Get reference particle.
+ * @brief Get real particle tracked by a virtual site.
  *
- * @param vs_rel Parameters to get the reference particle for.
- * @return Pointer to reference particle.
+ * @param p Virtual site.
+ * @return Pointer to real particle if @p p is a virtual site,
+ *         otherwise nullptr. If the real particle is not found,
+ *         queue a runtime error.
  */
-Particle &get_reference_particle(
-    ParticleProperties::VirtualSitesRelativeParameters const &vs_rel) {
+static Particle *get_reference_particle(Particle const &p) {
+  if (!p.is_virtual()) {
+    return nullptr;
+  }
+  auto const &vs_rel = p.vs_relative();
   auto p_ref_ptr = cell_structure.get_local_particle(vs_rel.to_particle_id);
   if (!p_ref_ptr) {
-    throw std::runtime_error("No real particle associated with virtual site.");
+    runtimeErrorMsg() << "No real particle with id " << vs_rel.to_particle_id
+                      << " for virtual site with id " << p.identity();
   }
-  return *p_ref_ptr;
+  return p_ref_ptr;
 }
 
 /**
  * @brief Constraint force to hold the particles at its prescribed position.
  *
- * @param f Force on the virtual site.
  * @param p_ref Reference particle.
- * @param vs_rel Parameters.
+ * @param p_vs Virtual site.
  * @return Constraint force.
  */
-auto constraint_stress(
-    const Utils::Vector3d &f, const Particle &p_ref,
-    const ParticleProperties::VirtualSitesRelativeParameters &vs_rel) {
+static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) {
   /* The constraint force is minus the force on the particle, make it force
    * free. The counter force is translated by the connection vector to the
    * real particle, hence the virial stress is */
-  return tensor_product(-f, connection_vector(p_ref, vs_rel));
+  return tensor_product(-p_vs.force(), connection_vector(p_ref, p_vs));
 }
-} // namespace
 
 void VirtualSitesRelative::update() const {
   cell_structure.ghosts_update(Cells::DATA_PART_POSITION |
@@ -114,12 +109,12 @@ void VirtualSitesRelative::update() const {
 
   auto const particles = cell_structure.local_particles();
   for (auto &p : particles) {
-    if (!p.is_virtual())
+    auto const *p_ref_ptr = get_reference_particle(p);
+    if (!p_ref_ptr)
       continue;
 
-    auto const &p_ref = get_reference_particle(p.vs_relative());
-
-    auto new_pos = p_ref.pos() + connection_vector(p_ref, p.vs_relative());
+    auto const &p_ref = *p_ref_ptr;
+    auto new_pos = p_ref.pos() + connection_vector(p_ref, p);
     /* The shift has to respect periodic boundaries: if the reference
      * particles is not in the same image box, we potentially avoid shifting
      * to the other side of the box. */
@@ -129,7 +124,7 @@ void VirtualSitesRelative::update() const {
     fold_position(shift, image_shift, box_geo);
     p.image_box() = p_ref.image_box() - image_shift;
 
-    p.v() = velocity(p_ref, p.vs_relative());
+    p.v() = velocity(p_ref, p);
 
     if (box_geo.type() == BoxType::LEES_EDWARDS) {
       auto const &shear_dir = box_geo.clees_edwards_bc().shear_direction;
@@ -158,16 +153,15 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const {
 
   // Iterate over all the particles in the local cells
   for (auto &p : cell_structure.local_particles()) {
-    // We only care about virtual particles
-    if (!p.is_virtual())
+    auto *p_ref_ptr = get_reference_particle(p);
+    if (!p_ref_ptr)
       continue;
-    auto &p_ref = get_reference_particle(p.vs_relative());
 
     // Add forces and torques
+    auto &p_ref = *p_ref_ptr;
     p_ref.force() += p.force();
     p_ref.torque() +=
-        vector_product(connection_vector(p_ref, p.vs_relative()), p.force()) +
-        p.torque();
+        vector_product(connection_vector(p_ref, p), p.force()) + p.torque();
   }
 }
 
@@ -175,13 +169,11 @@ void VirtualSitesRelative::back_transfer_forces_and_torques() const {
 Utils::Matrix<double, 3, 3> VirtualSitesRelative::pressure_tensor() const {
   Utils::Matrix<double, 3, 3> pressure_tensor = {};
 
-  for (auto &p : cell_structure.local_particles()) {
-    if (!p.is_virtual())
-      continue;
-
-    auto const &p_ref = get_reference_particle(p.vs_relative());
-
-    pressure_tensor += constraint_stress(p.force(), p_ref, p.vs_relative());
+  for (auto const &p : cell_structure.local_particles()) {
+    auto const *p_ref_ptr = get_reference_particle(p);
+    if (p_ref_ptr) {
+      pressure_tensor += constraint_stress(*p_ref_ptr, p);
+    }
   }
 
   return pressure_tensor;
diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py
index d48fb5f16eb..14657efbac1 100644
--- a/testsuite/python/virtual_sites_relative.py
+++ b/testsuite/python/virtual_sites_relative.py
@@ -39,6 +39,7 @@ def tearDown(self):
         self.system.integrator.set_vv()
         self.system.non_bonded_inter[0, 0].lennard_jones.set_params(
             epsilon=0., sigma=0., cutoff=0., shift=0.)
+        self.system.virtual_sites = espressomd.virtual_sites.VirtualSitesOff()
 
     def multiply_quaternions(self, a, b):
         return np.array(
@@ -135,9 +136,20 @@ def test_vs_quat(self):
         self.system.integrator.run(1)
         self.assertAlmostEqual(np.dot(p1.director, p2.director), 0.0)
 
-        # Check exceptions.
+    def test_vs_exceptions(self):
+        system = self.system
+        system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
+        system.time_step = 0.01
+        p2 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=[1, 1, 1], id=2)
+        p3 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=[1, 1, 1], id=3)
+        # relating to anything else other than a particle or id is not allowed
         with self.assertRaisesRegex(ValueError, "Argument of vs_auto_relate_to has to be of type ParticleHandle or int"):
             p2.vs_auto_relate_to('0')
+        # relating to a deleted particle is not allowed
+        with self.assertRaisesRegex(Exception, "No real particle with id 3 for virtual site with id 2"):
+            p2.vs_auto_relate_to(p3)
+            p3.remove()
+            system.integrator.run(0, recalc_forces=True)
 
     def test_pos_vel_forces(self):
         system = self.system

From 6e3cf3235be1b591348c211e8a570c8931c8f4e0 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Thu, 17 Mar 2022 22:31:15 +0100
Subject: [PATCH 05/10] Improve exception mechanism

Add missing calls to handle_errors() and remove superfluous ones.
When calling a script interface method, call handle_errors() with a
clear error message. Catch null quaternions in the script interface
to avoid triggering a fatal error in boost::qvm. Prevent relative
virtual sites from tracking themselves.
---
 src/core/integrate.cpp                        | 12 ++---
 src/core/io/writer/h5md_core.hpp              |  5 +-
 .../reaction_methods/ReactionAlgorithm.hpp    | 27 ++++++-----
 .../tests/particle_tracking_test.cpp          |  3 ++
 src/core/virtual_sites.cpp                    | 14 +++---
 .../virtual_sites/VirtualSitesRelative.cpp    |  4 +-
 src/python/espressomd/cellsystem.pyx          |  2 +-
 src/python/espressomd/integrate.pyx           |  3 +-
 src/python/espressomd/interactions.pyx        | 15 +++---
 src/python/espressomd/particle_data.pxd       |  2 +-
 src/python/espressomd/particle_data.pyx       | 15 ++++--
 src/python/espressomd/script_interface.pyx    |  2 +-
 src/python/espressomd/system.pyx              |  4 +-
 .../reaction_methods/ReactionAlgorithm.hpp    |  4 +-
 testsuite/python/CMakeLists.txt               |  1 +
 testsuite/python/field_test.py                | 17 +++++++
 testsuite/python/h5md.py                      | 31 +++++++++++++
 testsuite/python/integrator_exceptions.py     | 45 +++++++++++++++++-
 .../python/interactions_bonded_interface.py   | 46 +++++++++++++++++++
 .../interactions_non-bonded_interface.py      | 13 ++++++
 testsuite/python/particle.py                  | 10 ++++
 testsuite/python/reaction_methods.py          | 18 ++++++--
 testsuite/python/virtual_sites_relative.py    |  9 ++++
 23 files changed, 247 insertions(+), 55 deletions(-)

diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp
index d3777995dec..b56f36c8b3a 100644
--- a/src/core/integrate.cpp
+++ b/src/core/integrate.cpp
@@ -61,6 +61,7 @@
 #include <boost/range/algorithm/min_element.hpp>
 
 #include <algorithm>
+#include <cassert>
 #include <cmath>
 #include <csignal>
 #include <functional>
@@ -419,6 +420,9 @@ int integrate(int n_steps, int reuse_forces) {
 
 int python_integrate(int n_steps, bool recalc_forces_par,
                      bool reuse_forces_par) {
+
+  assert(n_steps >= 0);
+
   // Override the signal handler so that the integrator obeys Ctrl+C
   SignalHandler sa(SIGINT, [](int) { ctrl_C = 1; });
 
@@ -431,12 +435,6 @@ int python_integrate(int n_steps, bool recalc_forces_par,
     reuse_forces = -1;
   }
 
-  /* go on with integrate <n_steps> */
-  if (n_steps < 0) {
-    runtimeErrorMsg() << "illegal number of steps (must be >0)";
-    return ES_ERROR;
-  }
-
   /* if skin wasn't set, do an educated guess now */
   if (!skin_set) {
     auto const max_cut = maximal_cutoff();
@@ -552,7 +550,7 @@ REGISTER_CALLBACK(mpi_set_time_step_local)
 
 void mpi_set_time_step(double time_s) {
   if (time_s <= 0.)
-    throw std::invalid_argument("time_step must be > 0.");
+    throw std::domain_error("time_step must be > 0.");
   if (lb_lbfluid_get_lattice_switch() != ActiveLB::NONE)
     check_tau_time_step_consistency(lb_lbfluid_get_tau(), time_s);
   mpi_call_all(mpi_set_time_step_local, time_s);
diff --git a/src/core/io/writer/h5md_core.hpp b/src/core/io/writer/h5md_core.hpp
index 30f4e4b15b4..31c08497375 100644
--- a/src/core/io/writer/h5md_core.hpp
+++ b/src/core/io/writer/h5md_core.hpp
@@ -217,9 +217,8 @@ struct incompatible_h5mdfile : public std::exception {
 
 struct left_backupfile : public std::exception {
   const char *what() const noexcept override {
-    return "A backup of the .h5 file exists. This usually means \
-that either you forgot to call the 'close' method or your simulation \
-crashed.";
+    return "A backup of the .h5 file exists. This usually means that either "
+           "you forgot to call the 'close' method or your simulation crashed.";
   }
 };
 
diff --git a/src/core/reaction_methods/ReactionAlgorithm.hpp b/src/core/reaction_methods/ReactionAlgorithm.hpp
index c99a1109ef2..2444fe6112e 100644
--- a/src/core/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/core/reaction_methods/ReactionAlgorithm.hpp
@@ -61,18 +61,7 @@ class ReactionAlgorithm {
     }
     this->kT = kT;
     this->exclusion_range = exclusion_range;
-    for (auto const &item : exclusion_radius_per_type) {
-      auto type = item.first;
-      auto exclusion_radius = item.second;
-      if (exclusion_radius < 0) {
-        std::stringstream ss;
-        ss << "Invalid excluded_radius value for type " << type
-           << " value: " << exclusion_radius;
-        std::string error_message = ss.str();
-        throw std::domain_error(error_message);
-      }
-    }
-    this->exclusion_radius_per_type = exclusion_radius_per_type;
+    set_exclusion_radius_per_type(exclusion_radius_per_type);
     update_volume();
   }
 
@@ -109,6 +98,20 @@ class ReactionAlgorithm {
     volume = new_volume;
   }
   void update_volume();
+  void
+  set_exclusion_radius_per_type(std::unordered_map<int, double> const &map) {
+    for (auto const &item : map) {
+      auto const type = item.first;
+      auto const exclusion_radius = item.second;
+      if (exclusion_radius < 0.) {
+        throw std::domain_error("Invalid excluded_radius value for type " +
+                                std::to_string(type) + ": radius " +
+                                std::to_string(exclusion_radius));
+      }
+    }
+    exclusion_radius_per_type = map;
+  }
+
   virtual int do_reaction(int reaction_steps);
   void check_reaction_method() const;
   void remove_constraint() { m_reaction_constraint = ReactionConstraint::NONE; }
diff --git a/src/core/reaction_methods/tests/particle_tracking_test.cpp b/src/core/reaction_methods/tests/particle_tracking_test.cpp
index df14f1e5f32..b792b9e1761 100644
--- a/src/core/reaction_methods/tests/particle_tracking_test.cpp
+++ b/src/core/reaction_methods/tests/particle_tracking_test.cpp
@@ -57,6 +57,9 @@ BOOST_FIXTURE_TEST_CASE(particle_type_map_test, ParticleFactory) {
   // exception for random index that exceeds the number of particles
   BOOST_CHECK_THROW(get_random_p_id(type, 10), std::runtime_error);
 
+  // exception for untracked particle types
+  BOOST_CHECK_THROW(get_random_p_id(type + 1, 0), std::runtime_error);
+
   // check particle selection
   BOOST_CHECK_EQUAL(get_random_p_id(type, 0), pid);
 }
diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp
index e591277f0a2..629468a9b00 100644
--- a/src/core/virtual_sites.cpp
+++ b/src/core/virtual_sites.cpp
@@ -38,6 +38,7 @@
 
 #include <cmath>
 #include <cstdio>
+#include <stdexcept>
 #include <tuple>
 
 namespace {
@@ -66,12 +67,10 @@ calculate_vs_relate_to_params(Particle const &p_current,
   if (dist > min_global_cut && n_nodes > 1) {
     runtimeErrorMsg()
         << "Warning: The distance between virtual and non-virtual particle ("
-        << dist << ") is\nlarger than the minimum global cutoff ("
-        << min_global_cut
-        << "). This may lead to incorrect simulations\nunder certain "
-           "conditions. Set the \"System()\" "
-           "class property \"min_global_cut\" to\nincrease the minimum "
-           "cutoff.\n";
+        << dist << ") is larger than the minimum global cutoff ("
+        << min_global_cut << "). This may lead to incorrect simulations under "
+        << "certain conditions. Adjust the property system.min_global_cut to "
+        << "increase the minimum cutoff.";
   }
 
   // Now, calculate the quaternions which specify the angle between
@@ -136,6 +135,9 @@ void local_vs_relate_to(Particle &p_current, Particle const &p_relate_to) {
 }
 
 void vs_relate_to(int part_num, int relate_to) {
+  if (part_num == relate_to) {
+    throw std::invalid_argument("A virtual site cannot relate to itself");
+  }
   // Get the data for the particle we act on and the one we want to relate
   // it to.
   auto const &p_current = get_particle_data(part_num);
diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp
index b80996bbef9..943c1c8f25a 100644
--- a/src/core/virtual_sites/VirtualSitesRelative.cpp
+++ b/src/core/virtual_sites/VirtualSitesRelative.cpp
@@ -72,9 +72,7 @@ static Utils::Vector3d velocity(Particle const &p_ref, Particle const &p_vs) {
  * @brief Get real particle tracked by a virtual site.
  *
  * @param p Virtual site.
- * @return Pointer to real particle if @p p is a virtual site,
- *         otherwise nullptr. If the real particle is not found,
- *         queue a runtime error.
+ * @return Pointer to real particle.
  */
 static Particle *get_reference_particle(Particle const &p) {
   if (!p.is_virtual()) {
diff --git a/src/python/espressomd/cellsystem.pyx b/src/python/espressomd/cellsystem.pyx
index 3c023ea39f7..8a76ea2f759 100644
--- a/src/python/espressomd/cellsystem.pyx
+++ b/src/python/espressomd/cellsystem.pyx
@@ -126,7 +126,7 @@ cdef class CellSystem:
             pairs = mpi_get_pairs(distance)
         else:
             pairs = self._get_pairs_of_types(distance, types)
-        handle_errors("")
+        handle_errors("Error in get_pairs()")
         return pairs
 
     def _get_pairs_of_types(self, distance, types):
diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx
index 9851626ab2f..7cad1ef1861 100644
--- a/src/python/espressomd/integrate.pyx
+++ b/src/python/espressomd/integrate.pyx
@@ -210,11 +210,12 @@ cdef class Integrator:
 
         """
         check_type_or_throw_except(steps, 1, int, "steps must be an int")
-        assert steps >= 0, "steps has to be positive"
         check_type_or_throw_except(
             recalc_forces, 1, bool, "recalc_forces has to be a bool")
         check_type_or_throw_except(
             reuse_forces, 1, bool, "reuse_forces has to be a bool")
+        if steps < 0:
+            raise ValueError("steps must be positive")
 
         _integrate(steps, recalc_forces, reuse_forces)
 
diff --git a/src/python/espressomd/interactions.pyx b/src/python/espressomd/interactions.pyx
index 9b82e08cd61..2a411679c22 100644
--- a/src/python/espressomd/interactions.pyx
+++ b/src/python/espressomd/interactions.pyx
@@ -24,7 +24,7 @@ import collections
 include "myconfig.pxi"
 from . import utils
 from .utils import is_valid_type
-from .utils cimport check_type_or_throw_except
+from .utils cimport check_type_or_throw_except, handle_errors
 from .script_interface import ScriptObjectRegistry, ScriptInterfaceHelper, script_interface_register
 
 
@@ -73,9 +73,7 @@ cdef class NonBondedInteraction:
 
         """
         temp_params = self._get_params_from_es_core()
-        if self._params != temp_params:
-            return False
-        return True
+        return self._params == temp_params
 
     def get_params(self):
         """Get interaction parameters.
@@ -116,14 +114,15 @@ cdef class NonBondedInteraction:
 
         # If this instance refers to an interaction defined in the ESPResSo core,
         # load the parameters from there
+        is_valid_ia = self._part_types[0] >= 0 and self._part_types[1] >= 0
 
-        if self._part_types[0] >= 0 and self._part_types[1] >= 0:
+        if is_valid_ia:
             self._params = self._get_params_from_es_core()
 
         # Put in values given by the user
         self._params.update(p)
 
-        if self._part_types[0] >= 0 and self._part_types[1] >= 0:
+        if is_valid_ia:
             self._set_params_in_es_core()
 
         # update interaction dict when user sets interaction
@@ -137,6 +136,10 @@ cdef class NonBondedInteraction:
         self.user_interactions[self._part_types[0]][
             self._part_types[1]]['type_name'] = self.type_name()
 
+        # defer exception (core and interface must always agree on parameters)
+        if is_valid_ia:
+            handle_errors(f'setting {self.type_name()} raised an error')
+
     def validate_params(self):
         """Check that parameters are valid.
 
diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd
index 99b076af16f..b5c3f9d17bc 100644
--- a/src/python/espressomd/particle_data.pxd
+++ b/src/python/espressomd/particle_data.pxd
@@ -175,7 +175,7 @@ cdef extern from "particle_data.hpp":
 
 cdef extern from "virtual_sites.hpp":
     IF VIRTUAL_SITES_RELATIVE == 1:
-        void vs_relate_to(int part_num, int relate_to)
+        void vs_relate_to(int part_num, int relate_to) except +
 
 cdef extern from "rotation.hpp":
     Vector3d convert_vector_body_to_space(const particle & p, const Vector3d & v)
diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index 623ef1e439f..de82bebd3ba 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -29,7 +29,7 @@ from .analyze cimport max_seen_particle_type
 from copy import copy
 import collections
 import functools
-from .utils import nesting_level, array_locked, is_valid_type
+from .utils import nesting_level, array_locked, is_valid_type, handle_errors
 from .utils cimport make_array_locked, make_const_span, check_type_or_throw_except
 from .utils cimport Vector3i, Vector3d, Vector4d
 from .utils cimport make_Vector3d
@@ -373,8 +373,8 @@ cdef class ParticleHandle:
                     _mass, 1, float, "Mass has to be 1 float")
                 set_particle_mass(self._id, _mass)
             ELSE:
-                raise AttributeError("You are trying to set the particle mass \
-                                     but the mass feature is not compiled in.")
+                raise AttributeError("You are trying to set the particle mass "
+                                     "but the MASS feature is not compiled in.")
 
         def __get__(self):
             self.update_particle_data()
@@ -427,6 +427,8 @@ cdef class ParticleHandle:
                 cdef Quaternion[double] q
                 check_type_or_throw_except(
                     _q, 4, float, "Quaternions has to be 4 floats.")
+                if np.linalg.norm(_q) == 0.:
+                    raise ValueError("quaternion is zero")
                 for i in range(4):
                     q[i] = _q[i]
                 set_particle_quat(self._id, q)
@@ -642,6 +644,8 @@ cdef class ParticleHandle:
             def __set__(self, q):
                 check_type_or_throw_except(
                     q, 4, float, "vs_quat has to be an array-like of length 4")
+                if np.linalg.norm(q) == 0.:
+                    raise ValueError("quaternion is zero")
                 cdef Quaternion[double] _q
                 for i in range(4):
                     _q[i] = q[i]
@@ -680,6 +684,8 @@ cdef class ParticleHandle:
                     dist, 1, float, "The distance has to be given as a float.")
                 check_type_or_throw_except(
                     quat, 4, float, "The quaternion has to be given as a tuple of 4 floats.")
+                if np.linalg.norm(quat) == 0.:
+                    raise ValueError("quaternion is zero")
                 cdef Quaternion[double] q
                 for i in range(4):
                     q[i] = quat[i]
@@ -698,7 +704,7 @@ cdef class ParticleHandle:
         def vs_auto_relate_to(self, rel_to):
             """
             Setup this particle as virtual site relative to the particle
-            in argument ``rel_to``.
+            in argument ``rel_to``. A particle cannot relate to itself.
 
             Parameters
             -----------
@@ -713,6 +719,7 @@ cdef class ParticleHandle:
             check_type_or_throw_except(
                 rel_to, 1, int, "Argument of vs_auto_relate_to has to be of type ParticleHandle or int.")
             vs_relate_to(self._id, rel_to)
+            handle_errors('vs_auto_relate_to')
 
     IF DIPOLES:
         property dip:
diff --git a/src/python/espressomd/script_interface.pyx b/src/python/espressomd/script_interface.pyx
index 531639c8f1c..42746780e47 100644
--- a/src/python/espressomd/script_interface.pyx
+++ b/src/python/espressomd/script_interface.pyx
@@ -150,7 +150,7 @@ cdef class PScriptInterface:
 
         res = variant_to_python_object(
             self.sip.get().call_method(to_char_pointer(method), parameters))
-        handle_errors("")
+        handle_errors(f'while calling method {method}()')
         return res
 
     def name(self):
diff --git a/src/python/espressomd/system.pyx b/src/python/espressomd/system.pyx
index c25c8bd279c..e4a46691ea5 100644
--- a/src/python/espressomd/system.pyx
+++ b/src/python/espressomd/system.pyx
@@ -511,6 +511,4 @@ cdef class System:
 
         """
         check_type_or_throw_except(type, 1, int, "type must be 1 int")
-        number = number_of_particles_with_type(type)
-        handle_errors("")
-        return int(number)
+        return number_of_particles_with_type(type)
diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
index 1b5868fe461..a3e5fcd8b86 100644
--- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
@@ -70,8 +70,8 @@ class ReactionAlgorithm : public AutoParameters<ReactionAlgorithm> {
           [this]() { return RE()->get_exclusion_range(); }},
          {"exclusion_radius_per_type",
           [this](Variant const &v) {
-            RE()->exclusion_radius_per_type = get_map<int, double>(
-                get_value<std::unordered_map<int, Variant>>(v));
+            RE()->set_exclusion_radius_per_type(get_map<int, double>(
+                get_value<std::unordered_map<int, Variant>>(v)));
           },
           [this]() { return make_map(RE()->exclusion_radius_per_type); }}});
   }
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index 38a9b15527f..c5197445c75 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -264,6 +264,7 @@ python_test(FILE elc_vs_analytic.py MAX_NUM_PROC 2)
 python_test(FILE rotation.py MAX_NUM_PROC 1)
 python_test(FILE shapes.py MAX_NUM_PROC 1)
 python_test(FILE h5md.py MAX_NUM_PROC 2)
+python_test(FILE h5md.py MAX_NUM_PROC 1 SUFFIX 1_core)
 python_test(FILE mdanalysis.py MAX_NUM_PROC 2)
 python_test(FILE p3m_fft.py MAX_NUM_PROC 6)
 if(${TEST_NP} GREATER_EQUAL 8)
diff --git a/testsuite/python/field_test.py b/testsuite/python/field_test.py
index 27b5aa13e11..7f80073756c 100644
--- a/testsuite/python/field_test.py
+++ b/testsuite/python/field_test.py
@@ -267,6 +267,23 @@ def test_flow_field(self):
             np.testing.assert_allclose(
                 -gamma * (p.v - f_val), np.copy(p.f), atol=1e-12)
 
+    def test_invalid_interpolated_field(self):
+        h = np.array([1., 1., 1.])
+        Field = espressomd.constraints.FlowField
+
+        grid = espressomd.constraints.FlowField.field_from_fn(
+            self.system.box_l, h, self.force)
+
+        err_msg = "Constraint not compatible with box size"
+        with self.assertRaisesRegex(RuntimeError, err_msg):
+            field = Field(field=grid[:, :, :-2], grid_spacing=h, gamma=1.)
+            self.system.constraints.add(field)
+        with self.assertRaisesRegex(RuntimeError, err_msg):
+            field = Field(field=grid[:, 2:, :], grid_spacing=h, gamma=1.)
+            self.system.constraints.add(field)
+
+        self.assertEqual(len(self.system.constraints), 0)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/python/h5md.py b/testsuite/python/h5md.py
index 57ede3d2ff2..8bf74639585 100644
--- a/testsuite/python/h5md.py
+++ b/testsuite/python/h5md.py
@@ -72,6 +72,7 @@ class H5mdTests(ut.TestCase):
 
     system.integrator.run(steps=0)
     system.time = 12.3
+    n_nodes = system.cell_system.get_state()["n_nodes"]
 
     @classmethod
     def setUpClass(cls):
@@ -85,6 +86,7 @@ def setUpClass(cls):
         h5.write()
         h5.flush()
         h5.close()
+        cls.h5_params = h5.get_params()
         cls.py_file = h5py.File(cls.temp_file, 'r')
         cls.py_pos = cls.py_file['particles/atoms/position/value'][1]
         cls.py_img = cls.py_file['particles/atoms/image/value'][1]
@@ -106,6 +108,23 @@ def test_opening(self):
         h5 = espressomd.io.writer.h5md.H5md(file_path=self.temp_file)
         h5.close()
 
+    @ut.skipIf(n_nodes > 1, "only runs for 1 MPI rank")
+    def test_exceptions(self):
+        h5md = espressomd.io.writer.h5md
+        h5_units = h5md.UnitSystem(time='ps', mass='u', length='m', charge='e')
+        temp_file = os.path.join(self.temp_dir.name, 'exceptions.h5')
+        # write a non-compliant file
+        with open(temp_file, 'w') as _:
+            pass
+        with self.assertRaisesRegex(RuntimeError, 'not a valid HDF5 file'):
+            h5md.H5md(file_path=temp_file, unit_system=h5_units)
+        # write a leftover backup file
+        os.remove(temp_file)
+        with open(temp_file + '.bak', 'w') as _:
+            pass
+        with self.assertRaisesRegex(RuntimeError, 'A backup of the .h5 file exists'):
+            h5md.H5md(file_path=temp_file, unit_system=h5_units)
+
     def test_box(self):
         np.testing.assert_allclose(self.py_box, self.box_l)
 
@@ -194,6 +213,18 @@ def get_unit(path):
         self.assertEqual(get_unit('particles/atoms/force/value'), b'm u ps-2')
         self.assertEqual(get_unit('particles/atoms/velocity/value'), b'm ps-1')
 
+    def test_getters(self):
+        self.assertEqual(self.h5_params['file_path'], self.temp_file)
+        self.assertEqual(os.path.abspath(self.h5_params['script_path']),
+                         os.path.abspath(__file__))
+        self.assertEqual(self.h5_params['time_unit'], 'ps')
+        if espressomd.has_features(['ELECTROSTATICS']):
+            self.assertEqual(self.h5_params['charge_unit'], 'e')
+        if espressomd.has_features(['MASS']):
+            self.assertEqual(self.h5_params['mass_unit'], 'u')
+        self.assertEqual(self.h5_params['force_unit'], 'm u ps-2')
+        self.assertEqual(self.h5_params['velocity_unit'], 'm ps-1')
+
     def test_links(self):
         time_ref = self.py_id_time
         step_ref = self.py_id_step
diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py
index 08360150939..1f697854bfe 100644
--- a/testsuite/python/integrator_exceptions.py
+++ b/testsuite/python/integrator_exceptions.py
@@ -17,6 +17,9 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 import espressomd
+import espressomd.interactions
+import espressomd.lees_edwards
+import espressomd.shapes
 import unittest as ut
 import unittest_decorators as utx
 
@@ -25,7 +28,6 @@ class Integrator_test(ut.TestCase):
 
     system = espressomd.System(box_l=[1., 1., 1.])
     system.time_step = 0.01
-    system.cell_system.skin = 0.4
     msg = 'Encountered errors during integrate: ERROR: '
 
     def setUp(self):
@@ -36,27 +38,67 @@ def setUp(self):
     def tearDown(self):
         self.system.thermostat.turn_off()
         self.system.part.clear()
+        self.system.constraints.clear()
+        self.system.lees_edwards.protocol = None
+
+    def test_00_common_interface(self):
+        self.system.integrator.set_vv()
+        with self.assertRaisesRegex(ValueError, 'time_step must be > 0.'):
+            self.system.time_step = -2.
+        with self.assertRaisesRegex(ValueError, 'time_step must be > 0.'):
+            self.system.time_step = 0.
+        with self.assertRaisesRegex(ValueError, 'steps must be positive'):
+            self.system.integrator.run(steps=-1)
+        with self.assertRaisesRegex(Exception, self.msg + 'cannot automatically determine skin, please set it manually'):
+            self.system.integrator.run()
+        self.system.cell_system.skin = 0.4
+        with self.assertRaisesRegex(Exception, self.msg + 'cannot reuse old forces and recalculate forces'):
+            self.system.integrator.run(recalc_forces=True, reuse_forces=True)
+        if espressomd.has_features("WCA"):
+            wca = self.system.non_bonded_inter[0, 0].wca
+            wca.set_params(epsilon=1., sigma=0.01)
+            wall = espressomd.shapes.Wall(normal=[0., 0., 1.], dist=100.)
+            self.system.constraints.add(
+                shape=wall, particle_type=0, penetrable=False)
+            with self.assertRaisesRegex(Exception, self.msg + 'Constraint violated by particle 0 dist -100'):
+                self.system.integrator.run(0)
+            with self.assertRaisesRegex(Exception, 'Constraint violated by particle 0'):
+                self.system.analysis.energy()
+            wca.set_params(epsilon=0., sigma=0.)
 
     def test_vv_integrator(self):
+        self.system.cell_system.skin = 0.4
         self.system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=42)
         self.system.integrator.set_vv()
         with self.assertRaisesRegex(Exception, self.msg + 'The VV integrator is incompatible with the currently active combination of thermostats'):
             self.system.integrator.run(0)
 
     def test_brownian_integrator(self):
+        self.system.cell_system.skin = 0.4
         self.system.integrator.set_brownian_dynamics()
         with self.assertRaisesRegex(Exception, self.msg + 'The BD integrator requires the BD thermostat'):
             self.system.integrator.run(0)
 
     @utx.skipIfMissingFeatures("NPT")
     def test_npt_integrator(self):
+        self.system.cell_system.skin = 0.4
         self.system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=42)
         self.system.integrator.set_isotropic_npt(ext_pressure=1.0, piston=1.0)
         with self.assertRaisesRegex(Exception, self.msg + 'The NpT integrator requires the NpT thermostat'):
             self.system.integrator.run(0)
+        self.system.thermostat.turn_off()
+        self.system.thermostat.set_npt(kT=1.0, gamma0=2., gammav=0.04, seed=42)
+        with self.assertRaisesRegex(Exception, self.msg + 'The NpT integrator cannot use Lees-Edwards'):
+            self.system.lees_edwards.protocol = espressomd.lees_edwards.LinearShear(
+                shear_velocity=1., initial_pos_offset=0., time_0=0.)
+            self.system.integrator.run(0)
+        with self.assertRaisesRegex(Exception, self.msg + 'The NpT integrator cannot use Lees-Edwards'):
+            self.system.lees_edwards.protocol = espressomd.lees_edwards.Off()
+            self.system.integrator.run(0)
 
     @utx.skipIfMissingFeatures("STOKESIAN_DYNAMICS")
     def test_stokesian_integrator(self):
+        self.system.cell_system.skin = 0.4
         self.system.periodicity = 3 * [False]
         self.system.thermostat.set_langevin(kT=1.0, gamma=1.0, seed=42)
         self.system.integrator.set_stokesian_dynamics(
@@ -65,6 +107,7 @@ def test_stokesian_integrator(self):
             self.system.integrator.run(0)
 
     def test_steepest_descent_integrator(self):
+        self.system.cell_system.skin = 0.4
         with self.assertRaisesRegex(ValueError, r"The following keys have to be given as keyword arguments: "
                                                 r"\['f_max', 'gamma', 'max_displacement'\], got "
                                                 r"\['f_max', 'gamma', 'max_d'\] \(missing \['max_displacement'\]\)"):
diff --git a/testsuite/python/interactions_bonded_interface.py b/testsuite/python/interactions_bonded_interface.py
index dda1e79f3f4..d01069639e3 100644
--- a/testsuite/python/interactions_bonded_interface.py
+++ b/testsuite/python/interactions_bonded_interface.py
@@ -279,6 +279,52 @@ def test_exceptions(self):
         with self.assertRaisesRegex(ValueError, error_msg_not_yet_defined):
             self.system.bonded_inter[0]
 
+        # check cutoff exceptions
+        skin = self.system.cell_system.skin
+        box_l = self.system.box_l
+        node_grid = self.system.cell_system.node_grid
+        n_nodes = self.system.cell_system.get_state()['n_nodes']
+        max_ia_cutoff = min(box_l / node_grid) - skin * (n_nodes > 1)
+        err_msg = r"while calling method insert\(\): ERROR: "
+        if n_nodes == 1:
+            safe_cut = 0.5 * max_ia_cutoff
+            err_msg += r"number of cells 1 is smaller than minimum 8"
+        else:
+            safe_cut = 1.0 * max_ia_cutoff
+            err_msg += r"interaction range .+ in direction [0-2] is larger than the local box size"
+
+        # a bond with exactly the right cutoff can be added
+        h0 = espressomd.interactions.HarmonicBond(k=1., r_0=0., r_cut=safe_cut)
+        self.system.bonded_inter.add(h0)
+        self.assertEqual(len(self.system.bonded_inter), 1)
+        self.system.bonded_inter.clear()
+
+        # a bond with a cutoff larger than the box can be added, but throws
+        big_cut = 1.001 * safe_cut
+        h1 = espressomd.interactions.HarmonicBond(k=1., r_0=0., r_cut=big_cut)
+        with self.assertRaisesRegex(Exception, err_msg):
+            self.system.bonded_inter.add(h1)
+        self.assertEqual(len(self.system.bonded_inter), 1)
+        self.system.bonded_inter.clear()
+
+        # a dihedral halves the cutoff
+        safe_cut /= 2.
+        h2 = espressomd.interactions.HarmonicBond(k=1., r_0=0., r_cut=safe_cut)
+        self.system.bonded_inter.add(h2)
+        dihe = espressomd.interactions.Dihedral(bend=1., mult=1, phase=0.)
+        self.system.bonded_inter.add(dihe)
+        self.assertEqual(len(self.system.bonded_inter), 2)
+        self.system.bonded_inter.clear()
+
+        # a dihedral halves the cutoff, safe cutoffs become unsafe
+        half_cut = big_cut / 2.
+        h3 = espressomd.interactions.HarmonicBond(k=1., r_0=0., r_cut=half_cut)
+        self.system.bonded_inter.add(h3)
+        with self.assertRaisesRegex(Exception, err_msg):
+            self.system.bonded_inter.add(dihe)
+        self.assertEqual(len(self.system.bonded_inter), 2)
+        self.system.bonded_inter.clear()
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/python/interactions_non-bonded_interface.py b/testsuite/python/interactions_non-bonded_interface.py
index b842fb8394b..59e9ebdb14d 100644
--- a/testsuite/python/interactions_non-bonded_interface.py
+++ b/testsuite/python/interactions_non-bonded_interface.py
@@ -183,6 +183,19 @@ def test_exceptions(self):
             self.system.non_bonded_inter[0, 0].lennard_jones.set_params(
                 epsilon=1., sigma=2., cutoff=3., shift=4., unknown=5.)
 
+        skin = self.system.cell_system.skin
+        box_l = self.system.box_l
+        node_grid = self.system.cell_system.node_grid
+        n_nodes = self.system.cell_system.get_state()['n_nodes']
+        max_ia_cutoff = min(box_l / node_grid) - skin * (n_nodes > 1)
+        wrong_cutoff = 1.01 * max_ia_cutoff
+        lennard_jones = self.system.non_bonded_inter[0, 0].lennard_jones
+        with self.assertRaisesRegex(Exception, "setting LennardJones raised an error: ERROR: interaction range .+ in direction [0-2] is larger than the local box size"):
+            lennard_jones.set_params(
+                epsilon=1., sigma=1., cutoff=wrong_cutoff, shift="auto")
+        self.assertAlmostEqual(
+            lennard_jones.get_params()['cutoff'], wrong_cutoff, delta=1e-10)
+
 
 if __name__ == "__main__":
     ut.main()
diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py
index 43f704c105e..bf182147949 100644
--- a/testsuite/python/particle.py
+++ b/testsuite/python/particle.py
@@ -205,6 +205,10 @@ def test_vs_relative(self):
             p1.vs_relative = (0, '5', (1, 0, 0, 0))
         with self.assertRaisesRegex(ValueError, "quaternion has to be given as a tuple of 4 floats"):
             p1.vs_relative = (0, 5.0, (1, 0, 0))
+        with self.assertRaisesRegex(ValueError, "quaternion is zero"):
+            p1.vs_relative = (0, 5.0, (0, 0, 0, 0))
+        with self.assertRaisesRegex(ValueError, "quaternion is zero"):
+            p1.vs_quat = [0, 0, 0, 0]
 
     @utx.skipIfMissingFeatures("DIPOLES")
     def test_contradicting_properties_dip_dipm(self):
@@ -217,6 +221,12 @@ def test_contradicting_properties_dip_quat(self):
             self.system.part.add(pos=[0, 0, 0], dip=[1, 1, 1],
                                  quat=[1.0, 1.0, 1.0, 1.0])
 
+    @utx.skipIfMissingFeatures(["ROTATION"])
+    def test_invalid_quat(self):
+        system = self.system
+        with self.assertRaisesRegex(ValueError, "quaternion is zero"):
+            system.part.add(pos=[0., 0., 0.], quat=[0., 0., 0., 0.])
+
     @utx.skipIfMissingFeatures("ELECTROSTATICS")
     def test_particle_selection(self):
         self.system.part.clear()
diff --git a/testsuite/python/reaction_methods.py b/testsuite/python/reaction_methods.py
index d82a5119698..f9076579d53 100644
--- a/testsuite/python/reaction_methods.py
+++ b/testsuite/python/reaction_methods.py
@@ -269,19 +269,29 @@ def test_exceptions(self):
                 x=1, **single_reaction_params)
         with self.assertRaisesRegex(ValueError, err_msg):
             espressomd.reaction_ensemble.ReactionEnsemble(
-                kT=1., exclusion_range=1., seed=12, x=1, exclusion_radius_per_type={1: 0.1})
+                kT=1., exclusion_range=1., seed=12, x=1)
         with self.assertRaisesRegex(ValueError, err_msg):
             espressomd.reaction_ensemble.ConstantpHEnsemble(
-                kT=1., exclusion_range=1., seed=12, x=1, constant_pH=2, exclusion_radius_per_type={1: 0.1})
+                kT=1., exclusion_range=1., seed=12, x=1, constant_pH=2)
         with self.assertRaisesRegex(ValueError, err_msg):
             espressomd.reaction_ensemble.WidomInsertion(
                 kT=1., seed=12, x=1)
         with self.assertRaisesRegex(ValueError, "Invalid value for 'kT'"):
             espressomd.reaction_ensemble.ReactionEnsemble(
-                kT=-1., exclusion_range=1., seed=12, exclusion_radius_per_type={1: 0.1})
+                kT=-1., exclusion_range=1., seed=12)
+
+        # check invalid exclusion ranges and radii
         with self.assertRaisesRegex(ValueError, "Invalid value for 'exclusion_range'"):
             espressomd.reaction_ensemble.ReactionEnsemble(
-                kT=1., exclusion_range=-1., seed=12, exclusion_radius_per_type={1: 0.1})
+                kT=1., seed=12, exclusion_range=-1.)
+        with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 1: radius -0.10"):
+            espressomd.reaction_ensemble.ReactionEnsemble(
+                kT=1., seed=12, exclusion_range=1., exclusion_radius_per_type={1: -0.1})
+        method = espressomd.reaction_ensemble.ReactionEnsemble(
+            kT=1., exclusion_range=1., seed=12, exclusion_radius_per_type={1: 0.1})
+        with self.assertRaisesRegex(ValueError, "Invalid excluded_radius value for type 2: radius -0.10"):
+            method.exclusion_radius_per_type = {2: -0.1}
+        self.assertEqual(list(method.exclusion_radius_per_type.keys()), [1])
 
 
 if __name__ == "__main__":
diff --git a/testsuite/python/virtual_sites_relative.py b/testsuite/python/virtual_sites_relative.py
index 14657efbac1..1d2dd0c64d8 100644
--- a/testsuite/python/virtual_sites_relative.py
+++ b/testsuite/python/virtual_sites_relative.py
@@ -140,16 +140,25 @@ def test_vs_exceptions(self):
         system = self.system
         system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative()
         system.time_step = 0.01
+        system.cell_system.skin = 0.1
+        system.min_global_cut = 0.1
+        p1 = system.part.add(pos=[0.0, 0.0, 0.0], rotation=[1, 1, 1], id=1)
         p2 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=[1, 1, 1], id=2)
         p3 = system.part.add(pos=[1.0, 1.0, 1.0], rotation=[1, 1, 1], id=3)
         # relating to anything else other than a particle or id is not allowed
         with self.assertRaisesRegex(ValueError, "Argument of vs_auto_relate_to has to be of type ParticleHandle or int"):
             p2.vs_auto_relate_to('0')
+        # relating to itself is not allowed
+        with self.assertRaisesRegex(ValueError, "A virtual site cannot relate to itself"):
+            p2.vs_auto_relate_to(p2)
         # relating to a deleted particle is not allowed
         with self.assertRaisesRegex(Exception, "No real particle with id 3 for virtual site with id 2"):
             p2.vs_auto_relate_to(p3)
             p3.remove()
             system.integrator.run(0, recalc_forces=True)
+        if system.cell_system.get_state()["n_nodes"] > 1:
+            with self.assertRaisesRegex(Exception, r"The distance between virtual and non-virtual particle \([0-9\.]+\) is larger than the minimum global cutoff"):
+                p2.vs_auto_relate_to(p1)
 
     def test_pos_vel_forces(self):
         system = self.system

From 61e025abbac2637a85bd8f5b3ca0d19617ee15d7 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 18 Mar 2022 02:32:26 +0100
Subject: [PATCH 06/10] script_interface: Refactor get_value() framework

Rename `make_map()` to `make_unordered_map_of_variants()` to better
reflect its role. Remove `get_map()` and instead implement getter
`template <class T> get_value<std::unordered_map<int, T>>(Variant)`.
Make Variant runtime error messages human-readable by substituting
the Variant demangled name by the string "ScriptInterface::Variant".
Reduce code duplication.
---
 src/script_interface/ObjectMap.hpp            |   7 +-
 .../constraints/couplings.hpp                 |  10 +-
 src/script_interface/get_value.hpp            | 146 ++++++++++++------
 .../reaction_methods/ConstantpHEnsemble.hpp   |   5 +-
 .../reaction_methods/ReactionAlgorithm.hpp    |  10 +-
 .../reaction_methods/ReactionEnsemble.hpp     |   6 +-
 .../reaction_methods/WidomInsertion.hpp       |   3 +-
 src/script_interface/tests/get_value_test.cpp | 130 +++++++++++-----
 8 files changed, 204 insertions(+), 113 deletions(-)

diff --git a/src/script_interface/ObjectMap.hpp b/src/script_interface/ObjectMap.hpp
index d6af7da67df..a41d7ddc902 100644
--- a/src/script_interface/ObjectMap.hpp
+++ b/src/script_interface/ObjectMap.hpp
@@ -130,12 +130,7 @@ class ObjectMap : public BaseType {
     }
 
     if (method == "get_map") {
-      std::unordered_map<KeyType, Variant> ret;
-
-      for (auto const &kv : m_elements)
-        ret[kv.first] = kv.second;
-
-      return ret;
+      return make_unordered_map_of_variants(m_elements);
     }
 
     if (method == "keys") {
diff --git a/src/script_interface/constraints/couplings.hpp b/src/script_interface/constraints/couplings.hpp
index c1ceab0de49..b96831bd883 100644
--- a/src/script_interface/constraints/couplings.hpp
+++ b/src/script_interface/constraints/couplings.hpp
@@ -74,8 +74,9 @@ template <> struct coupling_parameters_impl<Scaled> {
                 AutoParameter::read_only,
                 [this_]() { return this_().default_scale(); },
             },
-            {"particle_scales", AutoParameter::read_only,
-             [this_]() { return make_map(this_().particle_scales()); }}};
+            {"particle_scales", AutoParameter::read_only, [this_]() {
+               return make_unordered_map_of_variants(this_().particle_scales());
+             }}};
   }
 };
 
@@ -90,9 +91,8 @@ template <> inline Viscous make_coupling<Viscous>(const VariantMap &params) {
 }
 
 template <> inline Scaled make_coupling<Scaled>(const VariantMap &params) {
-  auto const particle_scales = get_value_or<std::unordered_map<int, Variant>>(
-      params, "particle_scales", {});
-  return Scaled{get_map<int, double>(particle_scales),
+  return Scaled{get_value_or<std::unordered_map<int, double>>(
+                    params, "particle_scales", {}),
                 get_value<double>(params, "default_scale")};
 }
 } // namespace detail
diff --git a/src/script_interface/get_value.hpp b/src/script_interface/get_value.hpp
index acff18edf23..904230b8079 100644
--- a/src/script_interface/get_value.hpp
+++ b/src/script_interface/get_value.hpp
@@ -38,16 +38,64 @@
 
 namespace ScriptInterface {
 namespace detail {
+
+/**
+ * @brief Demangle the symbol of a container @c value_type or @c mapped_type.
+ * Only for recursive variant container types used in @ref Variant.
+ * Other container types and non-container objects return an empty string.
+ */
+template <typename Container> struct demangle_container_value_type {
+private:
+  template <typename T> struct is_std_vector : std::false_type {};
+  template <typename... Args>
+  struct is_std_vector<std::vector<Args...>> : std::true_type {};
+
+  template <typename T> struct is_std_unordered_map : std::false_type {};
+  template <typename... Args>
+  struct is_std_unordered_map<std::unordered_map<Args...>> : std::true_type {};
+
+public:
+  template <typename T = Container,
+            std::enable_if_t<!is_std_vector<T>::value and
+                                 !is_std_unordered_map<T>::value,
+                             bool> = true>
+  auto operator()() const {
+    return std::string("");
+  }
+
+  template <typename T = Container,
+            std::enable_if_t<is_std_vector<T>::value, bool> = true>
+  auto operator()() const {
+    return Utils::demangle<typename T::value_type>();
+  }
+
+  template <typename T = Container,
+            std::enable_if_t<is_std_unordered_map<T>::value, bool> = true>
+  auto operator()() const {
+    return Utils::demangle<typename T::mapped_type>();
+  }
+};
+
+/**
+ * @brief Demangle the data type of an object wrapped inside a @ref Variant.
+ * When the data type involves a recursive variant (e.g. containers), the
+ * demangled symbol name is processed to replace all occurences of the
+ * variant symbol name by a human-readable string "ScriptInterface::Variant".
+ */
 struct type_label_visitor : boost::static_visitor<std::string> {
   template <class T> std::string operator()(const T &) const {
-    return Utils::demangle<T>();
+    auto const symbol_for_variant = Utils::demangle<Variant>();
+    auto const name_for_variant = std::string("ScriptInterface::Variant");
+    auto symbol = Utils::demangle<T>();
+    for (std::string::size_type pos{};
+         (pos = symbol.find(symbol_for_variant, pos)) != symbol.npos;
+         pos += name_for_variant.length()) {
+      symbol.replace(pos, symbol_for_variant.length(), name_for_variant);
+    }
+    return symbol;
   }
 };
 
-inline std::string type_label(const Variant &v) {
-  return boost::apply_visitor(type_label_visitor{}, v);
-}
-
 /*
  * Allows
  * T -> T,
@@ -137,6 +185,9 @@ struct GetVectorOrEmpty : boost::static_visitor<std::vector<T>> {
 
   /* Standard case, correct type */
   std::vector<T> operator()(std::vector<T> const &v) const { return v; }
+
+  template <typename V = T,
+            std::enable_if_t<!std::is_same<V, Variant>::value, bool> = true>
   std::vector<T> operator()(std::vector<Variant> const &vv) const {
     std::vector<T> ret(vv.size());
 
@@ -147,17 +198,10 @@ struct GetVectorOrEmpty : boost::static_visitor<std::vector<T>> {
   }
 };
 
-/* std::vector cases
- * We implicitly transform an empty vector<Variant> into an empty vector<T>. */
-template <> struct get_value_helper<std::vector<int>, void> {
-  std::vector<int> operator()(Variant const &v) const {
-    return boost::apply_visitor(GetVectorOrEmpty<int>{}, v);
-  }
-};
-
-template <> struct get_value_helper<std::vector<double>, void> {
-  std::vector<double> operator()(Variant const &v) const {
-    return boost::apply_visitor(GetVectorOrEmpty<double>{}, v);
+/* std::vector cases */
+template <typename T> struct get_value_helper<std::vector<T>, void> {
+  std::vector<T> operator()(Variant const &v) const {
+    return boost::apply_visitor(GetVectorOrEmpty<T>{}, v);
   }
 };
 
@@ -172,12 +216,24 @@ struct GetMapOrEmpty : boost::static_visitor<std::unordered_map<K, T>> {
   std::unordered_map<K, T> operator()(std::unordered_map<K, T> const &v) const {
     return v;
   }
+
+  template <typename V = T,
+            std::enable_if_t<!std::is_same<V, Variant>::value, bool> = true>
+  std::unordered_map<K, T>
+  operator()(std::unordered_map<K, Variant> const &v) const {
+    std::unordered_map<K, T> ret;
+    for (auto it = v.begin(); it != v.end(); ++it) {
+      ret.insert({it->first, get_value_helper<T>{}(it->second)});
+    }
+    return ret;
+  }
 };
 
 /* std::unordered_map cases */
-template <> struct get_value_helper<std::unordered_map<int, Variant>, void> {
-  std::unordered_map<int, Variant> operator()(Variant const &v) const {
-    return boost::apply_visitor(GetMapOrEmpty<int, Variant>{}, v);
+template <typename T>
+struct get_value_helper<std::unordered_map<int, T>, void> {
+  std::unordered_map<int, T> operator()(Variant const &v) const {
+    return boost::apply_visitor(GetMapOrEmpty<int, T>{}, v);
   }
 };
 
@@ -213,24 +269,27 @@ struct get_value_helper<
  * @brief Re-throw a @c boost::bad_get exception wrapped in an @ref Exception.
  * Write a custom error message for invalid conversions due to type mismatch
  * and due to nullptr values, possibly with context information if the variant
- * is in a container.
- * @tparam T     Type which the variant was supposed to convert to
- * @tparam U     Type of the container the variant is contained in
+ * is a container.
+ * @tparam T     Which type the variant was supposed to convert to
  */
-template <typename T, typename U = void>
-inline void handle_bad_get(Variant const &v) {
-  using Utils::demangle;
-  auto const what = "Provided argument of type " + type_label(v);
-  auto const context =
-      (std::is_void<U>::value)
-          ? ""
-          : " (raised during the creation of a " + demangle<U>() + ")";
+template <typename T> inline void handle_bad_get(Variant const &v) {
+  auto const object_symbol_name = boost::apply_visitor(type_label_visitor{}, v);
+  auto const element_symbol_name = demangle_container_value_type<T>{}();
+  auto const is_container = !element_symbol_name.empty();
+  auto const what = "Provided argument of type " + object_symbol_name;
   try {
     throw;
   } catch (bad_get_nullptr const &) {
-    throw Exception(what + " is a null pointer" + context);
+    auto const item_error = (is_container) ? " contains a value that" : "";
+    throw Exception(what + item_error + " is a null pointer");
   } catch (boost::bad_get const &) {
-    throw Exception(what + " is not convertible to " + demangle<T>() + context);
+    auto const non_convertible = " is not convertible to ";
+    auto item_error = std::string("");
+    if (is_container) {
+      item_error += " because it contains a value that";
+      item_error += non_convertible + element_symbol_name;
+    }
+    throw Exception(what + non_convertible + Utils::demangle<T>() + item_error);
   }
 }
 
@@ -241,9 +300,9 @@ inline void handle_bad_get(Variant const &v) {
  *
  * This is a wrapper around boost::get that allows us to
  * customize the behavior for different types. This is
- * needed e.g. to deal with Vector types that are not
- * explicitly contained in Variant, but can easily
- * be converted.
+ * needed e.g. to deal with containers whose elements
+ * have mixed types that are implicitly convertible
+ * to a requested type.
  */
 template <typename T> T get_value(Variant const &v) {
   try {
@@ -255,22 +314,7 @@ template <typename T> T get_value(Variant const &v) {
 }
 
 template <typename K, typename V>
-std::unordered_map<K, V> get_map(std::unordered_map<K, Variant> const &v) {
-  std::unordered_map<K, V> ret;
-  auto it = v.begin();
-  try {
-    for (; it != v.end(); ++it) {
-      ret.insert({it->first, detail::get_value_helper<V>{}(it->second)});
-    }
-  } catch (...) {
-    detail::handle_bad_get<V, std::unordered_map<K, V>>(v);
-    throw;
-  }
-  return ret;
-}
-
-template <typename K, typename V>
-std::unordered_map<K, Variant> make_map(std::unordered_map<K, V> const &v) {
+auto make_unordered_map_of_variants(std::unordered_map<K, V> const &v) {
   std::unordered_map<K, Variant> ret;
   for (auto const &it : v) {
     ret.insert({it.first, Variant(it.second)});
diff --git a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
index 834bf651062..47c0a672e3f 100644
--- a/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ConstantpHEnsemble.hpp
@@ -28,6 +28,7 @@
 #include "core/reaction_methods/ReactionAlgorithm.hpp"
 
 #include <memory>
+#include <unordered_map>
 
 namespace ScriptInterface {
 namespace ReactionMethods {
@@ -53,8 +54,8 @@ class ConstantpHEnsemble : public ReactionAlgorithm {
         get_value<int>(params, "seed"), get_value<double>(params, "kT"),
         get_value<double>(params, "exclusion_range"),
         get_value<double>(params, "constant_pH"),
-        get_map<int, double>(get_value_or<std::unordered_map<int, Variant>>(
-            params, "exclusion_radius_per_type", {})));
+        get_value_or<std::unordered_map<int, double>>(
+            params, "exclusion_radius_per_type", {}));
   }
 
 private:
diff --git a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
index a3e5fcd8b86..13f990207ea 100644
--- a/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
+++ b/src/script_interface/reaction_methods/ReactionAlgorithm.hpp
@@ -29,6 +29,7 @@
 #include <memory>
 #include <stdexcept>
 #include <string>
+#include <unordered_map>
 #include <vector>
 
 namespace ScriptInterface {
@@ -70,10 +71,13 @@ class ReactionAlgorithm : public AutoParameters<ReactionAlgorithm> {
           [this]() { return RE()->get_exclusion_range(); }},
          {"exclusion_radius_per_type",
           [this](Variant const &v) {
-            RE()->set_exclusion_radius_per_type(get_map<int, double>(
-                get_value<std::unordered_map<int, Variant>>(v)));
+            RE()->set_exclusion_radius_per_type(
+                get_value<std::unordered_map<int, double>>(v));
           },
-          [this]() { return make_map(RE()->exclusion_radius_per_type); }}});
+          [this]() {
+            return make_unordered_map_of_variants(
+                RE()->exclusion_radius_per_type);
+          }}});
   }
 
   Variant do_call_method(std::string const &name,
diff --git a/src/script_interface/reaction_methods/ReactionEnsemble.hpp b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
index 03b66128420..20478e1e8b6 100644
--- a/src/script_interface/reaction_methods/ReactionEnsemble.hpp
+++ b/src/script_interface/reaction_methods/ReactionEnsemble.hpp
@@ -28,6 +28,7 @@
 #include "core/reaction_methods/ReactionEnsemble.hpp"
 
 #include <memory>
+#include <unordered_map>
 
 namespace ScriptInterface {
 namespace ReactionMethods {
@@ -39,12 +40,11 @@ class ReactionEnsemble : public ReactionAlgorithm {
   }
 
   void do_construct(VariantMap const &params) override {
-
     m_re = std::make_shared<::ReactionMethods::ReactionEnsemble>(
         get_value<int>(params, "seed"), get_value<double>(params, "kT"),
         get_value<double>(params, "exclusion_range"),
-        get_map<int, double>(get_value_or<std::unordered_map<int, Variant>>(
-            params, "exclusion_radius_per_type", {})));
+        get_value_or<std::unordered_map<int, double>>(
+            params, "exclusion_radius_per_type", {}));
   }
 
 private:
diff --git a/src/script_interface/reaction_methods/WidomInsertion.hpp b/src/script_interface/reaction_methods/WidomInsertion.hpp
index f563047bd9e..343679a0184 100644
--- a/src/script_interface/reaction_methods/WidomInsertion.hpp
+++ b/src/script_interface/reaction_methods/WidomInsertion.hpp
@@ -41,10 +41,9 @@ class WidomInsertion : public ReactionAlgorithm {
   }
 
   void do_construct(VariantMap const &params) override {
-    std::unordered_map<int, double> exclusion_radius_per_type;
     m_re = std::make_shared<::ReactionMethods::WidomInsertion>(
         get_value<int>(params, "seed"), get_value<double>(params, "kT"), 0.,
-        exclusion_radius_per_type);
+        std::unordered_map<int, double>{});
   }
 
   Variant do_call_method(std::string const &name,
diff --git a/src/script_interface/tests/get_value_test.cpp b/src/script_interface/tests/get_value_test.cpp
index d6d346ace08..3410e3fd672 100644
--- a/src/script_interface/tests/get_value_test.cpp
+++ b/src/script_interface/tests/get_value_test.cpp
@@ -24,6 +24,7 @@
 #include "script_interface/ObjectHandle.hpp"
 #include "script_interface/get_value.hpp"
 
+#include <cassert>
 #include <memory>
 #include <regex>
 #include <stdexcept>
@@ -127,67 +128,112 @@ BOOST_AUTO_TEST_CASE(get_value_from_map) {
   BOOST_CHECK_THROW((get_value<int>(map, "unknown")), std::exception);
 }
 
-BOOST_AUTO_TEST_CASE(get_map_value) {
-  using ScriptInterface::get_map;
-  using ScriptInterface::Variant;
-
-  std::unordered_map<int, Variant> const map_variant{{1, 1.5}, {2, 2.5}};
-  std::unordered_map<int, double> const map = get_map<int, double>(map_variant);
-  BOOST_CHECK_EQUAL(map.at(1), 1.5);
-  BOOST_CHECK_EQUAL(map.at(2), 2.5);
-}
-
-BOOST_AUTO_TEST_CASE(get_unordered_map) {
+BOOST_AUTO_TEST_CASE(unordered_map) {
   using ScriptInterface::get_value;
+  using ScriptInterface::make_unordered_map_of_variants;
   using ScriptInterface::Variant;
 
   auto const var = Variant{std::unordered_map<int, Variant>{{1, 1}, {2, 2.5}}};
-  auto const map = get_value<std::unordered_map<int, Variant>>(var);
-  BOOST_CHECK_EQUAL(get_value<int>(map.at(1)), 1);
-  BOOST_CHECK_EQUAL(get_value<double>(map.at(2)), 2.5);
+  {
+    auto const map_var = get_value<std::unordered_map<int, Variant>>(var);
+    BOOST_CHECK_EQUAL(get_value<int>(map_var.at(1)), 1);
+    BOOST_CHECK_EQUAL(get_value<double>(map_var.at(2)), 2.5);
+  }
+  {
+    auto const map_dbl = get_value<std::unordered_map<int, double>>(var);
+    BOOST_CHECK_EQUAL(map_dbl.at(1), 1.0);
+    BOOST_CHECK_EQUAL(map_dbl.at(2), 2.5);
+  }
+  {
+    auto const map_dbl_input = get_value<std::unordered_map<int, double>>(var);
+    auto const map_var = make_unordered_map_of_variants(map_dbl_input);
+    auto const map_dbl =
+        get_value<std::unordered_map<int, double>>(Variant{map_var});
+    BOOST_CHECK_EQUAL(map_dbl.at(1), 1.0);
+    BOOST_CHECK_EQUAL(map_dbl.at(2), 2.5);
+  }
+}
+
+auto exception_message_predicate(std::string const &pattern) {
+  return [=](std::exception const &ex) {
+    boost::test_tools::predicate_result result = true;
+    std::string const what = ex.what();
+    std::smatch match;
+    if (!std::regex_search(what, match, std::regex(pattern))) {
+      result = false;
+      result.message() << "Error message \"" << what << "\" "
+                       << "doesn't match pattern \"" << pattern << "\"";
+    }
+    return result;
+  };
 }
 
-BOOST_AUTO_TEST_CASE(exceptions) {
-  using ScriptInterface::get_map;
+BOOST_AUTO_TEST_CASE(check_exceptions) {
   using ScriptInterface::get_value;
   using ScriptInterface::Variant;
 
+  assert(!!exception_message_predicate("A")(std::runtime_error("A")));
+  assert(!exception_message_predicate("A")(std::runtime_error("B")));
+
   using so_ptr_t = std::shared_ptr<ScriptInterface::ObjectHandle>;
 
   auto const so_obj = so_ptr_t();
-  auto const so_ptr_tn = Utils::demangle<so_ptr_t>();
+  auto const msg_prefix = std::string("Provided argument of type ");
+  auto const variant_name = std::string("ScriptInterface::Variant");
 
   {
-    auto const predicate = [](std::string const &type, std::string const &why) {
-      auto const message = "Provided argument of type " + type + " is " + why;
-      return [=](std::exception const &ex) { return ex.what() == message; };
-    };
-    auto const so_variant = Variant(so_obj);
-    BOOST_CHECK_EXCEPTION((get_value<so_ptr_t>(so_variant)), std::exception,
-                          predicate(so_ptr_tn, "a null pointer"));
-    BOOST_CHECK_EXCEPTION((get_value<int>(so_variant)), std::exception,
-                          predicate(so_ptr_tn, "not convertible to int"));
+    // basic types
+    auto const obj_variant = Variant{so_obj};
+    auto const obj_variant_pattern = Utils::demangle<so_ptr_t>();
+    auto const what = msg_prefix + obj_variant_pattern;
+    auto const predicate_nullptr =
+        exception_message_predicate(what + " is a null pointer");
+    auto const predicate_conversion =
+        exception_message_predicate(what + " is not convertible to int");
+    BOOST_CHECK_EXCEPTION(get_value<so_ptr_t>(obj_variant), std::exception,
+                          predicate_nullptr);
+    BOOST_CHECK_EXCEPTION(get_value<int>(obj_variant), std::exception,
+                          predicate_conversion);
+  }
+  {
+    // vectors
+    auto const vec_variant = Variant{std::vector<Variant>{{so_obj}}};
+    auto const vec_variant_pattern =
+        "std::(__1::)?vector<" + variant_name + ", .*?>";
+    auto const what = msg_prefix + vec_variant_pattern;
+    auto const predicate_nullptr = exception_message_predicate(
+        what + " contains a value that is a null pointer");
+    auto const predicate_conversion = exception_message_predicate(
+        what + " is not convertible to std::(__1::)?vector<int, .*?> because "
+               "it contains a value that is not convertible to int");
+    BOOST_CHECK_EXCEPTION(get_value<std::vector<so_ptr_t>>(vec_variant),
+                          std::exception, predicate_nullptr);
+    BOOST_CHECK_EXCEPTION(get_value<std::vector<int>>(vec_variant),
+                          std::exception, predicate_conversion);
   }
   {
-    auto const predicate = [](std::string const &why, std::string const &type) {
-      auto const msg = "during the creation of a std::(__1::)?unordered_map";
-      auto const pattern = why + " \\(raised " + msg + "<int, " + type + ", ";
-      return [=](std::exception const &ex) {
-        std::string const what = ex.what();
-        std::smatch match;
-        return std::regex_search(what, match, std::regex(pattern));
-      };
-    };
-    auto const so_map = std::unordered_map<int, Variant>{{1, so_obj}};
-    BOOST_CHECK_EXCEPTION((get_map<int, so_ptr_t>(so_map)), std::exception,
-                          predicate("is a null pointer", so_ptr_tn));
-    BOOST_CHECK_EXCEPTION((get_map<int, int>(so_map)), std::exception,
-                          predicate("is not convertible to int", "int"));
+    // unordered maps
+    auto const map_variant =
+        Variant{std::unordered_map<int, Variant>{{1, so_obj}}};
+    auto const map_variant_pattern =
+        "std::(__1::)?unordered_map<int, " + variant_name + ", .*?>";
+    auto const what = msg_prefix + map_variant_pattern;
+    auto const predicate_nullptr = exception_message_predicate(
+        what + " contains a value that is a null pointer");
+    auto const predicate_conversion = exception_message_predicate(
+        what +
+        " is not convertible to std::(__1::)?unordered_map<int, int, .*?> "
+        "because it contains a value that is not convertible to int");
+    BOOST_CHECK_EXCEPTION(
+        (get_value<std::unordered_map<int, so_ptr_t>>(map_variant)),
+        std::exception, predicate_nullptr);
+    BOOST_CHECK_EXCEPTION(
+        (get_value<std::unordered_map<int, int>>(map_variant)), std::exception,
+        predicate_conversion);
   }
   {
     using Utils::Vector3d;
     std::unordered_map<int, Variant> const mixed{{1, 1}, {2, std::string("2")}};
-    BOOST_CHECK_THROW((get_map<int, double>(mixed)), std::exception);
     std::vector<Variant> const v_var = {1., 2.};
     std::vector<double> const v_dbl = {1., 2.};
     BOOST_CHECK_THROW((get_value<Vector3d>(Variant{v_var})), std::exception);
@@ -199,5 +245,7 @@ BOOST_AUTO_TEST_CASE(exceptions) {
     BOOST_CHECK_THROW((get_value<int>(Variant{v_dbl})), std::exception);
     BOOST_CHECK_THROW((get_value<std::unordered_map<int, Variant>>(Variant{})),
                       std::exception);
+    BOOST_CHECK_THROW((get_value<std::unordered_map<int, int>>(Variant{mixed})),
+                      std::exception);
   }
 }

From 78a22db3104a881f4df7b6318fbc5514c32da18d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 18 Mar 2022 17:15:08 +0100
Subject: [PATCH 07/10] core: Fix undefined behavior

When LB_BOUNDARIES is not compiled in, the return type changes and
the function pointer stored in the MpiCallbacks framework accesses
the wrong data type.
---
 src/core/grid_based_algorithms/lb_collective_interface.cpp | 4 ++--
 testsuite/python/lb_vtk.py                                 | 2 ++
 2 files changed, 4 insertions(+), 2 deletions(-)

diff --git a/src/core/grid_based_algorithms/lb_collective_interface.cpp b/src/core/grid_based_algorithms/lb_collective_interface.cpp
index 4ea9ddcd5ef..2638f091602 100644
--- a/src/core/grid_based_algorithms/lb_collective_interface.cpp
+++ b/src/core/grid_based_algorithms/lb_collective_interface.cpp
@@ -109,14 +109,14 @@ auto mpi_lb_get_populations(Utils::Vector3i const &index) {
 
 REGISTER_CALLBACK_ONE_RANK(mpi_lb_get_populations)
 
-auto mpi_lb_get_boundary_flag(Utils::Vector3i const &index) {
+boost::optional<int> mpi_lb_get_boundary_flag(Utils::Vector3i const &index) {
   return detail::lb_calc(index, [&](auto index) {
 #ifdef LB_BOUNDARIES
     auto const linear_index =
         get_linear_index(lblattice.local_index(index), lblattice.halo_grid);
     return lbfields[linear_index].boundary;
 #else
-    return false;
+    return 0;
 #endif
   });
 }
diff --git a/testsuite/python/lb_vtk.py b/testsuite/python/lb_vtk.py
index 8cd8e065c02..2e4e28f8088 100644
--- a/testsuite/python/lb_vtk.py
+++ b/testsuite/python/lb_vtk.py
@@ -142,6 +142,8 @@ def test_vtk(self):
         vtk_boundary = self.parse_vtk(
             'vtk_out/boundary.vtk', 'boundary', shape)
         np.testing.assert_equal(vtk_boundary, node_boundary.astype(int))
+        if self.system.lbboundaries is None:
+            np.testing.assert_equal(np.sum(node_boundary), 0.)
 
     def test_print(self):
         '''

From 9ad55b9003e8a113b02c990e3382b56b5f631ec5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 18 Mar 2022 22:58:00 +0100
Subject: [PATCH 08/10] python: Fix OpenGL visualization

---
 samples/lj-demo.py                            |  2 +-
 src/python/espressomd/visualization_opengl.py | 23 ++++++++-----------
 2 files changed, 10 insertions(+), 15 deletions(-)

diff --git a/samples/lj-demo.py b/samples/lj-demo.py
index 46173f31cbf..eeab45f217a 100644
--- a/samples/lj-demo.py
+++ b/samples/lj-demo.py
@@ -454,7 +454,7 @@ def main_loop():
         new_box = np.ones(3) * controls.volume**(1. / 3.)
         if np.any(np.array(system.box_l) != new_box):
             for p in system.part:
-                p.pos *= new_box / system.box_l[0]
+                p.pos = p.pos * new_box / system.box_l[0]
             print("volume changed")
             system.force_cap = lj_cap
         system.box_l = new_box
diff --git a/src/python/espressomd/visualization_opengl.py b/src/python/espressomd/visualization_opengl.py
index 3984349bc66..9b4b7588fe4 100644
--- a/src/python/espressomd/visualization_opengl.py
+++ b/src/python/espressomd/visualization_opengl.py
@@ -1587,22 +1587,17 @@ def _init_espresso_visualization(self):
         self.depth = 0
 
         # LOOK FOR LB ACTOR
-        if self.specs['LB_draw_velocity_plane'] or \
-                self.specs['LB_draw_nodes'] or \
-                self.specs['LB_draw_node_boundaries']:
-            lb_types = [espressomd.lb.LBFluid]
-            if espressomd.has_features('CUDA'):
-                lb_types.append(espressomd.lb.LBFluidGPU)
-            for a in self.system.actors:
-                if isinstance(a, tuple(lb_types)):
-                    # if 'agrid' in pa:
-                    self.lb_params = a.get_params()
-                    self.lb = a
-                    self.lb_is_cpu = isinstance(a, espressomd.lb.LBFluid)
-                    break
+        lb_types = [espressomd.lb.LBFluid]
+        if espressomd.has_features('CUDA'):
+            lb_types.append(espressomd.lb.LBFluidGPU)
+        for actor in self.system.actors:
+            if isinstance(actor, tuple(lb_types)):
+                self.lb_params = actor.get_params()
+                self.lb = actor
+                self.lb_is_cpu = isinstance(actor, espressomd.lb.LBFluid)
+                break
 
         if self.specs['LB_draw_velocity_plane']:
-
             if self.specs['LB_plane_axis'] == 0:
                 pn = [1.0, 0.0, 0.0]
                 self.lb_plane_b1 = [0.0, 1.0, 0.0]

From 4d319c37584ae4499a5a9d9d79295a03218dec95 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 18 Mar 2022 23:04:40 +0100
Subject: [PATCH 09/10] testsuite: Fix regressions

---
 src/core/unit_tests/rotation_test.cpp | 10 +++++++---
 testsuite/python/CMakeLists.txt       |  1 +
 2 files changed, 8 insertions(+), 3 deletions(-)

diff --git a/src/core/unit_tests/rotation_test.cpp b/src/core/unit_tests/rotation_test.cpp
index 5ff5d98cb6c..336f16cd75f 100644
--- a/src/core/unit_tests/rotation_test.cpp
+++ b/src/core/unit_tests/rotation_test.cpp
@@ -19,10 +19,14 @@
 
 #define BOOST_TEST_MODULE rotation test
 #define BOOST_TEST_DYN_LINK
+
+#include "config.hpp"
+
+#ifdef ROTATION
+
 #include <boost/test/unit_test.hpp>
 
 #include "Particle.hpp"
-#include "config.hpp"
 #include "rotation.hpp"
 
 #include <utils/Vector.hpp>
@@ -33,8 +37,6 @@
 #include <stdexcept>
 #include <tuple>
 
-#ifdef ROTATION
-
 auto constexpr tol = 5. * 100. * std::numeric_limits<double>::epsilon();
 
 namespace Testing {
@@ -216,4 +218,6 @@ BOOST_AUTO_TEST_CASE(convert_dip_to_quat_test) {
   }
 }
 #endif // DIPOLES
+#else  // ROTATION
+int main(int argc, char **argv) {}
 #endif // ROTATION
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index c5197445c75..6b91021da26 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -127,6 +127,7 @@ python_test(FILE mass-and-rinertia_per_particle.py MAX_NUM_PROC 2 LABELS long)
 python_test(FILE integrate.py MAX_NUM_PROC 4)
 python_test(FILE interactions_bond_angle.py MAX_NUM_PROC 4)
 python_test(FILE interactions_bonded_interface.py MAX_NUM_PROC 4)
+python_test(FILE interactions_bonded_interface.py MAX_NUM_PROC 1 SUFFIX 1_core)
 python_test(FILE interactions_bonded.py MAX_NUM_PROC 2)
 python_test(FILE interactions_dihedral.py MAX_NUM_PROC 4)
 python_test(FILE interactions_non-bonded_interface.py MAX_NUM_PROC 4)

From db253de9126b060ccf00cae7bfd2f9127a0df5af Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 18 Mar 2022 23:37:15 +0100
Subject: [PATCH 10/10] Improve particle exception handling

Convert fatal errors into ValueError for invalid particle ids.
Make python code for particle creation more readable.
---
 src/core/particle_data.cpp              |  8 ++-
 src/python/espressomd/particle_data.pxd |  2 +-
 src/python/espressomd/particle_data.pyx | 65 ++++++++++++-------------
 testsuite/python/particle.py            | 14 +++++-
 4 files changed, 51 insertions(+), 38 deletions(-)

diff --git a/src/core/particle_data.cpp b/src/core/particle_data.cpp
index ff011db2825..17d8864a258 100644
--- a/src/core/particle_data.cpp
+++ b/src/core/particle_data.cpp
@@ -487,8 +487,9 @@ void build_particle_node() { mpi_who_has(); }
  *  @brief Get the mpi rank which owns the particle with id.
  */
 int get_particle_node(int id) {
-  if (id < 0)
-    throw std::runtime_error("Invalid particle id!");
+  if (id < 0) {
+    throw std::domain_error("Invalid particle id: " + std::to_string(id));
+  }
 
   if (particle_node.empty())
     build_particle_node();
@@ -724,6 +725,9 @@ void mpi_place_particle(int node, int p_id, const Utils::Vector3d &pos) {
 }
 
 int place_particle(int p_id, Utils::Vector3d const &pos) {
+  if (p_id < 0) {
+    throw std::domain_error("Invalid particle id: " + std::to_string(p_id));
+  }
   if (particle_exists(p_id)) {
     mpi_place_particle(get_particle_node(p_id), p_id, pos);
 
diff --git a/src/python/espressomd/particle_data.pxd b/src/python/espressomd/particle_data.pxd
index b5c3f9d17bc..695047b4090 100644
--- a/src/python/espressomd/particle_data.pxd
+++ b/src/python/espressomd/particle_data.pxd
@@ -74,7 +74,7 @@ cdef extern from "particle_data.hpp":
     # Setter/getter/modifier functions functions
     void prefetch_particle_data(vector[int] ids)
 
-    int place_particle(int part, const Vector3d & p)
+    int place_particle(int part, const Vector3d & p) except +
 
     void set_particle_v(int part, const Vector3d & v)
 
diff --git a/src/python/espressomd/particle_data.pyx b/src/python/espressomd/particle_data.pyx
index de82bebd3ba..28b8afc9b56 100644
--- a/src/python/espressomd/particle_data.pyx
+++ b/src/python/espressomd/particle_data.pyx
@@ -1747,84 +1747,83 @@ cdef class ParticleList:
         # Did we get a dictionary
         if len(args) == 1:
             if hasattr(args[0], "__getitem__"):
-                P = args[0]
+                particles_dict = args[0]
         else:
-            if len(args) == 0 and len(kwargs.keys()) != 0:
-                P = kwargs
+            if len(args) == 0 and len(kwargs) != 0:
+                particles_dict = kwargs
             else:
                 raise ValueError(
                     "add() takes either a dictionary or a bunch of keyword args.")
 
         # Check for presence of pos attribute
-        if "pos" not in P:
+        if "pos" not in particles_dict:
             raise ValueError(
                 "pos attribute must be specified for new particle")
 
-        if len(np.array(P["pos"]).shape) == 2:
-            return self._place_new_particles(P)
+        if len(np.array(particles_dict["pos"]).shape) == 2:
+            return self._place_new_particles(particles_dict)
         else:
-            return self._place_new_particle(P)
+            return self._place_new_particle(particles_dict)
 
-    def _place_new_particle(self, P):
+    def _place_new_particle(self, p_dict):
         # Handling of particle id
-        if "id" not in P:
+        if "id" not in p_dict:
             # Generate particle id
-            P["id"] = get_maximal_particle_id() + 1
+            p_dict["id"] = get_maximal_particle_id() + 1
         else:
-            if particle_exists(P["id"]):
-                raise Exception(f"Particle {P['id']} already exists.")
+            if particle_exists(p_dict["id"]):
+                raise Exception(f"Particle {p_dict['id']} already exists.")
 
         # Prevent setting of contradicting attributes
         IF DIPOLES:
-            if 'dip' in P and 'dipm' in P:
+            if 'dip' in p_dict and 'dipm' in p_dict:
                 raise ValueError("Contradicting attributes: dip and dipm. Setting \
 dip is sufficient as the length of the vector defines the scalar dipole moment.")
             IF ROTATION:
-                if 'dip' in P and 'quat' in P:
+                if 'dip' in p_dict and 'quat' in p_dict:
                     raise ValueError("Contradicting attributes: dip and quat. \
 Setting dip overwrites the rotation of the particle around the dipole axis. \
 Set quat and scalar dipole moment (dipm) instead.")
 
         # The ParticleList can not be used yet, as the particle
         # doesn't yet exist. Hence, the setting of position has to be
-        # done here. the code is from the pos:property of ParticleHandle
+        # done here.
         check_type_or_throw_except(
-            P["pos"], 3, float, "Position must be 3 floats.")
-        if place_particle(P["id"], make_Vector3d(P["pos"])) == -1:
+            p_dict["pos"], 3, float, "Position must be 3 floats.")
+        error_code = place_particle(p_dict["id"], make_Vector3d(p_dict["pos"]))
+        if error_code == -1:
             raise Exception("particle could not be set.")
 
-        # Pos is taken care of
-        del P["pos"]
-        pid = P["id"]
-        del P["id"]
+        # position is taken care of
+        del p_dict["pos"]
+        pid = p_dict.pop("id")
 
-        if P != {}:
-            self.by_id(pid).update(P)
+        if p_dict != {}:
+            self.by_id(pid).update(p_dict)
 
         return self.by_id(pid)
 
-    def _place_new_particles(self, Ps):
+    def _place_new_particles(self, p_list_dict):
         # Check if all entries have the same length
-        n_parts = len(Ps["pos"])
-        if not all(np.shape(Ps[k]) and len(Ps[k]) == n_parts for k in Ps):
+        n_parts = len(p_list_dict["pos"])
+        if not all(np.shape(v) and len(v) ==
+                   n_parts for v in p_list_dict.values()):
             raise ValueError(
                 "When adding several particles at once, all lists of attributes have to have the same size")
 
         # If particle ids haven't been provided, use free ones
         # beyond the highest existing one
-        if not "id" in Ps:
+        if "id" not in p_list_dict:
             first_id = get_maximal_particle_id() + 1
-            Ps["id"] = range(first_id, first_id + n_parts)
+            p_list_dict["id"] = range(first_id, first_id + n_parts)
 
         # Place the particles
         for i in range(n_parts):
-            P = {}
-            for k in Ps:
-                P[k] = Ps[k][i]
-            self._place_new_particle(P)
+            p_dict = {k: v[i] for k, v in p_list_dict.items()}
+            self._place_new_particle(p_dict)
 
         # Return slice of added particles
-        return self.by_ids(Ps["id"])
+        return self.by_ids(p_list_dict["id"])
 
     # Iteration over all existing particles
     def __iter__(self):
diff --git a/testsuite/python/particle.py b/testsuite/python/particle.py
index bf182147949..e812848bcac 100644
--- a/testsuite/python/particle.py
+++ b/testsuite/python/particle.py
@@ -274,11 +274,21 @@ def test_image_box(self):
 
         np.testing.assert_equal(np.copy(p.image_box), [1, 1, 1])
 
-    def test_accessing_invalid_id_raises(self):
+    def test_invalid_particle_ids_exceptions(self):
         self.system.part.clear()
         handle_to_non_existing_particle = self.system.part.by_id(42)
-        with self.assertRaises(RuntimeError):
+        with self.assertRaisesRegex(RuntimeError, "Particle node for id 42 not found"):
             handle_to_non_existing_particle.id
+        p = self.system.part.add(pos=[0., 0., 0.], id=0)
+        with self.assertRaisesRegex(RuntimeError, "Particle node for id 42 not found"):
+            p._id = 42
+            p.node
+        for i in range(1, 10):
+            with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"):
+                p._id = -i
+                p.node
+            with self.assertRaisesRegex(ValueError, f"Invalid particle id: {-i}"):
+                self.system.part.add(pos=[0., 0., 0.], id=-i)
 
     def test_parallel_property_setters(self):
         s = self.system