Disable box length reset when particles are present

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

@@ -105,12 +105,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/ek/EKNone.hpp

@@ -34,6 +34,7 @@ struct EKNone {
   void veto_time_step(double) const { throw NoEKActive{}; }
   void veto_kT(double) const { throw NoEKActive{}; }
   void sanity_checks(System::System const &) const { throw NoEKActive{}; }
+  void veto_boxl_change() const { throw NoEKActive{}; }
   void on_cell_structure_change() const { throw NoEKActive{}; }
   void on_boxl_change() const { throw NoEKActive{}; }
   void on_node_grid_change() const { throw NoEKActive{}; }

src/core/ek/EKWalberla.hpp

@@ -62,9 +62,10 @@ struct EKWalberla {
   void perform_reactions();
 
   void on_cell_structure_change() const {}
-  void on_boxl_change() const {
+  void veto_boxl_change() const {
     throw std::runtime_error("MD cell geometry change not supported by EK");
   }
+  void on_boxl_change() const { veto_boxl_change(); }
   void on_node_grid_change() const {
     throw std::runtime_error("MPI topology change not supported by EK");
   }

src/core/ek/Solver.cpp

@@ -90,6 +90,12 @@ void Solver::veto_kT(double kT) const {
   }
 }
 
+void Solver::veto_boxl_change() const {
+  if (impl->solver) {
+    std::visit([](auto const &ptr) { ptr->veto_boxl_change(); }, *impl->solver);
+  }
+}
+
 void Solver::on_cell_structure_change() {
   if (impl->solver) {
     auto &solver = *impl->solver;

src/core/ek/Solver.hpp

@@ -95,6 +95,7 @@ struct Solver : public System::Leaf<Solver> {
   void on_cell_structure_change();
   void on_timestep_change();
   void on_temperature_change();
+  void veto_boxl_change() const;
 
 private:
   /** @brief Pointer-to-implementation. */

src/core/lb/LBNone.hpp

@@ -52,6 +52,7 @@ struct LBNone {
   void lebc_sanity_checks(unsigned int, unsigned int) const {
     throw NoLBActive{};
   }
+  void veto_boxl_change() const { throw NoLBActive{}; }
   void on_cell_structure_change() const { throw NoLBActive{}; }
   void on_boxl_change() const { throw NoLBActive{}; }
   void on_node_grid_change() const { throw NoLBActive{}; }

src/core/lb/LBWalberla.hpp

@@ -74,9 +74,10 @@ struct LBWalberla {
                           unsigned int shear_plane_normal) const;
 
   void on_cell_structure_change() const {}
-  void on_boxl_change() const {
+  void veto_boxl_change() const {
     throw std::runtime_error("MD cell geometry change not supported by LB");
   }
+  void on_boxl_change() const { veto_boxl_change(); }
   void on_node_grid_change() const {
     throw std::runtime_error("MPI topology change not supported by LB");
   }

src/core/lb/Solver.cpp

@@ -107,6 +107,12 @@ void Solver::on_cell_structure_change() {
   }
 }
 
+void Solver::veto_boxl_change() const {
+  if (impl->solver) {
+    std::visit([](auto const &ptr) { ptr->veto_boxl_change(); }, *impl->solver);
+  }
+}
+
 void Solver::on_boxl_change() {
   if (impl->solver) {
     std::visit([](auto &ptr) { ptr->on_boxl_change(); }, *impl->solver);

src/core/lb/Solver.hpp

@@ -154,6 +154,7 @@ struct Solver : public System::Leaf<Solver> {
   void on_cell_structure_change();
   void on_timestep_change();
   void on_temperature_change();
+  void veto_boxl_change() const;
 
 private:
   /** @brief Pointer-to-implementation. */

src/core/system/System.cpp

@@ -167,6 +167,20 @@ void System::on_boxl_change(bool skip_method_adaption) {
   Constraints::constraints.on_boxl_change();
 }
 
+void System::veto_boxl_change(bool skip_particle_checks) const {
+  if (not skip_particle_checks) {
+    auto const n_part = boost::mpi::all_reduce(
+        ::comm_cart, cell_structure->local_particles().size(), std::plus<>());
+    if (n_part > 0ul) {
+      throw std::runtime_error(
+          "Cannot reset the box length when particles are present");
+    }
+  }
+  Constraints::constraints.veto_boxl_change();
+  lb.veto_boxl_change();
+  ek.veto_boxl_change();
+}
+
 void System::on_node_grid_change() {
   update_local_geo();
   lb.on_node_grid_change();

src/core/system/System.hpp

@@ -241,6 +241,7 @@ class System : public std::enable_shared_from_this<System> {
    *  initialized immediately (P3M etc.).
    */
   void on_observable_calc();
+  void veto_boxl_change(bool skip_particle_checks = false) const;
   /**@}*/
 
   /**

src/core/unit_tests/ek_interface_test.cpp

@@ -188,6 +188,7 @@ BOOST_AUTO_TEST_CASE(ek_interface_walberla) {
       espresso::ek_container->add(ek_species);
       BOOST_CHECK_THROW(ek.veto_kT(params.kT + 1.), std::runtime_error);
       BOOST_CHECK_THROW(ek.on_boxl_change(), std::runtime_error);
+      BOOST_CHECK_THROW(ek.veto_boxl_change(), std::runtime_error);
       BOOST_CHECK_THROW(ek.veto_time_step(ek.get_tau() * 2.),
                         std::invalid_argument);
       BOOST_CHECK_THROW(ek.veto_time_step(ek.get_tau() / 2.5),
@@ -215,6 +216,7 @@ BOOST_AUTO_TEST_CASE(ek_interface_none) {
     BOOST_CHECK_THROW(ek.propagate(), NoEKActive);
     BOOST_CHECK_THROW(ek.get_tau(), NoEKActive);
     BOOST_CHECK_THROW(ek.sanity_checks(), NoEKActive);
+    BOOST_CHECK_THROW(ek.veto_boxl_change(), NoEKActive);
     BOOST_CHECK_THROW(ek.veto_time_step(0.), NoEKActive);
     BOOST_CHECK_THROW(ek.veto_kT(0.), NoEKActive);
     BOOST_CHECK_THROW(ek.on_cell_structure_change(), NoEKActive);

src/core/unit_tests/lb_particle_coupling_test.cpp

@@ -520,6 +520,7 @@ BOOST_AUTO_TEST_CASE(runtime_exceptions) {
   // LB prevents changing most of the system state
   {
     BOOST_CHECK_THROW(lb.on_boxl_change(), std::runtime_error);
+    BOOST_CHECK_THROW(lb.veto_boxl_change(), std::runtime_error);
     BOOST_CHECK_THROW(lb.veto_time_step(lb.get_tau() * 2.),
                       std::invalid_argument);
     BOOST_CHECK_THROW(lb.veto_time_step(lb.get_tau() / 2.5),
@@ -618,6 +619,7 @@ BOOST_AUTO_TEST_CASE(lb_exceptions) {
     BOOST_CHECK_THROW(lb.get_pressure_tensor(), NoLBActive);
     BOOST_CHECK_THROW(lb.get_momentum(), NoLBActive);
     BOOST_CHECK_THROW(lb.sanity_checks(), NoLBActive);
+    BOOST_CHECK_THROW(lb.veto_boxl_change(), NoLBActive);
     BOOST_CHECK_THROW(lb.veto_time_step(0.), NoLBActive);
     BOOST_CHECK_THROW(lb.veto_kT(0.), NoLBActive);
     BOOST_CHECK_THROW(lb.lebc_sanity_checks(0u, 1u), NoLBActive);

src/script_interface/system/System.cpp

@@ -111,6 +111,7 @@ System::System() : m_instance{}, m_leaves{std::make_shared<Leaves>()} {
            if (not(new_value > Utils::Vector3d::broadcast(0.))) {
              throw std::domain_error("Attribute 'box_l' must be > 0");
            }
+           m_instance->veto_boxl_change();
            m_instance->box_geo->set_length(new_value);
            m_instance->on_boxl_change();
          });
@@ -256,12 +257,15 @@ Variant System::do_call_method(std::string const &name,
     auto &box_geo = *m_instance->box_geo;
     auto const coord = get_value<int>(parameters, "coord");
     auto const length = get_value<double>(parameters, "length");
+    assert(coord != 3 or ((box_geo.length()[0] == box_geo.length()[1]) and
+                          (box_geo.length()[1] == box_geo.length()[2])));
     auto const scale = (coord == 3) ? length * box_geo.length_inv()[0]
                                     : length * box_geo.length_inv()[coord];
     context()->parallel_try_catch([&]() {
       if (length <= 0.) {
         throw std::domain_error("Parameter 'd_new' must be > 0");
       }
+      m_instance->veto_boxl_change(true);
     });
     auto new_value = Utils::Vector3d{};
     if (coord == 3) {

testsuite/python/analyze_chains.py

@@ -126,7 +126,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,14 +102,15 @@ 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):
         if espressomd.has_features("WALBERLA"):
             self.system.lb = None
+        self.system.part.clear()
 
     def assert_copy_is_writable(self, array):
         cpy = np.copy(array)

testsuite/python/constraint_shape_based.py

@@ -1072,14 +1072,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(RuntimeError, "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

@@ -253,12 +253,12 @@ def test_elc_p3m_exceptions(self):
             self.assertEqual(
                 list(self.system.cell_system.node_grid),
                 list(self.original_node_grid))
-        with self.assertRaisesRegex(Exception, "while setting parameter 'box_l': 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.]
+        with self.assertRaisesRegex(Exception, "ERROR: ELC gap size .+ larger than box length in z-direction"):
+            self.system.change_volume_and_rescale_particles(2.5, "z")
+        self.system.change_volume_and_rescale_particles(10., "z")
         self.system.electrostatics.solver = None
         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.,
@@ -271,7 +271,7 @@ def test_elc_p3m_exceptions(self):
             )
             self.system.electrostatics.solver = elc
         self.assertIsNone(self.system.electrostatics.solver)
-        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, periodicity_err_msg):
             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 = [True, True, False]
         self.system.cell_system.set_regular_decomposition()
+        self.setup_particles()
 
         scafacos = espressomd.electrostatics.Scafacos(
             prefactor=1,

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

@@ -80,7 +80,7 @@ def test(self):
                 prefactor=pf_dds_gpu)
             system.magnetostatics.solver = 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)
@@ -97,7 +97,7 @@ def test(self):
                 prefactor=pf_bh_gpu, epssq=200.0, itolsq=8.0)
             system.magnetostatics.solver = 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

@@ -49,7 +49,8 @@ def dds_gpu_data(self):
         dds_cpu = espressomd.magnetostatics.DipolarDirectSumGpu(prefactor=1.2)
         system.magnetostatics.solver = 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
 
@@ -68,7 +69,8 @@ def dds_data(self):
         dds_cpu = espressomd.magnetostatics.DipolarDirectSumCpu(prefactor=1.2)
         system.magnetostatics.solver = 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

@@ -112,7 +112,7 @@ def test_exceptions_non_p3m(self):
             self.system.magnetostatics.solver = mdlc
         self.assertIsNone(self.system.magnetostatics.solver)
         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.magnetostatics.solver = mdlc
@@ -128,14 +128,14 @@ def test_exceptions_non_p3m(self):
             self.system.magnetostatics.clear()
         # 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_replicas=1)
             mdlc = MDLC(gap_size=1., maxPWerror=1e-5, actor=ddsr)
             self.system.magnetostatics.solver = mdlc
-            self.system.box_l = [9., 10., 10.]
+            self.system.change_volume_and_rescale_particles(9., "x")
         self.system.magnetostatics.clear()
-        self.system.box_l = [10., 10., 10.]
+        self.system.change_volume_and_rescale_particles(10., "x")
 
     @utx.skipIfMissingFeatures(["DP3M"])
     def test_exceptions_p3m(self):