analysis: accept p_type=-1 for center_of_mass and angular_momentum (bug-sweep #50)

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
         -------

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<int>(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<int>(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();
   }

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 * \