diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 252296760c..f7fd6f5abb 100644 --- a/src/core/Observable_stat.cpp +++ b/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(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 -Observable_stat::non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const { +Observable_stat::get_non_bonded_contribution(Utils::Span base_pointer, + int type1, int type2) const { auto const offset = static_cast(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}; diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index b34157004b..b84bcfdda0 100644 --- a/src/core/Observable_stat.hpp +++ b/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 non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const; + Utils::Span + get_non_bonded_contribution(Utils::Span 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 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 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 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. */ diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 8a3a1b8554..9b4e2fe210 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -164,7 +164,9 @@ void ShapeBasedConstraint::add_energy(const Particle &p, } } // NOLINTNEXTLINE(clang-analyzer-cplusplus.NewDeleteLeaks) - if (part_rep.type() >= 0) - obs_energy.add_non_bonded_contribution(p.type(), part_rep.type(), energy); + 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 diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index eaa24054fd..f4c09b93a7 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -183,7 +183,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)); diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index a1f6354eb8..7a6ef6589a 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -70,10 +70,8 @@ inline void add_non_bonded_pair_virials( .f + calc_non_central_force(p1, p2, ia_params, d, dist).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 diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index abd19fb80f..d5d23a3857 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -114,6 +114,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { BOOST_REQUIRE_GE(get_particle_node_parallel(pid2), rank ? -1 : 1); BOOST_REQUIRE_GE(get_particle_node_parallel(pid3), rank ? -1 : 1); } + set_particle_property(pid1, &Particle::mol_id, type_a); + set_particle_property(pid2, &Particle::mol_id, type_b); + set_particle_property(pid3, &Particle::mol_id, type_b); auto const reset_particle_positions = [&start_positions]() { for (auto const &kv : start_positions) { @@ -225,10 +228,10 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { on_non_bonded_ia_change(); // 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); @@ -236,6 +239,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory) { auto const obs_energy = calculate_energy(); if (rank == 0) { 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); diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py index cbe2e829e1..9b1eeff806 100644 --- a/testsuite/python/analyze_energy.py +++ b/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() @@ -86,12 +86,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) @@ -101,6 +103,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)