Disable box length reset when particles are present (espressomd#4901)

Fixes espressomd#4834

Description of changes:
- Resizing the box via `system.box_l = new_box_l` now throws a runtime error if there are particles present, because particle position folding cannot be guaranteed to be correct; use `system.change_volume_and_rescale_particles()` instead, which properly rescales particle positions.

samples/h5md_trajectory.py

@@ -66,7 +66,8 @@
 xyz_folded.append(system.part.all().pos_folded[:])
 xyz_unfolded.append(system.part.all().pos[:])
 # resize box (simulates NpT)
-system.box_l = system.box_l + 1.
+for i in range(3):
+    system.change_volume_and_rescale_particles(system.box_l[i] + 1., "xyz"[i])
 system.integrator.run(10)
 h5.write()
 xyz_folded.append(system.part.all().pos_folded[:])

src/core/constraints/Constraints.hpp

@@ -101,12 +101,14 @@ template <class ParticleRange, class Constraint> class Constraints {
     }
   }
 
-  void on_boxl_change() const {
+  void veto_boxl_change() const {
     if (not m_constraints.empty()) {
       throw std::runtime_error("The box size can not be changed because there "
                                "are active constraints.");
     }
   }
+
+  void on_boxl_change() const { veto_boxl_change(); }
 };
 } // namespace Constraints
 

src/core/event.cpp

@@ -269,6 +269,16 @@ void on_lbboundary_change() {
 #endif
 }
 
