Update type map when clearing the particle list (#4648)

Fixes #4644
Fixes #4645
Pre-requisite to #4629

Description of changes:
- bugfix: the type map is now properly updated when clearing particles
- add missing feature guards in the testsuite
- add subtests to provide more context when a test fails

src/core/particle_node.cpp

@@ -322,6 +322,12 @@ int get_particle_node(int p_id) {
 
 void clear_particle_node() { particle_node.clear(); }
 
+static void clear_particle_type_map() {
+  for (auto &kv : ::particle_type_map) {
+    kv.second.clear();
+  }
+}
+
 /**
  * @brief Calculate the largest particle id.
  * Traversing the @ref particle_node to find the largest particle id
@@ -404,6 +410,7 @@ REGISTER_CALLBACK(mpi_remove_all_particles_local)
 void remove_all_particles() {
   mpi_call_all(mpi_remove_all_particles_local);
   clear_particle_node();
+  clear_particle_type_map();
 }
 
 void remove_particle(int p_id) {

testsuite/python/ek_eof_one_species.py

@@ -421,7 +421,9 @@ def test_vtk(self):
             ek.write_vtk_velocity(str(path_vtk_velocity))
             ek.write_vtk_potential(str(path_vtk_potential))
             ek.write_vtk_density(str(path_vtk_lbdensity))
-            ek.write_vtk_lbforce(str(path_vtk_lbforce))
+            if espressomd.has_features('EK_DEBUG') or espressomd.has_features(
+                    'VIRTUAL_SITES_INERTIALESS_TRACERS'):
+                ek.write_vtk_lbforce(str(path_vtk_lbforce))
             counterions.write_vtk_density(str(path_vtk_density))
             counterions.write_vtk_flux(str(path_vtk_flux))
             if espressomd.has_features('EK_DEBUG'):
@@ -436,7 +438,9 @@ def test_vtk(self):
             vtk_potential = get_vtk(path_vtk_potential, "potential", grid_dims)
             vtk_lbdensity = get_vtk(
                 path_vtk_lbdensity, "density_lb", grid_dims)
-            get_vtk(path_vtk_lbforce, "lbforce", grid_dims + [3])
+            if espressomd.has_features('EK_DEBUG') or espressomd.has_features(
+                    'VIRTUAL_SITES_INERTIALESS_TRACERS'):
+                get_vtk(path_vtk_lbforce, "lbforce", grid_dims + [3])
             vtk_density = get_vtk(path_vtk_density, "density_1", grid_dims)
             vtk_flux = get_vtk(path_vtk_flux, "flux_1", grid_dims + [3])
             if espressomd.has_features('EK_DEBUG'):

testsuite/python/particle.py

@@ -290,13 +290,9 @@ def test_particle_selection(self):
         np.testing.assert_equal(sorted(res.id), np.arange(0, 16, 2, dtype=int))
 
     def test_image_box(self):
-        s = self.system
-        s.part.clear()
-
-        pos = 1.5 * s.box_l
-
-        p = s.part.add(pos=pos)
-
+        self.system.part.clear()
+        pos = 1.5 * self.system.box_l
+        p = self.system.part.add(pos=pos)
         np.testing.assert_equal(np.copy(p.image_box), [1, 1, 1])
 
     def test_particle_numbering(self):

testsuite/python/reaction_bookkeeping.py

@@ -94,10 +94,8 @@ def test_reaction_bookeeping(self):
         for _ in range(50):
             self.widom.calculate_particle_insertion_potential_energy(
                 reaction_id=0)
-            charge = 0.
-            for p in self.system.part:
-                charge += p.q
-            self.assertEqual(charge, 0.)
+            total_charge = sum(self.system.part.all().q)
+            self.assertEqual(total_charge, 0.)
 
 
 if __name__ == "__main__":

testsuite/python/reaction_methods_interface.py

@@ -45,6 +45,9 @@ def check_reaction_parameters(reactions, parameters):
                     else:
                         self.assertEqual(getattr(reaction, key), params[key])
 
+        def count_by_type(types):
+            return [self.system.number_of_particles(type=x) for x in types]
+
         reaction_forward = {
             'gamma': gamma,
             'reactant_types': [5],
@@ -156,11 +159,15 @@ def check_reaction_parameters(reactions, parameters):
             potential_energy = method.calculate_particle_insertion_potential_energy(
                 reaction_id=0)
             self.assertEqual(potential_energy, 0.)
+        self.assertEqual(count_by_type([5, 2, 3, 0]), [1, 1, 1, 0])
         method.delete_particle(p_id=p3.id)
+        self.assertEqual(count_by_type([5, 2, 3, 0]), [1, 1, 0, 0])
         self.assertEqual(len(self.system.part), 2)
-        method.delete_particle(p_id=p1.id)
+        p1.remove()
+        self.assertEqual(count_by_type([5, 2, 3, 0]), [0, 1, 0, 0])
         self.assertEqual(len(self.system.part), 1)
         self.system.part.clear()
+        self.assertEqual(count_by_type([5, 2, 3, 0]), [0, 0, 0, 0])
 
         # check reaction deletion
         method.delete_reaction(reaction_id=0)
@@ -171,21 +178,25 @@ def test_interface(self):
                   'exclusion_radius_per_type': {1: 0.1}}
 
         # reaction ensemble
-        method = espressomd.reaction_methods.ReactionEnsemble(
-            kT=1.4, seed=12, search_algorithm="order_n", **params)
-        self.check_interface(method, kT=1.4, gamma=1.2,
-                             search_algorithm="order_n", **params)
-
-        # constant pH ensemble
-        method = espressomd.reaction_methods.ConstantpHEnsemble(
-            kT=1.5, seed=14, search_algorithm="parallel", constant_pH=10., **params)
-        self.check_interface(method, kT=1.5, gamma=1.2,
-                             search_algorithm="parallel", **params)
-
-        # Widom insertion
-        method = espressomd.reaction_methods.WidomInsertion(kT=1.6, seed=16)
-        self.check_interface(method, kT=1.6, gamma=1., exclusion_range=0.,
-                             exclusion_radius_per_type={}, search_algorithm=None)
+        with self.subTest(msg="reaction ensemble"):
+            method = espressomd.reaction_methods.ReactionEnsemble(
+                kT=1.4, seed=12, search_algorithm="order_n", **params)
+            self.check_interface(method, kT=1.4, gamma=1.2,
+                                 search_algorithm="order_n", **params)
+
+        with self.subTest(msg="constant pH ensemble"):
+            method = espressomd.reaction_methods.ConstantpHEnsemble(
+                kT=1.5, seed=14, search_algorithm="parallel", constant_pH=10.,
+                **params)
+            self.check_interface(method, kT=1.5, gamma=1.2,
+                                 search_algorithm="parallel", **params)
+
+        with self.subTest(msg="Widom insertion"):
+            method = espressomd.reaction_methods.WidomInsertion(
+                kT=1.6, seed=16)
+            self.check_interface(method, kT=1.6, gamma=1., exclusion_range=0.,
+                                 exclusion_radius_per_type={},
+                                 search_algorithm=None)
 
     def test_exceptions(self):
         single_reaction_params = {

testsuite/python/virtual_sites_relative.py

@@ -32,6 +32,8 @@ class VirtualSites(ut.TestCase):
 
     def setUp(self):
         self.system.box_l = [10.0, 10.0, 10.0]
+        self.system.cell_system.set_regular_decomposition(
+            use_verlet_lists=True)
 
     def tearDown(self):
         self.system.part.clear()
@@ -266,10 +268,7 @@ def run_test_lj(self):
         # Setup
         system.box_l = [l, l, l]
         system.min_global_cut = 0.501
-        system.part.clear()
-
         system.time_step = 0.01
-        system.thermostat.turn_off()
 
         # Dumbbells consist of 2 virtual lj spheres + central particle
         # w/o interactions. For n spheres, n/2 dumbbells.
@@ -300,8 +299,12 @@ def run_test_lj(self):
         # Remove overlap
         system.integrator.set_steepest_descent(
             f_max=0, gamma=0.1, max_displacement=0.1)
-        while system.analysis.energy()["total"] > 10 * n:
+        n_loops = 0
+        n_max = 10
+        while system.analysis.energy()["total"] > 10 * n and n_loops < n_max:
             system.integrator.run(20)
+            n_loops += 1
+        assert n_loops < n_max, "Steepest descent didn't converge"
         # Integrate
         system.integrator.set_vv()
         for i in range(10):
@@ -312,7 +315,7 @@ def run_test_lj(self):
             # Constant energy to get rid of thermostat forces in the
             # verification
             system.integrator.run(2)
-            # Check the virtual sites config,pos and vel of the lj spheres
+            # Check the virtual sites config, pos and vel of the lj spheres
             for j in range(int(n / 2)):
                 self.verify_vs(system.part.by_id(3 * j + 1))
                 self.verify_vs(system.part.by_id(3 * j + 2))
@@ -323,26 +326,33 @@ def run_test_lj(self):
             tests_common.verify_lj_forces(system, 1E-10, 3 * np.arange(n // 2))
 
         # Test applying changes
-        enegry_pre_change = system.analysis.energy()['total']
+        energy_pre_change = system.analysis.energy()['total']
         pressure_pre_change = system.analysis.pressure()['total']
         p0 = system.part.by_id(0)
         p0.pos = p0.pos + (2.2, -1.4, 4.2)
-        enegry_post_change = system.analysis.energy()['total']
+        energy_post_change = system.analysis.energy()['total']
         pressure_post_change = system.analysis.pressure()['total']
-        self.assertNotAlmostEqual(enegry_pre_change, enegry_post_change)
+        self.assertNotAlmostEqual(energy_pre_change, energy_post_change)
         self.assertNotAlmostEqual(pressure_pre_change, pressure_post_change)
 
     def test_lj(self):
         """Run LJ fluid test for different cell systems."""
-        system = self.system
 
-        system.cell_system.skin = 0.4
-        system.cell_system.set_n_square(use_verlet_lists=True)
-        self.run_test_lj()
-        system.cell_system.set_regular_decomposition(use_verlet_lists=True)
-        self.run_test_lj()
-        system.cell_system.set_regular_decomposition(use_verlet_lists=False)
-        self.run_test_lj()
+        self.system.cell_system.skin = 0.4
+        with self.subTest(msg='N-square cell system with Verlet lists'):
+            self.system.cell_system.set_n_square(use_verlet_lists=True)
+            self.run_test_lj()
+            self.tearDown()
+        with self.subTest(msg='regular decomposition cell system with Verlet lists'):
+            self.system.cell_system.set_regular_decomposition(
+                use_verlet_lists=True)
+            self.run_test_lj()
+            self.tearDown()
+        with self.subTest(msg='regular decomposition cell system without Verlet lists'):
+            self.system.cell_system.set_regular_decomposition(
+                use_verlet_lists=False)
+            self.run_test_lj()
+            self.tearDown()
 
     @utx.skipIfMissingFeatures("EXTERNAL_FORCES")
     def test_zz_pressure_tensor(self):