tests: Improve testing of bond angle forces

samples/electrophoresis.py

@@ -33,9 +33,6 @@
 
 logging.basicConfig(level=logging.INFO)
 
-# Use a fixed int for a deterministic behavior
-np.random.seed()
-
 required_features = ["P3M", "EXTERNAL_FORCES", "WCA"]
 espressomd.assert_features(required_features)
 
@@ -133,7 +130,8 @@
 
 # activate electrostatics
 #############################################################
-p3m = espressomd.electrostatics.P3M(prefactor=1.0, accuracy=1e-2)
+p3m_params = {"prefactor": 1., "accuracy": 1e-2}
+p3m = espressomd.electrostatics.P3M(**p3m_params)
 system.electrostatics.solver = p3m
 
 # apply external force (external electric field)
@@ -160,17 +158,15 @@
     obs=obs_bond_length, delta_N=1)
 system.auto_update_accumulators.add(acc_bond_length)
 
-# data storage for python analysis
-#############################################################
-pos = np.full((N_SAMPLES, N_MONOMERS, 3), np.nan)
+obs_pos = espressomd.observables.ParticlePositions(ids=range(N_MONOMERS))
+acc_pos = espressomd.accumulators.TimeSeries(obs=obs_pos, delta_N=100)
+system.auto_update_accumulators.add(acc_pos)
 
 # sampling Loop
 #############################################################
-for i in range(N_SAMPLES):
-    if i % 100 == 0:
-        logging.info(f"\rsampling: {i:4d}")
-    system.integrator.run(N_INT_STEPS)
-    pos[i] = system.part.by_ids(range(N_MONOMERS)).pos
+for i in range(N_SAMPLES // 100):
+    logging.info(f"\rsampling: {100 * i:4d}")
+    system.integrator.run(100 * N_INT_STEPS)
 
 logging.info("\nsampling finished!\n")
 
@@ -179,6 +175,7 @@
 
 # calculate center of mass (COM) and its velocity
 #############################################################
+pos = acc_pos.time_series()
 COM = pos.sum(axis=1) / N_MONOMERS
 COM_v = (COM[1:] - COM[:-1]) / (N_INT_STEPS * system.time_step)
 
@@ -220,7 +217,7 @@ def exponential(x, lp):
 # ...second by using observables
 
 
-def persistence_length_obs(
+def persistence_length_from_obs(
         acc_bond_length, acc_persistence_angles, exponential):
     bond_lengths_obs = np.array(acc_bond_length.mean())
     sampling_positions_obs = np.insert(
@@ -233,7 +230,7 @@ def persistence_length_obs(
     return sampling_positions_obs, cos_thetas_obs, opt_obs[0]
 
 
-sampling_positions_obs, cos_thetas_obs, persistence_length_obs = persistence_length_obs(
+sampling_positions_obs, cos_thetas_obs, persistence_length_obs = persistence_length_from_obs(
     acc_bond_length, acc_persistence_angles, exponential)
 logging.info(f"persistence length (observables): {persistence_length_obs}")
 

testsuite/python/interactions_bond_angle.py

@@ -93,17 +93,18 @@ def run_test(self, bond_instance, force_func, energy_func, phi0):
 
         d_phi = np.pi / self.N
         for i in range(1, self.N):  # avoid corner cases at phi = 0 or phi = pi
+            phi = i * d_phi
             self.p2.pos = self.start_pos + \
-                self.rotate_vector(self.rel_pos, self.axis, i * d_phi)
+                self.rotate_vector(self.rel_pos, self.axis, phi)
             self.system.integrator.run(recalc_forces=True, steps=0)
 
             # Calculate energies
             E_sim = self.system.analysis.energy()["bonded"]
-            E_ref = energy_func(i * d_phi)
+            E_ref = energy_func(phi)
             # Check that energies match
             np.testing.assert_almost_equal(E_sim, E_ref, decimal=4)
 
-            f_ref = force_func(i * d_phi)
+            f_ref = force_func(phi)
             phi_diff = i * d_phi - phi0
             for p, sign in [[self.p1, -1], [self.p2, +1]]:
                 # Check that force is perpendicular
@@ -123,6 +124,18 @@ def run_test(self, bond_instance, force_func, energy_func, phi0):
                         np.sign(sign * np.dot(force_axis, self.axis)),
                         msg="The force moves particles in the wrong direction")
 
+            # Check that force vectors are correct
+            v1 = np.copy(self.system.distance_vec(self.p1, self.p0))
+            v2 = np.copy(self.system.distance_vec(self.p2, self.p0))
+            d1 = np.linalg.norm(v1)
+            d2 = np.linalg.norm(v2)
+            v1 /= d1
+            v2 /= d2
+            ref_f1 = (f_ref / d1 / np.sin(phi)) * (v1 * np.cos(phi) - v2)
+            ref_f2 = (f_ref / d2 / np.sin(phi)) * (v2 * np.cos(phi) - v1)
+            np.testing.assert_allclose(np.copy(self.p1.f), ref_f1, atol=1E-8)
+            np.testing.assert_allclose(np.copy(self.p2.f), ref_f2, atol=1E-8)
+
             # Total force =0?
             np.testing.assert_allclose(
                 np.sum(np.copy(self.system.part.all().f), 0), [0, 0, 0], atol=1E-12)

testsuite/scripts/samples/test_electrophoresis.py

@@ -19,28 +19,12 @@
 import importlib_wrapper
 import numpy as np
 
-
-def set_seed(code):
-    np_seed = "np.random.seed()"
-    assert np_seed in code
-    return code.replace(np_seed, "np.random.seed(42)")
-
-
-def set_p3m_params(code):
-    p3m = "electrostatics.P3M(prefactor=1.0, accuracy=1e-2)"
-    assert p3m in code
-    return code.replace(
-        p3m, "electrostatics.P3M(prefactor=1, mesh=[16, 16, 16], cao=1, accuracy=1e-2, r_cut=3.7, alpha=0.2, tune=False)")
-
-
-def make_deterministic(code):
-    code = set_seed(code)
-    code = set_p3m_params(code)
-    return code
-
+np.random.seed(seed=42)
 
 sample, skipIfMissingFeatures = importlib_wrapper.configure_and_import(
-    "@SAMPLES_DIR@/electrophoresis.py", N_SAMPLES=400, substitutions=make_deterministic)
+    "@SAMPLES_DIR@/electrophoresis.py", N_SAMPLES=400,
+    p3m_params={"prefactor": 1., "mesh": [14, 14, 14], "cao": 1, "tune": False,
+                "accuracy": 1e-2, "r_cut": 5.3, "alpha": 0.16})
 
 
 @skipIfMissingFeatures
@@ -51,7 +35,7 @@ def test_persistence_length(self):
         # These two values differ due to undersampling, they converge
         # to the same value around N_SAMPLES=1000
         self.assertAlmostEqual(sample.persistence_length, 38.2, delta=1)
-        self.assertAlmostEqual(sample.persistence_length_obs, 48.3, delta=1)
+        self.assertAlmostEqual(sample.persistence_length_obs, 48.6, delta=1)
 
     def test_mobility(self):
         self.assertAlmostEqual(sample.mu, 1.05, delta=0.02)