+void veto_boxl_change(bool skip_particle_checks) {
+  if (not skip_particle_checks) {
+    if (get_n_part() > 0) {
+      throw std::runtime_error(
+          "Cannot reset the box length when particles are present");
+    }
+  }
+  Constraints::constraints.veto_boxl_change();
+}
+
 void on_boxl_change(bool skip_method_adaption) {
   grid_changed_box_l(box_geo);
   /* Electrostatics cutoffs mostly depend on the system size,

src/core/event.hpp

@@ -82,6 +82,13 @@ void on_short_range_ia_change();
 /** called every time a constraint is changed. */
 void on_constraint_change();
 
+/**
+ * @brief Called before the box length is changed.
+ *
+ * @param skip_particle_checks skip the particle checks
+ */
+void veto_boxl_change(bool skip_particle_checks);
+
 /**
  * @brief Called when the box length has changed. This routine is relatively
  * fast, and changing the box length every time step as for example necessary

src/core/grid.cpp

@@ -107,10 +107,26 @@ void grid_changed_n_nodes() {
   grid_changed_box_l(box_geo);
 }
 
+void mpi_set_box_length_local(const Utils::Vector3d &length) {
+  box_geo.set_length(length);
+  try {
+    on_boxl_change();
+  } catch (std::exception const &err) {
+    runtimeErrorMsg() << err.what();
+  }
+}
+
+REGISTER_CALLBACK(mpi_set_box_length_local)
+
 void rescale_boxl(int dir, double d_new) {
+  assert(dir != 3 or ((box_geo.length()[0] == box_geo.length()[1]) and
+                      (box_geo.length()[1] == box_geo.length()[2])));
   double scale = (dir - 3) ? d_new * box_geo.length_inv()[dir]
                            : d_new * box_geo.length_inv()[0];
 
+  assert(scale > 0.);
+  veto_boxl_change(true);
+
   /* If shrinking, rescale the particles first. */
   if (scale <= 1.) {
     mpi_rescale_particles(dir, scale);
@@ -119,32 +135,26 @@ void rescale_boxl(int dir, double d_new) {
   if (dir < 3) {
     auto box_l = box_geo.length();
     box_l[dir] = d_new;
-    mpi_set_box_length(box_l);
+    mpi_call_all(mpi_set_box_length_local, box_l);
   } else {
-    mpi_set_box_length({d_new, d_new, d_new});
+    mpi_call_all(mpi_set_box_length_local, Utils::Vector3d::broadcast(d_new));
   }
 
   if (scale > 1.) {
     mpi_rescale_particles(dir, scale);
   }
 }
 
-void mpi_set_box_length_local(const Utils::Vector3d &length) {
-  box_geo.set_length(length);
-  try {
-    on_boxl_change();
-  } catch (std::exception const &err) {
-    runtimeErrorMsg() << err.what();
-  }
-}
-
-REGISTER_CALLBACK(mpi_set_box_length_local)
-
-void mpi_set_box_length(const Utils::Vector3d &length) {
+static void check_new_length(const Utils::Vector3d &length) {
   if (boost::algorithm::any_of(length,
                                [](double value) { return value <= 0; })) {
     throw std::domain_error("Box length must be >0");
   }
+  veto_boxl_change(false);
+}
+
+void mpi_set_box_length(const Utils::Vector3d &length) {
+  check_new_length(length);
 
   mpi_call_all(mpi_set_box_length_local, length);
 }

src/python/espressomd/system.pxd

@@ -24,7 +24,7 @@ from .utils cimport Vector3d
 cdef extern from "grid.hpp":
     void mpi_set_box_length(Vector3d length) except +
     void mpi_set_periodicity(bool x, bool y, bool z)
-    void rescale_boxl(int dir, double d_new)
+    void rescale_boxl(int dir, double d_new) except +
 
 cdef extern from "rotate_system.hpp":
     void mpi_rotate_system(double phi, double theta, double alpha)

src/python/espressomd/system.pyx

@@ -362,8 +362,8 @@ cdef class System:
 
         """
 
-        if d_new < 0:
-            raise ValueError("No negative lengths")
+        if d_new <= 0.:
+            raise ValueError("Parameter 'd_new' must be > 0")
         if dir == "xyz":
             rescale_boxl(3, d_new)
         elif dir == "x" or dir == 0:
@@ -375,6 +375,7 @@ cdef class System:
         else:
             raise ValueError(
                 'Usage: change_volume_and_rescale_particles(<L_new>, [{ "x" | "y" | "z" | "xyz" }])')
+        utils.handle_errors("Exception while updating the box length")
 
     def volume(self):
         """Return box volume of the cuboid box.

testsuite/python/analyze_chains.py

@@ -125,7 +125,8 @@ def test_observables_no_pbc(self):
         # increase PBC to remove interactions with periodic images
         all_partcls = self.system.part.all()
         old_pos = all_partcls.pos.copy()
-        self.system.box_l = self.system.box_l * 2.
+        self.system.change_volume_and_rescale_particles(
+            2. * self.system.box_l[0], "xyz")
         all_partcls.pos = old_pos
         # compare calc_re()
         core_re = self.system.analysis.calc_re(chain_start=0,

testsuite/python/array_properties.py

@@ -102,13 +102,14 @@ class ArrayPropertyTest(ArrayCommon):
     system.box_l = [12.0, 12.0, 12.0]
     system.time_step = 0.01
     system.cell_system.skin = 0.01
-    partcl = system.part.add(pos=[0, 0, 0])
 
     def setUp(self):
         self.system.box_l = [12.0, 12.0, 12.0]
+        self.partcl = self.system.part.add(pos=[0, 0, 0])
 
     def tearDown(self):
         self.system.actors.clear()
+        self.system.part.clear()
 
     def assert_copy_is_writable(self, array):
         cpy = np.copy(array)

testsuite/python/constraint_shape_based.py

@@ -1073,14 +1073,26 @@ def test_torus(self):
 
     def test_exceptions(self):
         system = self.system
+        box_l = self.box_l
         wall = espressomd.shapes.Wall(normal=[0., 1., 0.], dist=0.)
         constraint = espressomd.constraints.ShapeBasedConstraint(
             shape=wall, particle_type=1)
         system.constraints.add(constraint)
         with self.assertRaisesRegex(Exception, "there are active constraints"):
             system.box_l = 0.5 * system.box_l
+        np.testing.assert_allclose(np.copy(system.box_l), box_l, atol=1e-7)
+        with self.assertRaisesRegex(RuntimeError, "there are active constraints"):
+            system.change_volume_and_rescale_particles(
+                0.5 * system.box_l[0], "xyz")
+        np.testing.assert_allclose(np.copy(system.box_l), box_l, atol=1e-7)
         system.constraints.remove(constraint)
         system.box_l = 0.75 * system.box_l
+        np.testing.assert_allclose(
+            np.copy(system.box_l), 0.75 * box_l, atol=1e-7)
+        system.change_volume_and_rescale_particles(
+            0.5 * system.box_l[0], "xyz")
+        np.testing.assert_allclose(
+            np.copy(system.box_l), 0.75 * 0.5 * box_l, atol=1e-7)
 
 
 if __name__ == "__main__":

testsuite/python/coulomb_interface.py

@@ -220,11 +220,11 @@ def test_elc_p3m_exceptions(self):
                 self.system.cell_system.node_grid,
                 self.original_node_grid)
         with self.assertRaisesRegex(Exception, "Exception while updating the box length: ERROR: ELC gap size .+ larger than box length in z-direction"):
-            self.system.box_l = [10., 10., 2.5]
-        self.system.box_l = [10., 10., 10.]
+            self.system.change_volume_and_rescale_particles(2.5, "z")
+        self.system.change_volume_and_rescale_particles(10., "z")
         self.system.actors.clear()
         with self.assertRaisesRegex(RuntimeError, "P3M real-space cutoff too large for ELC w/ dielectric contrast"):
-            self.system.box_l = [10., 10., 5.]
+            self.system.change_volume_and_rescale_particles(5., "z")
             elc = espressomd.electrostatics.ELC(
                 actor=p3m,
                 gap_size=1.,
@@ -237,7 +237,7 @@ def test_elc_p3m_exceptions(self):
             )
             self.system.actors.add(elc)
         self.assertEqual(len(self.system.actors), 0)
-        self.system.box_l = [10., 10., 10.]
+        self.system.change_volume_and_rescale_particles(10., "z")
         self.system.periodicity = [True, True, False]
         with self.assertRaisesRegex(RuntimeError, r"ELC: requires periodicity \(1 1 1\)"):
             elc = espressomd.electrostatics.ELC(

testsuite/python/coulomb_mixed_periodicity.py

@@ -41,6 +41,7 @@ def setUp(self):
         self.system.time_step = 0.01
         self.system.cell_system.skin = 0.
 
+    def setup_particles(self):
         # Add particles to system and store reference forces in hash
         # Input format: id pos q f
         self.system.part.add(pos=self.data[:, 1:4], q=self.data[:, 4])
@@ -64,17 +65,21 @@ def compare(self, method_name, force_tol, energy_tol):
         # triggering a solver re-initialization via a box resize
         # should not affect the forces nor the energies
         original_box_l = np.copy(self.system.box_l)
-        self.system.box_l = original_box_l * 1.1
-        self.system.box_l = original_box_l
+        for i in range(3):
+            self.system.change_volume_and_rescale_particles(
+                1.05 * original_box_l[i], "xyz"[i])
+        for i in range(3):
+            self.system.change_volume_and_rescale_particles(
+                original_box_l[i], "xyz"[i])
         self.system.integrator.run(0)
         forces_step2 = np.copy(self.system.part.all().f)
         energy_step2 = self.system.analysis.energy()["total"]
 
         err_msg = f"method {method_name} deviates after cells reinitialization"
-        np.testing.assert_allclose(forces_step1, forces_step2, atol=1e-12,
-                                   err_msg=f"Force {err_msg}")
-        np.testing.assert_allclose(energy_step2, energy_step1, rtol=1e-12,
-                                   err_msg=f"Energy {err_msg}")
+        np.testing.assert_allclose(forces_step1, forces_step2, rtol=1e-9,
+                                   atol=0., err_msg=f"Force {err_msg}")
+        np.testing.assert_allclose(energy_step2, energy_step1, rtol=1e-9,
+                                   atol=0., err_msg=f"Energy {err_msg}")
 
     def setup_elc_system(self):
         # Make sure, the data satisfies the gap
@@ -89,6 +94,7 @@ def setup_elc_system(self):
     @utx.skipIfMissingFeatures(["P3M"])
     def test_elc_cpu(self):
         self.system.box_l = [10., 10., 12.]
+        self.setup_particles()
         self.setup_elc_system()
 
         p3m = espressomd.electrostatics.P3M(
@@ -103,6 +109,7 @@ def test_elc_cpu(self):
     @utx.skipIfMissingFeatures(["P3M"])
     def test_elc_gpu(self):
         self.system.box_l = [10., 10., 12.]
+        self.setup_particles()
         self.setup_elc_system()
 
         p3m = espressomd.electrostatics.P3M(
@@ -119,6 +126,7 @@ def test_scafacos_p2nfft(self):
         self.system.box_l = [10., 10., 10.]
         self.system.periodicity = [1, 1, 0]
         self.system.cell_system.set_regular_decomposition()
+        self.setup_particles()
 
         scafacos = espressomd.electrostatics.Scafacos(
             prefactor=1,

testsuite/python/dds-and-bh-gpu.py

@@ -83,7 +83,7 @@ def test(self):
                 prefactor=pf_dds_gpu)
             system.actors.add(dds_gpu)
             # check MD cell reset has no impact
-            system.box_l = system.box_l
+            system.change_volume_and_rescale_particles(system.box_l[0], "x")
             system.periodicity = system.periodicity
             system.cell_system.node_grid = system.cell_system.node_grid
             system.integrator.run(steps=0, recalc_forces=True)
@@ -100,7 +100,7 @@ def test(self):
                 prefactor=pf_bh_gpu, epssq=200.0, itolsq=8.0)
             system.actors.add(bh_gpu)
             # check MD cell reset has no impact
-            system.box_l = system.box_l
+            system.change_volume_and_rescale_particles(system.box_l[0], "x")
             system.periodicity = system.periodicity
             system.cell_system.node_grid = system.cell_system.node_grid
             system.integrator.run(steps=0, recalc_forces=True)

testsuite/python/dipolar_direct_summation.py

@@ -48,7 +48,8 @@ def dds_gpu_data(self):
         dds_cpu = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1.2)
         system.actors.add(dds_cpu)
         # check MD cell reset has no impact
-        self.system.box_l = self.system.box_l
+        self.system.change_volume_and_rescale_particles(
+            self.system.box_l[0], "x")
         self.system.periodicity = self.system.periodicity
         self.system.cell_system.node_grid = self.system.cell_system.node_grid
 
@@ -67,7 +68,8 @@ def dds_data(self):
         dds_cpu = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1.2)
         system.actors.add(dds_cpu)
         # check MD cell reset has no impact
-        self.system.box_l = self.system.box_l
+        self.system.change_volume_and_rescale_particles(
+            self.system.box_l[0], "x")
         self.system.periodicity = self.system.periodicity
         self.system.cell_system.node_grid = self.system.cell_system.node_grid
 
@@ -87,7 +89,8 @@ def dds_replica_data(self):
             prefactor=1.2, n_replica=0)
         system.actors.add(dds_cpu)
         # check MD cell reset has no impact
-        self.system.box_l = self.system.box_l
+        self.system.change_volume_and_rescale_particles(
+            self.system.box_l[0], "x")
         self.system.periodicity = self.system.periodicity
         self.system.cell_system.node_grid = self.system.cell_system.node_grid
 

testsuite/python/dipolar_interface.py

@@ -141,21 +141,21 @@ def test_exceptions_serial(self):
             self.system.actors.add(mdlc)
         self.assertEqual(len(self.system.actors), 0)
         self.system.periodicity = [True, True, True]
-        self.system.box_l = [10., 10. + 2e-3, 10.]
+        self.system.change_volume_and_rescale_particles(10. + 2e-3, "y")
         with self.assertRaisesRegex(Exception, "box size in x direction is different from y direction"):
             mdlc = MDLC(gap_size=1., maxPWerror=1e-5, actor=ddsr)
             self.system.actors.add(mdlc)
         self.assertEqual(len(self.system.actors), 0)
         # check it's safe to resize the box, i.e. there are no currently
         # active sanity check in the core
-        self.system.box_l = [10., 10., 10.]
+        self.system.change_volume_and_rescale_particles(10., "y")
         with self.assertRaisesRegex(Exception, "box size in x direction is different from y direction"):
             ddsr = DDSR(prefactor=1., n_replica=1)
             mdlc = MDLC(gap_size=1., maxPWerror=1e-5, actor=ddsr)
             self.system.actors.add(mdlc)
-            self.system.box_l = [9., 10., 10.]
+            self.system.change_volume_and_rescale_particles(9., "x")
         self.system.actors.clear()
-        self.system.box_l = [10., 10., 10.]
+        self.system.change_volume_and_rescale_particles(10., "x")
 
     @utx.skipIfMissingFeatures(["DP3M"])
     def test_exceptions_parallel(self):

testsuite/python/dipolar_p3m.py

@@ -92,7 +92,8 @@ def test_dp3m(self):
         ref_mdlc_torque_metallic = np.copy(partcls.torque_lab)
 
         # actors should remain in a valid state after a cell system reset
-        self.system.box_l = self.system.box_l
+        self.system.change_volume_and_rescale_particles(
+            self.system.box_l[0], "x")
         self.system.periodicity = self.system.periodicity
         self.system.cell_system.node_grid = self.system.cell_system.node_grid
         self.system.integrator.run(0, recalc_forces=True)

testsuite/python/elc_vs_analytic.py

@@ -53,14 +53,13 @@ def test_elc(self):
         simulation box with dielectric contrast on the bottom of the box,
         which can be calculated analytically with image charges.
         """
-        self.system.part.add(pos=self.system.box_l / 2., q=self.q[0])
-        self.system.part.add(pos=self.system.box_l / 2. + [0, 0, self.distance],
-                             q=-self.q[0])
-
         self.system.box_l = [self.box_l, self.box_l, self.box_l + self.elc_gap]
         self.system.cell_system.set_regular_decomposition(
             use_verlet_lists=True)
         self.system.periodicity = [1, 1, 1]
+        self.system.part.add(pos=self.system.box_l / 2., q=self.q[0])
+        self.system.part.add(pos=self.system.box_l / 2. + [0, 0, self.distance],
+                             q=-self.q[0])
         prefactor = 2.0
         p3m = self.p3m_class(prefactor=prefactor, accuracy=self.accuracy,
                              mesh=[58, 58, 70], cao=4)