From 3986556e33ca5ddecc8fe83ad9fab7a1eef8325d Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Mon, 29 Jun 2026 16:49:14 +0200 Subject: [PATCH 1/6] Implemented virial contiribution for rigid bond --- src/core/Particle.hpp | 12 +- src/core/ghosts.cpp | 6 +- src/core/integrate.cpp | 9 ++ src/core/integrators/Propagation.hpp | 9 ++ src/core/pressure.cpp | 30 +++++ src/core/rattle.cpp | 10 ++ src/core/unit_tests/Particle_test.cpp | 56 ++++++++ testsuite/python/CMakeLists.txt | 1 + testsuite/python/integrator_exceptions.py | 34 +++++ testsuite/python/rigid_bond_virial.py | 148 ++++++++++++++++++++++ 10 files changed, 310 insertions(+), 5 deletions(-) create mode 100644 testsuite/python/rigid_bond_virial.py diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 7c2ac3df91..4db73fb80b 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -415,18 +415,26 @@ struct ParticleLocal { struct ParticleRattle { /** position/velocity correction */ Utils::Vector3d correction = {0., 0., 0.}; + /** total over all RATTLE iterations */ + Utils::Vector3d accumulated_correction = {0., 0., 0.}; friend ParticleRattle operator+(ParticleRattle const &lhs, ParticleRattle const &rhs) { - return {lhs.correction + rhs.correction}; + return {lhs.correction + rhs.correction, + lhs.accumulated_correction + rhs.accumulated_correction}; } + // Intentionally excludes accumulated_correction: ghost reduction must not + // overwrite the local particle's accumulated value with ghost data. ParticleRattle &operator+=(ParticleRattle const &rhs) { - return *this = *this + rhs; + correction += rhs.correction; + return *this; } template void serialize(Archive &ar, long int /* version */) { ar & correction; + // accumulated_correction is reset at the start of every SHAKE call; + // serializing it would waste checkpoint space with a transient value. } }; #endif diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index ded254f9e5..14d2415fdd 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -354,9 +354,9 @@ add_rattle_correction_from_recv_buffer(CommBuf &recv_buffer, auto archiver = Utils::MemcpyIArchive{recv_buffer.make_span()}; for (auto &part_list : ghost_comm.part_lists) { for (Particle &part : *part_list) { - ParticleRattle pr; - archiver >> pr; - part.rattle_params() += pr; + Utils::Vector3d correction; + archiver >> correction; + part.rattle_correction() += correction; } } } diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index d955de7ce7..b6fd5609e6 100644 --- a/src/core/integrate.cpp +++ b/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()) { diff --git a/src/core/integrators/Propagation.hpp b/src/core/integrators/Propagation.hpp index 76de5ee327..db62b2e465 100644 --- a/src/core/integrators/Propagation.hpp +++ b/src/core/integrators/Propagation.hpp @@ -48,4 +48,13 @@ 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; + } }; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 2130adf849..fdfff9b72c 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -151,6 +151,36 @@ std::shared_ptr System::calculate_pressure() { } #endif +#ifdef ESPRESSO_BOND_CONSTRAINT + if (propagation->is_inertial() and bonded_ias->get_n_rigid_bonds() >= 1) { + auto const dt = get_time_step(); + auto const sq_dt = dt * dt; + cell_structure->bond_loop( + [&obs_pressure, &bonded_ias = *bonded_ias, &box_geo = *box_geo, + sq_dt](Particle &p1, int bond_id, std::span partners) { + if (std::holds_alternative(*bonded_ias.at(bond_id))) { + auto const &p2 = *partners[0]; + // Recover pre-SHAKE bond vector from local accumulated_correction + // only: accum2 = -(m1/m2)*accum1, so + // r12_pre = r12_post - accum1 + accum2 + // = r12_post - (m1+m2)/m2 * accum1 + // This avoids reading accumulated_correction from ghost p2. + auto const &accum = p1.rattle_params().accumulated_correction; + auto const r12 = + box_geo.get_mi_vector(p1.pos(), p2.pos()) - + ((p1.mass() + p2.mass()) / p2.mass()) * accum; + auto const F_c = 2.0 * accum * p1.mass() / sq_dt; + auto const stress = + Utils::flatten(Utils::tensor_product(F_c, r12)); + auto dest = obs_pressure.bonded_contribution(bond_id); + for (std::size_t k = 0; k < 9u; ++k) + dest[k] += stress[k]; + } + return false; + }); + } +#endif + obs_pressure.rescale(volume); obs_pressure.mpi_reduce(); diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index 4a738807e9..0171d29d6a 100644 --- a/src/core/rattle.cpp +++ b/src/core/rattle.cpp @@ -162,6 +162,11 @@ void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, auto particles = cs.local_particles(); auto ghost_particles = cs.ghost_particles(); + // Reset the accumulator for this timestep + boost::for_each(particles, [](Particle &p) { + p.rattle_params().accumulated_correction.fill(0.); + }); + int cnt; for (cnt = 0; cnt < shake_max_iterations; ++cnt) { init_correction_vector(particles, ghost_particles); @@ -176,6 +181,11 @@ void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, cs.ghosts_reduce_rattle_correction(); + // Accumulate AFTER ghost reduction so cross-rank bonds are included + boost::for_each(particles, [](Particle &p) { + p.rattle_params().accumulated_correction += p.rattle_params().correction; + }); + apply_positional_correction(particles); cs.ghosts_update(Cells::DATA_PART_POSITION | Cells::DATA_PART_MOMENTUM); } diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp index 8f5948af09..a7706e6a41 100644 --- a/src/core/unit_tests/Particle_test.cpp +++ b/src/core/unit_tests/Particle_test.cpp @@ -201,6 +201,8 @@ void check_particle_rattle(ParticleRattle const &out, ParticleRattle const &ref) { BOOST_TEST(out.correction == ref.correction, boost::test_tools::per_element()); + BOOST_TEST(out.accumulated_correction == ref.accumulated_correction, + boost::test_tools::per_element()); } BOOST_AUTO_TEST_CASE(rattle_serialization) { @@ -232,6 +234,27 @@ BOOST_AUTO_TEST_CASE(rattle_serialization) { } } +BOOST_AUTO_TEST_CASE(rattle_checkpoint_serialization) { + // accumulated_correction is transient; it must NOT survive a checkpoint + // round-trip so that restored simulations start each SHAKE call clean. + auto pr = ParticleRattle{{1., 2., 3.}, {4., 5., 6.}}; + + std::stringstream stream; + { + boost::archive::text_oarchive oa(stream); + oa << pr; + } + + boost::archive::text_iarchive ia(stream); + ParticleRattle out; + ia >> out; + + BOOST_TEST(out.correction == (Utils::Vector3d{1., 2., 3.}), + boost::test_tools::per_element()); + BOOST_TEST(out.accumulated_correction == Utils::Vector3d{}, + boost::test_tools::per_element()); +} + BOOST_AUTO_TEST_CASE(rattle_constructors) { auto pr = ParticleRattle{{1, 2, 3}}; @@ -248,6 +271,39 @@ BOOST_AUTO_TEST_CASE(rattle_constructors) { check_particle_rattle(out, pr); } } + +BOOST_AUTO_TEST_CASE(rattle_accumulated_correction) { + // default-initialized accumulated_correction must be zero + ParticleRattle pr{}; + auto const zero = Utils::Vector3d{}; + BOOST_TEST(pr.accumulated_correction == zero, + boost::test_tools::per_element()); + + // operator+ must add both fields independently + auto const corr_a = Utils::Vector3d{1., 0., 0.}; + auto const accum_a = Utils::Vector3d{0., 1., 0.}; + auto const corr_b = Utils::Vector3d{0., 0., 1.}; + auto const accum_b = Utils::Vector3d{1., 0., 0.}; + ParticleRattle a{corr_a, accum_a}; + ParticleRattle b{corr_b, accum_b}; + + auto const expected_corr = corr_a + corr_b; // {1,0,1} + auto const expected_accum = accum_a + accum_b; // {1,1,0} + ParticleRattle expected_sum{expected_corr, expected_accum}; + check_particle_rattle(a + b, expected_sum); + + // operator+= must add only correction, NOT accumulated_correction + // (ghost reduction must not bleed stale accumulated values into locals) + auto const dst_accum = Utils::Vector3d{5., 5., 5.}; + ParticleRattle dst{Utils::Vector3d{2., 0., 0.}, dst_accum}; + dst += b; + auto const expected_dst_corr = Utils::Vector3d{2., 0., 1.}; + BOOST_TEST(dst.correction == expected_dst_corr, + boost::test_tools::per_element()); + // accumulated_correction must be unchanged after += + BOOST_TEST(dst.accumulated_correction == dst_accum, + boost::test_tools::per_element()); +} #endif // ESPRESSO_BOND_CONSTRAINT #ifdef ESPRESSO_THERMAL_STONER_WOHLFARTH diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 951e2c19ec..016c1d7a6a 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -264,6 +264,7 @@ 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 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) diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py index a7a9974552..729506e01e 100644 --- a/testsuite/python/integrator_exceptions.py +++ b/testsuite/python/integrator_exceptions.py @@ -278,6 +278,40 @@ 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"), diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py new file mode 100644 index 0000000000..ee35a2b52c --- /dev/null +++ b/testsuite/python/rigid_bond_virial.py @@ -0,0 +1,148 @@ +# +# Copyright (C) 2013-2026 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import unittest as ut +import unittest_decorators as utx +import espressomd +import espressomd.interactions +import numpy as np + +@utx.skipIfMissingFeatures("BOND_CONSTRAINT") +class RigidBondVirialTest(ut.TestCase): + + system = espressomd.System(box_l=[10., 10., 10.]) + system.cell_system.skin = 0.4 + + def setUp(self): + self.system.time_step = 1e-2 + + def tearDown(self): + self.system.non_bonded_inter.reset() + self.system.part.clear() + self.system.thermostat.turn_off() + + def _make_dimer(self, m1=1.0, m2=1.0, v=1.0): + x1 = 5. - m2 / (m1 + m2) + x2 = 5. + m1 / (m1 + m2) + v1 = (0., v, 0.) + v2 = (0., -(m1 / m2) * v, 0.) + bond = espressomd.interactions.RigidBond(r=1.0, ptol=1e-4, vtol=1e-4) + self.system.bonded_inter.add(bond) + p1 = self.system.part.add(pos=[x1, 5.0, 5.0], v=v1, mass=m1) + p2 = self.system.part.add(pos=[x2, 5.0, 5.0], v=v2, mass=m2) + p2.add_bond((bond, p1)) + return p1, p2 + + # Rotation test: after one step, constraint virial = F_centripetal·d + # F_centripetal = m v²/r = 1·1/0.5 = 2, W = -virial/(3V) = -2/(3V) + def _virial_in_rotation(self, set_integrator): + V = self.system.volume() + set_integrator() + self._make_dimer() + self.system.integrator.run(1) + v_p = np.trace(self.system.analysis.pressure_tensor()['bonded']) / 3. + v_theory = -2.0 / (3. * V) + self.assertAlmostEqual(v_p, v_theory, delta=0.01 * abs(v_theory)) + v_xx = self.system.analysis.pressure_tensor()['bonded'][0,0] + v_yy = self.system.analysis.pressure_tensor()['bonded'][1,1] + v_zz = self.system.analysis.pressure_tensor()['bonded'][2,2] + self.assertAlmostEqual(v_xx, -2.0 / V, delta=0.01 * abs(2.0 / V)) + self.assertAlmostEqual(v_yy, 0.0, delta=1e-8) + self.assertAlmostEqual(v_zz, 0.0, delta=1e-8) + + def test_virial_in_rotation_vv(self): + """VV: constraint virial of a rotating rigid dimer matches centripetal theory.""" + self._virial_in_rotation(self.system.integrator.set_vv) + + def test_virial_in_rotation_se(self): + """SE: constraint virial of a rotating rigid dimer matches centripetal theory.""" + self._virial_in_rotation(self.system.integrator.set_symplectic_euler) + + def _virial_unequal_masses(self, set_integrator): + m1 = 1.0 + m2 = 2.0 + d = 1.0 + v = 1.0 # velocity + V = self.system.volume() + + set_integrator() + self._make_dimer(m1=m1, m2=m2, v=v) + self.system.integrator.run(1) + + mu = m1 * m2 / (m1 + m2) # reduced mass = 2/3 + omega = v * (m1 + m2) / (m2 * d) # |v_rel| / d = 3/2 + v_theory = -mu * omega**2 * d**2 / (3. * V) # = -1/(2V) + + pt = self.system.analysis.pressure_tensor()['bonded'] + v_p = np.trace(pt) / 3. + self.assertAlmostEqual(v_p, v_theory, delta=0.01 * abs(v_theory)) + # Bond is along x: all virial goes into xx, none into yy or zz + self.assertAlmostEqual(pt[0, 0], 3. * v_theory, delta=0.01 * abs(3. * v_theory)) + self.assertAlmostEqual(pt[1, 1], 0., delta=1e-8) + self.assertAlmostEqual(pt[2, 2], 0., delta=1e-8) + + def test_virial_unequal_masses_vv(self): + """VV: constraint virial with m1!=m2 matches centripetal theory.""" + self._virial_unequal_masses(self.system.integrator.set_vv) + + def test_virial_unequal_masses_se(self): + """SE: constraint virial with m1!=m2 matches centripetal theory.""" + self._virial_unequal_masses(self.system.integrator.set_symplectic_euler) + + # Langevin consistency: mean of (P_bond + P_kin) = kT/V, + # std matches analytic fluctuation formula. + def _virial_consistency(self, set_integrator, noise_prefactor): + kT = 1.0 + gamma = 1.0 + mass = 1. + V = self.system.volume() + dt = self.system.time_step + self.system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42) + set_integrator() + + std_theory = ((6 * kT**2 + noise_prefactor * gamma * mass * kT / dt) / (9 * V**2))**0.5 + + self._make_dimer(v=0.0) + + self.system.integrator.run(1000) # equilibrate + n_loop = 2000 + n_steps = 100 + virial = [] + for _ in range(n_loop): + self.system.integrator.run(n_steps) + v_p = np.trace(self.system.analysis.pressure_tensor()['bonded']) / 3. + v_k = np.trace(self.system.analysis.pressure_tensor()['kinetic']) / 3. + virial.append(v_p + v_k) + + rigid_p = np.mean(virial) + rigid_std = np.std(virial) + self.assertAlmostEqual(rigid_p, 1. / V, delta=2.*std_theory/n_loop**0.5) + self.assertAlmostEqual(rigid_std, std_theory, delta=0.02*std_theory) + + def test_virial_consistency_vv(self): + """VV+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" + self._virial_consistency(self.system.integrator.set_vv, 1) + + def test_virial_consistency_se(self): + """SE+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" + self._virial_consistency(self.system.integrator.set_symplectic_euler, 4) + + +if __name__ == "__main__": + ut.main() + From ce4e20638858129fb000119a63c220d8a374132c Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Mon, 29 Jun 2026 17:37:34 +0200 Subject: [PATCH 2/6] Fixed style --- src/core/integrators/Propagation.hpp | 3 +-- src/core/pressure.cpp | 8 +++--- testsuite/python/integrator_exceptions.py | 3 ++- testsuite/python/rigid_bond_virial.py | 33 ++++++++++++++--------- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/core/integrators/Propagation.hpp b/src/core/integrators/Propagation.hpp index db62b2e465..0a39c4b721 100644 --- a/src/core/integrators/Propagation.hpp +++ b/src/core/integrators/Propagation.hpp @@ -54,7 +54,6 @@ class Propagation { * 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; + integ_switch != INTEG_METHOD_BD && integ_switch != INTEG_METHOD_SD; } }; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index fdfff9b72c..c07feff2bf 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -166,12 +166,10 @@ std::shared_ptr System::calculate_pressure() { // = r12_post - (m1+m2)/m2 * accum1 // This avoids reading accumulated_correction from ghost p2. auto const &accum = p1.rattle_params().accumulated_correction; - auto const r12 = - box_geo.get_mi_vector(p1.pos(), p2.pos()) - - ((p1.mass() + p2.mass()) / p2.mass()) * accum; + auto const r12 = box_geo.get_mi_vector(p1.pos(), p2.pos()) - + ((p1.mass() + p2.mass()) / p2.mass()) * accum; auto const F_c = 2.0 * accum * p1.mass() / sq_dt; - auto const stress = - Utils::flatten(Utils::tensor_product(F_c, r12)); + auto const stress = Utils::flatten(Utils::tensor_product(F_c, r12)); auto dest = obs_pressure.bonded_contribution(bond_id); for (std::size_t k = 0; k < 9u; ++k) dest[k] += stress[k]; diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py index 729506e01e..bbb99281fd 100644 --- a/testsuite/python/integrator_exceptions.py +++ b/testsuite/python/integrator_exceptions.py @@ -290,7 +290,8 @@ def test_rigid_bond_incompatible_integrators(self): 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' + 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) diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py index ee35a2b52c..3f54b442c7 100644 --- a/testsuite/python/rigid_bond_virial.py +++ b/testsuite/python/rigid_bond_virial.py @@ -22,6 +22,7 @@ import espressomd.interactions import numpy as np + @utx.skipIfMissingFeatures("BOND_CONSTRAINT") class RigidBondVirialTest(ut.TestCase): @@ -58,9 +59,9 @@ def _virial_in_rotation(self, set_integrator): v_p = np.trace(self.system.analysis.pressure_tensor()['bonded']) / 3. v_theory = -2.0 / (3. * V) self.assertAlmostEqual(v_p, v_theory, delta=0.01 * abs(v_theory)) - v_xx = self.system.analysis.pressure_tensor()['bonded'][0,0] - v_yy = self.system.analysis.pressure_tensor()['bonded'][1,1] - v_zz = self.system.analysis.pressure_tensor()['bonded'][2,2] + v_xx = self.system.analysis.pressure_tensor()['bonded'][0, 0] + v_yy = self.system.analysis.pressure_tensor()['bonded'][1, 1] + v_zz = self.system.analysis.pressure_tensor()['bonded'][2, 2] self.assertAlmostEqual(v_xx, -2.0 / V, delta=0.01 * abs(2.0 / V)) self.assertAlmostEqual(v_yy, 0.0, delta=1e-8) self.assertAlmostEqual(v_zz, 0.0, delta=1e-8) @@ -77,7 +78,7 @@ def _virial_unequal_masses(self, set_integrator): m1 = 1.0 m2 = 2.0 d = 1.0 - v = 1.0 # velocity + v = 1.0 # velocity V = self.system.volume() set_integrator() @@ -92,7 +93,8 @@ def _virial_unequal_masses(self, set_integrator): v_p = np.trace(pt) / 3. self.assertAlmostEqual(v_p, v_theory, delta=0.01 * abs(v_theory)) # Bond is along x: all virial goes into xx, none into yy or zz - self.assertAlmostEqual(pt[0, 0], 3. * v_theory, delta=0.01 * abs(3. * v_theory)) + self.assertAlmostEqual(pt[0, 0], 3. * v_theory, + delta=0.01 * abs(3. * v_theory)) self.assertAlmostEqual(pt[1, 1], 0., delta=1e-8) self.assertAlmostEqual(pt[2, 2], 0., delta=1e-8) @@ -102,7 +104,8 @@ def test_virial_unequal_masses_vv(self): def test_virial_unequal_masses_se(self): """SE: constraint virial with m1!=m2 matches centripetal theory.""" - self._virial_unequal_masses(self.system.integrator.set_symplectic_euler) + self._virial_unequal_masses( + self.system.integrator.set_symplectic_euler) # Langevin consistency: mean of (P_bond + P_kin) = kT/V, # std matches analytic fluctuation formula. @@ -115,7 +118,8 @@ def _virial_consistency(self, set_integrator, noise_prefactor): self.system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42) set_integrator() - std_theory = ((6 * kT**2 + noise_prefactor * gamma * mass * kT / dt) / (9 * V**2))**0.5 + std_theory = ((6 * kT**2 + noise_prefactor * gamma * + mass * kT / dt) / (9 * V**2))**0.5 self._make_dimer(v=0.0) @@ -125,14 +129,17 @@ def _virial_consistency(self, set_integrator, noise_prefactor): virial = [] for _ in range(n_loop): self.system.integrator.run(n_steps) - v_p = np.trace(self.system.analysis.pressure_tensor()['bonded']) / 3. - v_k = np.trace(self.system.analysis.pressure_tensor()['kinetic']) / 3. + v_p = np.trace( + self.system.analysis.pressure_tensor()['bonded']) / 3. + v_k = np.trace(self.system.analysis.pressure_tensor()[ + 'kinetic']) / 3. virial.append(v_p + v_k) rigid_p = np.mean(virial) rigid_std = np.std(virial) - self.assertAlmostEqual(rigid_p, 1. / V, delta=2.*std_theory/n_loop**0.5) - self.assertAlmostEqual(rigid_std, std_theory, delta=0.02*std_theory) + self.assertAlmostEqual( + rigid_p, 1. / V, delta=2. * std_theory / n_loop**0.5) + self.assertAlmostEqual(rigid_std, std_theory, delta=0.02 * std_theory) def test_virial_consistency_vv(self): """VV+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" @@ -140,9 +147,9 @@ def test_virial_consistency_vv(self): def test_virial_consistency_se(self): """SE+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" - self._virial_consistency(self.system.integrator.set_symplectic_euler, 4) + self._virial_consistency( + self.system.integrator.set_symplectic_euler, 4) if __name__ == "__main__": ut.main() - From b26d5f5e840ca24857dc689e365af210d92a2509 Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Tue, 30 Jun 2026 16:47:29 +0200 Subject: [PATCH 3/6] Reshape test script --- testsuite/python/CMakeLists.txt | 3 +- testsuite/python/rigid_bond_virial.py | 46 +--------- testsuite/python/rigid_bond_virial_stats.py | 93 +++++++++++++++++++++ 3 files changed, 96 insertions(+), 46 deletions(-) create mode 100644 testsuite/python/rigid_bond_virial_stats.py diff --git a/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 016c1d7a6a..a53809e090 100644 --- a/testsuite/python/CMakeLists.txt +++ b/testsuite/python/CMakeLists.txt @@ -264,7 +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 LABELS long) +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) diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py index 3f54b442c7..09bc15bcd1 100644 --- a/testsuite/python/rigid_bond_virial.py +++ b/testsuite/python/rigid_bond_virial.py @@ -77,7 +77,7 @@ def test_virial_in_rotation_se(self): def _virial_unequal_masses(self, set_integrator): m1 = 1.0 m2 = 2.0 - d = 1.0 + d = 1.0 # distance between particle v = 1.0 # velocity V = self.system.volume() @@ -107,49 +107,5 @@ def test_virial_unequal_masses_se(self): self._virial_unequal_masses( self.system.integrator.set_symplectic_euler) - # Langevin consistency: mean of (P_bond + P_kin) = kT/V, - # std matches analytic fluctuation formula. - def _virial_consistency(self, set_integrator, noise_prefactor): - kT = 1.0 - gamma = 1.0 - mass = 1. - V = self.system.volume() - dt = self.system.time_step - self.system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42) - set_integrator() - - std_theory = ((6 * kT**2 + noise_prefactor * gamma * - mass * kT / dt) / (9 * V**2))**0.5 - - self._make_dimer(v=0.0) - - self.system.integrator.run(1000) # equilibrate - n_loop = 2000 - n_steps = 100 - virial = [] - for _ in range(n_loop): - self.system.integrator.run(n_steps) - v_p = np.trace( - self.system.analysis.pressure_tensor()['bonded']) / 3. - v_k = np.trace(self.system.analysis.pressure_tensor()[ - 'kinetic']) / 3. - virial.append(v_p + v_k) - - rigid_p = np.mean(virial) - rigid_std = np.std(virial) - self.assertAlmostEqual( - rigid_p, 1. / V, delta=2. * std_theory / n_loop**0.5) - self.assertAlmostEqual(rigid_std, std_theory, delta=0.02 * std_theory) - - def test_virial_consistency_vv(self): - """VV+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" - self._virial_consistency(self.system.integrator.set_vv, 1) - - def test_virial_consistency_se(self): - """SE+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" - self._virial_consistency( - self.system.integrator.set_symplectic_euler, 4) - - if __name__ == "__main__": ut.main() diff --git a/testsuite/python/rigid_bond_virial_stats.py b/testsuite/python/rigid_bond_virial_stats.py new file mode 100644 index 0000000000..506cdf7422 --- /dev/null +++ b/testsuite/python/rigid_bond_virial_stats.py @@ -0,0 +1,93 @@ +# +# Copyright (C) 2013-2026 The ESPResSo project +# +# This file is part of ESPResSo. +# +# ESPResSo is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# ESPResSo is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +import unittest as ut +import unittest_decorators as utx +import espressomd +import espressomd.interactions +import numpy as np + + +@utx.skipIfMissingFeatures("BOND_CONSTRAINT") +class RigidBondVirialStatsTest(ut.TestCase): + + system = espressomd.System(box_l=[10., 10., 10.]) + system.cell_system.skin = 0.4 + + def setUp(self): + self.system.time_step = 1e-2 + + def tearDown(self): + self.system.non_bonded_inter.reset() + self.system.part.clear() + self.system.thermostat.turn_off() + + def _make_dimer(self): + bond = espressomd.interactions.RigidBond(r=1.0, ptol=1e-4, vtol=1e-4) + self.system.bonded_inter.add(bond) + p1 = self.system.part.add(pos=[4.5, 5.0, 5.0], v=[0., 0., 0.], mass=1.0) + p2 = self.system.part.add(pos=[5.5, 5.0, 5.0], v=[0., 0., 0.], mass=1.0) + p2.add_bond((bond, p1)) + return p1, p2 + + # Langevin consistency: mean of (P_bond + P_kin) = kT/V, + # std matches analytic fluctuation formula. + def _virial_consistency(self, set_integrator, noise_prefactor): + kT = 1.0 + gamma = 1.0 + mass = 1. + V = self.system.volume() + dt = self.system.time_step + self.system.thermostat.set_langevin(kT=kT, gamma=gamma, seed=42) + set_integrator() + + std_theory = ((6 * kT**2 + noise_prefactor * gamma * + mass * kT / dt) / (9 * V**2))**0.5 + + self._make_dimer() + + self.system.integrator.run(1000) # equilibrate + n_loop = 2000 + n_steps = 100 + virial = [] + for _ in range(n_loop): + self.system.integrator.run(n_steps) + v_p = np.trace( + self.system.analysis.pressure_tensor()['bonded']) / 3. + v_k = np.trace(self.system.analysis.pressure_tensor()[ + 'kinetic']) / 3. + virial.append(v_p + v_k) + + rigid_p = np.mean(virial) + rigid_std = np.std(virial) + self.assertAlmostEqual( + rigid_p, 1. / V, delta=2. * std_theory / n_loop**0.5) + self.assertAlmostEqual(rigid_std, std_theory, delta=0.02 * std_theory) + + def test_virial_consistency_vv(self): + """VV+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" + self._virial_consistency(self.system.integrator.set_vv, 1) + + def test_virial_consistency_se(self): + """SE+Langevin: rigid bond virial satisfies equipartition and fluctuation formula.""" + self._virial_consistency( + self.system.integrator.set_symplectic_euler, 4) + + +if __name__ == "__main__": + ut.main() From 252c23ce6d2992ce16a2eb978e7e66e2de2c3663 Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Tue, 30 Jun 2026 16:48:45 +0200 Subject: [PATCH 4/6] Fixed style --- testsuite/python/rigid_bond_virial.py | 1 + testsuite/python/rigid_bond_virial_stats.py | 6 ++++-- 2 files changed, 5 insertions(+), 2 deletions(-) diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py index 09bc15bcd1..b5a19b1e67 100644 --- a/testsuite/python/rigid_bond_virial.py +++ b/testsuite/python/rigid_bond_virial.py @@ -107,5 +107,6 @@ def test_virial_unequal_masses_se(self): self._virial_unequal_masses( self.system.integrator.set_symplectic_euler) + if __name__ == "__main__": ut.main() diff --git a/testsuite/python/rigid_bond_virial_stats.py b/testsuite/python/rigid_bond_virial_stats.py index 506cdf7422..c1e9f0aeff 100644 --- a/testsuite/python/rigid_bond_virial_stats.py +++ b/testsuite/python/rigid_bond_virial_stats.py @@ -40,8 +40,10 @@ def tearDown(self): def _make_dimer(self): bond = espressomd.interactions.RigidBond(r=1.0, ptol=1e-4, vtol=1e-4) self.system.bonded_inter.add(bond) - p1 = self.system.part.add(pos=[4.5, 5.0, 5.0], v=[0., 0., 0.], mass=1.0) - p2 = self.system.part.add(pos=[5.5, 5.0, 5.0], v=[0., 0., 0.], mass=1.0) + p1 = self.system.part.add(pos=[4.5, 5.0, 5.0], v=[ + 0., 0., 0.], mass=1.0) + p2 = self.system.part.add(pos=[5.5, 5.0, 5.0], v=[ + 0., 0., 0.], mass=1.0) p2.add_bond((bond, p1)) return p1, p2 From 59b0880809453172413aebd1462a98a7582e88fe Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Tue, 7 Jul 2026 17:51:21 +0200 Subject: [PATCH 5/6] Fixed bug --- src/core/Particle.hpp | 12 +--- .../bonded_interaction_data.hpp | 14 ++++ src/core/pressure.cpp | 37 ++++------ src/core/rattle.cpp | 57 ++++++++++----- src/core/rattle.hpp | 5 +- src/core/unit_tests/Particle_test.cpp | 55 -------------- testsuite/python/rigid_bond_virial.py | 71 ++++++++++++++++++- 7 files changed, 145 insertions(+), 106 deletions(-) diff --git a/src/core/Particle.hpp b/src/core/Particle.hpp index 4db73fb80b..7c2ac3df91 100644 --- a/src/core/Particle.hpp +++ b/src/core/Particle.hpp @@ -415,26 +415,18 @@ struct ParticleLocal { struct ParticleRattle { /** position/velocity correction */ Utils::Vector3d correction = {0., 0., 0.}; - /** total over all RATTLE iterations */ - Utils::Vector3d accumulated_correction = {0., 0., 0.}; friend ParticleRattle operator+(ParticleRattle const &lhs, ParticleRattle const &rhs) { - return {lhs.correction + rhs.correction, - lhs.accumulated_correction + rhs.accumulated_correction}; + return {lhs.correction + rhs.correction}; } - // Intentionally excludes accumulated_correction: ghost reduction must not - // overwrite the local particle's accumulated value with ghost data. ParticleRattle &operator+=(ParticleRattle const &rhs) { - correction += rhs.correction; - return *this; + return *this = *this + rhs; } template void serialize(Archive &ar, long int /* version */) { ar & correction; - // accumulated_correction is reset at the start of every SHAKE call; - // serializing it would waste checkpoint space with a transient value. } }; #endif diff --git a/src/core/bonded_interactions/bonded_interaction_data.hpp b/src/core/bonded_interactions/bonded_interaction_data.hpp index 5824eda3ad..d31dcceeb5 100644 --- a/src/core/bonded_interactions/bonded_interaction_data.hpp +++ b/src/core/bonded_interactions/bonded_interaction_data.hpp @@ -48,6 +48,8 @@ #include "Particle.hpp" #include "system/Leaf.hpp" +#include + #include #include #include @@ -55,6 +57,7 @@ #include #include #include +#include /** Interaction type for unused bonded interaction slots */ struct NoneBond { @@ -144,6 +147,17 @@ class BondedInteractionsMap : public System::Leaf { 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 rigid_bond_virial; #endif std::optional find_bond_id(mapped_type const &target_bond) const { for (auto const &[bond_id, bond] : m_params) { diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index c07feff2bf..c71f8961af 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -40,6 +40,7 @@ #include "virtual_sites/relative.hpp" #include +#include #include #include @@ -153,29 +154,19 @@ std::shared_ptr System::calculate_pressure() { #ifdef ESPRESSO_BOND_CONSTRAINT if (propagation->is_inertial() and bonded_ias->get_n_rigid_bonds() >= 1) { - auto const dt = get_time_step(); - auto const sq_dt = dt * dt; - cell_structure->bond_loop( - [&obs_pressure, &bonded_ias = *bonded_ias, &box_geo = *box_geo, - sq_dt](Particle &p1, int bond_id, std::span partners) { - if (std::holds_alternative(*bonded_ias.at(bond_id))) { - auto const &p2 = *partners[0]; - // Recover pre-SHAKE bond vector from local accumulated_correction - // only: accum2 = -(m1/m2)*accum1, so - // r12_pre = r12_post - accum1 + accum2 - // = r12_post - (m1+m2)/m2 * accum1 - // This avoids reading accumulated_correction from ghost p2. - auto const &accum = p1.rattle_params().accumulated_correction; - auto const r12 = box_geo.get_mi_vector(p1.pos(), p2.pos()) - - ((p1.mass() + p2.mass()) / p2.mass()) * accum; - auto const F_c = 2.0 * accum * p1.mass() / sq_dt; - auto const stress = Utils::flatten(Utils::tensor_product(F_c, r12)); - auto dest = obs_pressure.bonded_contribution(bond_id); - for (std::size_t k = 0; k < 9u; ++k) - dest[k] += stress[k]; - } - return false; - }); + // 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(bond_id)); + for (std::size_t k = 0; k < 9u; ++k) + dest[k] += stress[k]; + } } #endif diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index 0171d29d6a..8636be19e0 100644 --- a/src/core/rattle.cpp +++ b/src/core/rattle.cpp @@ -32,13 +32,19 @@ #include "communication.hpp" #include "errorhandling.hpp" +#include +#include +#include + #include #include #include +#include #include #include #include +#include /** 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 &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(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(&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,23 +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 accumulator for this timestep - boost::for_each(particles, [](Particle &p) { - p.rattle_params().accumulated_correction.fill(0.); - }); + // Reset the per-bond-type constraint virial for this timestep + bonded_ias.rigid_bond_virial.assign( + static_cast(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()); @@ -181,11 +207,6 @@ void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo, cs.ghosts_reduce_rattle_correction(); - // Accumulate AFTER ghost reduction so cross-rank bonds are included - boost::for_each(particles, [](Particle &p) { - p.rattle_params().accumulated_correction += p.rattle_params().correction; - }); - apply_positional_correction(particles); cs.ghosts_update(Cells::DATA_PART_POSITION | Cells::DATA_PART_MOMENTUM); } @@ -247,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()); diff --git a/src/core/rattle.hpp b/src/core/rattle.hpp index e3bbf2cdb3..71a6d5060c 100644 --- a/src/core/rattle.hpp +++ b/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. diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp index a7706e6a41..e735774a02 100644 --- a/src/core/unit_tests/Particle_test.cpp +++ b/src/core/unit_tests/Particle_test.cpp @@ -201,8 +201,6 @@ void check_particle_rattle(ParticleRattle const &out, ParticleRattle const &ref) { BOOST_TEST(out.correction == ref.correction, boost::test_tools::per_element()); - BOOST_TEST(out.accumulated_correction == ref.accumulated_correction, - boost::test_tools::per_element()); } BOOST_AUTO_TEST_CASE(rattle_serialization) { @@ -234,27 +232,6 @@ BOOST_AUTO_TEST_CASE(rattle_serialization) { } } -BOOST_AUTO_TEST_CASE(rattle_checkpoint_serialization) { - // accumulated_correction is transient; it must NOT survive a checkpoint - // round-trip so that restored simulations start each SHAKE call clean. - auto pr = ParticleRattle{{1., 2., 3.}, {4., 5., 6.}}; - - std::stringstream stream; - { - boost::archive::text_oarchive oa(stream); - oa << pr; - } - - boost::archive::text_iarchive ia(stream); - ParticleRattle out; - ia >> out; - - BOOST_TEST(out.correction == (Utils::Vector3d{1., 2., 3.}), - boost::test_tools::per_element()); - BOOST_TEST(out.accumulated_correction == Utils::Vector3d{}, - boost::test_tools::per_element()); -} - BOOST_AUTO_TEST_CASE(rattle_constructors) { auto pr = ParticleRattle{{1, 2, 3}}; @@ -272,38 +249,6 @@ BOOST_AUTO_TEST_CASE(rattle_constructors) { } } -BOOST_AUTO_TEST_CASE(rattle_accumulated_correction) { - // default-initialized accumulated_correction must be zero - ParticleRattle pr{}; - auto const zero = Utils::Vector3d{}; - BOOST_TEST(pr.accumulated_correction == zero, - boost::test_tools::per_element()); - - // operator+ must add both fields independently - auto const corr_a = Utils::Vector3d{1., 0., 0.}; - auto const accum_a = Utils::Vector3d{0., 1., 0.}; - auto const corr_b = Utils::Vector3d{0., 0., 1.}; - auto const accum_b = Utils::Vector3d{1., 0., 0.}; - ParticleRattle a{corr_a, accum_a}; - ParticleRattle b{corr_b, accum_b}; - - auto const expected_corr = corr_a + corr_b; // {1,0,1} - auto const expected_accum = accum_a + accum_b; // {1,1,0} - ParticleRattle expected_sum{expected_corr, expected_accum}; - check_particle_rattle(a + b, expected_sum); - - // operator+= must add only correction, NOT accumulated_correction - // (ghost reduction must not bleed stale accumulated values into locals) - auto const dst_accum = Utils::Vector3d{5., 5., 5.}; - ParticleRattle dst{Utils::Vector3d{2., 0., 0.}, dst_accum}; - dst += b; - auto const expected_dst_corr = Utils::Vector3d{2., 0., 1.}; - BOOST_TEST(dst.correction == expected_dst_corr, - boost::test_tools::per_element()); - // accumulated_correction must be unchanged after += - BOOST_TEST(dst.accumulated_correction == dst_accum, - boost::test_tools::per_element()); -} #endif // ESPRESSO_BOND_CONSTRAINT #ifdef ESPRESSO_THERMAL_STONER_WOHLFARTH diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py index b5a19b1e67..6608b9ab24 100644 --- a/testsuite/python/rigid_bond_virial.py +++ b/testsuite/python/rigid_bond_virial.py @@ -50,7 +50,7 @@ def _make_dimer(self, m1=1.0, m2=1.0, v=1.0): return p1, p2 # Rotation test: after one step, constraint virial = F_centripetal·d - # F_centripetal = m v²/r = 1·1/0.5 = 2, W = -virial/(3V) = -2/(3V) + # F_centripetal = m v²/r = 1·1/0.5 = 2, W = -virial/(3V) = -2/(3V) def _virial_in_rotation(self, set_integrator): V = self.system.volume() set_integrator() @@ -107,6 +107,75 @@ def test_virial_unequal_masses_se(self): self._virial_unequal_masses( self.system.integrator.set_symplectic_euler) + def _make_chain(self, masses, bond_lengths, omega=1.0): + """ + Build a rigid, collinear 3-particle chain 0-1-2 (two RigidBonds, + particle 1 shared by both) set up in rigid rotation about its + center of mass, so the net constraint force on each particle is + exactly the centripetal force F_i = -m_i * omega^2 * x_i. + """ + l1, l2 = bond_lengths + x1 = 0.0 + x0 = x1 - l1 + x2 = x1 + l2 + com = (masses[0] * x0 + masses[1] * x1 + masses[2] * x2) / sum(masses) + x = [x0 - com, x1 - com, x2 - com] + + parts = [] + for m, xi in zip(masses, x): + v = (0., omega * xi, 0.) + p = self.system.part.add(pos=[5. + xi, 5.0, 5.0], v=v, mass=m) + parts.append(p) + + bond01 = espressomd.interactions.RigidBond(r=l1, ptol=1e-9, vtol=1e-9) + bond12 = espressomd.interactions.RigidBond(r=l2, ptol=1e-9, vtol=1e-9) + self.system.bonded_inter.add(bond01) + self.system.bonded_inter.add(bond12) + parts[1].add_bond((bond01, parts[0])) + parts[2].add_bond((bond12, parts[1])) + return parts, x + + def _virial_chain(self, set_integrator, masses, bond_lengths, omega=1.0): + V = self.system.volume() + set_integrator() + _, x = self._make_chain(masses, bond_lengths, omega=omega) + self.system.integrator.run(1) + + # Centripetal theory, generalized from the dimer case: each particle's + # net constraint force is F_i = -m_i*omega^2*x_i, so the constraint + # virial is W_xx = sum_i x_i * F_i (bond is along x: all virial is + # in xx, none in yy or zz). + w_xx = sum(-m * omega**2 * xi**2 for m, xi in zip(masses, x)) + v_theory = w_xx / (3. * V) + + pt = self.system.analysis.pressure_tensor()['bonded'] + v_p = np.trace(pt) / 3. + self.assertAlmostEqual(v_p, v_theory, delta=0.01 * abs(v_theory)) + self.assertAlmostEqual(pt[0, 0], w_xx / V, + delta=0.01 * abs(w_xx / V)) + self.assertAlmostEqual(pt[1, 1], 0., delta=1e-8) + self.assertAlmostEqual(pt[2, 2], 0., delta=1e-8) + + def test_virial_chain_symmetric_vv(self): + """VV: constraint virial of a rotating rigid 3-particle chain (equal masses/bond lengths, shared middle bond) matches centripetal theory.""" + self._virial_chain(self.system.integrator.set_vv, + masses=[1.0, 1.0, 1.0], bond_lengths=(1.0, 1.0)) + + def test_virial_chain_symmetric_se(self): + """SE: constraint virial of a rotating rigid 3-particle chain (equal masses/bond lengths, shared middle bond) matches centripetal theory.""" + self._virial_chain(self.system.integrator.set_symplectic_euler, + masses=[1.0, 1.0, 1.0], bond_lengths=(1.0, 1.0)) + + def test_virial_chain_unequal_masses_vv(self): + """VV: constraint virial of a rigid 3-particle chain with unequal masses and bond lengths matches centripetal theory.""" + self._virial_chain(self.system.integrator.set_vv, + masses=[2.0, 1.0, 3.0], bond_lengths=(1.0, 1.5)) + + def test_virial_chain_unequal_masses_se(self): + """SE: constraint virial of a rigid 3-particle chain with unequal masses and bond lengths matches centripetal theory.""" + self._virial_chain(self.system.integrator.set_symplectic_euler, + masses=[2.0, 1.0, 3.0], bond_lengths=(1.0, 1.5)) + if __name__ == "__main__": ut.main() From d9bbd89bd1cc86ee4023f4809a7afff59fef7274 Mon Sep 17 00:00:00 2001 From: Hideki Kobayashi Date: Tue, 7 Jul 2026 19:51:53 +0200 Subject: [PATCH 6/6] Restored unnecessary modification --- src/core/ghosts.cpp | 6 +++--- src/core/unit_tests/Particle_test.cpp | 1 - 2 files changed, 3 insertions(+), 4 deletions(-) diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 14d2415fdd..ded254f9e5 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -354,9 +354,9 @@ add_rattle_correction_from_recv_buffer(CommBuf &recv_buffer, auto archiver = Utils::MemcpyIArchive{recv_buffer.make_span()}; for (auto &part_list : ghost_comm.part_lists) { for (Particle &part : *part_list) { - Utils::Vector3d correction; - archiver >> correction; - part.rattle_correction() += correction; + ParticleRattle pr; + archiver >> pr; + part.rattle_params() += pr; } } } diff --git a/src/core/unit_tests/Particle_test.cpp b/src/core/unit_tests/Particle_test.cpp index e735774a02..8f5948af09 100644 --- a/src/core/unit_tests/Particle_test.cpp +++ b/src/core/unit_tests/Particle_test.cpp @@ -248,7 +248,6 @@ BOOST_AUTO_TEST_CASE(rattle_constructors) { check_particle_rattle(out, pr); } } - #endif // ESPRESSO_BOND_CONSTRAINT #ifdef ESPRESSO_THERMAL_STONER_WOHLFARTH