Skip to content

Commit 8977e50

Browse files
jngradjngrad
committed
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).
1 parent 0500410 commit 8977e50

File tree

4 files changed

+160
-24
lines changed

4 files changed

+160
-24
lines changed

‎src/core/integrators/steepest_descent.cpp‎

Lines changed: 3 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -103,15 +103,15 @@ bool steepest_descent_step(const ParticleRange &particles) {
103103
return sqrt(f_max_global) < params.f_max;
104104
}
105105

106-
void mpi_bcast_steepest_descent_worker(int, int) {
106+
void mpi_bcast_steepest_descent_worker() {
107107
boost::mpi::broadcast(comm_cart, params, 0);
108108
}
109109

110110
REGISTER_CALLBACK(mpi_bcast_steepest_descent_worker)
111111

112112
/** Broadcast steepest descent parameters */
113113
void mpi_bcast_steepest_descent() {
114-
mpi_call_all(mpi_bcast_steepest_descent_worker, -1, 0);
114+
mpi_call_all(mpi_bcast_steepest_descent_worker);
115115
}
116116

117117
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,
126126
throw std::runtime_error("The maximal displacement must be positive.");
127127
}
128128

129-
params.f_max = f_max;
130-
params.gamma = gamma;
131-
params.max_displacement = max_displacement;
129+
params = SteepestDescentParameters{f_max, gamma, max_displacement};
132130

133131
mpi_bcast_steepest_descent();
134132
}

‎src/core/npt.cpp‎

