diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index aaedbdc1eb0..8a24cfe27f0 100644
--- a/src/core/CMakeLists.txt
+++ b/src/core/CMakeLists.txt
@@ -41,7 +41,6 @@ set(EspressoCore_SRC
     CellStructure.cpp
     PartCfg.cpp
     AtomDecomposition.cpp
-    reduce_observable_stat.cpp
     EspressoSystemStandAlone.cpp
     DomainDecomposition.cpp)
 
diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp
index 6822174b292..70a2c33f2b1 100644
--- a/src/core/Observable_stat.cpp
+++ b/src/core/Observable_stat.cpp
@@ -24,32 +24,42 @@
 #include "config.hpp"
 
 #include "bonded_interactions/bonded_interaction_data.hpp"
+#include "communication.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
 #include <utils/Span.hpp>
+#include <utils/index.hpp>
+
+#include <boost/mpi/collectives/reduce.hpp>
 
 #include <cassert>
 #include <cstddef>
+#include <functional>
 #include <vector>
 
+inline std::size_t max_non_bonded_pairs() {
+  auto const n = static_cast<std::size_t>(max_seen_particle_type);
+  return (n * (n + 1)) / 2;
+}
+
 Observable_stat::Observable_stat(std::size_t chunk_size)
     : m_chunk_size(chunk_size) {
   // number of chunks for different interaction types
-  auto constexpr n_coulomb = 2;
-  auto constexpr n_dipolar = 2;
+  constexpr std::size_t n_coulomb = 2;
+  constexpr std::size_t n_dipolar = 2;
 #ifdef VIRTUAL_SITES
-  auto constexpr n_vs = 1;
+  constexpr std::size_t n_vs = 1;
 #else
-  auto constexpr n_vs = 0;
+  constexpr std::size_t n_vs = 0;
 #endif
   auto const n_bonded = bonded_ia_params.size();
   auto const n_non_bonded = max_non_bonded_pairs();
   constexpr std::size_t n_ext_fields = 1; // reduction over all fields
   constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions
 
-  // resize vector
-  auto const total = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb +
-                     n_dipolar + n_vs + n_ext_fields;
-  m_data = std::vector<double>(m_chunk_size * total);
+  auto const n_elements = n_kinetic + n_bonded + 2 * 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
   kinetic = Utils::Span<double>(m_data.data(), m_chunk_size);
@@ -65,3 +75,20 @@ Observable_stat::Observable_stat(std::size_t chunk_size)
       Utils::Span<double>(non_bonded_intra.end(), n_non_bonded * m_chunk_size);
   assert(non_bonded_inter.end() == (m_data.data() + m_data.size()));
 }
+
+Utils::Span<double>
+Observable_stat::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};
+}
+
+void Observable_stat::mpi_reduce() {
+  if (comm_cart.rank() == 0) {
+    std::vector<double> temp(m_data);
+    boost::mpi::reduce(comm_cart, temp, m_data, std::plus<>{}, 0);
+  } else {
+    boost::mpi::reduce(comm_cart, m_data, std::plus<>{}, 0);
+  }
+}
diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp
index a0d44e5cc15..72109bcf892 100644
--- a/src/core/Observable_stat.hpp
+++ b/src/core/Observable_stat.hpp
@@ -23,7 +23,6 @@
 #include <boost/range/numeric.hpp>
 
 #include <utils/Span.hpp>
-#include <utils/index.hpp>
 
 #include <algorithm>
 #include <cassert>
@@ -38,32 +37,14 @@ class Observable_stat {
   /** Number of doubles per data item */
   std::size_t m_chunk_size;
 
-  /** Calculate the maximal number of non-bonded interaction pairs in the
-   *  system.
-   */
-  static std::size_t max_non_bonded_pairs() {
-    extern int max_seen_particle_type;
-    return static_cast<std::size_t>(
-        (max_seen_particle_type * (max_seen_particle_type + 1)) / 2);
-  }
-
   /** Get contribution from a non-bonded interaction */
   Utils::Span<double> non_bonded_contribution(Utils::Span<double> base_pointer,
-                                              int type1, int type2) const {
-    extern int max_seen_particle_type;
-    int offset = 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};
-  }
+                                              int type1, int type2) const;
 
 public:
   explicit Observable_stat(std::size_t chunk_size);
 
-  auto chunk_size() const { return m_chunk_size; }
-  Utils::Span<double> data_() { return {m_data.data(), m_data.size()}; }
-  Utils::Span<const double> data_() const {
-    return {m_data.data(), m_data.size()};
-  }
+  auto get_chunk_size() const { return m_chunk_size; }
 
   /** Accumulate values.
    *  @param acc    Initial value for the accumulator.
@@ -97,7 +78,7 @@ class Observable_stat {
   Utils::Span<double> coulomb;
   /** Contribution(s) from dipolar interactions. */
   Utils::Span<double> dipolar;
