ScaFaCoS maintenance (#3784)

Description of changes:
- fix a bug in the ScaFaCoS interface when calling the `get_params()` method of a ScaFaCos-based actor
- add a test for the ScaFaCos `get_params()` method
- remove unused ScaFaCoS function
- document and test FCS weights attribute
- improve testsuite readability

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();
 

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"

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
@@ -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]
 

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))

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(np.copy(acc.args), w**2)
+        w_squared = np.array([4, 5, 6])**2
+        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
         obs = espressomd.observables.ParticleVelocities(ids=(0,))

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')

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)

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):

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,23 +43,21 @@ 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[:]
         g = GalileiTransform()
 
         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)
 

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)
 

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)
 
 

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):

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):

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)

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)

testsuite/python/test_checkpoint.py

@@ -310,33 +310,33 @@ 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
         self.assertEqual(len(system.constraints),
                          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)

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)