Implemented virial contiribution for rigid bond

src/core/bonded_interactions/bonded_interaction_data.hpp

@@ -48,13 +48,16 @@
 #include "Particle.hpp"
 #include "system/Leaf.hpp"
 
+#include <utils/Vector.hpp>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <memory>
 #include <optional>
 #include <unordered_map>
 #include <variant>
+#include <vector>
 
 /** Interaction type for unused bonded interaction slots */
 struct NoneBond {
@@ -144,6 +147,17 @@ class BondedInteractionsMap : public System::Leaf<BondedInteractionsMap> {
     assert(n_rigid_bonds >= 0);
     return n_rigid_bonds;
   }
+  /**
+   * @brief Per-bond-type RATTLE constraint virial.
+   *
+   * Accumulated bond-by-bond inside correct_position_shake() (indexed by
+   * bond id, same convention as Observable_stat::bonded_contribution()),
+   * where the pairwise correction/mass/bond-vector are all unambiguous,
+   * regardless of how many rigid bonds a particle participates in. Reset
+   * at the start of every SHAKE call; consumed by
+   * System::calculate_pressure().
+   */
+  std::vector<Utils::Vector9d> rigid_bond_virial;
 #endif
   std::optional<key_type> find_bond_id(mapped_type const &target_bond) const {
     for (auto const &[bond_id, bond] : m_params) {

src/core/integrate.cpp

@@ -257,6 +257,15 @@ void System::System::integrator_sanity_checks() const {
     runtimeErrorMsg()
         << "Thermalized bonds require the thermalized_bond thermostat";
   }
+#ifdef ESPRESSO_BOND_CONSTRAINT
+  if (bonded_ias->get_n_rigid_bonds() >= 1) {
+    if (not propagation->is_inertial()) {
+      runtimeErrorMsg()
+          << "Rigid bonds (RATTLE) require an inertial integrator "
+             "(VV or symplectic Euler); BD and SD are not supported";
+    }
+  }
+#endif
 
 #ifdef ESPRESSO_ROTATION
   for (auto const &p : cell_structure->local_particles()) {

src/core/integrators/Propagation.hpp

@@ -48,4 +48,12 @@ class Propagation {
     recalc_forces = true;
     recalc_used_propagations = true;
   }
+
+  /** True for all integrators that use inertial dynamics (VV, SE, NPT).
+   *  False for non-inertial ones (BD, SD, steepest descent), which cannot
+   *  use RATTLE rigid bonds. */
+  bool is_inertial() const {
+    return integ_switch != INTEG_METHOD_STEEPEST_DESCENT &&
+           integ_switch != INTEG_METHOD_BD && integ_switch != INTEG_METHOD_SD;
+  }
 };

src/core/pressure.cpp

@@ -40,6 +40,7 @@
 #include "virtual_sites/relative.hpp"
 
 #include <utils/Vector.hpp>
+#include <utils/math/sqr.hpp>
 #include <utils/matrix.hpp>
 
 #include <algorithm>
@@ -151,6 +152,24 @@ std::shared_ptr<Observable_stat> System::calculate_pressure() {
   }
 #endif
 
+#ifdef ESPRESSO_BOND_CONSTRAINT
+  if (propagation->is_inertial() and bonded_ias->get_n_rigid_bonds() >= 1) {
+    // rigid_bond_virial was accumulated bond-by-bond inside
+    // correct_position_shake() during the last integration step, so it is
+    // already correct regardless of how many rigid bonds a particle
+    // participates in; only the deferred 1/dt^2 factor is applied here.
+    auto const sq_dt = Utils::sqr(get_time_step());
+    auto const &rigid_bond_virial = bonded_ias->rigid_bond_virial;
+    for (std::size_t bond_id = 0; bond_id < rigid_bond_virial.size();
+         ++bond_id) {
+      auto const stress = rigid_bond_virial[bond_id] / sq_dt;
+      auto dest = obs_pressure.bonded_contribution(static_cast<int>(bond_id));
+      for (std::size_t k = 0; k < 9u; ++k)
+        dest[k] += stress[k];
+    }
+  }
+#endif
+
   obs_pressure.rescale(volume);
 
   obs_pressure.mpi_reduce();

src/core/rattle.cpp

@@ -32,13 +32,19 @@
 #include "communication.hpp"
 #include "errorhandling.hpp"
 
+#include <utils/Vector.hpp>
+#include <utils/math/tensor_product.hpp>
+#include <utils/matrix.hpp>
+
 #include <boost/mpi/collectives/all_reduce.hpp>
 #include <boost/range/algorithm.hpp>
 
 #include <cmath>
+#include <cstddef>
 #include <functional>
 #include <span>
 #include <variant>
+#include <vector>
 
 /** Maximal number of iterations before the RATTLE algorithm bails out. */
 static constexpr auto shake_max_iterations = 1000;
@@ -85,11 +91,15 @@ static void init_correction_vector(const ParticleRange &particles,
  * @param box_geo Box geometry.
  * @param p1 First particle.
  * @param p2 Second particle.
+ * @param bond_id Bonded interaction id of this specific bond.
+ * @param rigid_bond_virial Per-bond-type RATTLE constraint virial
+ *                          accumulator, indexed by @p bond_id.
  * @return True if there was a correction.
  */
-static bool calculate_positional_correction(RigidBond const &ia_params,
-                                            BoxGeometry const &box_geo,
-                                            Particle &p1, Particle &p2) {
+static bool calculate_positional_correction(
+    RigidBond const &ia_params, BoxGeometry const &box_geo, Particle &p1,
+    Particle &p2, int bond_id,
+    std::vector<Utils::Vector9d> &rigid_bond_virial) {
   auto const r_ij = box_geo.get_mi_vector(p1.pos(), p2.pos());
   auto const r_ij2 = r_ij.norm2();
 
@@ -104,6 +114,17 @@ static bool calculate_positional_correction(RigidBond const &ia_params,
     p1.rattle_params().correction += pos_corr * p2.mass();
     p2.rattle_params().correction -= pos_corr * p1.mass();
 
+    // Constraint force implied by this bond alone during this iteration:
+    // the correction just applied to p1 is Δr1 = pos_corr*m2 = (1/2)*a1*dt²,
+    // so F1 = m1*a1 = 2*m1*Δr1/dt² = 2*m1*m2*pos_corr/dt². Division by dt²
+    // is deferred to System::calculate_pressure(), where the timestep is
+    // available. This uses r_ij_t, the bond vector at the start of the MD
+    // step (fixed across all SHAKE iterations of this step), so the
+    // contribution is exact for this bond alone, regardless of how many
+    // other rigid bonds p1 or p2 participate in.
+    rigid_bond_virial[static_cast<std::size_t>(bond_id)] += Utils::flatten(
+        Utils::tensor_product(2.0 * p1.mass() * p2.mass() * pos_corr, r_ij_t));
+
     return true;
   }
 
@@ -130,7 +151,7 @@ static bool compute_correction_vector(CellStructure &cs,
     auto const &iaparams = *bonded_ias.at(bond_id);
 
     if (auto const *bond = std::get_if<RigidBond>(&iaparams)) {
-      auto const corrected = kernel(*bond, box_geo, p1, *partners[0]);
+      auto const corrected = kernel(*bond, box_geo, p1, *partners[0], bond_id);
       if (corrected)
         correction = true;
     }
@@ -155,18 +176,28 @@ static void apply_positional_correction(const ParticleRange &particles) {
 }
 
 void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo,
-                            BondedInteractionsMap const &bonded_ias) {
+                            BondedInteractionsMap &bonded_ias) {
   unsigned const flag = Cells::DATA_PART_POSITION | Cells::DATA_PART_PROPERTIES;
   cs.update_ghosts_and_resort_particle(flag);
 
   auto particles = cs.local_particles();
   auto ghost_particles = cs.ghost_particles();
 
+  // Reset the per-bond-type constraint virial for this timestep
+  bonded_ias.rigid_bond_virial.assign(
+      static_cast<std::size_t>(bonded_ias.get_next_key()),
+      Utils::Vector9d::broadcast(0.));
+
   int cnt;
   for (cnt = 0; cnt < shake_max_iterations; ++cnt) {
     init_correction_vector(particles, ghost_particles);
     bool const repeat_ = compute_correction_vector(
-        cs, box_geo, bonded_ias, calculate_positional_correction);
+        cs, box_geo, bonded_ias,
+        [&bonded_ias](RigidBond const &bond, BoxGeometry const &box_geo_,
+                      Particle &p1, Particle &p2, int bond_id) {
+          return calculate_positional_correction(
+              bond, box_geo_, p1, p2, bond_id, bonded_ias.rigid_bond_virial);
+        });
     bool const repeat =
         boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or<bool>());
 
@@ -237,7 +268,11 @@ void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo,
   for (cnt = 0; cnt < shake_max_iterations; ++cnt) {
     init_correction_vector(particles, ghost_particles);
     bool const repeat_ = compute_correction_vector(
-        cs, box_geo, bonded_ias, calculate_velocity_correction);
+        cs, box_geo, bonded_ias,
+        [](RigidBond const &bond, BoxGeometry const &box_geo_, Particle &p1,
+           Particle &p2, int /* bond_id */) {
+          return calculate_velocity_correction(bond, box_geo_, p1, p2);
+        });
     bool const repeat =
         boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or<bool>());
 

src/core/rattle.hpp

@@ -46,9 +46,12 @@ void save_old_position(ParticleRange const &particles,
 /**
  * @brief Propagate velocity and position while using SHAKE algorithm for bond
  * constraint.
+ *
+ * Also accumulates the per-bond-type RATTLE constraint virial into
+ * @ref BondedInteractionsMap::rigid_bond_virial "bonded_ias.rigid_bond_virial".
  */
 void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo,
-                            BondedInteractionsMap const &bonded_ias);
+                            BondedInteractionsMap &bonded_ias);
 
 /**
  * @brief Correction of current velocities using RATTLE algorithm.

testsuite/python/CMakeLists.txt

@@ -264,6 +264,8 @@ python_test(FILE long_range_actors.py MAX_NUM_PROC 4 GPU_SLOTS 2)
 python_test(FILE tabulated.py MAX_NUM_PROC 2)
 python_test(FILE particle_slice.py MAX_NUM_PROC 4)
 python_test(FILE rigid_bond.py MAX_NUM_PROC 4)
+python_test(FILE rigid_bond_virial.py MAX_NUM_PROC 4)
+python_test(FILE rigid_bond_virial_stats.py MAX_NUM_PROC 4 LABELS long)
 python_test(FILE rotation_per_particle.py MAX_NUM_PROC 4)
 python_test(FILE rotational_inertia.py MAX_NUM_PROC 2)
 python_test(FILE rotational-diffusion-aniso.py MAX_NUM_PROC 1 LABELS long)

testsuite/python/integrator_exceptions.py

@@ -278,6 +278,41 @@ def test_temperature_change(self):
         with self.assertRaisesRegex(RuntimeError, f"Parameter 'kT' is read-only"):
             self.system.thermostat.kT = 2.
 
+    @utx.skipIfMissingFeatures("BOND_CONSTRAINT")
+    def test_rigid_bond_incompatible_integrators(self):
+        system = self.system
+        system.part.clear()
+        system.box_l = [4., 4., 4.]
+        system.cell_system.skin = 0.4
+        system.time_step = 0.01
+        bond = espressomd.interactions.RigidBond(r=1.0, ptol=1e-4, vtol=1e-4)
+        system.bonded_inter.add(bond)
+        p1 = system.part.add(pos=[1.5, 2.0, 2.0])
+        p2 = system.part.add(pos=[2.5, 2.0, 2.0])
+        p2.add_bond((bond, p1))
+        error_msg = self.msg + \
+            'Rigid bonds \\(RATTLE\\) require an inertial integrator'
+
+        system.integrator.set_brownian_dynamics()
+        system.thermostat.set_brownian(kT=1.0, gamma=1.0, seed=42)
+        with self.assertRaisesRegex(Exception, error_msg):
+            system.integrator.run(0)
+
+        system.thermostat.turn_off()
+        system.integrator.set_steepest_descent(
+            f_max=0., gamma=1., max_displacement=0.1)
+        with self.assertRaisesRegex(Exception, error_msg):
+            system.integrator.run(0)
+
+        if espressomd.has_features("STOKESIAN_DYNAMICS"):
+            system.thermostat.turn_off()
+            system.periodicity = 3 * [False]
+            system.thermostat.set_stokesian(kT=0)
+            system.integrator.set_stokesian_dynamics(
+                viscosity=1.0, radii={0: 1.0})
+            with self.assertRaisesRegex(Exception, error_msg):
+                system.integrator.run(0)
+
     def test_missing_features(self):
         thermostats = {
             "STOKESIAN_DYNAMICS": ("set_stokesian", "Stokesian"),