-  /** Contribution(s) from virtual sites (accumulated). */
+  /** Contribution from virtual sites (accumulated). */
   Utils::Span<double> virtual_sites;
   /** Contribution from external fields (accumulated). */
   Utils::Span<double> external_fields;
@@ -107,17 +88,16 @@ class Observable_stat {
   Utils::Span<double> non_bonded_inter;
 
   /** Get contribution from a bonded interaction */
-  Utils::Span<double> bonded_contribution(int bond_id) {
-    return Utils::Span<double>(bonded.data() + m_chunk_size * bond_id,
-                               m_chunk_size);
+  Utils::Span<double> bonded_contribution(int bond_id) const {
+    auto const offset = m_chunk_size * static_cast<std::size_t>(bond_id);
+    return Utils::Span<double>(bonded.data() + offset, m_chunk_size);
   }
 
   void add_non_bonded_contribution(int type1, int type2,
                                    Utils::Span<const double> data) {
-    auto const dest =
-        (type1 == type2)
-            ? non_bonded_contribution(non_bonded_intra, type1, type2)
-            : non_bonded_contribution(non_bonded_inter, type1, type2);
+    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);
 
     boost::transform(dest, data, dest.begin(), std::plus<>{});
   }
@@ -137,6 +117,9 @@ class Observable_stat {
                                                     int type2) const {
     return non_bonded_contribution(non_bonded_inter, type1, type2);
   }
+
+  /** MPI reduction. */
+  void mpi_reduce();
 };
 
 #endif // ESPRESSO_OBSERVABLE_STAT_HPP
diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp
index 22b1db2f30a..f94cd656b08 100644
--- a/src/core/constraints/ShapeBasedConstraint.cpp
+++ b/src/core/constraints/ShapeBasedConstraint.cpp
@@ -20,6 +20,7 @@
 #include "ShapeBasedConstraint.hpp"
 
 #include "BoxGeometry.hpp"
+#include "Observable_stat.hpp"
 #include "communication.hpp"
 #include "config.hpp"
 #include "dpd.hpp"
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index 0c4e20d495a..e5b67539bb8 100644
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -33,25 +33,24 @@
 #include "forces.hpp"
 #include "integrate.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "reduce_observable_stat.hpp"
 
 #include "short_range_loop.hpp"
 
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 
+#include <memory>
+
 ActorList energyActors;
 
