Remove Observable_stat external linkage

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

@@ -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/electrostatics_magnetostatics/coulomb.cpp

@@ -42,8 +42,7 @@
 Coulomb_parameters coulomb;
 
 namespace Coulomb {
-void calc_pressure_long_range(Observable_stat &virials,
-                              Observable_stat &p_tensor,
+void calc_pressure_long_range(Observable_stat &p_tensor,
                               const ParticleRange &particles) {
   switch (coulomb.method) {
 #ifdef P3M
@@ -60,7 +59,6 @@ void calc_pressure_long_range(Observable_stat &virials,
     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;
   }

src/core/electrostatics_magnetostatics/coulomb.hpp

@@ -60,8 +60,7 @@ struct Coulomb_parameters {
 extern Coulomb_parameters coulomb;
 
 namespace Coulomb {
-void calc_pressure_long_range(Observable_stat &virials,
-                              Observable_stat &p_tensor,
+void calc_pressure_long_range(Observable_stat &p_tensor,
                               const ParticleRange &particles);
 
 void sanity_checks(int &state);

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,11 +37,6 @@
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 
-/** Scalar pressure of the system */
-Observable_stat obs_scalar_pressure{1};
-/** Pressure tensor of the system */
-Observable_stat obs_pressure_tensor{9};
-
 nptiso_struct nptiso = {0.0,
                         0.0,
                         0.0,
@@ -56,70 +51,58 @@ 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
   /* calculate k-space part of electrostatic interaction. */
-  Coulomb::calc_pressure_long_range(obs_scalar_pressure, obs_pressure_tensor,
-                                    particles);
+  Coulomb::calc_pressure_long_range(obs_pressure, particles);
 #endif
 #ifdef DIPOLES
   /* calculate k-space part of magnetostatic interaction. */
   Dipole::calc_pressure_long_range();
 #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();
 
   if (!interactions_sanity_checks())
     return;
 
-  obs_scalar_pressure = Observable_stat{1};
-  obs_pressure_tensor = Observable_stat{9};
+  obs_pressure = Observable_stat{9};
 
   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());
 
 #ifdef VIRTUAL_SITES
-  if (!obs_scalar_pressure.virtual_sites.empty()) {
-    auto const vs_pressure_tensor = virtual_sites()->pressure_tensor();
-
-    obs_scalar_pressure.virtual_sites[0] += trace(vs_pressure_tensor);
-    boost::copy(flatten(vs_pressure_tensor),
-                obs_pressure_tensor.virtual_sites.begin());
+  if (!obs_pressure.virtual_sites.empty()) {
+    auto const vs_pressure = virtual_sites()->pressure_tensor();
+    boost::copy(flatten(vs_pressure), obs_pressure.virtual_sites.begin());
   }
 #endif
 
-  /* rescale kinetic energy (=ideal contribution) */
-  obs_scalar_pressure.rescale(3.0 * volume);
-
-  obs_pressure_tensor.rescale(volume);
+  obs_pressure.rescale(volume);
 
   /* gather data */
-  auto obs_scalar_pressure_res = reduce(comm_cart, obs_scalar_pressure);
-  if (obs_scalar_pressure_res) {
-    std::swap(obs_scalar_pressure, *obs_scalar_pressure_res);
-  }
-
-  auto obs_pressure_tensor_res = reduce(comm_cart, obs_pressure_tensor);
-  if (obs_pressure_tensor_res) {
-    std::swap(obs_pressure_tensor, *obs_pressure_tensor_res);
+  auto obs_pressure_res = reduce(comm_cart, obs_pressure);
+  if (obs_pressure_res) {
+    std::swap(obs_pressure, *obs_pressure_res);
   }
 }
 
@@ -129,7 +112,7 @@ Utils::Vector9d observable_compute_pressure_tensor() {
   update_pressure();
   Utils::Vector9d pressure_tensor{};
   for (size_t j = 0; j < 9; j++) {
-    pressure_tensor[j] = obs_pressure_tensor.accumulate(0, j);
+    pressure_tensor[j] = obs_pressure.accumulate(0, j);
   }
   return pressure_tensor;
 }

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