Make particle id less prominent in the testsuite (#4180)

Partial fix for #4029

Description of changes:
- use `ParticleHandle` objects instead of `system.part[pid]` in the python tests
- improve tests clarity (batch particle creation, utility functions from `itertools` and `tests_common.py`, etc.)

testsuite/python/accumulator_correlator.py

@@ -82,9 +82,9 @@ def check_pickling(self, acc):
     def test_square_distance_componentwise(self):
         s = self.system
         v = np.array([1, 2, 3])
-        s.part.add(id=0, pos=(0, 0, 0), v=v)
+        p = s.part.add(pos=(0, 0, 0), v=v)
 
-        obs = espressomd.observables.ParticlePositions(ids=(0,))
+        obs = espressomd.observables.ParticlePositions(ids=(p.id,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=10, tau_max=2, delta_N=1,
             corr_operation="square_distance_componentwise")
@@ -105,9 +105,9 @@ def test_square_distance_componentwise(self):
     def test_tensor_product(self):
         s = self.system
         v = np.array([1, 2, 3])
-        s.part.add(id=0, pos=(0, 0, 0), v=v)
+        p = s.part.add(pos=(0, 0, 0), v=v)
 
-        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        obs = espressomd.observables.ParticleVelocities(ids=(p.id,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=12, tau_max=2, delta_N=1,
             corr_operation="tensor_product")
@@ -128,9 +128,9 @@ def test_tensor_product(self):
     def test_componentwise_product(self):
         s = self.system
         v = np.array([1, 2, 3])
-        s.part.add(id=0, pos=(0, 0, 0), v=v)
+        p = s.part.add(pos=(0, 0, 0), v=v)
 
-        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        obs = espressomd.observables.ParticleVelocities(ids=(p.id,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=10, tau_max=2, delta_N=1,
             corr_operation="componentwise_product")
@@ -150,9 +150,9 @@ def test_componentwise_product(self):
     def test_scalar_product(self):
         s = self.system
         v = np.array([1, 2, 3])
-        s.part.add(id=0, pos=(0, 0, 0), v=v)
+        p = s.part.add(pos=(0, 0, 0), v=v)
 
-        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        obs = espressomd.observables.ParticleVelocities(ids=(p.id,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=10, tau_max=2, delta_N=1,
             corr_operation="scalar_product")
@@ -172,10 +172,10 @@ def test_scalar_product(self):
     def test_fcs(self):
         s = self.system
         v = np.array([1, 2, 3])
-        s.part.add(id=0, pos=(0, 0, 0), v=v)
+        p = s.part.add(pos=(0, 0, 0), v=v)
 
         w = np.array([3, 2, 1])
-        obs = espressomd.observables.ParticlePositions(ids=(0,))
+        obs = espressomd.observables.ParticlePositions(ids=(p.id,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=10, tau_max=9.9 * self.system.time_step,
             delta_N=1, corr_operation="fcs_acf", args=w)
@@ -202,7 +202,7 @@ def test_fcs(self):
 
     def test_correlator_interface(self):
         # test setters and getters
-        obs = espressomd.observables.ParticleVelocities(ids=(0,))
+        obs = espressomd.observables.ParticleVelocities(ids=(123,))
         acc = espressomd.accumulators.Correlator(
             obs1=obs, tau_lin=10, tau_max=12.0, delta_N=1,
             corr_operation="scalar_product")

testsuite/python/analyze_distribution.py

@@ -32,16 +32,14 @@ def setUpClass(cls):
         # start with a small box
         cls.system.box_l = np.array([box_l, box_l, box_l])
         cls.system.cell_system.set_n_square(use_verlet_lists=False)
-        for p in range(cls.num_part):
-            cls.system.part.add(
-                id=p,
-                pos=np.random.random() * cls.system.box_l)
+        cls.system.part.add(
+            pos=np.outer(np.random.random(cls.num_part), cls.system.box_l))
 
     def calc_min_distribution(self, bins):
         dist = []
-        for i in range(self.num_part):
-            dist.append(np.min([self.system.distance(
-                self.system.part[i], p.pos) for p in self.system.part if p.id != i]))
+        for p1 in self.system.part:
+            dist.append(min(self.system.distance(p1, p2.pos)
+                            for p2 in self.system.part if p1.id != p2.id))
         hist = np.histogram(dist, bins=bins, density=False)[0]
         return hist / (float(np.sum(hist)))
 

testsuite/python/analyze_energy.py

@@ -45,42 +45,44 @@ def setUpClass(cls):
         cls.system.bonded_inter.add(cls.harmonic)
 
     def setUp(self):
+        self.system.part.add(pos=[1, 2, 2], type=0)
+        self.system.part.add(pos=[5, 2, 2], type=0)
+
+    def tearDown(self):
         self.system.part.clear()
-        self.system.part.add(id=0, pos=[1, 2, 2], type=0)
-        self.system.part.add(id=1, pos=[5, 2, 2], type=0)
 
     def test_kinetic(self):
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [5, 2, 2]
-        self.system.part[0].v = [3, 4, 5]
-        self.system.part[1].v = [0, 0, 0]
+        p0, p1 = self.system.part[:]
+        p0.pos = [1, 2, 2]
+        p1.pos = [5, 2, 2]
         # single moving particle
+        p0.v = [3, 4, 5]
+        p1.v = [0, 0, 0]
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 25., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 25., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7)
         # two moving particles
-        self.system.part[1].v = [3, 4, 5]
+        p1.v = [3, 4, 5]
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 50., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7)
-        self.system.part[0].v = [0, 0, 0]
-        self.system.part[1].v = [0, 0, 0]
 
     def test_non_bonded(self):
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [2, 2, 2]
+        p0, p1 = self.system.part[:]
+        p0.pos = [1, 2, 2]
+        p1.pos = [2, 2, 2]
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 1., delta=1e-5)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
         # add another pair of particles
-        self.system.part.add(id=2, pos=[3, 2, 2], type=1)
-        self.system.part.add(id=3, pos=[4, 2, 2], type=1)
+        self.system.part.add(pos=[3, 2, 2], type=1)
+        self.system.part.add(pos=[4, 2, 2], type=1)
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 3., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
@@ -90,78 +92,74 @@ def test_non_bonded(self):
         self.assertAlmostEqual(energy["non_bonded", 0, 0]
                                + energy["non_bonded", 0, 1]
                                + energy["non_bonded", 1, 1], energy["total"], delta=1e-7)
-        self.system.part[2].remove()
-        self.system.part[3].remove()
 
     def test_bonded(self):
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [3, 2, 2]
-        self.system.part[0].v = [0, 0, 0]
-        self.system.part[1].v = [0, 0, 0]
+        p0, p1 = self.system.part[:]
+        p0.pos = [1, 2, 2]
+        p1.pos = [3, 2, 2]
         # single bond
-        self.system.part[0].add_bond((self.harmonic, 1))
+        p0.add_bond((self.harmonic, p1))
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 6, delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 6, delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7)
         # two bonds
-        self.system.part[1].add_bond((self.harmonic, 0))
+        p1.add_bond((self.harmonic, p0))
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 12, delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 12, delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7)
         # bonds deleted
-        self.system.part[0].delete_all_bonds()
-        self.system.part[1].delete_all_bonds()
+        p0.delete_all_bonds()
+        p1.delete_all_bonds()
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0, delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0., delta=1e-7)
 
     def test_all(self):
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [2, 2, 2]
-        self.system.part[0].v = [3, 4, 5]
-        self.system.part[1].v = [3, 4, 5]
+        p0, p1 = self.system.part[:]
+        p0.pos = [1, 2, 2]
+        p1.pos = [2, 2, 2]
+        p0.v = [3, 4, 5]
+        p1.v = [3, 4, 5]
         # single bond
