core: Make NpT and Steepest Descent init failsafe

Setting the NpT and Steepest Descent integrators with incorrect
parameters will raise a RuntimeError but won't leave the system
in an undefined state (the old integrator remains untouched).

src/core/integrators/steepest_descent.cpp

@@ -103,15 +103,15 @@ bool steepest_descent_step(const ParticleRange &particles) {
   return sqrt(f_max_global) < params.f_max;
 }
 
-void mpi_bcast_steepest_descent_worker(int, int) {
+void mpi_bcast_steepest_descent_worker() {
   boost::mpi::broadcast(comm_cart, params, 0);
 }
 
 REGISTER_CALLBACK(mpi_bcast_steepest_descent_worker)
 
 /** Broadcast steepest descent parameters */
 void mpi_bcast_steepest_descent() {
-  mpi_call_all(mpi_bcast_steepest_descent_worker, -1, 0);
+  mpi_call_all(mpi_bcast_steepest_descent_worker);
 }
 
 void steepest_descent_init(const double f_max, const double gamma,
@@ -126,9 +126,7 @@ void steepest_descent_init(const double f_max, const double gamma,
     throw std::runtime_error("The maximal displacement must be positive.");
   }
 
-  params.f_max = f_max;
-  params.gamma = gamma;
-  params.max_displacement = max_displacement;
+  params = SteepestDescentParameters{f_max, gamma, max_displacement};
 
   mpi_bcast_steepest_descent();
 }

src/core/npt.cpp

@@ -75,57 +75,68 @@ void mpi_bcast_nptiso_geom_barostat() {
 
 void nptiso_init(double ext_pressure, double piston, bool xdir_rescale,
                  bool ydir_rescale, bool zdir_rescale, bool cubic_box) {
-  nptiso.cubic_box = cubic_box;
-  nptiso.p_ext = ext_pressure;
-  nptiso.piston = piston;
 
   if (ext_pressure < 0.0) {
     throw std::runtime_error("The external pressure must be positive.");
   }
   if (piston <= 0.0) {
     throw std::runtime_error("The piston mass must be positive.");
   }
+
+  nptiso_struct new_nptiso = {piston,
+                              nptiso.inv_piston,
+                              nptiso.volume,
+                              ext_pressure,
+                              nptiso.p_inst,
+                              nptiso.p_diff,
+                              nptiso.p_vir,
+                              nptiso.p_vel,
+                              0,
+                              nptiso.nptgeom_dir,
+                              0,
+                              cubic_box,
+                              -1};
+
   /* set the NpT geometry */
-  nptiso.geometry = 0;
-  nptiso.dimension = 0;
-  nptiso.non_const_dim = -1;
   if (xdir_rescale) {
-    nptiso.geometry |= NPTGEOM_XDIR;
-    nptiso.dimension += 1;
-    nptiso.non_const_dim = 0;
+    new_nptiso.geometry |= NPTGEOM_XDIR;
+    new_nptiso.dimension += 1;
+    new_nptiso.non_const_dim = 0;
   }
   if (ydir_rescale) {
-    nptiso.geometry |= NPTGEOM_YDIR;
-    nptiso.dimension += 1;
-    nptiso.non_const_dim = 1;
+    new_nptiso.geometry |= NPTGEOM_YDIR;
+    new_nptiso.dimension += 1;
+    new_nptiso.non_const_dim = 1;
   }
   if (zdir_rescale) {
-    nptiso.geometry |= NPTGEOM_ZDIR;
-    nptiso.dimension += 1;
-    nptiso.non_const_dim = 2;
+    new_nptiso.geometry |= NPTGEOM_ZDIR;
+    new_nptiso.dimension += 1;
+    new_nptiso.non_const_dim = 2;
   }
 
   /* Sanity Checks */
 #ifdef ELECTROSTATICS
-  if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.prefactor > 0) {
+  if (new_nptiso.dimension < 3 && !cubic_box && coulomb.prefactor > 0) {
     throw std::runtime_error("If electrostatics is being used you must "
                              "use the cubic box npt.");
   }
 #endif
 
 #ifdef DIPOLES
-  if (nptiso.dimension < 3 && !nptiso.cubic_box && dipole.prefactor > 0) {
+  if (new_nptiso.dimension < 3 && !cubic_box && dipole.prefactor > 0) {
     throw std::runtime_error("If magnetostatics is being used you must "
                              "use the cubic box npt.");
   }
 #endif
 
-  if (nptiso.dimension == 0 || nptiso.non_const_dim == -1) {
+  if (new_nptiso.dimension == 0 || new_nptiso.non_const_dim == -1) {
     throw std::runtime_error(
         "You must enable at least one of the x y z components "
         "as fluctuating dimension(s) for box length motion!");
   }
 
+  nptiso = new_nptiso;
+
   mpi_bcast_nptiso_geom_barostat();
   on_parameter_change(FIELD_NPTISO_PISTON);
 }

testsuite/python/integrator_npt.py

@@ -18,7 +18,9 @@
 #
 import unittest as ut
 import unittest_decorators as utx
+import numpy as np
 import espressomd
+import espressomd.integrate
 
 
 @utx.skipIfMissingFeatures(["NPT", "LENNARD_JONES"])
@@ -52,6 +54,66 @@ def test_integrator_exceptions(self):
             self.system.integrator.set_isotropic_npt(
                 ext_pressure=1, piston=1, direction=[1, 0])
 
+    def test_integrator_recovery(self):
+        # the system is still in a valid state after a failure
+        system = self.system
+        np.random.seed(42)
+        npt_params = {'ext_pressure': 0.001, 'piston': 0.001}
+        system.box_l = [6] * 3
+        system.part.add(pos=np.random.uniform(0, system.box_l[0], (11, 3)))
+        system.non_bonded_inter[0, 0].lennard_jones.set_params(
+            epsilon=1, sigma=1, cutoff=2**(1 / 6), shift=0.25)
+        system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
+        system.integrator.set_isotropic_npt(**npt_params)
+
+        # get the equilibrium box length for the chosen NpT parameters
+        system.integrator.run(2000)
+        box_l_ref = system.box_l[0]
+
+        # resetting the NpT integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old parameters aren't overwritten)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_isotropic_npt(ext_pressure=-1, piston=100)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_isotropic_npt(ext_pressure=100, piston=-1)
+        # the core state is unchanged
+        system.integrator.run(500)
+        self.assertAlmostEqual(system.box_l[0], box_l_ref, delta=0.15)
+
+        # setting another integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old integrator is still active)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_steepest_descent(
+                f_max=-10, gamma=0, max_displacement=0.1)
+        # the interface state is unchanged
+        self.assertIsInstance(system.integrator.get_state(),
+                              espressomd.integrate.VelocityVerletIsotropicNPT)
+        params = system.integrator.get_state().get_params()
+        self.assertEqual(params['ext_pressure'], npt_params['ext_pressure'])
+        self.assertEqual(params['piston'], npt_params['piston'])
+        # the core state is unchanged
+        system.integrator.run(500)
+        self.assertAlmostEqual(system.box_l[0], box_l_ref, delta=0.15)
+
+        # setting the NpT integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old integrator is still active)
+        system.thermostat.turn_off()
+        system.integrator.set_vv()
+        system.part.clear()
+        system.box_l = [5] * 3
+        positions_start = np.array([[0, 0, 0], [1., 0, 0]])
+        system.part.add(pos=positions_start)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_isotropic_npt(ext_pressure=-1, piston=100)
+        # the interface state is unchanged
+        self.assertIsInstance(system.integrator.get_state(),
+                              espressomd.integrate.VelocityVerlet)
+        # the core state is unchanged
+        system.integrator.run(1)
+        np.testing.assert_allclose(
+            np.copy(system.part[:].pos),
+            positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))
+
 
 if __name__ == "__main__":
     ut.main()

testsuite/python/integrator_steepest_descent.py

@@ -38,7 +38,7 @@ class IntegratorSteepestDescent(ut.TestCase):
 
     lj_eps = 1.0
     lj_sig = 1.0
-    lj_cut = 1.12246
+    lj_cut = 2**(1 / 6)
 
     def setUp(self):
         self.system.box_l = 3 * [self.box_l]
@@ -150,6 +150,71 @@ def test_integrator_exceptions(self):
             self.system.integrator.set_steepest_descent(
                 f_max=0, gamma=1, max_displacement=-1)
 
+    def test_integrator_recovery(self):
+        # the system is still in a valid state after a failure
+        system = self.system
+        sd_params = {"f_max": 0, "gamma": 1, "max_displacement": 0.01}
+        positions_start = np.array([[0, 0, 0], [1., 0, 0]])
+        positions_ref = np.array([[-0.01, 0, 0], [1.01, 0, 0]])
+        system.part.add(pos=positions_start)
+        system.integrator.set_steepest_descent(**sd_params)
+
+        # get the positions after one step with the chosen parameters
+        system.integrator.run(1)
+        positions_ref = np.copy(system.part[:].pos)
+
+        # resetting the SD integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old parameters aren't overwritten)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_steepest_descent(
+                f_max=-10, gamma=1, max_displacement=0.01)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_steepest_descent(
+                f_max=0, gamma=-1, max_displacement=0.01)
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_steepest_descent(
+                f_max=0, gamma=1, max_displacement=-1)
+        # the core state is unchanged
+        system.part[:].pos = positions_start
+        system.integrator.run(1)
+        np.testing.assert_allclose(np.copy(system.part[:].pos), positions_ref)
+
+        # setting another integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old integrator is still active)
+        if espressomd.has_features("NPT"):
+            with self.assertRaises(RuntimeError):
+                system.integrator.set_isotropic_npt(ext_pressure=1, piston=-1)
+            # the interface state is unchanged
+            self.assertIsInstance(system.integrator.get_state(),
+                                  espressomd.integrate.SteepestDescent)
+            params = system.integrator.get_state().get_params()
+            self.assertEqual(params["f_max"], sd_params["f_max"])
+            self.assertEqual(params["gamma"], sd_params["gamma"])
+            self.assertEqual(
+                params["max_displacement"],
+                sd_params["max_displacement"])
+            # the core state is unchanged
+            system.part[:].pos = positions_start
+            system.integrator.run(1)
+            np.testing.assert_allclose(
+                np.copy(system.part[:].pos), positions_ref)
+
+        # setting the SD integrator with incorrect values doesn't leave the
+        # system in an undefined state (the old integrator is still active)
+        system.integrator.set_vv()
+        system.part[:].pos = positions_start
+        with self.assertRaises(RuntimeError):
+            system.integrator.set_steepest_descent(
+                f_max=0, gamma=1, max_displacement=-1)
+        # the interface state is unchanged
+        self.assertIsInstance(system.integrator.get_state(),
+                              espressomd.integrate.VelocityVerlet)
+        # the core state is unchanged
+        system.integrator.run(1)
+        np.testing.assert_allclose(
+            np.copy(system.part[:].pos),
+            positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))
+
 
 if __name__ == "__main__":
     ut.main()