-/** Energy of the system */
-Observable_stat obs_energy{1};
+static std::shared_ptr<Observable_stat> calculate_energy_local() {
 
-Observable_stat const &get_obs_energy() { return obs_energy; }
+  auto obs_energy_ptr = std::make_shared<Observable_stat>(1);
 
-void energy_calc(const double time) {
   if (!interactions_sanity_checks())
-    return;
+    return obs_energy_ptr;
 
-  obs_energy = Observable_stat{1};
+  auto &obs_energy = *obs_energy_ptr;
 
 #ifdef CUDA
   clear_energy_on_GPU();
@@ -65,12 +64,15 @@ void energy_calc(const double time) {
 
   on_observable_calc();
 
-  for (auto const &p : cell_structure.local_particles()) {
+  auto const local_parts = cell_structure.local_particles();
+
+  for (auto const &p : local_parts) {
     obs_energy.kinetic[0] += calc_kinetic_energy(p);
   }
 
   short_range_loop(
-      [](Particle &p1, int bond_id, Utils::Span<Particle *> partners) {
+      [&obs_energy](Particle const &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) {
@@ -79,16 +81,23 @@ void energy_calc(const double time) {
         }
         return true;
       },
-      [](Particle const &p1, Particle const &p2, Distance const &d) {
+      [&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,
                                    obs_energy);
       },
       maximal_cutoff(), maximal_cutoff_bonded());
 
-  calc_long_range_energies(cell_structure.local_particles());
+#ifdef ELECTROSTATICS
+  /* calculate k-space part of electrostatic interaction. */
+  obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(local_parts);
+#endif
+
+#ifdef DIPOLES
+  /* calculate k-space part of magnetostatic interaction. */
+  obs_energy.dipolar[1] = Dipole::calc_energy_long_range(local_parts);
+#endif
 
-  auto local_parts = cell_structure.local_particles();
-  Constraints::constraints.add_energy(local_parts, time, obs_energy);
+  Constraints::constraints.add_energy(local_parts, get_sim_time(), obs_energy);
 
 #ifdef CUDA
   auto const energy_host = copy_energy_from_GPU();
@@ -98,37 +107,22 @@ void energy_calc(const double time) {
     obs_energy.dipolar[1] += energy_host.dipolar;
 #endif
 
-  /* gather data */
-  auto obs_energy_res = reduce(comm_cart, obs_energy);
-  if (obs_energy_res) {
-    std::swap(obs_energy, *obs_energy_res);
-  }
+  obs_energy.mpi_reduce();
+  return obs_energy_ptr;
 }
 
-void update_energy_local(int, int) { energy_calc(get_sim_time()); }
+REGISTER_CALLBACK_MASTER_RANK(calculate_energy_local)
 
-REGISTER_CALLBACK(update_energy_local)
-
-void update_energy() { mpi_call_all(update_energy_local, -1, -1); }
-
-void calc_long_range_energies(const ParticleRange &particles) {
-#ifdef ELECTROSTATICS
-  /* calculate k-space part of electrostatic interaction. */
-  obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(particles);
-#endif
-
-#ifdef DIPOLES
-  /* calculate k-space part of magnetostatic interaction. */
-  obs_energy.dipolar[1] = Dipole::calc_energy_long_range(particles);
-#endif
+std::shared_ptr<Observable_stat> calculate_energy() {
+  return mpi_call(Communication::Result::master_rank, calculate_energy_local);
 }
 
 double calculate_current_potential_energy_of_system() {
-  update_energy();
-  return obs_energy.accumulate(-obs_energy.kinetic[0]);
+  auto const obs_energy = calculate_energy();
+  return obs_energy->accumulate(-obs_energy->kinetic[0]);
 }
 
 double observable_compute_energy() {
-  update_energy();
-  return obs_energy.accumulate(0);
+  auto const obs_energy = calculate_energy();
+  return obs_energy->accumulate(0);
 }
diff --git a/src/core/energy.hpp b/src/core/energy.hpp
index 00c439118ac..d9ad254c7a2 100644
--- a/src/core/energy.hpp
+++ b/src/core/energy.hpp
@@ -28,22 +28,14 @@
 #define _ENERGY_H
 
 #include "Observable_stat.hpp"
-#include "ParticleRange.hpp"
 #include "actor/ActorList.hpp"
 
+#include <memory>
+
 extern ActorList energyActors;
 
 /** 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);
+std::shared_ptr<Observable_stat> calculate_energy();
 
 /** Calculate the total energy of the system. */
 double calculate_current_potential_energy_of_system();
diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp
index cf6822297b6..0b7680a1683 100644
--- a/src/core/pressure.cpp
+++ b/src/core/pressure.cpp
@@ -34,7 +34,6 @@
 #include "grid.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 #include "pressure_inline.hpp"
-#include "reduce_observable_stat.hpp"
 #include "virtual_sites.hpp"
 
 #include "short_range_loop.hpp"
@@ -50,42 +49,30 @@
 #include <algorithm>
 #include <cmath>
 #include <cstddef>
+#include <memory>
 #include <utility>
 
-/** Pressure tensor of the system */
-Observable_stat obs_pressure{9};
+static std::shared_ptr<Observable_stat> calculate_pressure_local() {
 
-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. */
-  auto const coulomb_pressure = Coulomb::calc_pressure_long_range(particles);
-  boost::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9);
-#endif
-#ifdef DIPOLES
-  /* calculate k-space part of magnetostatic interaction. */
-  Dipole::calc_pressure_long_range();
-#endif
-}
-
-void pressure_calc() {
-  auto const volume = box_geo.volume();
+  auto obs_pressure_ptr = std::make_shared<Observable_stat>(9);
 
   if (!interactions_sanity_checks())
-    return;
+    return obs_pressure_ptr;
 
-  obs_pressure = Observable_stat{9};
+  auto &obs_pressure = *obs_pressure_ptr;
 
   on_observable_calc();
 
-  for (auto const &p : cell_structure.local_particles()) {
+  auto const volume = box_geo.volume();
+  auto const local_parts = cell_structure.local_particles();
+
+  for (auto const &p : local_parts) {
     add_kinetic_virials(p, obs_pressure);
   }
 
   short_range_loop(
-      [](Particle &p1, int bond_id, Utils::Span<Particle *> partners) {
+      [&obs_pressure](Particle const &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) {
@@ -100,13 +87,22 @@ void pressure_calc() {
         }
         return true;
       },
-      [](Particle &p1, Particle &p2, Distance const &d) {
+      [&obs_pressure](Particle const &p1, Particle const &p2,
+                      Distance const &d) {
         add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2),
                                     obs_pressure);
       },
       maximal_cutoff(), maximal_cutoff_bonded());
 
-  calc_long_range_virials(cell_structure.local_particles());
+#ifdef ELECTROSTATICS
+  /* calculate k-space part of electrostatic interaction. */
+  auto const coulomb_pressure = Coulomb::calc_pressure_long_range(local_parts);
+  boost::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9);
+#endif
+#ifdef DIPOLES
+  /* calculate k-space part of magnetostatic interaction. */
+  Dipole::calc_pressure_long_range();
+#endif
 
 #ifdef VIRTUAL_SITES
   if (!obs_pressure.virtual_sites.empty()) {
@@ -117,24 +113,21 @@ void pressure_calc() {
 
   obs_pressure.rescale(volume);
 
-  /* gather data */
-  auto obs_pressure_res = reduce(comm_cart, obs_pressure);
-  if (obs_pressure_res) {
-    std::swap(obs_pressure, *obs_pressure_res);
-  }
+  obs_pressure.mpi_reduce();
+  return obs_pressure_ptr;
 }
 
-void update_pressure_local(int, int) { pressure_calc(); }
-
-REGISTER_CALLBACK(update_pressure_local)
+REGISTER_CALLBACK_MASTER_RANK(calculate_pressure_local)
 
-void update_pressure() { mpi_call_all(update_pressure_local, -1, -1); }
+std::shared_ptr<Observable_stat> calculate_pressure() {
+  return mpi_call(Communication::Result::master_rank, calculate_pressure_local);
+}
 
 Utils::Vector9d observable_compute_pressure_tensor() {
-  update_pressure();
+  auto const obs_pressure = calculate_pressure();
   Utils::Vector9d pressure_tensor{};
   for (std::size_t j = 0; j < 9; j++) {
-    pressure_tensor[j] = obs_pressure.accumulate(0, j);
+    pressure_tensor[j] = obs_pressure->accumulate(0, j);
   }
   return pressure_tensor;
 }
diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp
index 9b08b22b340..426381ceb9e 100644
--- a/src/core/pressure.hpp
+++ b/src/core/pressure.hpp
@@ -29,14 +29,10 @@
 
 #include <utils/Vector.hpp>
 
-/** Parallel pressure calculation from a virial expansion. */
-void pressure_calc();
-
-/** Run @ref pressure_calc in parallel. */
-void update_pressure();
+#include <memory>
 
-/** Return the pressure observable. */
-Observable_stat const &get_obs_pressure();
+/** Parallel pressure calculation from a virial expansion. */
+std::shared_ptr<Observable_stat> calculate_pressure();
 
 /** Helper function for @ref Observables::PressureTensor. */
 Utils::Vector9d observable_compute_pressure_tensor();
diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp
index 54349d0cc79..eef2ff4718f 100644
--- a/src/core/pressure_inline.hpp
+++ b/src/core/pressure_inline.hpp
@@ -138,7 +138,8 @@ calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams,
 }
 
 inline boost::optional<Utils::Matrix<double, 3, 3>>
-calc_bonded_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle &p1,
+calc_bonded_pressure_tensor(Bonded_IA_Parameters const &iaparams,
+                            Particle const &p1,
                             Utils::Span<Particle *> partners) {
   switch (number_of_partners(iaparams)) {
   case 1:
diff --git a/src/core/reduce_observable_stat.cpp b/src/core/reduce_observable_stat.cpp
deleted file mode 100644
index aa43c4dbff4..00000000000
--- a/src/core/reduce_observable_stat.cpp
+++ /dev/null
@@ -1,45 +0,0 @@
-/*
- * Copyright (C) 2010-2020 The ESPResSo project
- * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
- *   Max-Planck-Institute for Polymer Research, Theory Group
- *
- * This file is part of ESPResSo.
- *
- * ESPResSo is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * ESPResSo is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#include "reduce_observable_stat.hpp"
-
-#include <boost/mpi/collectives/reduce.hpp>
-#include <boost/optional.hpp>
-
-#include <functional>
-
-/** Reduce contributions from all MPI ranks. */
-boost::optional<Observable_stat> reduce(boost::mpi::communicator const &comm,
-                                        Observable_stat const &obs) {
-  if (comm.rank() == 0) {
-    Observable_stat res{obs.chunk_size()};
-
-    boost::mpi::reduce(comm, obs.data_().data(), obs.data_().size(),
-                       res.data_().data(), std::plus<>{}, 0);
-
-    return res;
-  }
-
-  boost::mpi::reduce(comm, obs.data_().data(), obs.data_().size(),
-                     std::plus<>{}, 0);
-
-  return {};
-}
diff --git a/src/core/reduce_observable_stat.hpp b/src/core/reduce_observable_stat.hpp
deleted file mode 100644
index e54c02f3b32..00000000000
--- a/src/core/reduce_observable_stat.hpp
+++ /dev/null
@@ -1,33 +0,0 @@
-/*
- * Copyright (C) 2010-2020 The ESPResSo project
- * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
- *   Max-Planck-Institute for Polymer Research, Theory Group
- *
- * This file is part of ESPResSo.
- *
- * ESPResSo is free software: you can redistribute it and/or modify
- * it under the terms of the GNU General Public License as published by
- * the Free Software Foundation, either version 3 of the License, or
- * (at your option) any later version.
- *
- * ESPResSo is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU General Public License for more details.
- *
- * You should have received a copy of the GNU General Public License
- * along with this program.  If not, see <http://www.gnu.org/licenses/>.
- */
-
-#ifndef ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP
-#define ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP
-
-#include "Observable_stat.hpp"
-
-#include <boost/mpi/communicator.hpp>
-#include <boost/optional.hpp>
-
-/** Reduce contributions from all MPI ranks. */
-boost::optional<Observable_stat> reduce(boost::mpi::communicator const &comm,
-                                        Observable_stat const &obs);
-#endif // ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP
diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
index 2ec3a6725a5..ddd295170b5 100644
--- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
+++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp
@@ -79,7 +79,6 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
   auto const box_l = 8.;
   auto const box_center = box_l / 2.;
   espresso::system->set_box_l(Utils::Vector3d::broadcast(box_l));
-  auto const &obs_energy = get_obs_energy();
 
   // particle properties
   auto const pid1 = 9;
@@ -148,8 +147,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
       set_particle_v(pid2, {static_cast<double>(i), 0., 0.});
       auto const &p = get_particle_data(pid2);
       auto const kinetic_energy = 0.5 * p.p.mass * p.m.v.norm2();
-      update_energy();
-      BOOST_CHECK_CLOSE(obs_energy.kinetic[0], kinetic_energy, tol);
+      auto const obs_energy = calculate_energy();
+      BOOST_CHECK_CLOSE(obs_energy->kinetic[0], kinetic_energy, tol);
       BOOST_CHECK_CLOSE(observable_compute_energy(), kinetic_energy, tol);
     }
   }
@@ -179,12 +178,12 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
     auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift);
 
     // measure energies
-    update_energy();
+    auto const obs_energy = calculate_energy();
     for (int i = 0; i < n_pairs; ++i) {
       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, 500. * tol);
-      BOOST_CHECK_CLOSE(obs_energy.non_bonded_intra[i], ref_intra, 500. * tol);
+      BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 500. * tol);
+      BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 500. * tol);
     }
   }
 #endif // LENNARD_JONES
@@ -207,7 +206,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
     add_particle_bond(pid2, std::vector<int>{fene_bond_id, pid3});
 
     // measure energies
-    update_energy();
+    auto const obs_energy = calculate_energy();
     auto const &p1 = get_particle_data(pid1);
     auto const &p2 = get_particle_data(pid2);
     auto const &p3 = get_particle_data(pid3);
@@ -216,9 +215,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
     auto const fene_energy =
         -0.5 * fene_bond.k * Utils::sqr(fene_bond.drmax) *
         std::log(1.0 - Utils::sqr((dist - fene_bond.r0) / fene_bond.drmax));
-    BOOST_CHECK_CLOSE(obs_energy.bonded[none_bond_id], none_energy, 0.0);
-    BOOST_CHECK_CLOSE(obs_energy.bonded[harm_bond_id], harm_energy, 40. * tol);
-    BOOST_CHECK_CLOSE(obs_energy.bonded[fene_bond_id], fene_energy, 40. * tol);
+    BOOST_CHECK_CLOSE(obs_energy->bonded[none_bond_id], none_energy, 0.0);
+    BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 40. * tol);
+    BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 40. * tol);
   }
 
   // check electrostatics
