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/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..0a39c4b721 100644 --- a/src/core/integrators/Propagation.hpp +++ b/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; + } }; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 2130adf849..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 @@ -151,6 +152,24 @@ std::shared_ptr 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(bond_id)); + for (std::size_t k = 0; k < 9u; ++k) + dest[k] += stress[k]; + } + } +#endif + obs_pressure.rescale(volume); obs_pressure.mpi_reduce(); diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index 4a738807e9..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,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(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()); @@ -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()); 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/testsuite/python/CMakeLists.txt b/testsuite/python/CMakeLists.txt index 951e2c19ec..a53809e090 100644 --- a/testsuite/python/CMakeLists.txt +++ b/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) diff --git a/testsuite/python/integrator_exceptions.py b/testsuite/python/integrator_exceptions.py index a7a9974552..bbb99281fd 100644 --- a/testsuite/python/integrator_exceptions.py +++ b/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"), diff --git a/testsuite/python/rigid_bond_virial.py b/testsuite/python/rigid_bond_virial.py new file mode 100644 index 0000000000..6608b9ab24 --- /dev/null +++ b/testsuite/python/rigid_bond_virial.py @@ -0,0 +1,181 @@ +# +# 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 # distance between particle + 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) + + 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() diff --git a/testsuite/python/rigid_bond_virial_stats.py b/testsuite/python/rigid_bond_virial_stats.py new file mode 100644 index 0000000000..c1e9f0aeff --- /dev/null +++ b/testsuite/python/rigid_bond_virial_stats.py @@ -0,0 +1,95 @@ +# +# 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()