Fix energy/pressure observables (espressomd#4788)

Fixes espressomd#4787

Description of changes:
- the energy and pressure observables now correctly tally non-bonded contributions

src/core/Observable_stat.cpp

@@ -58,8 +58,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions
 
-  auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
-                          n_dipolar + n_vs + n_ext_fields;
+  auto const n_elements = n_kinetic + n_bonded + 2ul * n_non_bonded +
+                          n_coulomb + n_dipolar + n_vs + n_ext_fields;
   m_data = std::vector<double>(m_chunk_size * n_elements);
 
   // spans for the different contributions
@@ -78,8 +78,8 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
 }
 
 Utils::Span<double>
-Observable_stat::non_bonded_contribution(Utils::Span<double> base_pointer,
-                                         int type1, int type2) const {
+Observable_stat::get_non_bonded_contribution(Utils::Span<double> base_pointer,
+                                             int type1, int type2) const {
   auto const offset = static_cast<std::size_t>(Utils::upper_triangular(
       std::min(type1, type2), std::max(type1, type2), max_seen_particle_type));
   return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size};

src/core/Observable_stat.hpp

@@ -38,8 +38,9 @@ class Observable_stat {
   std::size_t m_chunk_size;
 
   /** Get contribution from a non-bonded interaction */
-  Utils::Span<double> non_bonded_contribution(Utils::Span<double> base_pointer,
-                                              int type1, int type2) const;
+  Utils::Span<double>
+  get_non_bonded_contribution(Utils::Span<double> base_pointer, int type1,
+                              int type2) const;
 
 public:
   explicit Observable_stat(std::size_t chunk_size);
@@ -91,29 +92,30 @@ class Observable_stat {
     return {bonded.data() + offset, m_chunk_size};
   }
 
-  void add_non_bonded_contribution(int type1, int type2,
+  void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2,
                                    Utils::Span<const double> data) {
     assert(data.size() == m_chunk_size);
-    auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter;
-    auto const dest = non_bonded_contribution(source, type1, type2);
+    auto const span = (molid1 == molid2) ? non_bonded_intra : non_bonded_inter;
+    auto const dest = get_non_bonded_contribution(span, type1, type2);
 
     boost::transform(dest, data, dest.begin(), std::plus<>{});
   }
 
-  void add_non_bonded_contribution(int type1, int type2, double data) {
-    add_non_bonded_contribution(type1, type2, {&data, 1});
+  void add_non_bonded_contribution(int type1, int type2, int molid1, int molid2,
+                                   double data) {
+    add_non_bonded_contribution(type1, type2, molid1, molid2, {&data, 1});
   }
 
   /** Get contribution from a non-bonded intramolecular interaction */
   Utils::Span<double> non_bonded_intra_contribution(int type1,
                                                     int type2) const {
-    return non_bonded_contribution(non_bonded_intra, type1, type2);
+    return get_non_bonded_contribution(non_bonded_intra, type1, type2);
   }
 
   /** Get contribution from a non-bonded intermolecular interaction */
   Utils::Span<double> non_bonded_inter_contribution(int type1,
                                                     int type2) const {
-    return non_bonded_contribution(non_bonded_inter, type1, type2);
+    return get_non_bonded_contribution(non_bonded_inter, type1, type2);
   }
 
   /** MPI reduction. */

src/core/constraints/ShapeBasedConstraint.cpp

@@ -154,7 +154,10 @@ void ShapeBasedConstraint::add_energy(const Particle &p,
       runtimeErrorMsg() << "Constraint violated by particle " << p.id();
     }
   }
-  if (part_rep.type() >= 0)
-    obs_energy.add_non_bonded_contribution(p.type(), part_rep.type(), energy);
+  // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks)
+  if (part_rep.type() >= 0) {
+    obs_energy.add_non_bonded_contribution(
+        p.type(), part_rep.type(), p.mol_id(), part_rep.mol_id(), energy);
+  }
 }
 } // namespace Constraints

src/core/energy_inline.hpp

@@ -184,7 +184,7 @@ inline void add_non_bonded_pair_energy(
   if (do_nonbonded(p1, p2))
 #endif
     obs_energy.add_non_bonded_contribution(
-        p1.type(), p2.type(),
+        p1.type(), p2.type(), p1.mol_id(), p2.mol_id(),
         calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist,
                                     coulomb_kernel));
 

