From 8977e50c597486d441690d4c592f0b6b6d388e80 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 4 Dec 2020 01:43:03 +0100 Subject: [PATCH] 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 | 8 +-- src/core/npt.cpp | 47 ++++++++----- testsuite/python/integrator_npt.py | 62 +++++++++++++++++ .../python/integrator_steepest_descent.py | 67 ++++++++++++++++++- 4 files changed, 160 insertions(+), 24 deletions(-) diff --git a/src/core/integrators/steepest_descent.cpp b/src/core/integrators/steepest_descent.cpp index c7808e1ce6e..affb6740207 100644 --- a/src/core/integrators/steepest_descent.cpp +++ b/src/core/integrators/steepest_descent.cpp @@ -103,7 +103,7 @@ 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); } @@ -111,7 +111,7 @@ 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(); } diff --git a/src/core/npt.cpp b/src/core/npt.cpp index 3f9e9420efc..eeb6fe429eb 100644 --- a/src/core/npt.cpp +++ b/src/core/npt.cpp @@ -75,9 +75,6 @@ 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."); @@ -85,47 +82,61 @@ void nptiso_init(double ext_pressure, double piston, bool xdir_rescale, 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); } diff --git a/testsuite/python/integrator_npt.py b/testsuite/python/integrator_npt.py index a2daf9a6439..2bfd296d486 100644 --- a/testsuite/python/integrator_npt.py +++ b/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() diff --git a/testsuite/python/integrator_steepest_descent.py b/testsuite/python/integrator_steepest_descent.py index d171307463c..f4520ea1867 100644 --- a/testsuite/python/integrator_steepest_descent.py +++ b/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()