Lines changed: 29 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -75,57 +75,68 @@ void mpi_bcast_nptiso_geom_barostat() {
7575

7676
void nptiso_init(double ext_pressure, double piston, bool xdir_rescale,
7777
bool ydir_rescale, bool zdir_rescale, bool cubic_box) {
78-
nptiso.cubic_box = cubic_box;
79-
nptiso.p_ext = ext_pressure;
80-
nptiso.piston = piston;
8178

8279
if (ext_pressure < 0.0) {
8380
throw std::runtime_error("The external pressure must be positive.");
8481
}
8582
if (piston <= 0.0) {
8683
throw std::runtime_error("The piston mass must be positive.");
8784
}
85+
86+
nptiso_struct new_nptiso = {piston,
87+
nptiso.inv_piston,
88+
nptiso.volume,
89+
ext_pressure,
90+
nptiso.p_inst,
91+
nptiso.p_diff,
92+
nptiso.p_vir,
93+
nptiso.p_vel,
94+
0,
95+
nptiso.nptgeom_dir,
96+
0,
97+
cubic_box,
98+
-1};
99+
88100
/* set the NpT geometry */
89-
nptiso.geometry = 0;
90-
nptiso.dimension = 0;
91-
nptiso.non_const_dim = -1;
92101
if (xdir_rescale) {
93-
nptiso.geometry |= NPTGEOM_XDIR;
94-
nptiso.dimension += 1;
95-
nptiso.non_const_dim = 0;
102+
new_nptiso.geometry |= NPTGEOM_XDIR;
103+
new_nptiso.dimension += 1;
104+
new_nptiso.non_const_dim = 0;
96105
}
97106
if (ydir_rescale) {
98-
nptiso.geometry |= NPTGEOM_YDIR;
99-
nptiso.dimension += 1;
100-
nptiso.non_const_dim = 1;
107+
new_nptiso.geometry |= NPTGEOM_YDIR;
108+
new_nptiso.dimension += 1;
109+
new_nptiso.non_const_dim = 1;
101110
}
102111
if (zdir_rescale) {
103-
nptiso.geometry |= NPTGEOM_ZDIR;
104-
nptiso.dimension += 1;
105-
nptiso.non_const_dim = 2;
112+
new_nptiso.geometry |= NPTGEOM_ZDIR;
113+
new_nptiso.dimension += 1;
114+
new_nptiso.non_const_dim = 2;
106115
}
107116

108117
/* Sanity Checks */
109118
#ifdef ELECTROSTATICS
110-
if (nptiso.dimension < 3 && !nptiso.cubic_box && coulomb.prefactor > 0) {
119+
if (new_nptiso.dimension < 3 && !cubic_box && coulomb.prefactor > 0) {
111120
throw std::runtime_error("If electrostatics is being used you must "
112121
"use the cubic box npt.");
113122
}
114123
#endif
115124

116125
#ifdef DIPOLES
117-
if (nptiso.dimension < 3 && !nptiso.cubic_box && dipole.prefactor > 0) {
126+
if (new_nptiso.dimension < 3 && !cubic_box && dipole.prefactor > 0) {
118127
throw std::runtime_error("If magnetostatics is being used you must "
119128
"use the cubic box npt.");
120129
}
121130
#endif
122131

123-
if (nptiso.dimension == 0 || nptiso.non_const_dim == -1) {
132+
if (new_nptiso.dimension == 0 || new_nptiso.non_const_dim == -1) {
124133
throw std::runtime_error(
125134
"You must enable at least one of the x y z components "
126135
"as fluctuating dimension(s) for box length motion!");
127136
}
128137

138+
nptiso = new_nptiso;
139+
129140
mpi_bcast_nptiso_geom_barostat();
130141
on_parameter_change(FIELD_NPTISO_PISTON);
131142
}

‎testsuite/python/integrator_npt.py‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,9 @@
1818
#
1919
import unittest as ut
2020
import unittest_decorators as utx
21+
import numpy as np
2122
import espressomd
23+
import espressomd.integrate
2224

2325

2426
@utx.skipIfMissingFeatures(["NPT", "LENNARD_JONES"])
@@ -52,6 +54,66 @@ def test_integrator_exceptions(self):
5254
self.system.integrator.set_isotropic_npt(
5355
ext_pressure=1, piston=1, direction=[1, 0])
5456

57+
def test_integrator_recovery(self):
58+
# the system is still in a valid state after a failure
59+
system = self.system
60+
np.random.seed(42)
61+
npt_params = {'ext_pressure': 0.001, 'piston': 0.001}
62+
system.box_l = [6] * 3
63+
system.part.add(pos=np.random.uniform(0, system.box_l[0], (11, 3)))
64+
system.non_bonded_inter[0, 0].lennard_jones.set_params(
65+
epsilon=1, sigma=1, cutoff=2**(1 / 6), shift=0.25)
66+
system.thermostat.set_npt(kT=1.0, gamma0=2, gammav=0.04, seed=42)
67+
system.integrator.set_isotropic_npt(**npt_params)
68+
69+
# get the equilibrium box length for the chosen NpT parameters
70+
system.integrator.run(2000)
71+
box_l_ref = system.box_l[0]
72+
73+
# resetting the NpT integrator with incorrect values doesn't leave the
74+
# system in an undefined state (the old parameters aren't overwritten)
75+
with self.assertRaises(RuntimeError):
76+
system.integrator.set_isotropic_npt(ext_pressure=-1, piston=100)
77+
with self.assertRaises(RuntimeError):
78+
system.integrator.set_isotropic_npt(ext_pressure=100, piston=-1)
79+
# the core state is unchanged
80+
system.integrator.run(500)
81+
self.assertAlmostEqual(system.box_l[0], box_l_ref, delta=0.15)
82+
83+
# setting another integrator with incorrect values doesn't leave the
84+
# system in an undefined state (the old integrator is still active)
85+
with self.assertRaises(RuntimeError):
86+
system.integrator.set_steepest_descent(
87+
f_max=-10, gamma=0, max_displacement=0.1)
88+
# the interface state is unchanged
89+
self.assertIsInstance(system.integrator.get_state(),
90+
espressomd.integrate.VelocityVerletIsotropicNPT)
91+
params = system.integrator.get_state().get_params()
92+
self.assertEqual(params['ext_pressure'], npt_params['ext_pressure'])
93+
self.assertEqual(params['piston'], npt_params['piston'])
94+
# the core state is unchanged
95+
system.integrator.run(500)
96+
self.assertAlmostEqual(system.box_l[0], box_l_ref, delta=0.15)
97+
98+
# setting the NpT integrator with incorrect values doesn't leave the
99+
# system in an undefined state (the old integrator is still active)
100+
system.thermostat.turn_off()
101+
system.integrator.set_vv()
102+
system.part.clear()
103+
system.box_l = [5] * 3
104+
positions_start = np.array([[0, 0, 0], [1., 0, 0]])
105+
system.part.add(pos=positions_start)
106+
with self.assertRaises(RuntimeError):
107+
system.integrator.set_isotropic_npt(ext_pressure=-1, piston=100)
108+
# the interface state is unchanged
109+
self.assertIsInstance(system.integrator.get_state(),
110+
espressomd.integrate.VelocityVerlet)
111+
# the core state is unchanged
112+
system.integrator.run(1)
113+
np.testing.assert_allclose(
114+
np.copy(system.part[:].pos),
115+
positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))
116+
55117

56118
if __name__ == "__main__":
57119
ut.main()

‎testsuite/python/integrator_steepest_descent.py‎

Lines changed: 66 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -38,7 +38,7 @@ class IntegratorSteepestDescent(ut.TestCase):
3838

3939
lj_eps = 1.0
4040
lj_sig = 1.0
41-
lj_cut = 1.12246
41+
lj_cut = 2**(1 / 6)
4242

4343
def setUp(self):
4444
self.system.box_l = 3 * [self.box_l]
@@ -150,6 +150,71 @@ def test_integrator_exceptions(self):
150150
self.system.integrator.set_steepest_descent(
151151
f_max=0, gamma=1, max_displacement=-1)
152152

153+
def test_integrator_recovery(self):
154+
# the system is still in a valid state after a failure
155+
system = self.system
156+
sd_params = {"f_max": 0, "gamma": 1, "max_displacement": 0.01}
157+
positions_start = np.array([[0, 0, 0], [1., 0, 0]])
158+
positions_ref = np.array([[-0.01, 0, 0], [1.01, 0, 0]])
159+
system.part.add(pos=positions_start)
160+
system.integrator.set_steepest_descent(**sd_params)
161+
162+
# get the positions after one step with the chosen parameters
163+
system.integrator.run(1)
164+
positions_ref = np.copy(system.part[:].pos)
165+
166+
# resetting the SD integrator with incorrect values doesn't leave the
167+
# system in an undefined state (the old parameters aren't overwritten)
168+
with self.assertRaises(RuntimeError):
169+
system.integrator.set_steepest_descent(
170+
f_max=-10, gamma=1, max_displacement=0.01)
171+
with self.assertRaises(RuntimeError):
172+
system.integrator.set_steepest_descent(
173+
f_max=0, gamma=-1, max_displacement=0.01)
174+
with self.assertRaises(RuntimeError):
175+
system.integrator.set_steepest_descent(
176+
f_max=0, gamma=1, max_displacement=-1)
177+
# the core state is unchanged
178+
system.part[:].pos = positions_start
179+
system.integrator.run(1)
180+
np.testing.assert_allclose(np.copy(system.part[:].pos), positions_ref)
181+
182+
# setting another integrator with incorrect values doesn't leave the
183+
# system in an undefined state (the old integrator is still active)
184+
if espressomd.has_features("NPT"):
185+
with self.assertRaises(RuntimeError):
186+
system.integrator.set_isotropic_npt(ext_pressure=1, piston=-1)
187+
# the interface state is unchanged
188+
self.assertIsInstance(system.integrator.get_state(),
189+
espressomd.integrate.SteepestDescent)
190+
params = system.integrator.get_state().get_params()
191+
self.assertEqual(params["f_max"], sd_params["f_max"])
192+
self.assertEqual(params["gamma"], sd_params["gamma"])
193+
self.assertEqual(
194+
params["max_displacement"],
195+
sd_params["max_displacement"])
196+
# the core state is unchanged
197+
system.part[:].pos = positions_start
198+
system.integrator.run(1)
199+
np.testing.assert_allclose(
200+
np.copy(system.part[:].pos), positions_ref)
201+
202+
# setting the SD integrator with incorrect values doesn't leave the
203+
# system in an undefined state (the old integrator is still active)
204+
system.integrator.set_vv()
205+
system.part[:].pos = positions_start
206+
with self.assertRaises(RuntimeError):
207+
system.integrator.set_steepest_descent(
208+
f_max=0, gamma=1, max_displacement=-1)
209+
# the interface state is unchanged
210+
self.assertIsInstance(system.integrator.get_state(),
211+
espressomd.integrate.VelocityVerlet)
212+
# the core state is unchanged
213+
system.integrator.run(1)
214+
np.testing.assert_allclose(
215+
np.copy(system.part[:].pos),
216+
positions_start + np.array([[-1.2e-3, 0, 0], [1.2e-3, 0, 0]]))
217+
153218

154219
if __name__ == "__main__":
155220
ut.main()

0 commit comments

Comments
 (0)