src/core/pressure_inline.hpp

@@ -66,10 +66,8 @@ inline void add_non_bonded_pair_virials(
     auto const force =
         calc_non_bonded_pair_force(p1, p2, ia_params, d, dist, kernel_forces).f;
     auto const stress = Utils::tensor_product(d, force);
-
-    auto const type1 = p1.mol_id();
-    auto const type2 = p2.mol_id();
-    obs_pressure.add_non_bonded_contribution(type1, type2, flatten(stress));
+    obs_pressure.add_non_bonded_contribution(p1.type(), p2.type(), p1.mol_id(),
+                                             p2.mol_id(), flatten(stress));
   }
 
 #ifdef ELECTROSTATICS

src/core/unit_tests/EspressoSystemStandAlone_test.cpp

@@ -150,6 +150,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
       BOOST_REQUIRE_GE(get_particle_node(pid3), 1);
     }
   }
+  set_particle_mol_id(pid1, type_a);
+  set_particle_mol_id(pid2, type_b);
+  set_particle_mol_id(pid3, type_b);
 
   auto const reset_particle_positions = [&start_positions]() {
     for (auto const &kv : start_positions) {
@@ -214,16 +217,17 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
     lennard_jones_set_params(type_b, type_b, eps, sig, cut, shift, offset, min);
 
     // matrix indices and reference energy value
-    auto const max_type = type_b + 1;
-    auto const n_pairs = Utils::upper_triangular(type_b, type_b, max_type) + 1;
-    auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, max_type);
-    auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, max_type);
+    auto const size = std::max(type_a, type_b) + 1;
+    auto const n_pairs = Utils::upper_triangular(type_b, type_b, size) + 1;
+    auto const lj_pair_ab = Utils::upper_triangular(type_a, type_b, size);
+    auto const lj_pair_bb = Utils::upper_triangular(type_b, type_b, size);
     auto const frac6 = Utils::int_pow<6>(sig / r_off);
     auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift);
 
     // measure energies
     auto const obs_energy = calculate_energy();
     for (int i = 0; i < n_pairs; ++i) {
+      // particles were set up with type == mol_id
       auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.;
       auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.;
       BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 1e-10);

testsuite/python/analyze_energy.py

@@ -49,8 +49,8 @@ def setUpClass(cls):
         cls.system.bonded_inter[5] = cls.harmonic
 
     def setUp(self):
-        self.system.part.add(pos=[1, 2, 2], type=0)
-        self.system.part.add(pos=[5, 2, 2], type=0)
+        self.system.part.add(pos=[1, 2, 2], type=0, mol_id=6)
+        self.system.part.add(pos=[5, 2, 2], type=0, mol_id=6)
 
     def tearDown(self):
         self.system.part.clear()
@@ -85,12 +85,14 @@ def test_non_bonded(self):
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 1., delta=1e-7)
+        self.assertAlmostEqual(energy["non_bonded_intra"], 1., delta=1e-7)
+        self.assertAlmostEqual(energy["non_bonded_inter"], 0., delta=1e-7)
         # Test the single particle energy function
         self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
             [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)
         # add another pair of particles
-        self.system.part.add(pos=[3, 2, 2], type=1)
-        self.system.part.add(pos=[4, 2, 2], type=1)
+        self.system.part.add(pos=[3, 2, 2], type=1, mol_id=7)
+        self.system.part.add(pos=[4, 2, 2], type=1, mol_id=7)
         energy = self.system.analysis.energy()
         self.assertAlmostEqual(energy["total"], 3., delta=1e-7)
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
@@ -100,6 +102,18 @@ def test_non_bonded(self):
         self.assertAlmostEqual(energy["non_bonded", 0, 0]
                                + energy["non_bonded", 0, 1]
                                + energy["non_bonded", 1, 1], energy["total"], delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_intra", 0, 0], 1., delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_intra", 1, 1], 1., delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_intra", 0, 1], 0., delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_inter", 0, 0], 0., delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_inter", 1, 1], 0., delta=1e-7)
+        self.assertAlmostEqual(
+            energy["non_bonded_inter", 0, 1], 1., delta=1e-7)
         # Test the single particle energy function
         self.assertAlmostEqual(energy["non_bonded"], 0.5 * sum(
             [self.system.analysis.particle_energy(p) for p in self.system.part.all()]), delta=1e-7)