-        self.system.part[0].add_bond((self.harmonic, 1))
+        p0.add_bond((self.harmonic, p1))
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 50. + 3. / 2. + 1., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 3. / 2., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
         # two bonds
-        self.system.part[1].add_bond((self.harmonic, 0))
+        p1.add_bond((self.harmonic, p0))
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 50. + 3 + 1., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 3., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
         # add another pair of particles
-        self.system.part.add(id=2, pos=[1, 5, 5], type=1)
-        self.system.part.add(id=3, pos=[2, 5, 5], type=1)
+        self.system.part.add(pos=[1, 5, 5], type=1)
+        self.system.part.add(pos=[2, 5, 5], type=1)
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(
             energy["total"], 50. + 3 + (1. + 1.), delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 50., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 3., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 1. + 1., delta=1e-7)
-        self.system.part[2].remove()
-        self.system.part[3].remove()
-        self.system.part[0].delete_all_bonds()
 
     @utx.skipIfMissingFeatures(["ELECTROSTATICS", "P3M"])
     def test_electrostatics(self):
 
         from espressomd import electrostatics
 
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [3, 2, 2]
-        self.system.part[0].q = 1
-        self.system.part[1].q = -1
+        p0, p1 = self.system.part[:]
+        p0.pos = [1, 2, 2]
+        p1.pos = [3, 2, 2]
+        p0.q = 1
+        p1.q = -1
         p3m = electrostatics.P3M(prefactor=1.0,
                                  accuracy=9.910945054074526e-08,
                                  mesh=[22, 22, 22],
@@ -180,10 +178,6 @@ def test_electrostatics(self):
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0, delta=1e-7)
         self.assertAlmostEqual(energy["coulomb"], u_p3m, delta=1e-5)
