From 644e8e02ebeea5756f7171c204f990b91cd80b3e Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= <jgrad@icp.uni-stuttgart.de>
Date: Fri, 7 Feb 2025 21:41:24 +0100
Subject: [PATCH] Bump Python requirements
---
CMakeLists.txt | 2 +-
doc/sphinx/appendix.rst | 4 +-
doc/sphinx/inter_non-bonded.rst | 2 +-
.../langevin_dynamics/langevin_dynamics.ipynb | 2 +-
testsuite/python/constraint_shape_based.py | 19 +-----
testsuite/python/integrator_langevin_stats.py | 3 +-
testsuite/python/interactions_bonded.py | 2 +
testsuite/python/interactions_non-bonded.py | 66 ++++++++++++-------
testsuite/python/lb_pressure_tensor.py | 5 +-
testsuite/python/p3m_madelung.py | 2 +-
testsuite/python/tests_common.py | 14 ++--
testsuite/python/thermalized_bond.py | 10 ++-
testsuite/python/thermostats_common.py | 3 +-
testsuite/python/widom_insertion.py | 3 +-
14 files changed, 77 insertions(+), 60 deletions(-)
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 6053603862..bd7a111f53 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -276,7 +276,7 @@ endif()
# Python interpreter and Cython interface library
if(ESPRESSO_BUILD_WITH_PYTHON)
- find_package(Python 3.10 REQUIRED COMPONENTS Interpreter Development NumPy)
+ find_package(Python 3.11 REQUIRED COMPONENTS Interpreter Development NumPy)
find_package(Cython 0.29.28...<3.0.12 REQUIRED)
find_program(IPYTHON_EXECUTABLE NAMES jupyter ipython3 ipython)
endif()
diff --git a/doc/sphinx/appendix.rst b/doc/sphinx/appendix.rst
index 43463b8e48..68cce711ea 100644
--- a/doc/sphinx/appendix.rst
+++ b/doc/sphinx/appendix.rst
@@ -510,7 +510,7 @@ prefactor, :math:`q` the electric charge and :math:`a` the lattice constant.
Likewise, the pressure per ion can be derived as :math:`MC\frac{q}{aV}`
with :math:`V` the simulation box volume. For details, see :cite:`ciftja19a`.
-For an infinite 2D or 3D NaCl crystal lattice, the Madelung constant can be
+For an infinite 2D or 3D ionic crystal lattice, the Madelung constant can be
obtained in a numerical simulation with the Evjen method :cite:`evjen32a` or
the Ewald method :cite:`ewald21a`.
@@ -527,6 +527,6 @@ with :math:`M` the orientation-dependent 1D Madelung constant,
:math:`C` the magnetostatics prefactor, :math:`\mu` the dipole moment and
:math:`a` the lattice constant :cite:`batle20a`.
-For an infinite 2D or 3D NaCl crystal lattice, the Madelung constant for
+For an infinite 2D or 3D magnetic crystal lattice, the Madelung constant for
the maximal energy and minimal energy dipole orientation can be estimated
in a numerical simulation :cite:`batle20a`.
diff --git a/doc/sphinx/inter_non-bonded.rst b/doc/sphinx/inter_non-bonded.rst
index 911ef599a6..2fec7a9231 100644
--- a/doc/sphinx/inter_non-bonded.rst
+++ b/doc/sphinx/inter_non-bonded.rst
@@ -192,7 +192,7 @@ the normal LJ potential is recovered for :math:`b_1=b_2=4`,
The optional ``LJGEN_SOFTCORE`` feature activates a softcore version of
the potential, where the following transformations apply:
:math:`\epsilon \rightarrow \lambda \epsilon` and
-:math:`r-r_\mathrm{off} \rightarrow \sqrt{(r-r_\mathrm{off})^2 +
+:math:`(r-r_\mathrm{off}) \rightarrow \sqrt{(r-r_\mathrm{off})^2 +
(1-\lambda) \delta \sigma^2}`. :math:`\lambda` allows to tune the strength of the
interaction, while :math:`\delta` varies how smoothly the potential goes to zero as
:math:`\lambda\rightarrow 0`. Such a feature allows one to perform
diff --git a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb
index 6c602341e1..65b47e525f 100644
--- a/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb
+++ b/doc/tutorials/langevin_dynamics/langevin_dynamics.ipynb
@@ -459,7 +459,7 @@
"source": [
"We find that the velocity-autocorrelation function quickly decays towards zero. However, owing to the relatively short overall sampling time, only the first part of the correlation function is well-sampled and a lot of noise is found in the tail of the autocorrelation function already early on. The obvious solution would be to increase the sampling time and in a production setting one would definitely have to do so in order to smoothly resolve at least several relaxation times. However, depending on a system's characteristics, under certain conditions it might still be necessary to replace a noisy long-time tail with an analytical expression, fitted to the short-time part of the autocorrelation function (again over at least several decay times; typically one would smoothly transition between numerical short-time data and the analytical tail-fit).\n",
"\n",
- "A perfect smoothly sampled autocorrelation function could be integrated numerically, using e.g. [<tt>numpy.trapz</tt>](https://numpy.org/doc/stable/reference/generated/numpy.trapz.html).\n",
+ "A perfect smoothly sampled autocorrelation function could be integrated numerically, using e.g. [<tt>scipy.integrate.trapezoid</tt>](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.trapezoid.html).\n",
"Here, however, we will use the initial part of the velocity-autocorrelation function to obtain a fully analytic data representation. For a Brownian particle the velocity-autocorrelation is expected to follow a simple exponential decay.\n",
"\n",
"Write a Python-function for the exponential decay. Fit your decay-function to the (short-time) correlation data and create a plot to visually verify that the analytic fits are indeed good representations of the data (the exponential decay should be a perfect match in the smooth region of the correlation function). You can copy and modify the plot script given above.\n",
diff --git a/testsuite/python/constraint_shape_based.py b/testsuite/python/constraint_shape_based.py
index d9b9379fa3..c62d1b081e 100644
--- a/testsuite/python/constraint_shape_based.py
+++ b/testsuite/python/constraint_shape_based.py
@@ -426,7 +426,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
-1.0 * outer_cylinder_wall.total_force()[1],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -445,7 +444,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
outer_cylinder_wall.total_normal_force(),
2 *
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -545,7 +543,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
-1.0 * outer_cylinder_constraint.total_force()[1],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -563,8 +560,8 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(outer_cylinder_constraint.total_force()[2], 0.0)
self.assertAlmostEqual(outer_cylinder_constraint.total_normal_force(),
2 * tests_common.lj_force(
- espressomd, cutoff=2.0, offset=0.,
- epsilon=1.0, sigma=1.0, r=dist_part2))
+ cutoff=2.0, offset=0., epsilon=1.0,
+ sigma=1.0, r=dist_part2))
# Reset
system.part.clear()
@@ -662,7 +659,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
p.f[1],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -672,7 +668,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
p.f[2],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.5,
@@ -684,7 +679,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
-1.0 * wall_xz.total_force()[1],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -694,7 +688,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
-1.0 * wall_xy.total_force()[2],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.5,
@@ -706,7 +699,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
wall_xy.total_normal_force(),
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.5,
@@ -818,7 +810,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
-p.f[2],
tests_common.lj_force(
- espressomd,
cutoff=2.,
offset=0.,
epsilon=1.,
@@ -828,7 +819,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
rhomboid_constraint.total_normal_force(),
tests_common.lj_force(
- espressomd,
cutoff=2.,
offset=0.,
epsilon=1.,
@@ -918,7 +908,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
rhomboid_constraint.total_normal_force(),
tests_common.lj_force(
- espressomd,
cutoff=2.,
offset=0.,
epsilon=1.,
@@ -933,7 +922,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
rhomboid_constraint.total_normal_force(),
tests_common.lj_force(
- espressomd,
cutoff=2.,
offset=0.,
epsilon=1.,
@@ -980,7 +968,6 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(
torus_wall.total_force()[1],
tests_common.lj_force(
- espressomd,
cutoff=2.0,
offset=0.,
epsilon=1.0,
@@ -996,7 +983,7 @@ class ShapeBasedConstraintTest(ut.TestCase):
self.assertAlmostEqual(torus_wall.total_force()[1], 0.0)
self.assertAlmostEqual(torus_wall.total_normal_force(), 2 * tests_common.lj_force(
- espressomd, cutoff=2.0, offset=0., epsilon=1.0, sigma=1.0,
+ cutoff=2.0, offset=0., epsilon=1.0, sigma=1.0,
r=radius - tube_radius - part_offset))
# Test the geometry of the shape directly
diff --git a/testsuite/python/integrator_langevin_stats.py b/testsuite/python/integrator_langevin_stats.py
index 882be1f709..a3220045b5 100644
--- a/testsuite/python/integrator_langevin_stats.py
+++ b/testsuite/python/integrator_langevin_stats.py
@@ -19,6 +19,7 @@
import unittest as ut
import unittest_decorators as utx
import numpy as np
+import scipy.integrate
import espressomd
import espressomd.accumulators
@@ -123,7 +124,7 @@ class LangevinThermostat(ut.TestCase, thermostats_common.ThermostatsCommon):
# Integrate with trapezoidal rule
for i in range(3):
- I = np.trapz(acf[:, p.id, i], tau)
+ I = scipy.integrate.trapezoid(acf[:, p.id, i], tau)
ratio = I / (kT / gamma[i])
self.assertAlmostEqual(ratio, 1., delta=0.07)
diff --git a/testsuite/python/interactions_bonded.py b/testsuite/python/interactions_bonded.py
index 6f8adb031f..03ea079ffb 100644
--- a/testsuite/python/interactions_bonded.py
+++ b/testsuite/python/interactions_bonded.py
@@ -256,6 +256,8 @@ class InteractionsBondedTest(ut.TestCase):
p2.pos = p1.pos + self.axis * cutoff * 1.01
with self.assertRaisesRegex(Exception, r"while calling method integrate\(\)"):
self.system.integrator.run(recalc_forces=True, steps=0)
+ with self.assertRaisesRegex(Exception, r"while calling method calculate_energy\(\)"):
+ self.system.analysis.energy()["total"]
if test_same_pos_exception:
p2.pos = p1.pos
with self.assertRaises(Exception):
diff --git a/testsuite/python/interactions_non-bonded.py b/testsuite/python/interactions_non-bonded.py
index ba55bfced4..aaa835c543 100644
--- a/testsuite/python/interactions_non-bonded.py
+++ b/testsuite/python/interactions_non-bonded.py
@@ -47,12 +47,12 @@ def lj_cos_potential(r, epsilon, sigma, cutoff, offset):
return V
-def lj_cos_force(espressomd, r, epsilon, sigma, cutoff, offset):
+def lj_cos_force(r, epsilon, sigma, cutoff, offset):
f = 0.
r_min = offset + np.power(2., 1. / 6.) * sigma
r_cut = cutoff + offset
if r < r_min:
- f = tests_common.lj_force(espressomd, r, epsilon=epsilon, sigma=sigma,
+ f = tests_common.lj_force(r, epsilon=epsilon, sigma=sigma,
cutoff=cutoff, offset=offset)
elif r < r_cut:
alpha = np.pi / \
@@ -78,12 +78,12 @@ def lj_cos2_potential(r, epsilon, sigma, offset, width):
return V
-def lj_cos2_force(espressomd, r, epsilon, sigma, offset, width):
+def lj_cos2_force(r, epsilon, sigma, offset, width):
f = 0.
r_min = offset + np.power(2., 1. / 6.) * sigma
r_cut = r_min + width
if r < r_min:
- f = tests_common.lj_force(espressomd, r, epsilon=epsilon, sigma=sigma,
+ f = tests_common.lj_force(r, epsilon=epsilon, sigma=sigma,
cutoff=r_cut, offset=offset)
elif r < r_cut:
f = - np.pi * epsilon * \
@@ -304,8 +304,7 @@ class InteractionsNonBondedTest(ut.TestCase):
params,
force_kernel=tests_common.lj_generic_force,
energy_kernel=tests_common.lj_generic_potential,
- n_steps=231,
- force_kernel_needs_espressomd=True)
+ n_steps=231)
params["shift"] = "auto"
obj = espressomd.interactions.GenericLennardJonesInteraction(**params)
@@ -317,6 +316,31 @@ class InteractionsNonBondedTest(ut.TestCase):
obj = espressomd.interactions.GenericLennardJonesInteraction(**params)
self.assertEqual(obj.shift, 0.)
+ def test_lj_generic_reference_potential(self):
+
+ class HasFeatures:
+ """Mock class to override ``espressomd.has_features()``."""
+
+ def __init__(self, feature_name, has_feature):
+ self.feature_name = feature_name
+ self.has_feature = has_feature
+
+ def has_features(self, feature_name):
+ assert isinstance(feature_name, str)
+ result = espressomd.has_features(feature_name)
+ if feature_name == self.feature_name:
+ result = self.has_feature
+ return result
+
+ params = {
+ "r": 1., "epsilon": 1., "sigma": 1., "cutoff": 5., "offset": 2.}
+ lj_gen_force_with_softcore = tests_common.lj_generic_force(
+ espressomd=HasFeatures("LJGEN_SOFTCORE", True), **params)
+ lj_gen_force_without_softcore = tests_common.lj_generic_force(
+ espressomd=HasFeatures("LJGEN_SOFTCORE", False), **params)
+ self.assertAlmostEqual(lj_gen_force_with_softcore, -24., delta=1e-7)
+ self.assertAlmostEqual(lj_gen_force_without_softcore, 24., delta=1e-7)
+
# Test WCA Potential
@utx.skipIfMissingFeatures("WCA")
def test_wca(self):
@@ -333,8 +357,7 @@ class InteractionsNonBondedTest(ut.TestCase):
espressomd, r, epsilon=epsilon, sigma=sigma, cutoff=wca_cutoff),
energy_kernel=lambda r, epsilon, sigma: tests_common.lj_generic_potential(
r, epsilon=epsilon, sigma=sigma, cutoff=wca_cutoff, shift=4. * wca_shift),
- n_steps=231,
- force_kernel_needs_espressomd=True)
+ n_steps=231)
# Test Generic Lennard-Jones Softcore Potential
@utx.skipIfMissingFeatures("LJGEN_SOFTCORE")
@@ -354,8 +377,7 @@ class InteractionsNonBondedTest(ut.TestCase):
"lam": 0.34},
force_kernel=tests_common.lj_generic_force,
energy_kernel=tests_common.lj_generic_potential,
- n_steps=231,
- force_kernel_needs_espressomd=True)
+ n_steps=231)
# Test Lennard-Jones Potential
@utx.skipIfMissingFeatures("LENNARD_JONES")
@@ -368,8 +390,7 @@ class InteractionsNonBondedTest(ut.TestCase):
"shift": 0.92},
force_kernel=tests_common.lj_force,
energy_kernel=tests_common.lj_potential,
- n_steps=113,
- force_kernel_needs_espressomd=True)
+ n_steps=113)
# Test Lennard-Jones Cosine Potential
@utx.skipIfMissingFeatures("LJCOS")
@@ -382,8 +403,7 @@ class InteractionsNonBondedTest(ut.TestCase):
"offset": 0.223},
force_kernel=lj_cos_force,
energy_kernel=lj_cos_potential,
- n_steps=175,
- force_kernel_needs_espressomd=True)
+ n_steps=175)
# Test Lennard-Jones Cosine^2 Potential
@utx.skipIfMissingFeatures("LJCOS2")
@@ -396,8 +416,7 @@ class InteractionsNonBondedTest(ut.TestCase):
"offset": 0.321},
force_kernel=lj_cos2_force,
energy_kernel=lj_cos2_potential,
- n_steps=267,
- force_kernel_needs_espressomd=True)
+ n_steps=267)
# Test Smooth-step Potential
@utx.skipIfMissingFeatures("SMOOTH_STEP")
@@ -456,8 +475,7 @@ class InteractionsNonBondedTest(ut.TestCase):
"shift": 0.133},
force_kernel=buckingham_force,
energy_kernel=buckingham_potential,
- n_steps=226,
- force_kernel_remove_shift=False)
+ n_steps=226)
# Test Soft-sphere Potential
@utx.skipIfMissingFeatures("SOFT_SPHERE")
@@ -640,9 +658,7 @@ class InteractionsNonBondedTest(ut.TestCase):
self.assertEqual(self.system.analysis.energy()["non_bonded"], 0.0)
def run_test(self, name, parameters, force_kernel,
- energy_kernel, n_steps, n_initial_steps=0,
- force_kernel_needs_espressomd=False,
- force_kernel_remove_shift=True):
+ energy_kernel, n_steps, n_initial_steps=0):
getattr(self.system.non_bonded_inter[0, 0], name).set_params(
**parameters)
@@ -650,9 +666,11 @@ class InteractionsNonBondedTest(ut.TestCase):
p1.pos = p0.pos + self.step * n_initial_steps
force_parameters = parameters.copy()
- if "shift" in force_parameters and force_kernel_remove_shift:
+ energy_parameters = parameters.copy()
+ force_kernel_varnames = force_kernel.__code__.co_varnames
+ if "shift" in parameters and "shift" not in force_kernel_varnames:
del force_parameters["shift"]
- if force_kernel_needs_espressomd:
+ if "espressomd" in force_kernel_varnames:
force_parameters["espressomd"] = espressomd
for _ in range(n_steps):
@@ -662,7 +680,7 @@ class InteractionsNonBondedTest(ut.TestCase):
# Calculate energies
E_sim = self.system.analysis.energy()["non_bonded"]
- E_ref = energy_kernel(r=d, **parameters)
+ E_ref = energy_kernel(r=d, **energy_parameters)
# Calculate forces
f0_sim = np.copy(p0.f)
diff --git a/testsuite/python/lb_pressure_tensor.py b/testsuite/python/lb_pressure_tensor.py
index 263ff9295f..f8e2cb4ad6 100644
--- a/testsuite/python/lb_pressure_tensor.py
+++ b/testsuite/python/lb_pressure_tensor.py
@@ -19,10 +19,11 @@
import unittest as ut
import unittest_decorators as utx
import numpy as np
+# import scipy.optimize
+# import scipy.integrate
import espressomd
import espressomd.lb
-# import scipy.optimize
N_CELLS = 12
@@ -195,7 +196,7 @@ class TestLBPressureTensorGPU(TestLBPressureTensor, ut.TestCase):
# integrate first part numerically, fit exponential to tail
t_max_fit = 50 * tau
ts = np.arange(0, t_max_fit, 2 * tau)
- numeric_integral = np.trapz(acf[:len(ts)], dx=2 * self.params["tau"])
+ numeric_integral = scipy.integrate.trapezoid(acf[:len(ts)], dx=2 * self.params["tau"])
# fit tail
def fit(x, a, b): return a * np.exp(-b * x)
diff --git a/testsuite/python/p3m_madelung.py b/testsuite/python/p3m_madelung.py
index af98348672..95f263910a 100644
--- a/testsuite/python/p3m_madelung.py
+++ b/testsuite/python/p3m_madelung.py
@@ -29,7 +29,7 @@ import espressomd.magnetostatics
class Test(ut.TestCase):
"""
Check all P3M algorithms against the Madelung constant of 1D, 2D and 3D
- NaCl lattices. See user guide sections :ref:`Madelung electrostatics` and
+ lattices. See user guide sections :ref:`Madelung electrostatics` and
:ref:`Madelung magnetostatics` for more details.
"""
diff --git a/testsuite/python/tests_common.py b/testsuite/python/tests_common.py
index 9b7f9b8722..042fece595 100644
--- a/testsuite/python/tests_common.py
+++ b/testsuite/python/tests_common.py
@@ -318,16 +318,14 @@ def lj_generic_potential(r, epsilon, sigma, cutoff, offset=0., shift=0.,
def lj_generic_force(espressomd, r, epsilon, sigma, cutoff, offset=0., e1=12,
e2=6, b1=4., b2=4., delta=0., lam=1., generic=True):
- f = 1.
- if r >= offset + cutoff:
- f = 0.
- else:
- has_ljgen_softcore = espressomd.has_features("LJGEN_SOFTCORE")
+ f = 0.
+ if r < offset + cutoff:
h = (r - offset)**2 + delta * (1. - lam) * sigma**2
f = (r - offset) * epsilon * lam * (
b1 * e1 * np.power(sigma / np.sqrt(h), e1) -
b2 * e2 * np.power(sigma / np.sqrt(h), e2)) / h
- f *= np.sign(r - offset) if generic and not has_ljgen_softcore else 1
+ if generic and not espressomd.has_features("LJGEN_SOFTCORE"):
+ f *= np.sign(r - offset)
return f
# Lennard-Jones
@@ -339,9 +337,9 @@ def lj_potential(r, epsilon, sigma, cutoff, shift, offset=0.):
return V
-def lj_force(espressomd, r, epsilon, sigma, cutoff, offset=0.):
+def lj_force(r, epsilon, sigma, cutoff, offset=0.):
f = lj_generic_force(
- espressomd, r, epsilon, sigma, cutoff, offset=offset, generic=False)
+ None, r, epsilon, sigma, cutoff, offset=offset, generic=False)
return f
diff --git a/testsuite/python/thermalized_bond.py b/testsuite/python/thermalized_bond.py
index 527bef5095..434d118bb4 100644
--- a/testsuite/python/thermalized_bond.py
+++ b/testsuite/python/thermalized_bond.py
@@ -45,6 +45,7 @@ class ThermalizedBond(ut.TestCase, thermostats_common.ThermostatsCommon):
def tearDown(self):
self.system.part.clear()
+ self.system.bonded_inter.clear()
self.system.thermostat.turn_off()
def test_com_langevin(self):
@@ -149,11 +150,18 @@ class ThermalizedBond(ut.TestCase, thermostats_common.ThermostatsCommon):
v_stored, v_minmax, bins, error_tol, t_dist)
def test_exceptions(self):
- espressomd.interactions.ThermalizedBond(
+ bond = espressomd.interactions.ThermalizedBond(
temp_com=1., gamma_com=1., temp_distance=1., gamma_distance=1.,
r_cut=3.)
+ self.system.bonded_inter.add(bond)
with self.assertRaisesRegex(Exception, "Thermalized bonds require the thermalized_bond thermostat"):
self.system.integrator.run(0)
+ self.system.thermostat.set_thermalized_bond(seed=51)
+ p1 = self.system.part.add(pos=[0., 0., 0.])
+ p2 = self.system.part.add(pos=[0., 0., bond.r_cut * 1.01])
+ p1.bonds = ((bond, p2),)
+ with self.assertRaisesRegex(Exception, r"while calling method integrate\(\)"):
+ self.system.integrator.run(steps=0, recalc_forces=True)
if __name__ == "__main__":
diff --git a/testsuite/python/thermostats_common.py b/testsuite/python/thermostats_common.py
index f0457d4d13..ce49112fb1 100644
--- a/testsuite/python/thermostats_common.py
+++ b/testsuite/python/thermostats_common.py
@@ -17,6 +17,7 @@
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
import numpy as np
+import scipy.integrate
import espressomd
import espressomd.accumulators
@@ -28,7 +29,7 @@ def single_component_maxwell(x1, x2, kT, mass=1.):
x = np.linspace(x1, x2, 1000)
maxwell_distr = np.exp(- mass * x**2 / (2. * kT)) * \
np.sqrt(mass / (2. * np.pi * kT))
- return np.trapz(maxwell_distr, x=x)
+ return scipy.integrate.trapezoid(maxwell_distr, x=x)
class ThermostatsCommon:
diff --git a/testsuite/python/widom_insertion.py b/testsuite/python/widom_insertion.py
index 6cbc32b0e7..60b4754fa8 100644
--- a/testsuite/python/widom_insertion.py
+++ b/testsuite/python/widom_insertion.py
@@ -22,6 +22,7 @@
import unittest as ut
import unittest_decorators as utx
import numpy as np
+import scipy.integrate
import espressomd
import espressomd.reaction_methods
import tests_common
@@ -50,7 +51,7 @@ class WidomInsertionTest(ut.TestCase):
radius = np.linspace(1e-10, LJ_CUT, 1000)
# numerical integration for radii smaller than the cut-off in spherical
# coordinates
- integrateUpToCutOff = 4 * np.pi * np.trapz(
+ integrateUpToCutOff = 4 * np.pi * scipy.integrate.trapezoid(
radius**2 * np.exp(-tests_common.lj_potential(radius,
LJ_EPS,
LJ_SIG,
--
GitLab