diff --git a/src/core/npt.cpp b/src/core/npt.cpp
index 1a202777d05..ac9b2d42316 100644
--- a/src/core/npt.cpp
+++ b/src/core/npt.cpp
@@ -16,11 +16,13 @@
* You should have received a copy of the GNU General Public License
* along with this program. If not, see .
*/
-#include "npt.hpp"
+#include
+
#include "communication.hpp"
#include "config.hpp"
#include "errorhandling.hpp"
#include "integrate.hpp"
+#include "npt.hpp"
#ifdef NPT
void synchronize_npt_state(int n_steps) {
@@ -28,8 +30,13 @@ void synchronize_npt_state(int n_steps) {
MPI_Bcast(&nptiso.p_inst, 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&nptiso.p_diff, 1, MPI_DOUBLE, 0, comm_cart);
MPI_Bcast(&nptiso.volume, 1, MPI_DOUBLE, 0, comm_cart);
- if (this_node == 0)
- nptiso.p_inst_av /= 1.0 * n_steps;
+ if (this_node == 0) {
+ if (n_steps == 0) {
+ nptiso.p_inst_av = std::nan("");
+ } else {
+ nptiso.p_inst_av /= 1.0 * n_steps;
+ }
+ }
MPI_Bcast(&nptiso.p_inst_av, 1, MPI_DOUBLE, 0, comm_cart);
}
diff --git a/src/python/espressomd/integrate.pxd b/src/python/espressomd/integrate.pxd
index 2f26462123c..1eff4454a4e 100644
--- a/src/python/espressomd/integrate.pxd
+++ b/src/python/espressomd/integrate.pxd
@@ -58,8 +58,3 @@ cdef extern from "integrators/steepest_descent.hpp":
void minimize_energy_init(const double f_max, const double gamma, const int max_steps, const double max_displacement)
cdef extern from "communication.hpp":
int mpi_minimize_energy()
-
-
-# cdef class Integrator:
-# cdef public _method
-# cdef public _steepest_descent_params
diff --git a/src/python/espressomd/integrate.pyx b/src/python/espressomd/integrate.pyx
index 3cb2f08523c..3aab04f97d1 100644
--- a/src/python/espressomd/integrate.pyx
+++ b/src/python/espressomd/integrate.pyx
@@ -39,6 +39,7 @@ cdef class Integrator:
self._steepest_descent_params = {}
self._isotropic_npt_params = {}
+ # __getstate__ and __setstate__ define the pickle interaction
def __getstate__(self):
state = {}
state['_method'] = self._method
@@ -47,15 +48,18 @@ cdef class Integrator:
return state
def __setstate__(self, state):
+ self.__init__()
self._method = state['_method']
if self._method == "STEEPEST_DESCENT":
- self.set_steepest_descent(state['_steepest_descent_params'])
+ self.set_steepest_descent(**state['_steepest_descent_params'])
elif self._method == "NVT":
self.set_nvt()
elif self._method == "NPT":
- npt_params = state['_isotropic_npt_params']
- self.set_isotropic_npt(npt_params['ext_pressure'], npt_params[
- 'piston'], direction=npt_params['direction'], cubic_box=npt_params['cubic_box'])
+ self.set_isotropic_npt(**state['_isotropic_npt_params'])
+
+ def get_state(self):
+ """Returns the integrator status."""
+ return self.__getstate__()
def run(self, steps=1, recalc_forces=False, reuse_forces=False):
"""
diff --git a/src/python/espressomd/thermostat.pyx b/src/python/espressomd/thermostat.pyx
index a49a1750c1e..68477ff6c3d 100644
--- a/src/python/espressomd/thermostat.pyx
+++ b/src/python/espressomd/thermostat.pyx
@@ -104,10 +104,11 @@ cdef class Thermostat:
self.set_lb(
LB_fluid=thmst["LB_fluid"],
act_on_virtual=thmst["act_on_virtual"],
+ gamma=thmst["gamma"],
seed=thmst["rng_counter_fluid"])
if thmst["type"] == "NPT_ISO":
- self.set_npt(kT=thmst["kT"], p_diff=thmst["p_diff"],
- piston=thmst["piston"])
+ self.set_npt(kT=thmst["kT"], gamma0=thmst["gamma0"],
+ gammav=thmst["gammav"])
if thmst["type"] == "DPD":
self.set_dpd(kT=thmst["kT"], seed=thmst["seed"])
@@ -157,16 +158,11 @@ cdef class Thermostat:
npt_dict = {}
npt_dict["type"] = "NPT_ISO"
npt_dict["kT"] = temperature
+ npt_dict["gamma0"] = nptiso_gamma0
+ npt_dict["gammav"] = nptiso_gammav
npt_dict.update(nptiso)
- # thermo_dict["gamma0"] = nptiso_gamma0
- # thermo_dict["gammav"] = nptiso_gammav
- # thermo_dict["p_ext"] = nptiso.p_ext
- # thermo_dict["p_inst"] = nptiso.p_inst
- # thermo_dict["p_inst_av"] = nptiso.p_inst_av
- # thermo_dict["piston"] = nptiso.piston
- # thermo_dict["p_diff"] = nptiso.p_diff
thermo_list.append(npt_dict)
- if (thermo_switch & THERMO_DPD):
+ if thermo_switch & THERMO_DPD:
IF DPD:
dpd_dict = {}
dpd_dict["type"] = "DPD"
@@ -221,7 +217,7 @@ cdef class Thermostat:
Thermal energy of the simulated heat bath.
gamma : :obj:`float`
Contains the friction coefficient of the bath. If the feature
- ``PARTICLE_ANISOTROPY`` is compiled in then ``gamma`` can be a list
+ ``PARTICLE_ANISOTROPY`` is compiled in, then ``gamma`` can be a list
of three positive floats, for the friction coefficient in each
cardinal direction.
gamma_rotation : :obj:`float`, optional
diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt
index f8cdcc7229e..407a98321d4 100644
--- a/testsuite/python/CMakeLists.txt
+++ b/testsuite/python/CMakeLists.txt
@@ -64,13 +64,18 @@ endfunction(PYTHON_TEST)
# Checkpointing tests: semicolon-separated list of mutually-compatible features.
# Separate features with hyphens, use a period to add an optional flag.
-foreach(TEST_COMBINATION lb.cpu-p3m.cpu-lj-lbtherm;lb.gpu-p3m.cpu-lj-lbtherm;ek.gpu)
+foreach(TEST_COMBINATION lb.cpu-p3m.cpu-lj-therm.lb;lb.gpu-p3m.cpu-lj-therm.lb;ek.gpu;lb.off-therm.npt-int.npt;lb.off-therm.langevin-int.nvt;lb.off-therm.dpd-int.sd)
if(${TEST_COMBINATION} MATCHES "\\.gpu")
set(TEST_LABELS "gpu")
else()
set(TEST_LABELS "")
endif()
- foreach(TEST_BINARY 1;0)
+ if(${TEST_COMBINATION} MATCHES "lb\\.off")
+ set(TEST_BINARY_LIST "1")
+ else()
+ set(TEST_BINARY_LIST "1;0")
+ endif()
+ foreach(TEST_BINARY ${TEST_BINARY_LIST})
python_test(FILE save_checkpoint.py MAX_NUM_PROC 4 LABELS ${TEST_LABELS} SUFFIX ${TEST_COMBINATION}_${TEST_BINARY})
python_test(FILE test_checkpoint.py MAX_NUM_PROC 4 LABELS ${TEST_LABELS} SUFFIX ${TEST_COMBINATION}_${TEST_BINARY}
DEPENDS save_checkpoint_${TEST_COMBINATION}_${TEST_BINARY})
diff --git a/testsuite/python/save_checkpoint.py b/testsuite/python/save_checkpoint.py
index 812bc471931..36ef3a8a226 100644
--- a/testsuite/python/save_checkpoint.py
+++ b/testsuite/python/save_checkpoint.py
@@ -56,7 +56,7 @@
if LB_implementation:
lbf = LB_implementation(agrid=0.5, visc=1.3, dens=1.5, tau=0.01)
system.actors.add(lbf)
- if 'LBTHERM' in modes:
+ if 'THERM.LB' in modes:
system.thermostat.set_lb(LB_fluid=lbf, seed=23, gamma=2.0)
if any(has_features(i) for i in ["LB_BOUNDARIES", "LB_BOUNDARIES_GPU"]):
if 'EK.GPU' not in modes:
@@ -133,9 +133,23 @@
system.constraints.add(constraints.ElectricPlaneWave(
E0=[1., -2., 3.], k=[-.1, .2, .3], omega=5., phi=1.4))
-
-if 'LBTHERM' not in modes:
- system.thermostat.set_langevin(kT=1.0, gamma=2.0, seed=42)
+if 'LB.OFF' in modes:
+ # set thermostat
+ if 'THERM.LANGEVIN' in modes:
+ system.thermostat.set_langevin(kT=1.0, gamma=2.0, seed=42)
+ elif 'THERM.NPT' in modes and has_features('NPT'):
+ system.thermostat.set_npt(kT=1.0, gamma0=2.0, gammav=0.1)
+ elif 'THERM.DPD' in modes and has_features('DPD'):
+ system.thermostat.set_dpd(kT=1.0, seed=42)
+ # set integrator
+ if 'INT.NPT' in modes and has_features('NPT'):
+ system.integrator.set_isotropic_npt(ext_pressure=2.0, piston=0.01,
+ direction=[1, 0, 0])
+ elif 'INT.SD' in modes:
+ system.integrator.set_steepest_descent(f_max=2.0, gamma=0.1,
+ max_displacement=0.01)
+ elif 'INT.NVT' in modes:
+ system.integrator.set_nvt()
if espressomd.has_features(['VIRTUAL_SITES', 'VIRTUAL_SITES_RELATIVE']):
system.virtual_sites = espressomd.virtual_sites.VirtualSitesRelative(
@@ -151,7 +165,7 @@
harmonic_bond = espressomd.interactions.HarmonicBond(r_0=0.0, k=1.0)
system.bonded_inter.add(harmonic_bond)
system.part[1].add_bond((harmonic_bond, 0))
-if 'LBTHERM' not in modes:
+if 'THERM.LB' not in modes:
thermalized_bond = espressomd.interactions.ThermalizedBond(
temp_com=0.0, gamma_com=0.0, temp_distance=0.2, gamma_distance=0.5,
r_cut=2, seed=51)
diff --git a/testsuite/python/test_checkpoint.py b/testsuite/python/test_checkpoint.py
index e1abc5e14c5..9656500f512 100644
--- a/testsuite/python/test_checkpoint.py
+++ b/testsuite/python/test_checkpoint.py
@@ -139,19 +139,71 @@ def test_part(self):
np.testing.assert_allclose(np.copy(system.part[0].f), particle_force0)
np.testing.assert_allclose(np.copy(system.part[1].f), particle_force1)
- @ut.skipIf('LBTHERM' not in modes, 'LB thermostat not in modes')
- def test_thermostat(self):
- self.assertEqual(system.thermostat.get_state()[0]['type'], 'LB')
- self.assertEqual(system.thermostat.get_state()[0]['seed'], 23)
- self.assertEqual(system.thermostat.get_state()[0]['gamma'], 2.0)
-
- @ut.skipIf('LBTHERM' in modes, 'Langevin incompatible with LB thermostat')
- def test_thermostat(self):
- self.assertEqual(system.thermostat.get_state()[0]['type'], 'LANGEVIN')
- self.assertEqual(system.thermostat.get_state()[0]['kT'], 1.0)
- self.assertEqual(system.thermostat.get_state()[0]['seed'], 42)
- np.testing.assert_array_equal(system.thermostat.get_state()[
- 0]['gamma'], np.array([2.0, 2.0, 2.0]))
+ @ut.skipIf('THERM.LB' not in modes, 'LB thermostat not in modes')
+ def test_thermostat_LB(self):
+ thmst = system.thermostat.get_state()[0]
+ if 'LB.GPU' in modes and not espressomd.gpu_available():
+ self.assertEqual(thmst['type'], 'OFF')
+ else:
+ self.assertEqual(thmst['type'], 'LB')
+ # rng_counter_fluid = seed, seed is 0 because kT=0
+ self.assertEqual(thmst['rng_counter_fluid'], 0)
+ self.assertEqual(thmst['gamma'], 2.0)
+
+ @ut.skipIf('THERM.LANGEVIN' not in modes,
+ 'Langevin thermostat not in modes')
+ def test_thermostat_Langevin(self):
+ thmst = system.thermostat.get_state()[0]
+ self.assertEqual(thmst['type'], 'LANGEVIN')
+ self.assertEqual(thmst['kT'], 1.0)
+ self.assertEqual(thmst['seed'], 42)
+ np.testing.assert_array_equal(thmst['gamma'], np.array(3 * [2.0]))
+
+ @utx.skipIfMissingFeatures('DPD')
+ @ut.skipIf('THERM.DPD' not in modes, 'DPD thermostat not in modes')
+ def test_thermostat_DPD(self):
+ thmst = system.thermostat.get_state()[0]
+ self.assertEqual(thmst['type'], 'DPD')
+ self.assertEqual(thmst['kT'], 1.0)
+ self.assertEqual(thmst['seed'], 42 + 6)
+
+ @utx.skipIfMissingFeatures('NPT')
+ @ut.skipIf('THERM.NPT' not in modes, 'NPT thermostat not in modes')
+ def test_thermostat_NPT(self):
+ thmst = system.thermostat.get_state()[0]
+ self.assertEqual(thmst['type'], 'NPT_ISO')
+ self.assertEqual(thmst['gamma0'], 2.0)
+ self.assertEqual(thmst['gammav'], 0.1)
+
+ @utx.skipIfMissingFeatures('NPT')
+ @ut.skipIf('INT.NPT' not in modes, 'NPT integrator not in modes')
+ def test_integrator_NPT(self):
+ integ = system.integrator.get_state()
+ self.assertEqual(integ['_method'], 'NPT')
+ params = integ['_isotropic_npt_params']
+ self.assertEqual(params['ext_pressure'], 2.0)
+ self.assertEqual(params['piston'], 0.01)
+ self.assertEqual(params['direction'], [1, 0, 0])
+ self.assertEqual(params['cubic_box'], False)
+
+ @ut.skipIf('INT.SD' not in modes, 'SD integrator not in modes')
+ def test_integrator_SD(self):
+ integ = system.integrator.get_state()
+ self.assertEqual(integ['_method'], 'STEEPEST_DESCENT')
+ params = integ['_steepest_descent_params']
+ self.assertEqual(params['f_max'], 2.0)
+ self.assertEqual(params['gamma'], 0.1)
+ self.assertEqual(params['max_displacement'], 0.01)
+
+ @ut.skipIf('INT.NVT' not in modes, 'NVT integrator not in modes')
+ def test_integrator_VV(self):
+ integ = system.integrator.get_state()
+ self.assertEqual(integ['_method'], 'NVT')
+
+ @ut.skipIf('INT' in modes, 'VV integrator not the default')
+ def test_integrator_VV(self):
+ integ = system.integrator.get_state()
+ self.assertEqual(integ['_method'], 'VV')
@utx.skipIfMissingFeatures('LENNARD_JONES')
@ut.skipIf('LJ' not in modes, "Skipping test due to missing mode.")
@@ -174,7 +226,7 @@ def test_bonded_inter(self):
reference = {'r_0': 0.0, 'k': 1.0}
for key in reference.keys():
self.assertAlmostEqual(state[key], reference[key], delta=1E-10)
- if 'LBTHERM' not in modes:
+ if 'THERM.LB' not in modes:
state = system.part[1].bonds[1][0].params
reference = {'temp_com': 0., 'gamma_com': 0., 'temp_distance': 0.2,
'gamma_distance': 0.5, 'r_cut': 2.0, 'seed': 51}