@@ -246,11 +245,11 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory,
       place_particle(pid2, pos2);
       auto const r = (pos2 - pos1).norm();
       // check P3M energy
-      update_energy();
+      auto const obs_energy = calculate_energy();
       // at very short distances, the real-space contribution to
       // the energy is much larger than the k-space contribution
       auto const energy_ref = -prefactor / r;
-      auto const energy_p3m = obs_energy.coulomb[0] + obs_energy.coulomb[1];
+      auto const energy_p3m = obs_energy->coulomb[0] + obs_energy->coulomb[1];
       BOOST_CHECK_CLOSE(energy_p3m, energy_ref, 0.01);
     }
   }
diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd
index 2ab4f40e24c..c42750b5fdc 100644
--- a/src/python/espressomd/analyze.pxd
+++ b/src/python/espressomd/analyze.pxd
@@ -25,6 +25,11 @@ from .utils cimport Vector3i, Vector3d, Vector9d, Span
 from .utils cimport create_nparray_from_double_span
 from libcpp.vector cimport vector  # import std::vector as vector
 from libcpp cimport bool as cbool
+from libcpp.memory cimport shared_ptr
+from .interactions cimport bonded_ia_params_is_type
+from .interactions cimport bonded_ia_params_size
+from .interactions cimport CoreNoneBond
+include "myconfig.pxi"
 
 cdef extern from "<array>" namespace "std" nogil:
     cdef cppclass array4 "std::array<double, 4>":
