diff --git a/src/python/espressomd/analyze.py b/src/python/espressomd/analyze.py index 693267c643..ca2d4c106c 100644 --- a/src/python/espressomd/analyze.py +++ b/src/python/espressomd/analyze.py @@ -111,7 +111,8 @@ class Analysis(ScriptInterfaceHelper): ---------- p_type : :obj:`int` Particle :attr:`~espressomd.particle_data.ParticleHandle.type` - for which to calculate the center of mass. + for which to calculate the center of mass, or ``-1`` for all + particles. Returns ------- diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index d6aa626084..5398690585 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -167,13 +167,23 @@ Variant Analysis::do_call_method(std::string const &name, } if (name == "center_of_mass") { auto const p_type = get_value(parameters, "p_type"); - context()->parallel_try_catch([&]() { check_particle_type(p_type); }); + // p_type == -1 is the documented sentinel for "all (non-virtual) + // particles", which the core function handles directly; only validate + // explicit (non-sentinel) particle types. + if (p_type != -1) { + context()->parallel_try_catch([&]() { check_particle_type(p_type); }); + } auto const local = center_of_mass(get_system(), p_type); return mpi_reduce_sum(context()->get_comm(), local).as_vector(); } if (name == "angular_momentum") { auto const p_type = get_value(parameters, "p_type"); - context()->parallel_try_catch([&]() { check_particle_type(p_type); }); + // p_type == -1 is the documented sentinel for "all (non-virtual) + // particles", which the core function handles directly; only validate + // explicit (non-sentinel) particle types. + if (p_type != -1) { + context()->parallel_try_catch([&]() { check_particle_type(p_type); }); + } auto const local = angular_momentum(get_system(), p_type); return mpi_reduce_sum(context()->get_comm(), local).as_vector(); } diff --git a/testsuite/python/analyze_mass_related.py b/testsuite/python/analyze_mass_related.py index 66c8714eb9..b66c6b0575 100644 --- a/testsuite/python/analyze_mass_related.py +++ b/testsuite/python/analyze_mass_related.py @@ -97,6 +97,18 @@ def test_center_of_mass(self): np.testing.assert_allclose(com, com_ref) + def test_center_of_mass_all_particles(self): + # p_type=-1 is the documented sentinel for "all (non-virtual) + # particles", and must be distinct from the type-0-only result. + no_virtual = self.system.part.select(lambda p: not p.is_virtual()) + com_ref = np.zeros(3) + for p in no_virtual: + com_ref += p.pos * p.mass + com_ref /= np.sum(no_virtual.mass) + com = self.system.analysis.center_of_mass(p_type=-1) + + np.testing.assert_allclose(com, com_ref) + def test_galilei_transform(self): no_virtual = self.system.part.select(lambda p: not p.is_virtual()) @@ -120,6 +132,18 @@ def test_angularmomentum(self): am, self.system.analysis.angular_momentum(p_type=0)) + def test_angularmomentum_all_particles(self): + # p_type=-1 is the documented sentinel for "all (non-virtual) + # particles", and must be distinct from the type-0-only result. + no_virtual = self.system.part.select(lambda p: not p.is_virtual()) + am_ref = np.zeros(3) + for p in no_virtual: + am_ref += p.mass * np.cross(p.pos, p.v) + + np.testing.assert_allclose( + am_ref, + self.system.analysis.angular_momentum(p_type=-1)) + def test_kinetic_energy(self): no_virtual = self.system.part.select(lambda p: not p.is_virtual()) E_kin = 0.5 * \