diff --git a/samples/electrophoresis.py b/samples/electrophoresis.py index 907fa07bbe..efecee935c 100644 --- a/samples/electrophoresis.py +++ b/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}") diff --git a/testsuite/python/interactions_bond_angle.py b/testsuite/python/interactions_bond_angle.py index b21d3043bc..ff94dad18f 100644 --- a/testsuite/python/interactions_bond_angle.py +++ b/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) diff --git a/testsuite/scripts/samples/test_electrophoresis.py b/testsuite/scripts/samples/test_electrophoresis.py index c8a2e79dea..b057a843dc 100644 --- a/testsuite/scripts/samples/test_electrophoresis.py +++ b/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)