@@ -58,7 +63,7 @@ cdef extern from "Observable_stat.hpp":
         Span[double] bonded_contribution(int bond_id)
         Span[double] non_bonded_intra_contribution(int type1, int type2)
         Span[double] non_bonded_inter_contribution(int type1, int type2)
-        size_t chunk_size()
+        size_t get_chunk_size()
 
 cdef extern from "statistics.hpp":
     cdef void calc_structurefactor(PartCfg & , const vector[int] & p_types, int order, vector[double] & wavevectors, vector[double] & intensities) except +
@@ -82,18 +87,16 @@ cdef extern from "statistics_chain.hpp":
     array2 calc_rh(int, int, int)
 
 cdef extern from "pressure.hpp":
-    cdef void update_pressure()
-    cdef const Observable_stat & get_obs_pressure()
+    cdef shared_ptr[Observable_stat] calculate_pressure()
 
 cdef extern from "energy.hpp":
-    cdef void update_energy()
-    cdef const Observable_stat & get_obs_energy()
+    cdef shared_ptr[Observable_stat] calculate_energy()
     double calculate_current_potential_energy_of_system()
 
 cdef extern from "dpd.hpp":
     Vector9d dpd_stress()
 
-cdef inline get_obs_contribs(Span[double] contributions, int size,
+cdef inline get_obs_contribs(Span[double] contributions,
                              cbool calc_scalar_pressure):
     """
     Convert an Observable_stat range of contributions into a correctly
@@ -112,7 +115,7 @@ cdef inline get_obs_contribs(Span[double] contributions, int size,
     """
     cdef np.ndarray value
     value = create_nparray_from_double_span(contributions)
-    if size == 9:
+    if contributions.size() == 9:
         if calc_scalar_pressure:
             return np.einsum('...ii', value.reshape((-1, 3, 3))) / 3
         else:
@@ -120,7 +123,7 @@ cdef inline get_obs_contribs(Span[double] contributions, int size,
     else:
         return value
 
-cdef inline get_obs_contrib(Span[double] contribution, int size,
+cdef inline get_obs_contrib(Span[double] contribution,
                             cbool calc_scalar_pressure):
     """
     Convert an Observable_stat contribution into a correctly
@@ -130,15 +133,137 @@ cdef inline get_obs_contrib(Span[double] contribution, int size,
     ----------
     contributions : (N,) array_like of :obj:`float`
         Flattened array of energy/pressure contributions from an observable.
-    size : :obj:`int`, \{1, 9\}
-        Dimensionality of the data.
     calc_scalar_pressure : :obj:`bool`
         Whether to calculate a scalar pressure (only relevant when
         ``contributions`` is a pressure tensor).
 
     """
     cdef np.ndarray value
-    value = get_obs_contribs(contribution, size, calc_scalar_pressure)
+    value = get_obs_contribs(contribution, calc_scalar_pressure)
     if value.shape[0] == 1:
         return value[0]
     return value
+
+cdef inline observable_stat_matrix(size_t size, cbool calc_sp):
+    if size == 9 and not calc_sp:
+        return np.zeros((3, 3), dtype=float)
+    else:
+        return 0.0
+
+cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp):
+    """Transform an ``Observable_stat`` object to a python dict.
+
+    Parameters
+    ----------
+    obs :
+        Core observable.
+    calc_sp : :obj:`bool`
+        Whether to calculate a scalar pressure (only relevant when
+        ``obs`` is a pressure tensor observable).
+
+    Returns
+    -------
+    :obj:`dict`
+        A dictionary with the following keys:
+
+        * ``"total"``: total contribution
+        * ``"kinetic"``: kinetic contribution
+        * ``"bonded"``: total bonded contribution
+        * ``"bonded", <bond_type>``: bonded contribution which arises from the given bond_type
+        * ``"nonbonded"``: total nonbonded contribution
+        * ``"nonbonded", <type_i>, <type_j>``: nonbonded contribution which arises from the interactions between type_i and type_j
+        * ``"nonbonded_intra", <type_i>, <type_j>``: nonbonded contribution between short ranged forces between type i and j and with the same mol_id
+        * ``"nonbonded_inter", <type_i>, <type_j>``: nonbonded contribution between short ranged forces between type i and j and different mol_ids
+        * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method
+        * ``"coulomb", <i>``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``)
+        * ``"dipolar"``: dipolar contribution, how it is calculated depends on the method
+        * ``"dipolar", <i>``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``)
+        * ``"virtual_sites"``: virtual sites contribution
+        * ``"virtual_sites", <i>``: contribution from virtual site i
+
+    """
+
+    cdef size_t i
+    cdef size_t j
+    cdef size_t obs_dim = obs.get_chunk_size()
+    cdef size_t n_bonded = bonded_ia_params_size()
+    cdef size_t n_nonbonded = <size_t > max_seen_particle_type
+    cdef double[9] total
+
+    # Dict to store the results
+    p = {}
+
+    # Total contribution
+
+    for i in range(obs_dim):
+        total[i] = obs.accumulate(0.0, i)
+    p["total"] = get_obs_contrib(Span[double](total, obs_dim), calc_sp)
+
+    # Kinetic
+    p["kinetic"] = get_obs_contrib(obs.kinetic, calc_sp)
+
+    # External
+    p["external_fields"] = get_obs_contrib(obs.external_fields, calc_sp)
+
+    # Bonded
+    total_bonded = observable_stat_matrix(obs_dim, calc_sp)
+    for i in range(n_bonded):
+        if not bonded_ia_params_is_type[CoreNoneBond](i):
+            val = get_obs_contrib(obs.bonded_contribution(i), calc_sp)
+            p["bonded", i] = val
+            total_bonded += val
+    p["bonded"] = total_bonded
+
+    # Non-Bonded interactions, total as well as intra and inter molecular
+    total_intra = observable_stat_matrix(obs_dim, calc_sp)
+    total_inter = observable_stat_matrix(obs_dim, calc_sp)
+    for i in range(n_nonbonded):
+        for j in range(i, n_nonbonded):
+            intra = get_obs_contrib(
+                obs.non_bonded_intra_contribution(i, j), calc_sp)
+            total_intra += intra
+            p["non_bonded_intra", i, j] = intra
+            inter = get_obs_contrib(
+                obs.non_bonded_inter_contribution(i, j), calc_sp)
+            total_inter += inter
+            p["non_bonded_inter", i, j] = inter
+            p["non_bonded", i, j] = intra + inter
+
+    p["non_bonded_intra"] = total_intra
+    p["non_bonded_inter"] = total_inter
+    p["non_bonded"] = total_intra + total_inter
+
+    # Electrostatics
+    IF ELECTROSTATICS == 1:
+        cdef np.ndarray coulomb
+        coulomb = get_obs_contribs(obs.coulomb, calc_sp)
+        for i in range(coulomb.shape[0]):
+            p["coulomb", i] = coulomb[i]
+        p["coulomb"] = np.sum(coulomb, axis=0)
+
+    # Dipoles
+    IF DIPOLES == 1:
+        cdef np.ndarray dipolar
+        dipolar = get_obs_contribs(obs.dipolar, calc_sp)
+        for i in range(dipolar.shape[0]):
+            p["dipolar", i] = dipolar[i]
+        p["dipolar"] = np.sum(dipolar, axis=0)
+
+    # virtual sites
+    IF VIRTUAL_SITES == 1:
+        cdef np.ndarray virtual_sites
+        virtual_sites = get_obs_contribs(obs.virtual_sites, calc_sp)
+        for i in range(virtual_sites.shape[0]):
+            p["virtual_sites", i] = virtual_sites[i]
+        p["virtual_sites"] = np.sum(virtual_sites, axis=0)
+
+    return p
+
+cdef inline get_scalar_pressure():
+    return Observable_stat_to_dict(calculate_pressure().get(), True)
+
+cdef inline get_pressure_tensor():
+    return Observable_stat_to_dict(calculate_pressure().get(), False)
+
+cdef inline get_energy():
+    return Observable_stat_to_dict(calculate_energy().get(), False)
diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx
index 312ed782a2c..50df235268a 100644
--- a/src/python/espressomd/analyze.pyx
+++ b/src/python/espressomd/analyze.pyx
@@ -20,16 +20,11 @@
 include "myconfig.pxi"
 from . cimport analyze
 from libcpp.vector cimport vector  # import std::vector as vector
-from libcpp cimport bool as cbool
-from .interactions cimport bonded_ia_params_is_type
-from .interactions cimport bonded_ia_params_size
-from .interactions cimport CoreNoneBond
 import numpy as np
 cimport numpy as np
 import scipy.signal
 from .grid cimport box_geo
 
-from collections import OrderedDict
 from .system import System
 from .utils import array_locked, is_valid_type, handle_errors
 from .utils cimport Vector3i, Vector3d, Vector9d
@@ -93,129 +88,6 @@ def autocorrelation(time_series):
                          f"are supported, got shape {time_series.shape}")
 
 
-cdef _Observable_stat_to_dict(Observable_stat obs, cbool calc_sp):
-    """Transform an Observable_stat struct to a python dict.
-
-    Parameters
-    ----------
-    obs :
-        Core observable.
-    calc_sp : :obj:`bool`
-        Whether to calculate a scalar pressure (only relevant when
-        ``obs`` is a pressure tensor observable).
-
-    Returns
-    -------
-    :obj:`dict`
-        A dictionary with the following keys:
-
-        * ``"total"``: total contribution
-        * ``"kinetic"``: kinetic contribution
-        * ``"bonded"``: total bonded contribution
-        * ``"bonded", <bond_type>``: bonded contribution which arises from the given bond_type
-        * ``"nonbonded"``: total nonbonded contribution
-        * ``"nonbonded", <type_i>, <type_j>``: nonbonded contribution which arises from the interactions between type_i and type_j
-        * ``"nonbonded_intra", <type_i>, <type_j>``: nonbonded contribution between short ranged forces between type i and j and with the same mol_id
-        * ``"nonbonded_inter", <type_i>, <type_j>``: nonbonded contribution between short ranged forces between type i and j and different mol_ids
-        * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method
-        * ``"coulomb", <i>``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``)
-        * ``"dipolar"``: dipolar contribution, how it is calculated depends on the method
-        * ``"dipolar", <i>``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``)
-        * ``"virtual_sites"``: virtual sites contribution
-        * ``"virtual_sites", <i>``: contribution from virtual site i
-
-    """
-
-    size = obs.chunk_size()
-
-    def set_initial():
-        if size == 9 and not calc_sp:
-            return np.zeros((3, 3), dtype=float)
-        else:
-            return 0.0
-
-    cdef int i
-    cdef int j
-
-    # Dict to store the results
-    p = OrderedDict()
-
-    # Total contribution
-    if size == 1:
-        p["total"] = obs.accumulate()
-    else:
-        total = np.zeros((9,), dtype=float)
-        for i in range(9):
-            total[i] = obs.accumulate(0.0, i)
-        total = total.reshape((3, 3))
-        if calc_sp:
-            p["total"] = np.einsum('ii', total) / 3
-        else:
-            p["total"] = total
-
-    # Kinetic
-    p["kinetic"] = get_obs_contrib(obs.kinetic, size, calc_sp)
-
-    # External
-    p["external_fields"] = get_obs_contrib(obs.external_fields, size, calc_sp)
-
-    # Bonded
-    total_bonded = set_initial()
-    for i in range(bonded_ia_params_size()):
-        if not bonded_ia_params_is_type[CoreNoneBond](i):
-            val = get_obs_contrib(obs.bonded_contribution(i), size, calc_sp)
-            p["bonded", i] = val
-            total_bonded += val
-    p["bonded"] = total_bonded
-
-    # Non-Bonded interactions, total as well as intra and inter molecular
-    total_intra = set_initial()
-    total_inter = set_initial()
-    for i in range(analyze.max_seen_particle_type):
-        for j in range(i, analyze.max_seen_particle_type):
-            intra = get_obs_contrib(
-                obs.non_bonded_intra_contribution(
-                    i, j), size, calc_sp)
-            total_intra += intra
-            p["non_bonded_intra", i, j] = intra
-            inter = get_obs_contrib(
-                obs.non_bonded_inter_contribution(
-                    i, j), size, calc_sp)
-            total_inter += inter
-            p["non_bonded_inter", i, j] = inter
-            p["non_bonded", i, j] = intra + inter
-
-    p["non_bonded_intra"] = total_intra
-    p["non_bonded_inter"] = total_inter
-    p["non_bonded"] = total_intra + total_inter
-
-    # Electrostatics
-    IF ELECTROSTATICS == 1:
-        cdef np.ndarray coulomb
-        coulomb = get_obs_contribs(obs.coulomb, size, calc_sp)
-        for i in range(coulomb.shape[0]):
-            p["coulomb", i] = coulomb[i]
-        p["coulomb"] = np.sum(coulomb, axis=0)
-
-    # Dipoles
-    IF DIPOLES == 1:
-        cdef np.ndarray dipolar
-        dipolar = get_obs_contribs(obs.dipolar, size, calc_sp)
-        for i in range(dipolar.shape[0]):
-            p["dipolar", i] = dipolar[i]
-        p["dipolar"] = np.sum(dipolar, axis=0)
-
-    # virtual sites
-    IF VIRTUAL_SITES == 1:
-        cdef np.ndarray virtual_sites
-        virtual_sites = get_obs_contribs(obs.virtual_sites, size, calc_sp)
-        for i in range(virtual_sites.shape[0]):
-            p["virtual_sites", i] = virtual_sites[i]
-        p["virtual_sites"] = np.sum(virtual_sites, axis=0)
-
-    return p
-
-
 class Analysis:
 
     def __init__(self, system):
@@ -385,10 +257,9 @@ class Analysis:
 
         """
 
-        # Update in ESPResSo core
-        analyze.update_pressure()
-
-        return _Observable_stat_to_dict(analyze.get_obs_pressure(), True)
+        obs = analyze.get_scalar_pressure()
+        handle_errors("calculate_pressure() failed")
+        return obs
 
     def pressure_tensor(self):
         """Calculate the instantaneous pressure_tensor (in parallel). This is
@@ -419,10 +290,9 @@ class Analysis:
 
         """
 
-        # Update in ESPResSo core
-        analyze.update_pressure()
-
-        return _Observable_stat_to_dict(analyze.get_obs_pressure(), False)
+        obs = analyze.get_pressure_tensor()
+        handle_errors("calculate_pressure() failed")
+        return obs
 
     IF DPD == 1:
         def dpd_stress(self):
@@ -470,10 +340,9 @@ class Analysis:
 
         """
 
-        analyze.update_energy()
-        handle_errors("calc_long_range_energies failed")
-
-        return _Observable_stat_to_dict(analyze.get_obs_energy(), False)
+        obs = analyze.get_energy()
+        handle_errors("calculate_energy() failed")
+        return obs
 
     def calc_re(self, chain_start=None, number_of_chains=None,
                 chain_length=None):