-        self.system.part[0].q = 0
-        self.system.part[1].q = 0
-        self.system.part[0].pos = [1, 2, 2]
-        self.system.part[1].pos = [5, 2, 2]
 
 
 if __name__ == "__main__":

testsuite/python/brownian_dynamics_stats.py

@@ -104,7 +104,7 @@ def test_msd_global_temp(self):
         dt = 0.5
 
         system = self.system
-        p = system.part.add(pos=(0, 0, 0), id=0)
+        p = system.part.add(pos=(0, 0, 0))
         system.time_step = dt
         system.thermostat.set_brownian(kT=kT, gamma=gamma, seed=41)
         system.cell_system.skin = 0.4
@@ -140,7 +140,7 @@ def test_08__noise_correlation(self):
         system.cell_system.skin = 0.1
         kT = 3.2
         system.thermostat.set_brownian(kT=kT, gamma=2.1, seed=17)
-        system.part.add(id=(0, 1), pos=np.zeros((2, 3)))
+        system.part.add(pos=np.zeros((2, 3)))
         steps = int(1e4)
         error_delta = 0.04
         self.check_noise_correlation(kT, steps, error_delta)

testsuite/python/constant_pH_stats.py

@@ -49,11 +49,9 @@ class ReactionEnsembleTest(ut.TestCase):
 
     @classmethod
     def setUpClass(cls):
-        for i in range(0, 2 * cls.N0, 2):
-            cls.system.part.add(id=i, pos=np.random.random(3) *
-                                cls.system.box_l, type=cls.type_A)
-            cls.system.part.add(id=i + 1, pos=np.random.random(3) *
-                                cls.system.box_l, type=cls.type_H)
+        cls.system.part.add(
+            pos=np.random.random((2 * cls.N0, 3)) * cls.system.box_l,
+            type=cls.N0 * [cls.type_A, cls.type_H])
 
         cls.RE.add_reaction(
             gamma=cls.Ka,

testsuite/python/constraint_homogeneous_magnetic_field.py

@@ -70,13 +70,13 @@ def test_add_energy_and_forces(self):
         self.assertEqual(self.S.analysis.energy()["dipolar"], 0.0)
 
         # check dipolar energy when adding dipole moments
-        self.S.part.add(id=0, pos=[0, 0, 0], dip=dip_mom0, rotation=(1, 1, 1))
+        p0 = self.S.part.add(pos=[0, 0, 0], dip=dip_mom0, rotation=(1, 1, 1))
         self.assertEqual(self.S.analysis.energy()["dipolar"],
                          -1.0 * np.dot(H_field, dip_mom0))
-        self.S.part.add(id=1, pos=[1, 1, 1], dip=dip_mom1, rotation=(1, 1, 1))
+        p1 = self.S.part.add(pos=[1, 1, 1], dip=dip_mom1, rotation=(1, 1, 1))
         self.assertEqual(self.S.analysis.energy()["dipolar"],
-                         (-1.0 * np.dot(H_field, dip_mom0)
-                          - 1.0 * np.dot(H_field, dip_mom1)))
+                         -(np.dot(H_field, dip_mom0) +
+                           np.dot(H_field, dip_mom1)))
 
         if espressomd.has_features(["ROTATION"]):
             # check that running the integrator leads to expected torques
@@ -85,11 +85,11 @@ def test_add_energy_and_forces(self):
             torque_expected1 = np.cross(dip_mom1, H_field)
             for i in range(3):
                 self.assertAlmostEqual(
-                    self.S.part[0].torque_lab[i],
+                    p0.torque_lab[i],
                     torque_expected0[i],
                     places=10)
                 self.assertAlmostEqual(
-                    self.S.part[1].torque_lab[i],
+                    p1.torque_lab[i],
                     torque_expected1[i],
                     places=10)