Remove Observable_stat external linkage (#3770)

Description of changes:
- re-use code in the cython interface to extract data from `Observable_stat` objects
- remove `obs_scalar_pressure` global variable
   - write-only variable in the core -> can be calculated from `obs_pressure_tensor` in cython
   - simplify pressure kernels
   - less MPI communication
- rename `obs_pressure_tensor` to `obs_pressure` for consistency with `obs_energy`
- remove external linkage for  `obs_pressure` and `obs_energy`
   - these globals are no longer visible outside of their translation unit
   - use a const reference getter for the script interface
   - partial fix for #2628
- narrow includes of `Observable_stat.hpp`
- add ELC/DLC energy corrections to the P3M/DP3M/DDS solver energy slots
   - the `analyze.energy()`, `analyze.pressure()`, `analyze.pressure_tensor()` functions now only reports 2 slots for Coulomb and dipolar contributions: particles contribution and electrostatics/magnetostatics solver contribution

src/core/Observable_stat.cpp

@@ -27,8 +27,8 @@
 
 Observable_stat::Observable_stat(size_t chunk_size) : m_chunk_size(chunk_size) {
   // number of chunks for different interaction types
-  auto constexpr n_coulomb = 3;
-  auto constexpr n_dipolar = 3;
+  auto constexpr n_coulomb = 2;
+  auto constexpr n_dipolar = 2;
 #ifdef VIRTUAL_SITES
   auto constexpr n_vs = 1;
 #else

src/core/Observable_stat.hpp

@@ -29,7 +29,7 @@
 #include <utility>
 #include <vector>
 
-/** Observable for the scalar pressure, pressure tensor and energy. */
+/** Observable for the pressure and energy. */
 class Observable_stat {
   /** Array for observables on each node. */
   std::vector<double> m_data;
@@ -68,7 +68,7 @@ class Observable_stat {
    *  @param column Which column to sum up (only relevant for multi-dimensional
    *                observables).
    */
-  double accumulate(double acc = 0.0, size_t column = 0) {
+  double accumulate(double acc = 0.0, size_t column = 0) const {
     assert(column < m_chunk_size);
     if (m_chunk_size == 1)
       return boost::accumulate(m_data, acc);

src/core/bonded_interactions/bonded_interactions.dox

@@ -126,25 +126,20 @@
  *        result = fene_pair_force(iaparams, dx);
  *        break;
  *      @endcode
- *    - in function @ref add_bonded_force(): add the new interaction to the
- *      switch statement if it is not a pairwise bond. For the harmonic angle,
- *      the code looks like this:
+ *    - in functions called by @ref add_bonded_force(): add the new interaction
+ *      to the switch statements if it is not a pairwise bond. For the harmonic
+ *      angle, the code looks like this:
  *      @code{.cpp}
  *      case BONDED_IA_ANGLE_HARMONIC:
- *        std::tie(force1, force2, force3) =
- *          angle_harmonic_force(p1->r.p, p2->r.p, p3->r.p, iaparams);
- *        break;
+ *        return angle_harmonic_force(p1.r.p, p2.r.p, p3.r.p, iaparams);
  *      @endcode
  *  * In energy_inline.hpp:
  *    - include the header file of the new interaction
- *    - in function @ref add_bonded_energy(): add the new interaction to the
+ *    - in function @ref calc_bonded_energy(): add the new interaction to the
  *      switch statement. For the FENE bond, e.g., the code looks like this:
  *      @code{.cpp}
- *      boost::optional<double> retval;
- *      // ...
  *      case BONDED_IA_FENE:
- *        retval = fene_pair_energy(iaparams, dx);
- *        break;
+ *        return fene_pair_energy(iaparams, dx);
  *      @endcode
  *  * Pressure tensor and virial calculation: if your bonded
  *    interaction is not a pair bond or modifies the particles involved,

src/core/constraints/ShapeBasedConstraint.cpp

@@ -113,7 +113,8 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p,
 
 void ShapeBasedConstraint::add_energy(const Particle &p,
                                       const Utils::Vector3d &folded_pos,
-                                      double t, Observable_stat &energy) const {
+                                      double t,
+                                      Observable_stat &obs_energy) const {
   double dist;
   double nonbonded_en = 0.0;
 

src/core/cuda_common_cuda.cu

@@ -22,7 +22,6 @@
 
 #include "debug.hpp"
 
-#include "Observable_stat.hpp"
 #include "ParticleRange.hpp"
 #include "cuda_init.hpp"
 #include "cuda_interface.hpp"
@@ -221,16 +220,12 @@ void clear_energy_on_GPU() {
   cuda_safe_mem(cudaMemset(energy_device, 0, sizeof(CUDA_energy)));
 }
 
-void copy_energy_from_GPU() {
-  extern Observable_stat obs_energy;
+CUDA_energy copy_energy_from_GPU() {
   if (!global_part_vars_host.communication_enabled)
-    return;
+    return {};
   cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy),
                            cudaMemcpyDeviceToHost));
-  if (!obs_energy.coulomb.empty())
-    obs_energy.coulomb[0] += energy_host.coulomb;
-  if (obs_energy.dipolar.size() >= 2)
-    obs_energy.dipolar[1] += energy_host.dipolar;
+  return energy_host;
 }
 
 void _cuda_safe_mem(cudaError_t CU_err, const char *file, unsigned int line) {

src/core/cuda_interface.hpp

@@ -107,7 +107,7 @@ typedef struct {
 } CUDA_global_part_vars;
 
 void copy_forces_from_GPU(ParticleRange &particles);
-void copy_energy_from_GPU();
+CUDA_energy copy_energy_from_GPU();
 void clear_energy_on_GPU();
 
 CUDA_global_part_vars *gpu_get_global_particle_vars_pointer_host();

src/core/electrostatics_magnetostatics/coulomb.cpp

@@ -42,9 +42,7 @@
 Coulomb_parameters coulomb;
 
 namespace Coulomb {
-void calc_pressure_long_range(Observable_stat &virials,
-                              Observable_stat &p_tensor,
-                              const ParticleRange &particles) {
+Utils::Vector9d calc_pressure_long_range(const ParticleRange &particles) {
   switch (coulomb.method) {
 #ifdef P3M
   case COULOMB_ELC_P3M:
@@ -58,11 +56,7 @@ void calc_pressure_long_range(Observable_stat &virials,
     break;
   case COULOMB_P3M: {
     p3m_charge_assign(particles);
-    auto const p3m_p_tensor = p3m_calc_kspace_pressure_tensor();
-    std::copy_n(p3m_p_tensor.data(), 9, p_tensor.coulomb.begin() + 9);
-    virials.coulomb[1] = p3m_p_tensor[0] + p3m_p_tensor[4] + p3m_p_tensor[8];
-
-    break;
+    return p3m_calc_kspace_pressure_tensor();
   }
 #endif
   case COULOMB_MMM1D:
@@ -74,6 +68,7 @@ void calc_pressure_long_range(Observable_stat &virials,
   default:
     break;
   }
+  return {};
 }
 
 void sanity_checks(int &state) {
@@ -299,8 +294,8 @@ void calc_long_range_force(const ParticleRange &particles) {
 #endif
 }
 
-void calc_energy_long_range(Observable_stat &energy,
-                            const ParticleRange &particles) {
+double calc_energy_long_range(const ParticleRange &particles) {
+  double energy = 0.0;
   switch (coulomb.method) {
 #ifdef P3M
   case COULOMB_P3M_GPU:
@@ -309,46 +304,47 @@ void calc_energy_long_range(Observable_stat &energy,
     break;
   case COULOMB_P3M:
     p3m_charge_assign(particles);
-    energy.coulomb[1] = p3m_calc_kspace_forces(false, true, particles);
+    energy = p3m_calc_kspace_forces(false, true, particles);
     break;
   case COULOMB_ELC_P3M:
     // assign the original charges first
     // they may not have been assigned yet
     p3m_charge_assign(particles);
     if (!elc_params.dielectric_contrast_on)
-      energy.coulomb[1] = p3m_calc_kspace_forces(false, true, particles);
+      energy = p3m_calc_kspace_forces(false, true, particles);
     else {
-      energy.coulomb[1] = 0.5 * p3m_calc_kspace_forces(false, true, particles);
-      energy.coulomb[1] += 0.5 * coulomb.prefactor *
-                           ELC_P3M_dielectric_layers_energy_self(particles);
+      energy = 0.5 * p3m_calc_kspace_forces(false, true, particles);
+      energy += 0.5 * coulomb.prefactor *
+                ELC_P3M_dielectric_layers_energy_self(particles);
 
       // assign both original and image charges now
       ELC_p3m_charge_assign_both(particles);
       ELC_P3M_modify_p3m_sums_both(particles);
 
-      energy.coulomb[1] += 0.5 * p3m_calc_kspace_forces(false, true, particles);
+      energy += 0.5 * p3m_calc_kspace_forces(false, true, particles);
 
       // assign only the image charges now
       ELC_p3m_charge_assign_image(particles);
       ELC_P3M_modify_p3m_sums_image(particles);
 
-      energy.coulomb[1] -= 0.5 * p3m_calc_kspace_forces(false, true, particles);
+      energy -= 0.5 * p3m_calc_kspace_forces(false, true, particles);
 
       // restore modified sums
       ELC_P3M_restore_p3m_sums(particles);
     }
-    energy.coulomb[2] = ELC_energy(particles);
+    energy += ELC_energy(particles);
     break;
 #endif
 #ifdef SCAFACOS
   case COULOMB_SCAFACOS:
     assert(!Scafacos::dipolar());
-    energy.coulomb[1] += Scafacos::long_range_energy();
+    energy += Scafacos::long_range_energy();
     break;
 #endif
   default:
     break;
   }
+  return energy;
 }
 
 int iccp3m_sanity_check() {

src/core/electrostatics_magnetostatics/coulomb.hpp

@@ -24,7 +24,6 @@
 #include <cstddef>
 
 #ifdef ELECTROSTATICS
-#include "Observable_stat.hpp"
 
 #include <ParticleRange.hpp>
 #include <utils/Vector.hpp>
@@ -60,9 +59,7 @@ struct Coulomb_parameters {
 extern Coulomb_parameters coulomb;
 
 namespace Coulomb {
-void calc_pressure_long_range(Observable_stat &virials,
-                              Observable_stat &p_tensor,
-                              const ParticleRange &particles);
+Utils::Vector9d calc_pressure_long_range(const ParticleRange &particles);
 
 void sanity_checks(int &state);
 double cutoff(const Utils::Vector3d &box_l);
@@ -75,8 +72,7 @@ void init();
 
 void calc_long_range_force(const ParticleRange &particles);
 
-void calc_energy_long_range(Observable_stat &energy,
-                            const ParticleRange &particles);
+double calc_energy_long_range(const ParticleRange &particles);
 
 int iccp3m_sanity_check();
 

src/core/electrostatics_magnetostatics/dipole.cpp

@@ -205,33 +205,31 @@ void calc_long_range_force(const ParticleRange &particles) {
   }
 }
 
-void calc_energy_long_range(Observable_stat &energy,
-                            const ParticleRange &particles) {
+double calc_energy_long_range(const ParticleRange &particles) {
+  double energy = 0.;
   switch (dipole.method) {
 #ifdef DP3M
   case DIPOLAR_P3M:
     dp3m_dipole_assign(particles);
-    energy.dipolar[1] = dp3m_calc_kspace_forces(false, true, particles);
+    energy = dp3m_calc_kspace_forces(false, true, particles);
     break;
   case DIPOLAR_MDLC_P3M:
     dp3m_dipole_assign(particles);
-    energy.dipolar[1] = dp3m_calc_kspace_forces(false, true, particles);
-    energy.dipolar[2] = add_mdlc_energy_corrections(particles);
+    energy = dp3m_calc_kspace_forces(false, true, particles);
+    energy += add_mdlc_energy_corrections(particles);
     break;
 #endif
   case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA:
-    energy.dipolar[1] = dawaanr_calculations(false, true, particles);
+    energy = dawaanr_calculations(false, true, particles);
     break;
 #ifdef DP3M
   case DIPOLAR_MDLC_DS:
-    energy.dipolar[1] =
-        magnetic_dipolar_direct_sum_calculations(false, true, particles);
-    energy.dipolar[2] = add_mdlc_energy_corrections(particles);
+    energy = magnetic_dipolar_direct_sum_calculations(false, true, particles);
+    energy += add_mdlc_energy_corrections(particles);
     break;
 #endif
   case DIPOLAR_DS:
-    energy.dipolar[1] =
-        magnetic_dipolar_direct_sum_calculations(false, true, particles);
+    energy = magnetic_dipolar_direct_sum_calculations(false, true, particles);
     break;
   case DIPOLAR_DS_GPU: // NOLINT(bugprone-branch-clone)
     // do nothing: it's an actor
@@ -244,7 +242,7 @@ void calc_energy_long_range(Observable_stat &energy,
 #ifdef SCAFACOS_DIPOLES
   case DIPOLAR_SCAFACOS:
     assert(Scafacos::dipolar());
-    energy.dipolar[1] = Scafacos::long_range_energy();
+    energy = Scafacos::long_range_energy();
 #endif
   case DIPOLAR_NONE:
     break;
@@ -253,6 +251,7 @@ void calc_energy_long_range(Observable_stat &energy,
         << "energy calculation not implemented for dipolar method.";
     break;
   }
+  return energy;
 }
 
 int set_mesh() {

src/core/electrostatics_magnetostatics/dipole.hpp

@@ -24,8 +24,6 @@
 #include <cstddef>
 
 #ifdef DIPOLES
-#include "Observable_stat.hpp"
-
 #include <utils/Vector.hpp>
 
 #include <ParticleRange.hpp>
@@ -80,8 +78,7 @@ void init();
 
 void calc_long_range_force(const ParticleRange &particles);
 
-void calc_energy_long_range(Observable_stat &energy,
-                            const ParticleRange &particles);
+double calc_energy_long_range(const ParticleRange &particles);
 
 int set_mesh();
 void bcast_params(const boost::mpi::communicator &comm);

src/core/energy.cpp

@@ -43,6 +43,8 @@ ActorList energyActors;
 /** Energy of the system */
 Observable_stat obs_energy{1};
 
+Observable_stat const &get_obs_energy() { return obs_energy; }
+
 /** Reduce the system energy from all MPI ranks. */
 void master_energy_calc() { mpi_gather_stats(GatherStats::energy); }
 
@@ -65,13 +67,14 @@ void energy_calc(const double time) {
   on_observable_calc();
 
   for (auto const &p : cell_structure.local_particles()) {
-    add_kinetic_energy(p);
+    obs_energy.kinetic[0] += calc_kinetic_energy(p);
   }
 
   short_range_loop(
-      [](Particle &p) { add_single_particle_energy(p); },
+      [](Particle &p) { add_bonded_energy(p, obs_energy); },
       [](Particle const &p1, Particle const &p2, Distance const &d) {
-        add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2);
+        add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2,
+                                   obs_energy);
       });
 
   calc_long_range_energies(cell_structure.local_particles());
@@ -80,7 +83,11 @@ void energy_calc(const double time) {
   Constraints::constraints.add_energy(local_parts, time, obs_energy);
 
 #ifdef CUDA
-  copy_energy_from_GPU();
+  auto const energy_host = copy_energy_from_GPU();
+  if (!obs_energy.coulomb.empty())
+    obs_energy.coulomb[1] += energy_host.coulomb;
+  if (!obs_energy.dipolar.empty())
+    obs_energy.dipolar[1] += energy_host.dipolar;
 #endif
 
   /* gather data */
@@ -95,12 +102,13 @@ void update_energy() { master_energy_calc(); }
 void calc_long_range_energies(const ParticleRange &particles) {
 #ifdef ELECTROSTATICS
   /* calculate k-space part of electrostatic interaction. */
-  Coulomb::calc_energy_long_range(obs_energy, particles);
-#endif /* ifdef ELECTROSTATICS */
+  obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(particles);
+#endif
 
 #ifdef DIPOLES
-  Dipole::calc_energy_long_range(obs_energy, particles);
-#endif /* ifdef DIPOLES */
+  /* calculate k-space part of magnetostatic interaction. */
+  obs_energy.dipolar[1] = Dipole::calc_energy_long_range(particles);
+#endif
 }
 
 double calculate_current_potential_energy_of_system() {

src/core/energy.hpp

@@ -27,17 +27,21 @@
 #ifndef _ENERGY_H
 #define _ENERGY_H
 
+#include "Observable_stat.hpp"
 #include "ParticleRange.hpp"
 #include "actor/ActorList.hpp"
 
 extern ActorList energyActors;
 
-/** Calculate energies. */
-void update_energy();
-
 /** Parallel energy calculation. */
 void energy_calc(double time);
 
+/** Run @ref energy_calc in parallel. */
+void update_energy();
+
+/** Return the energy observable. */
+Observable_stat const &get_obs_energy();
+
 /** Calculate long-range energies (P3M, ...). */
 void calc_long_range_energies(const ParticleRange &particles);