Fix thermostat and integrator checkpointing

src/core/npt.cpp

@@ -16,20 +16,27 @@
  * You should have received a copy of the GNU General Public License
  * along with this program. If not, see <http://www.gnu.org/licenses/>.
  */
-#include "npt.hpp"
+#include <cmath>
+
 #include "communication.hpp"
 #include "config.hpp"
 #include "errorhandling.hpp"
 #include "integrate.hpp"
+#include "npt.hpp"
 
 #ifdef NPT
 void synchronize_npt_state(int n_steps) {
  nptiso.invalidate_p_vel = false;
  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);
 }
 

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

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

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

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

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)

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}