core: Make Observable_stat globals local

Pass local variable obs_energy and obs_pressure to energy/pressure
kernels as mutable references to remove external linkage.

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

@@ -221,8 +221,7 @@ 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;
+void copy_energy_from_GPU(Observable_stat &obs_energy) {
   if (!global_part_vars_host.communication_enabled)
     return;
   cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy),

src/core/cuda_interface.hpp

@@ -24,6 +24,7 @@
 #ifdef CUDA
 
 #include "CudaHostAllocator.hpp"
+#include "Observable_stat.hpp"
 #include "ParticleRange.hpp"
 
 #include <utils/Span.hpp>
@@ -107,7 +108,7 @@ typedef struct {
 } CUDA_global_part_vars;
 
 void copy_forces_from_GPU(ParticleRange &particles);
-void copy_energy_from_GPU();
+void copy_energy_from_GPU(Observable_stat &obs_energy);
 void clear_energy_on_GPU();
 
 CUDA_global_part_vars *gpu_get_global_particle_vars_pointer_host();

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);
+    add_kinetic_energy(p, obs_energy);
   }
 
   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,7 @@ void energy_calc(const double time) {
   Constraints::constraints.add_energy(local_parts, time, obs_energy);
 
 #ifdef CUDA
-  copy_energy_from_GPU();
+  copy_energy_from_GPU(obs_energy);
 #endif
 
   /* gather data */

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

src/core/energy_inline.hpp

@@ -71,8 +71,6 @@
 #include "electrostatics_magnetostatics/dipole_inline.hpp"
 #endif
 
-extern Observable_stat obs_energy;
-
 /** Calculate non-bonded energies between a pair of particles.
  *  @param p1         particle 1.
  *  @param p2         particle 2.
@@ -179,10 +177,12 @@ inline double calc_non_bonded_pair_energy(Particle const &p1,
  *  @param d         vector between p1 and p2.
  *  @param dist      distance between p1 and p2.
  *  @param dist2     distance squared between p1 and p2.
+ *  @param[in,out] obs_energy   energy observable.
  */
 inline void add_non_bonded_pair_energy(Particle const &p1, Particle const &p2,
                                        Utils::Vector3d const &d,
-                                       double const dist, double const dist2) {
+                                       double const dist, double const dist2,
+                                       Observable_stat &obs_energy) {
   IA_parameters const &ia_params = *get_ia_param(p1.p.type, p2.p.type);
 
 #ifdef EXCLUSIONS
@@ -279,32 +279,12 @@ calc_bonded_energy(Bonded_ia_parameters const &iaparams, Particle const &p1,
   throw BondInvalidSizeError(n_partners);
 }
 
-/** Add bonded energies for one particle to the energy observable.
- *  @param[in] p1   particle for which to calculate energies
- *  @param[in] bond_id Numeric id of the bond
- *  @param[in] partners Bond partners of particle.
- *
- *  @return True if bond was broken, false otherwise.
- */
-inline bool add_bonded_energy(Particle &p1, int bond_id,
-                              Utils::Span<Particle *> partners) {
-  auto const &iaparams = bonded_ia_params[bond_id];
-
-  auto const result = calc_bonded_energy(iaparams, p1, partners);
-
-  if (result) {
-    obs_energy.bonded_contribution(bond_id)[0] += result.get();
-
-    return false;
-  }
-
-  return true;
-}
-
 /** Add kinetic energies for one particle to the energy observable.
  *  @param[in] p1   particle for which to calculate energies
+ *  @param[out]    obs_energy   energy observable
  */
-inline void add_kinetic_energy(Particle const &p1) {
+inline void add_kinetic_energy(Particle const &p1,
+                               Observable_stat &obs_energy) {
   if (p1.p.is_virtual)
     return;
 
@@ -323,11 +303,22 @@ inline void add_kinetic_energy(Particle const &p1) {
 #endif
 }
 
-/** Add kinetic and bonded energies for one particle to the energy observable.
+/** Add bonded energies for one particle to the energy observable.
  *  @param[in] p   particle for which to calculate energies
+ *  @param[out]    obs_energy   energy observable
  */
-inline void add_single_particle_energy(Particle &p) {
-  cell_structure.execute_bond_handler(p, add_bonded_energy);
+inline void add_bonded_energy(Particle &p, Observable_stat &obs_energy) {
+  cell_structure.execute_bond_handler(
+      p, [&obs_energy](Particle &p1, int bond_id,
+                       Utils::Span<Particle *> partners) {
+        auto const &iaparams = bonded_ia_params[bond_id];
+        auto const result = calc_bonded_energy(iaparams, p1, partners);
+        if (result) {
+          obs_energy.bonded_contribution(bond_id)[0] += result.get();
+          return false;
+        }
+        return true;
+      });
 }
 
 #endif // ENERGY_INLINE_HPP

src/core/pressure.cpp

@@ -37,9 +37,6 @@
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 
-/** Pressure tensor of the system */
-Observable_stat obs_pressure{9};
-
 nptiso_struct nptiso = {0.0,
                         0.0,
                         0.0,
@@ -54,6 +51,11 @@ nptiso_struct nptiso = {0.0,
                         false,
                         0};
 
+/** Pressure tensor of the system */
+Observable_stat obs_pressure{9};
+
+Observable_stat const &get_obs_pressure() { return obs_pressure; }
+
 /** Calculate long-range virials (P3M, ...). */
 void calc_long_range_virials(const ParticleRange &particles) {
 #ifdef ELECTROSTATICS
@@ -66,10 +68,6 @@ void calc_long_range_virials(const ParticleRange &particles) {
 #endif
 }
 
-static void add_single_particle_virials(Particle &p) {
-  cell_structure.execute_bond_handler(p, add_bonded_pressure_tensor);
-}
-
 void pressure_calc() {
   auto const volume = box_geo.volume();
 
@@ -81,13 +79,13 @@ void pressure_calc() {
   on_observable_calc();
 
   for (auto const &p : cell_structure.local_particles()) {
-    add_kinetic_virials(p);
+    add_kinetic_virials(p, obs_pressure);
   }
 
-  short_range_loop([](Particle &p) { add_single_particle_virials(p); },
+  short_range_loop([](Particle &p) { add_bonded_virials(p, obs_pressure); },
                    [](Particle &p1, Particle &p2, Distance const &d) {
-                     add_non_bonded_pair_virials(p1, p2, d.vec21,
-                                                 sqrt(d.dist2));
+                     add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2),
+                                                 obs_pressure);
                    });
 
   calc_long_range_virials(cell_structure.local_particles());

src/core/pressure.hpp

@@ -25,16 +25,20 @@
 #ifndef CORE_PRESSURE_HPP
 #define CORE_PRESSURE_HPP
 
+#include "Observable_stat.hpp"
+
 #include <utils/Vector.hpp>
 
-/** Parallel pressure calculation from a virial expansion.
- */
+/** Parallel pressure calculation from a virial expansion. */
 void pressure_calc();
 
+/** Run @ref pressure_calc in parallel. */
+void update_pressure();
+
+/** Return the pressure observable. */
+Observable_stat const &get_obs_pressure();
+
 /** Helper function for @ref Observables::PressureTensor. */
 Utils::Vector9d observable_compute_pressure_tensor();
 
-/** Calculate the scalar pressure and pressure tensor. */
-void update_pressure();
-
 #endif

src/core/pressure_inline.hpp

@@ -32,16 +32,16 @@
 
 #include <utils/math/tensor_product.hpp>
 
-extern Observable_stat obs_pressure;
-
 /** Calculate non bonded energies between a pair of particles.
  *  @param p1        pointer to particle 1.
  *  @param p2        pointer to particle 2.
  *  @param d         vector between p1 and p2.
  *  @param dist      distance between p1 and p2.
+ *  @param[in,out] obs_pressure   pressure observable.
  */
 inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2,
-                                        Utils::Vector3d const &d, double dist) {
+                                        Utils::Vector3d const &d, double dist,
+                                        Observable_stat &obs_pressure) {
 #ifdef EXCLUSIONS
   if (do_nonbonded(p1, p2))
 #endif
@@ -136,29 +136,12 @@ calc_bonded_pressure_tensor(Bonded_ia_parameters const &iaparams, Particle &p1,
   }
 }
 
-inline bool add_bonded_pressure_tensor(Particle &p1, int bond_id,
-                                       Utils::Span<Particle *> partners) {
-  auto const &iaparams = bonded_ia_params[bond_id];
-
-  auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners);
-  if (result) {
-    auto const &tensor = result.get();
-
-    /* pressure tensor part */
-    for (int k = 0; k < 3; k++)
-      for (int l = 0; l < 3; l++)
-        obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l];
-
-    return false;
-  }
-
-  return true;
-}
-
 /** Calculate kinetic pressure (aka energy) for one particle.
- *  @param p1 particle for which to calculate pressure
+ *  @param[in] p1   particle for which to calculate pressure
+ *  @param[out]    obs_energy  pressure observable
  */
-inline void add_kinetic_virials(Particle const &p1) {
+inline void add_kinetic_virials(Particle const &p1,
+                                Observable_stat &obs_pressure) {
   if (p1.p.is_virtual)
     return;
 
@@ -168,4 +151,27 @@ inline void add_kinetic_virials(Particle const &p1) {
       obs_pressure.kinetic[k * 3 + l] += p1.m.v[k] * p1.m.v[l] * p1.p.mass;
 }
 
+/** Add bonded energies for one particle to the energy observable.
+ *  @param[in] p   particle for which to calculate pressure
+ *  @param[out]    obs_pressure   pressure observable
+ */
+inline void add_bonded_virials(Particle &p, Observable_stat &obs_pressure) {
+  cell_structure.execute_bond_handler(p, [&obs_pressure](
+                                             Particle &p1, int bond_id,
+                                             Utils::Span<Particle *> partners) {
+    auto const &iaparams = bonded_ia_params[bond_id];
+    auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners);
+    if (result) {
+      auto const &tensor = result.get();
+      /* pressure tensor part */
+      for (int k = 0; k < 3; k++)
+        for (int l = 0; l < 3; l++)
+          obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l];
+
+      return false;
+    }
+    return true;
+  });
+}
+
 #endif

src/python/espressomd/analyze.pxd

@@ -81,17 +81,13 @@ cdef extern from "statistics_chain.hpp":
     array4 calc_rg(int, int, int) except +
     array2 calc_rh(int, int, int)
 
-cdef extern from "pressure_inline.hpp":
-    cdef Observable_stat obs_pressure
-
 cdef extern from "pressure.hpp":
     cdef void update_pressure()
-
-cdef extern from "energy_inline.hpp":
-    cdef Observable_stat obs_energy
+    cdef const Observable_stat & get_obs_pressure()
 
 cdef extern from "energy.hpp":
     cdef void update_energy()
+    cdef const Observable_stat & get_obs_energy()
     double calculate_current_potential_energy_of_system()
 
 cdef extern from "dpd.hpp":