From 95f2b67ad58d20125eaf3a17712fa4e61885c8f2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 17:27:57 +0200
Subject: [PATCH 1/8] correlators: Document FCS weight parameter

---
 src/python/espressomd/accumulators.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/src/python/espressomd/accumulators.py b/src/python/espressomd/accumulators.py
index ed0f0dd8b5e..2a1f84f2bc2 100644
--- a/src/python/espressomd/accumulators.py
+++ b/src/python/espressomd/accumulators.py
@@ -243,8 +243,11 @@ class Correlator(ScriptInterfaceHelper):
 
     args: :obj:`float` of length 3
         Three floats which are passed as arguments to the correlation
-        function.  Currently it is only used by ``"fcs_acf"``.
-        Other correlation operations will ignore these values.
+        function. Currently it is only used by ``"fcs_acf"``, which
+        will square these values in the core; if you later decide to
+        update these weights with ``obs.args = [...]``, you'll have to
+        provide already squared values! Other correlation operations
+        will ignore these values.
     """
 
     _so_name = "Accumulators::Correlator"

From 476a71f7a66a0cf0475e59040c7aaa84eb5d2e4b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 17:28:28 +0200
Subject: [PATCH 2/8] correlators: Test FCS weight parameter

---
 testsuite/python/correlation.py | 6 ++++++
 1 file changed, 6 insertions(+)

diff --git a/testsuite/python/correlation.py b/testsuite/python/correlation.py
index 65f24226896..0a89e226197 100644
--- a/testsuite/python/correlation.py
+++ b/testsuite/python/correlation.py
@@ -176,6 +176,12 @@ def test_fcs(self):
                 corr[i, 2:],
                 np.exp(-np.linalg.norm(v / w * tau[i])**2), decimal=10)
 
+        # check setter and getter
+        np.testing.assert_array_almost_equal(acc.args, w**2)
+        w_squared = np.array([4, 5, 6])**2
+        acc.args = w_squared
+        np.testing.assert_array_almost_equal(acc.args, w_squared)
+
     def test_correlator_interface(self):
         # test setters and getters
         obs = espressomd.observables.ParticleVelocities(ids=(0,))

From 0cecf5f584ef9366347dc95c0381f03d2c070b32 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 17:29:29 +0200
Subject: [PATCH 3/8] scafacos: Update Cython class

---
 src/python/espressomd/scafacos.pyx | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx
index 715adf7f1b2..d8069898da6 100644
--- a/src/python/espressomd/scafacos.pyx
+++ b/src/python/espressomd/scafacos.pyx
@@ -45,13 +45,13 @@ IF SCAFACOS == 1:
             self._set_params_in_es_core()
 
         def valid_keys(self):
-            tmp = ["method_name", "method_params", "prefactor"]
+            tmp = self.required_keys().copy()
             if not self.dipolar:
-                tmp.append("check_neutrality")
+                tmp.add("check_neutrality")
             return tmp
 
         def required_keys(self):
-            return "method_name", "method_params", "prefactor"
+            return {"method_name", "method_params", "prefactor"}
 
         def validate_params(self):
             return True

From 3be57258bbc580b4b5cb5235ca1b70331a03c2a3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 17:31:11 +0200
Subject: [PATCH 4/8] scafacos: Fix bytestring bug

---
 src/python/espressomd/scafacos.pyx | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/python/espressomd/scafacos.pyx b/src/python/espressomd/scafacos.pyx
index d8069898da6..da6133a81e0 100644
--- a/src/python/espressomd/scafacos.pyx
+++ b/src/python/espressomd/scafacos.pyx
@@ -60,7 +60,7 @@ IF SCAFACOS == 1:
             # Parameters are returned as strings
             # First word is method name, rest are key value pairs
             # which we convert to a dict
-            p = to_str(get_method_and_parameters().split(" "))
+            p = to_str(get_method_and_parameters()).split(" ")
             res = {}
             res["method_name"] = p[0]
 

From c760bc67ee29e45f7019a775d6ab6171a7892821 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 17:54:13 +0200
Subject: [PATCH 5/8] scafacos: Test scafacos tuning and interface

---
 testsuite/python/coulomb_cloud_wall.py | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/testsuite/python/coulomb_cloud_wall.py b/testsuite/python/coulomb_cloud_wall.py
index ab72e11b1b2..29c4e8675a3 100644
--- a/testsuite/python/coulomb_cloud_wall.py
+++ b/testsuite/python/coulomb_cloud_wall.py
@@ -131,6 +131,25 @@ def test_scafacos_p3m(self):
         self.S.integrator.run(0)
         self.compare("scafacos_p3m", energy=True, prefactor=0.5)
 
+    @ut.skipIf(not espressomd.has_features(["SCAFACOS"])
+               or 'p3m' not in scafacos.available_methods(),
+               'Skipping test: missing feature SCAFACOS or p3m method')
+    def test_scafacos_p3m_tuning(self):
+        # check that the tuning function can be called without throwing
+        # an exception or causing an MPI deadlock
+        self.S.actors.add(
+            espressomd.electrostatics.Scafacos(
+                prefactor=0.5,
+                method_name="p3m",
+                method_params={
+                    "p3m_r_cut": -1.5,
+                    "p3m_grid": 64,
+                    "p3m_cao": 7,
+                    "p3m_alpha": 2.70746}))
+        self.S.integrator.run(0)
+        # check the scafacos script interface
+        self.assertEqual(self.S.actors[-1].get_params()['prefactor'], 0.5)
+
     @ut.skipIf(not espressomd.has_features("SCAFACOS")
                or 'p2nfft' not in scafacos.available_methods(),
                'Skipping test: missing feature SCAFACOS or p2nfft method')

From 4420b0c86267ee071185c52a1b4b7ebf82c5e330 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 18:29:20 +0200
Subject: [PATCH 6/8] scafacos: Remove unused function

ScaFaCoS has its own exception mechanism for the case where
particles leave the box, e.g. in 2d or 1d system without walls.
The function is no longer used since 4.1.0.
---
 .../scafacos.cpp                              | 19 -------------------
 1 file changed, 19 deletions(-)

diff --git a/src/core/electrostatics_magnetostatics/scafacos.cpp b/src/core/electrostatics_magnetostatics/scafacos.cpp
index 835bc6bd086..d6db6c1f381 100644
--- a/src/core/electrostatics_magnetostatics/scafacos.cpp
+++ b/src/core/electrostatics_magnetostatics/scafacos.cpp
@@ -180,25 +180,6 @@ double pair_energy(double q1q2, double dist) {
   return 0.;
 }
 
-// Issues a runtime error if positions are outside the box domain
-// This is needed, because the scafacos grid sort produces an mpi deadlock
-// otherwise
-// Returns true if calculations can continue.
-bool check_position_validity(const std::vector<double> &pos) {
-  assert(pos.size() % 3 == 0);
-  for (int i = 0; i < pos.size() / 3; i++) {
-    for (int j = 0; j < 3; j++) {
-      if (pos[3 * i + j] < 0 || pos[3 * i + j] > box_geo.length()[j]) {
-        // Throwing exception rather than runtime error, because continuing will
-        // result in mpi deadlock
-        throw std::runtime_error("Particle position outside the box domain not "
-                                 "allowed for scafacos-based methods.");
-      }
-    }
-  }
-  return true;
-}
-
 void add_long_range_force() {
   particles.update_particle_data();
 

From 939b0f6e5a3cee3e9baf28bd2e8c37ec90f7be70 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 19:34:02 +0200
Subject: [PATCH 7/8] testsuite: Simplify assertions

Use the automatic casting feature of np.testing.assert_* functions.
Use specialized unittest assertions. Simplify numpy reshapes. Use
numpy.linalg.norm() and espressomd.system.System.volume().
---
 testsuite/python/analyze_chains.py          |  2 +-
 testsuite/python/dpd.py                     |  4 ++--
 testsuite/python/field_test.py              |  4 +---
 testsuite/python/galilei.py                 | 12 +++++-------
 testsuite/python/lb_electrohydrodynamics.py |  2 +-
 testsuite/python/lb_poiseuille.py           |  4 ++--
 testsuite/python/lb_poiseuille_cylinder.py  |  2 +-
 testsuite/python/pressure.py                | 12 ++++++------
 testsuite/python/reaction_ensemble.py       |  2 +-
 testsuite/python/scafacos_dipoles_1d_2d.py  |  8 ++++----
 testsuite/python/test_checkpoint.py         | 18 +++++++++---------
 testsuite/python/widom_insertion.py         |  2 +-
 12 files changed, 34 insertions(+), 38 deletions(-)

diff --git a/testsuite/python/analyze_chains.py b/testsuite/python/analyze_chains.py
index 39738edf709..6ebc0ed99b9 100644
--- a/testsuite/python/analyze_chains.py
+++ b/testsuite/python/analyze_chains.py
@@ -95,7 +95,7 @@ def calc_rh(self):
             # this generates indices for all i<j combinations
             ij = np.triu_indices(len(r), k=1)
             r_ij = r[ij[0]] - r[ij[1]]
-            dist = np.sqrt(np.sum(r_ij**2, axis=1))
+            dist = np.linalg.norm(r_ij, axis=1)
             # rh.append(self.num_mono*self.num_mono*0.5/(np.sum(1./dist)))
             # the other way do it, with the proper prefactor of N(N-1)
             rh.append(1. / np.mean(1. / dist))
diff --git a/testsuite/python/dpd.py b/testsuite/python/dpd.py
index f95f22acc31..5fffffefaf3 100644
--- a/testsuite/python/dpd.py
+++ b/testsuite/python/dpd.py
@@ -450,12 +450,12 @@ def calc_stress(dist, vel_diff):
                     vel_diff = pair[1].v - pair[0].v
                     stress += calc_stress(dist, vel_diff)
 
-            stress /= np.prod(np.copy(s.box_l))
+            stress /= s.volume()
 
             dpd_stress = s.analysis.dpd_stress()
 
             dpd_obs = DPDStress()
-            obs_stress = np.array(dpd_obs.calculate()).reshape((3, 3))
+            obs_stress = dpd_obs.calculate()
 
             np.testing.assert_array_almost_equal(np.copy(dpd_stress), stress)
             np.testing.assert_array_almost_equal(np.copy(obs_stress), stress)
diff --git a/testsuite/python/field_test.py b/testsuite/python/field_test.py
index 32fd93f0828..3cf32e30658 100644
--- a/testsuite/python/field_test.py
+++ b/testsuite/python/field_test.py
@@ -65,9 +65,7 @@ def test_gravity(self):
         if espressomd.has_features("VIRTUAL_SITES"):
             self.system.part[0].virtual = True
             self.system.integrator.run(0)
-            np.testing.assert_allclose(
-                np.zeros(3),
-                np.copy(self.system.part[0].f))
+            np.testing.assert_allclose(np.copy(self.system.part[0].f), 0)
 
     @utx.skipIfMissingFeatures("ELECTROSTATICS")
     def test_linear_electric_potential(self):
diff --git a/testsuite/python/galilei.py b/testsuite/python/galilei.py
index 45d7c8c47f8..40d2e97128a 100644
--- a/testsuite/python/galilei.py
+++ b/testsuite/python/galilei.py
@@ -23,13 +23,13 @@
 from espressomd.galilei import GalileiTransform
 
 BOX_L = np.array([10, 20, 30])
-N_PART = 500
 
 
 class Galilei(ut.TestCase):
     system = espressomd.System(box_l=BOX_L)
 
     def setUp(self):
+        N_PART = 500
         self.system.part.add(pos=BOX_L * np.random.random((N_PART, 3)),
                              v=-5. + 10. * np.random.random((N_PART, 3)),
                              f=np.random.random((N_PART, 3)))
@@ -43,15 +43,13 @@ def test_kill_particle_motion(self):
         g = GalileiTransform()
         g.kill_particle_motion()
 
-        np.testing.assert_array_equal(
-            np.copy(self.system.part[:].v), np.zeros((N_PART, 3)))
+        np.testing.assert_array_equal(np.copy(self.system.part[:].v), 0)
 
     def test_kill_particle_forces(self):
         g = GalileiTransform()
         g.kill_particle_forces()
 
-        np.testing.assert_array_equal(
-            np.copy(self.system.part[:].f), np.zeros((N_PART, 3)))
+        np.testing.assert_array_equal(np.copy(self.system.part[:].f), 0)
 
     def test_cms(self):
         parts = self.system.part[:]
@@ -59,7 +57,7 @@ def test_cms(self):
 
         total_mass = np.sum(parts.mass)
         com = np.sum(
-            np.multiply(parts.mass.reshape((N_PART, 1)), parts.pos), axis=0) / total_mass
+            np.multiply(parts.mass.reshape((-1, 1)), parts.pos), axis=0) / total_mass
 
         np.testing.assert_allclose(np.copy(g.system_CMS()), com)
 
@@ -68,7 +66,7 @@ def test_cms_velocity(self):
         g = GalileiTransform()
         total_mass = np.sum(parts.mass)
         com_v = np.sum(
-            np.multiply(parts.mass.reshape((N_PART, 1)), parts.v), axis=0) / total_mass
+            np.multiply(parts.mass.reshape((-1, 1)), parts.v), axis=0) / total_mass
 
         np.testing.assert_allclose(np.copy(g.system_CMS_velocity()), com_v)
 
diff --git a/testsuite/python/lb_electrohydrodynamics.py b/testsuite/python/lb_electrohydrodynamics.py
index 55b132b5338..a716e5a8d5c 100644
--- a/testsuite/python/lb_electrohydrodynamics.py
+++ b/testsuite/python/lb_electrohydrodynamics.py
@@ -64,7 +64,7 @@ def test(self):
         mu_E = np.array(self.params['muE'])
         # Terminal velocity is mu_E minus the momentum the fluid
         # got by accelerating the particle in the beginning.
-        v_term = (1. - 1. / (np.prod(s.box_l) * self.params['dens'])) * mu_E
+        v_term = (1. - 1. / (s.volume() * self.params['dens'])) * mu_E
 
         s.integrator.run(steps=500)
 
diff --git a/testsuite/python/lb_poiseuille.py b/testsuite/python/lb_poiseuille.py
index 598d0bd0763..8bd38bd9511 100644
--- a/testsuite/python/lb_poiseuille.py
+++ b/testsuite/python/lb_poiseuille.py
@@ -120,7 +120,7 @@ def test_profile(self):
                                      self.system.box_l[0] - 2.0 * AGRID,
                                      EXT_FORCE,
                                      VISC * DENS)
-        rmsd = np.sqrt(np.sum(np.square(v_expected - v_measured)))
+        rmsd = np.linalg.norm(v_expected - v_measured)
         self.assertLess(rmsd, 0.015 * AGRID / TIME_STEP)
 
 
@@ -180,7 +180,7 @@ def test_profile(self):
                                      self.system.box_l[0] - 2.0 * AGRID,
                                      EXT_FORCE,
                                      VISC * DENS)
-        rmsd = np.sqrt(np.sum(np.square(v_expected - velocities[:, 1])))
+        rmsd = np.linalg.norm(v_expected - velocities[:, 1])
         self.assertLess(rmsd, 0.02 * AGRID / TIME_STEP)
 
 
diff --git a/testsuite/python/lb_poiseuille_cylinder.py b/testsuite/python/lb_poiseuille_cylinder.py
index c553acf7437..88790bd992c 100644
--- a/testsuite/python/lb_poiseuille_cylinder.py
+++ b/testsuite/python/lb_poiseuille_cylinder.py
@@ -136,7 +136,7 @@ def compare_to_analytical(self):
             BOX_L / 2.0 - 1.0,
             EXT_FORCE,
             VISC * DENS)
-        rmsd = np.sqrt(np.sum(np.square(v_expected - v_measured)))
+        rmsd = np.linalg.norm(v_expected - v_measured)
         self.assertLess(rmsd, 0.02 * AGRID / TIME_STEP)
 
     def prepare_obs(self):
diff --git a/testsuite/python/pressure.py b/testsuite/python/pressure.py
index 0bf600a8e1d..80111af414a 100644
--- a/testsuite/python/pressure.py
+++ b/testsuite/python/pressure.py
@@ -34,7 +34,7 @@
 
 def pressure_tensor_kinetic(vel):
     '''Analytical result for convective pressure tensor'''
-    return np.einsum('ij,ik->jk', vel, vel) / np.prod(np.copy(system.box_l))
+    return np.einsum('ij,ik->jk', vel, vel) / system.volume()
 
 
 def pressure_tensor_bonded(pos):
@@ -43,7 +43,7 @@ def pressure_tensor_bonded(pos):
     for p1, p2 in zip(pos[0::2], pos[1::2]):
         r = p1 - p2
         f = -1.0e4 * r
-        tensor += np.einsum('i,j', f, r) / np.prod(np.copy(system.box_l))
+        tensor += np.einsum('i,j', f, r) / system.volume()
     return tensor
 
 
@@ -56,7 +56,7 @@ def pressure_tensor_nonbonded(particle_pairs):
             r = np.linalg.norm(d)
             r_hat = d / r
             f = (24.0 * 1.0 * (2.0 * 1.0**12 / r**13 - 1.0**6 / r**7)) * r_hat
-            tensor += np.einsum('i,j', f, d) / np.prod(np.copy(system.box_l))
+            tensor += np.einsum('i,j', f, d) / system.volume()
     return tensor
 
 
@@ -68,7 +68,7 @@ def pressure_tensor_nonbonded_inter(particle_pairs):
             d = np.linalg.norm(r)
             r_hat = r / d
             f = (24.0 * 1.0 * (2.0 * 1.0**12 / d**13 - 1.0**6 / d**7)) * r_hat
-            tensor += np.einsum('i,j', f, r) / np.prod(np.copy(system.box_l))
+            tensor += np.einsum('i,j', f, r) / system.volume()
     return tensor
 
 
@@ -80,7 +80,7 @@ def pressure_tensor_nonbonded_intra(particle_pairs):
             d = np.linalg.norm(r)
             r_hat = r / d
             f = (24.0 * 1.0 * (2.0 * 1.0**12 / d**13 - 1.0**6 / d**7)) * r_hat
-            tensor += np.einsum('i,j', f, r) / np.prod(np.copy(system.box_l))
+            tensor += np.einsum('i,j', f, r) / system.volume()
     return tensor
 
 
@@ -265,7 +265,7 @@ def get_anal_pressure_tensor_fene(self, pos_1, pos_2, k, d_r_max, r_0):
         tensor = np.zeros([3, 3])
         vec_r = pos_1 - pos_2
         f = -fene_force2(vec_r, k, d_r_max, r_0)
-        tensor += np.einsum('i,j', f, vec_r) / np.prod(np.copy(system.box_l))
+        tensor += np.einsum('i,j', f, vec_r) / system.volume()
         return tensor
 
     def test_fene(self):
diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py
index 4b4992ec93e..9567ef4b10b 100644
--- a/testsuite/python/reaction_ensemble.py
+++ b/testsuite/python/reaction_ensemble.py
@@ -54,7 +54,7 @@ class ReactionEnsembleTest(ut.TestCase):
     system = espressomd.System(box_l=np.ones(3) * (N0 / c0)**(1.0 / 3.0))
     np.random.seed(69)  # make reaction code fully deterministic
     system.cell_system.skin = 0.4
-    volume = np.prod(system.box_l)  # cuboid box
+    volume = system.volume()
     # Calculate gamma which should lead to target_alpha with given N0 and V
     # Requires N0>>1, otherwise discrete character of N changes the statistics (N>20 should suffice)
     # gamma = prod_i (N_i / V) = alpha^2 N0 / (1-alpha)*V**(-nubar)
diff --git a/testsuite/python/scafacos_dipoles_1d_2d.py b/testsuite/python/scafacos_dipoles_1d_2d.py
index 23ec027c79e..4442236e607 100644
--- a/testsuite/python/scafacos_dipoles_1d_2d.py
+++ b/testsuite/python/scafacos_dipoles_1d_2d.py
@@ -120,10 +120,10 @@ def test_scafacos(self):
 
             # Calculate errors
 
-            err_f = np.sum(np.sqrt(
-                np.sum((s.part[:].f - data[:, 7:10])**2, 1)), 0) / np.sqrt(data.shape[0])
-            err_t = np.sum(np.sqrt(np.sum(
-                (s.part[:].torque_lab - data[:, 10:13])**2, 1)), 0) / np.sqrt(data.shape[0])
+            err_f = np.sum(np.linalg.norm(
+                s.part[:].f - data[:, 7:10], axis=1)) / np.sqrt(data.shape[0])
+            err_t = np.sum(np.linalg.norm(
+                s.part[:].torque_lab - data[:, 10:13], axis=1)) / np.sqrt(data.shape[0])
             err_e = s.analysis.energy()["dipolar"] - ref_E
             print("Energy difference", err_e)
             print("Force difference", err_f)
diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py
index 43f82e0dd38..4faf7b9d1f9 100644
--- a/testsuite/python/test_checkpoint.py
+++ b/testsuite/python/test_checkpoint.py
@@ -310,7 +310,7 @@ def test_lb_boundaries(self):
         self.assertEqual(len(system.lbboundaries), 1)
         np.testing.assert_allclose(
             np.copy(system.lbboundaries[0].velocity), [1e-4, 1e-4, 0])
-        self.assertEqual(type(system.lbboundaries[0].shape), Wall)
+        self.assertIsInstance(system.lbboundaries[0].shape, Wall)
 
     def test_constraints(self):
         from espressomd import constraints
@@ -318,25 +318,25 @@ def test_constraints(self):
                          8 - int(not espressomd.has_features("ELECTROSTATICS")))
         c = system.constraints
 
-        self.assertEqual(type(c[0].shape), Sphere)
+        self.assertIsInstance(c[0].shape, Sphere)
         self.assertAlmostEqual(c[0].shape.radius, 0.1, delta=1E-10)
         self.assertEqual(c[0].particle_type, 17)
 
-        self.assertEqual(type(c[1].shape), Wall)
+        self.assertIsInstance(c[1].shape, Wall)
         np.testing.assert_allclose(np.copy(c[1].shape.normal),
                                    [1. / np.sqrt(3)] * 3)
 
-        self.assertEqual(type(c[2]), constraints.Gravity)
+        self.assertIsInstance(c[2], constraints.Gravity)
         np.testing.assert_allclose(np.copy(c[2].g), [1., 2., 3.])
 
-        self.assertEqual(type(c[3]), constraints.HomogeneousMagneticField)
+        self.assertIsInstance(c[3], constraints.HomogeneousMagneticField)
         np.testing.assert_allclose(np.copy(c[3].H), [1., 2., 3.])
 
-        self.assertEqual(type(c[4]), constraints.HomogeneousFlowField)
+        self.assertIsInstance(c[4], constraints.HomogeneousFlowField)
         np.testing.assert_allclose(np.copy(c[4].u), [1., 2., 3.])
         self.assertAlmostEqual(c[4].gamma, 2.3, delta=1E-10)
 
-        self.assertEqual(type(c[5]), constraints.PotentialField)
+        self.assertIsInstance(c[5], constraints.PotentialField)
         self.assertEqual(c[5].field.shape, (14, 16, 18, 1))
         self.assertAlmostEqual(c[5].default_scale, 1.6, delta=1E-10)
         np.testing.assert_allclose(np.copy(c[5].origin), [-0.5, -0.5, -0.5])
@@ -346,7 +346,7 @@ def test_constraints(self):
         np.testing.assert_allclose(np.copy(c[5].field), np.copy(ref_pot.field),
                                    atol=1e-10)
 
-        self.assertEqual(type(c[6]), constraints.ForceField)
+        self.assertIsInstance(c[6], constraints.ForceField)
         self.assertEqual(c[6].field.shape, (14, 16, 18, 3))
         self.assertAlmostEqual(c[6].default_scale, 1.4, delta=1E-10)
         np.testing.assert_allclose(np.copy(c[6].origin), [-0.5, -0.5, -0.5])
@@ -357,7 +357,7 @@ def test_constraints(self):
                                    atol=1e-10)
 
         if espressomd.has_features("ELECTROSTATICS"):
-            self.assertEqual(type(c[7]), constraints.ElectricPlaneWave)
+            self.assertIsInstance(c[7], constraints.ElectricPlaneWave)
             np.testing.assert_allclose(np.copy(c[7].E0), [1., -2., 3.])
             np.testing.assert_allclose(np.copy(c[7].k), [-.1, .2, .3])
             self.assertAlmostEqual(c[7].omega, 5., delta=1E-10)
diff --git a/testsuite/python/widom_insertion.py b/testsuite/python/widom_insertion.py
index 4730426c285..4947e37ffb9 100644
--- a/testsuite/python/widom_insertion.py
+++ b/testsuite/python/widom_insertion.py
@@ -69,7 +69,7 @@ class WidomInsertionTest(ut.TestCase):
     system.cell_system.set_n_square()
     np.random.seed(69)  # make reaction code fully deterministic
     system.cell_system.skin = 0.4
-    volume = np.prod(system.box_l)  # cuboid box
+    volume = system.volume()
 
     Widom = reaction_ensemble.WidomInsertion(
         temperature=TEMPERATURE, seed=1)

From 4bd11e1273c1577b3f755595646a50b7d949062e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Tue, 30 Jun 2020 20:43:06 +0200
Subject: [PATCH 8/8] testsuite: Add missing np.copy()

---
 testsuite/python/correlation.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/testsuite/python/correlation.py b/testsuite/python/correlation.py
index 0a89e226197..4b3e8f527ff 100644
--- a/testsuite/python/correlation.py
+++ b/testsuite/python/correlation.py
@@ -177,10 +177,10 @@ def test_fcs(self):
                 np.exp(-np.linalg.norm(v / w * tau[i])**2), decimal=10)
 
         # check setter and getter
-        np.testing.assert_array_almost_equal(acc.args, w**2)
+        np.testing.assert_array_almost_equal(np.copy(acc.args), w**2)
         w_squared = np.array([4, 5, 6])**2
         acc.args = w_squared
-        np.testing.assert_array_almost_equal(acc.args, w_squared)
+        np.testing.assert_array_almost_equal(np.copy(acc.args), w_squared)
 
     def test_correlator_interface(self):
         # test setters and getters