Create observables based on Observable_stat and remove state tracking (…

…#3723)

Fixes #2200, partial fix for #3718

Description of changes:
- removed state tracking from `Observable_stat`
- created `Observables` classes `PressureTensor`, `Pressure`, `Energy`
- renamed Stress Tensor to Pressure Tensor everywhere

doc/sphinx/analysis.rst

@@ -24,6 +24,7 @@ The direct analysis commands can be classified into two types:
 
     - :ref:`Energies`
     - :ref:`Pressure`
+    - :ref:`Pressure Tensor`
     - :ref:`Momentum of the System`
     - :ref:`Minimal distances between particles`
     - :ref:`Particles in the neighborhood`
@@ -32,7 +33,6 @@ The direct analysis commands can be classified into two types:
     - :ref:`Center of mass`
     - :ref:`Moment of inertia matrix`
     - :ref:`Gyration tensor`
-    - :ref:`Stress Tensor`
 
 .. _Energies:
 
@@ -41,7 +41,7 @@ Energies
 :meth:`espressomd.analyze.Analysis.energy`
 
 Returns the energies of the system.
-The different energetic contributions to the total energy can also be obtained (kinetic, bonded,non-bonded, Coulomb).
+The different energetic contributions to the total energy can also be obtained (kinetic, bonded, non-bonded, Coulomb).
 
 For example, ::
 
@@ -203,26 +203,26 @@ constraints of any kind are not currently accounted for in the pressure
 calculations. The pressure is no longer correct, e.g., when particles
 are confined to a plane.
 
-Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure formula. Later when the user calculates the ensemble average via e.g. :math:`\langle p \rangle \approx 1/N \sum_{i=1}^N p_i` however the ensemble average with all interactions present is performed. That means the contributions are not easy to interpret! Those are the contributions to the stress/pressure in a system where all interactions are present and therefore in a coupled system.
+Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure formula. Later when the user calculates the ensemble average via e.g. :math:`\langle p \rangle \approx 1/N \sum_{i=1}^N p_i` however the ensemble average with all interactions present is performed. That means the contributions are not easy to interpret! Those are the contributions to the pressure in a system where all interactions are present and therefore in a coupled system.
 
-.. _Stress Tensor:
+.. _Pressure Tensor:
 
-Stress Tensor
-~~~~~~~~~~~~~
-:meth:`espressomd.analyze.Analysis.stress_tensor`
+Pressure Tensor
+~~~~~~~~~~~~~~~
+:meth:`espressomd.analyze.Analysis.pressure_tensor`
 
-Computes the volume averaged instantaneous stress tensor of the system with options which are
-described by in :meth:`espressomd.analyze.Analysis.stress_tensor`. It is called a stress tensor but the sign convention follows that of a pressure tensor.
-In general do only use it for (on average) homogeneous systems. For inhomogeneous systems you need to use the local stress tensor.
+Computes the volume averaged instantaneous pressure tensor of the system with options which are
+described by in :meth:`espressomd.analyze.Analysis.pressure_tensor`.
+In general do only use it for (on average) homogeneous systems. For inhomogeneous systems you need to use the local pressure tensor.
 
-The instantaneous virial stress tensor is calculated by
+The instantaneous virial pressure tensor is calculated by
 
 .. math:: p_{(k,l)} = \frac{\sum_{i} {m_{i}v_{i}^{(k)}v_{i}^{(l)}}}{V} + \frac{\sum_{j>i}{F_{ij}^{(k)}r_{ij}^{(l)}}}{V}
 
 where the notation is the same as for the pressure. The superscripts :math:`k`
 and :math:`l` correspond to the components in the tensors and vectors.
 
-If electrostatic interactions are present then also the coulombic parts of the stress tensor need to be calculated. If P3M is present, then the instantaneous stress tensor is added to the above equation in accordance with :cite:`essmann95a` :
+If electrostatic interactions are present then also the coulombic parts of the pressure tensor need to be calculated. If P3M is present, then the instantaneous pressure tensor is added to the above equation in accordance with :cite:`essmann95a` :
 
 .. math :: p^\text{Coulomb, P3M}_{(k,l)} =p^\text{Coulomb, P3M, dir}_{(k,l)} + p^\text{Coulomb, P3M, rec}_{(k,l)},
 
@@ -239,12 +239,12 @@ The long ranged (k-space) part is given by:
 
 where :math:`S(\vec{k})` is the Fourier transformed charge density. Compared to Essmann we do not have the contribution :math:`p^\text{corr}_{k,l}` since we want to calculate the pressure that arises from all particles in the system.
 
-Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial stress tensor formula.
-Later when the user calculates the stress tensor via :math:`\langle p_{(k,l)}\rangle  \approx 1/N \sum_{i=1}^N p_{k,l}` however the ensemble average with all interactions present is performed.
-That means the contributions are not easy to interpret! Those are the contributions to the stress/pressure in a system where all interactions are present and therefore in a coupled system.
+Note: The different contributions which are returned are the summands that arise from force splitting :math:`\vec{F}_{i,j}={\vec{F}_{i,j}}_\text{bonded}+{\vec{F}_{i,j}}_\text{nonbonded}+...` in the virial pressure tensor formula.
+Later when the user calculates the pressure tensor via :math:`\langle p_{(k,l)}\rangle  \approx 1/N \sum_{i=1}^N p_{k,l}` however the ensemble average with all interactions present is performed.
+That means the contributions are not easy to interpret! Those are the contributions to the pressure in a system where all interactions are present and therefore in a coupled system.
 
 Note that the angular velocities of the particles are not included in
-the calculation of the stress tensor.
+the calculation of the pressure tensor.
 
 .. _Chains:
 
@@ -500,7 +500,11 @@ documentation for all available observables in :mod:`espressomd.observables`.
 
 - System-wide observables
 
-   - :class:`~espressomd.observables.StressTensor`: Total stress tensor (see :ref:`stress tensor`)
+   - :class:`~espressomd.observables.Energy`: Total energy (see :ref:`Energies`)
+
+   - :class:`~espressomd.observables.Pressure`: Total scalar pressure (see :ref:`Pressure`)
+
+   - :class:`~espressomd.observables.PressureTensor`: Total pressure tensor (see :ref:`Pressure Tensor`)
 
    - :class:`~espressomd.observables.DPDStress`
 

doc/sphinx/lb.rst

@@ -181,12 +181,12 @@ Reading and setting properties of single lattice nodes
 
 Appending three indices to the ``lb`` object returns an object that represents the selected LB grid node and allows one to access all of its properties::
 
-    lb[x, y, z].density     # fluid density (one scalar for LB and CUDA)
-    lb[x, y, z].velocity    # fluid velocity (a numpy array of three floats)
-    lb[x, y, z].stress      # fluid pressure tensor (a symmetric 3x3 numpy array of floats)
-    lb[x, y, z].stress_neq  # nonequilbrium part of the pressure tensor (as above)
-    lb[x, y, z].boundary    # flag indicating whether the node is fluid or boundary (fluid: boundary=0, boundary: boundary != 0)
-    lb[x, y, z].population  # 19 LB populations (a numpy array of 19 floats, check order from the source code)
+    lb[x, y, z].density              # fluid density (one scalar for LB and CUDA)
+    lb[x, y, z].velocity             # fluid velocity (a numpy array of three floats)
+    lb[x, y, z].pressure_tensor      # fluid pressure tensor (a symmetric 3x3 numpy array of floats)
+    lb[x, y, z].pressure_tensor_neq  # nonequilibrium part of the pressure tensor (as above)
+    lb[x, y, z].boundary             # flag indicating whether the node is fluid or boundary (fluid: boundary=0, boundary: boundary != 0)
+    lb[x, y, z].population           # 19 LB populations (a numpy array of 19 floats, check order from the source code)
 
 All of these properties can be read and used in further calculations. Only the property ``population`` can be modified. The indices ``x,y,z`` are integers and enumerate the LB nodes in the three directions, starts with 0. To modify ``boundary``, refer to :ref:`Setting up boundary conditions`.
 

doc/sphinx/particles.rst

@@ -398,7 +398,7 @@ Please note:
    needs to be adapted.
 
 -  The presence of rigid bodies constructed by means of virtual sites
-   adds a contribution to the pressure and stress tensor.
+   adds a contribution to the scalar pressure and pressure tensor.
 
 .. _Inertialess lattice-Boltzmann tracers:
 

doc/tutorials/04-lattice_boltzmann/04-lattice_boltzmann_part1.ipynb

@@ -300,8 +300,8 @@
     "The possible properties are:\n",
     "\n",
     "+ <tt>velocity</tt>: the fluid velocity (list of three floats)\n",
-    "+ <tt>stress</tt>: the pressure tensor (list of six floats: $\\Pi_{xx}$, $\\Pi_{xy}$, $\\Pi_{yy}$, $\\Pi_{xz}$, $\\Pi_{yz}$, and $\\Pi_{zz}$)\n",
-    "+ <tt>stress_neq</tt>: the nonequilibrium part of the pressure tensor, components as above.\n",
+    "+ <tt>pressure_tensor</tt>: the pressure tensor (3x3 matrix)\n",
+    "+ <tt>pressure_tensor_neq</tt>: the nonequilibrium part of the pressure tensor (3x3 matrix).\n",
     "+ <tt>population</tt>: the 19 populations of the D3Q19 lattice.\n",
     "+ <tt>boundary</tt>: the boundary flag.\n",
     "+ <tt>density</tt>: the local density."

src/core/Observable_stat.cpp

@@ -23,30 +23,13 @@
 
 #include "config.hpp"
 
+#include "communication.hpp"
+
 #include "bonded_interactions/bonded_interaction_data.hpp"
 #include "electrostatics_magnetostatics/coulomb.hpp"
 #include "electrostatics_magnetostatics/dipole.hpp"
 #include "nonbonded_interactions/nonbonded_interaction_data.hpp"
 
-#include <boost/mpi/communicator.hpp>
-
-extern boost::mpi::communicator comm_cart;
-
-/** Tracker of observables */
-std::vector<Observable_stat_wrapper *> &registered_observables() {
-  static std::vector<Observable_stat_wrapper *> s_registered_observables;
-  return s_registered_observables;
-}
-
-void Observable_stat_wrapper::register_obs() {
-  registered_observables().push_back(this);
-}
-
-void invalidate_obs() {
-  for (auto *obs : registered_observables())
-    obs->is_initialized = false;
-}
-
 void Observable_stat::resize() {
   // number of chunks for different interaction types
   auto const n_coulomb =

src/core/Observable_stat.hpp

@@ -28,9 +28,9 @@
 #include <utility>
 #include <vector>
 
-/** Cache for system statistics.
+/** Observable for system statistics.
  *  Store unidimensional (energy, scalar pressure) and multi-dimensional
- *  (stress tensors) properties of the system and provide accumulation and
+ *  (pressure tensor) properties of the system and provide accumulation and
  *  reduction functionality.
  */
 class Observable_stat_base {
@@ -106,9 +106,7 @@ class Observable_stat_base {
   }
 };
 
-/** Structure used to cache the results of the scalar pressure, stress tensor
- *  and energy calculations.
- */
+/** Observable for the scalar pressure, pressure tensor and energy. */
 class Observable_stat : public Observable_stat_base {
 private:
   /** Whether this observable is a pressure or energy observable */
@@ -151,8 +149,8 @@ class Observable_stat : public Observable_stat_base {
   }
 };
 
-/** Structure used only in the pressure and stress tensor calculation to
- *  distinguish non-bonded intra- and inter- molecular contributions.
+/** Structure used only in the pressure calculation to distinguish
+ *  non-bonded intra- and inter- molecular contributions.
  */
 class Observable_stat_non_bonded : public Observable_stat_base {
 public:
@@ -182,20 +180,10 @@ class Observable_stat_non_bonded : public Observable_stat_base {
   }
 };
 
-/** Invalidate observables.
- *  This function is called whenever the system has changed in such a way
- *  that the cached observables are no longer accurate or that the size of
- *  the cache is no longer suitable (e.g. after addition of a new actor or
- *  interaction).
- */
-void invalidate_obs();
-
 class Observable_stat_wrapper : public Observable_stat {
 public:
   /** Observed statistic for the current MPI rank. */
   Observable_stat local;
-  /** Flag to signal if the observable is initialized. */
-  bool is_initialized;
   /** Flag to signal if the observable measures instantaneous pressure, i.e.
    *  the pressure with velocity compensation (half a time step), instead of
    *  the conventional pressure. Only relevant for NpT simulations.
@@ -205,16 +193,10 @@ class Observable_stat_wrapper : public Observable_stat {
   explicit Observable_stat_wrapper(size_t chunk_size, bool pressure_obs = true)
       : Observable_stat{chunk_size, pressure_obs}, local{chunk_size,
                                                          pressure_obs},
-        is_initialized(false), v_comp(false) {
-    register_obs();
-  }
+        v_comp(false) {}
 
   /** Gather the contributions from all MPI ranks. */
   void reduce() { local.reduce(data.data()); }
-
-private:
-  /** Register this observable. */
-  void register_obs();
 };
 
 class Observable_stat_non_bonded_wrapper : public Observable_stat_non_bonded {

src/core/bonded_interactions/bonded_interactions.dox

@@ -31,7 +31,7 @@
  *   - write functions for calculating force and energy
  *   - add calls to the force and energy calculation functions to the force
  *     calculation in the integration loop as well as to energy and
- *     pressure/stress tensor analysis
+ *     pressure tensor analysis
  *  * Python interface:
  *   - import the definition of the interaction struct from the core
  *   - implement a class for the bonded interaction derived from the Python
@@ -146,7 +146,7 @@
  *        retval = fene_pair_energy(iaparams, dx);
  *        break;
  *      @endcode
- *  * Pressure, stress tensor and virial calculation: if your bonded
+ *  * Pressure tensor and virial calculation: if your bonded
  *    interaction is not a pair bond or modifies the particles involved,
  *    you have to implement a custom solution for virial calculation.
  *    The pressure occurs twice, once for the parallelized isotropic pressure

src/core/communication.hpp

@@ -220,7 +220,7 @@ void mpi_bcast_max_seen_particle_type(int s);
  *      \arg for \ref GatherStats::energy, calculate and reduce (sum up)
  *           energies, using \ref energy_calc.
  *      \arg for \ref GatherStats::pressure, calculate and reduce (sum up)
- *           pressure, stress tensor, using \ref pressure_calc.
+ *           pressure, using \ref pressure_calc.
  *      \arg for \ref GatherStats::pressure_v_comp, calculate and reduce
  *           (sum up) instantaneous pressure, using \ref pressure_calc.
  *      \arg for \ref GatherStats::lb_fluid_momentum, use

src/core/electrostatics_magnetostatics/coulomb.cpp

@@ -58,9 +58,9 @@ void calc_pressure_long_range(Observable_stat &virials,
     break;
   case COULOMB_P3M: {
     p3m_charge_assign(particles);
-    auto const p3m_stress = p3m_calc_kspace_stress();
-    std::copy_n(p3m_stress.data(), 9, p_tensor.coulomb.begin() + 9);
-    virials.coulomb[1] = p3m_stress[0] + p3m_stress[4] + p3m_stress[8];
+    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,7 +60,10 @@ struct Coulomb_parameters {
 extern Coulomb_parameters coulomb;
 
 namespace Coulomb {
-/** Number of electrostatic contributions to the system pressure calculation. */
+/** Number of electrostatic contributions to the system pressure calculation.
+ *  - slot 0: pressure from particle pairs
+ *  - slot 1: energies from electrostatics solvers
+ */
 inline size_t pressure_n() {
   switch (coulomb.method) {
   case COULOMB_NONE:
@@ -73,7 +76,11 @@ inline size_t pressure_n() {
   }
 }
 
-/** Number of electrostatic contributions to the system energy calculation. */
+/** Number of electrostatic contributions to the system energy calculation.
+ *  - slot 0: energies from particle pairs
+ *  - slot 1: energies from electrostatics solvers
+ *  - slot 2: energy corrections
+ */
 inline size_t energy_n() {
   switch (coulomb.method) {
   case COULOMB_NONE:

src/core/electrostatics_magnetostatics/dipole.hpp

@@ -71,7 +71,11 @@ namespace Dipole {
 /** Number of electrostatic contributions to the system pressure calculation. */
 constexpr size_t pressure_n() { return 0; }
 
-/** Number of electrostatic contributions to the system energy calculation. */
+/** Number of electrostatic contributions to the system energy calculation.
+ *  - slot 0: energies from particle pairs and magnetic field constraints
+ *  - slot 1: energies from magnetostatics solvers
+ *  - slot 2: energy corrections
+ */
 inline size_t energy_n() {
   switch (dipole.method) {
   case DIPOLAR_NONE:

src/core/electrostatics_magnetostatics/p3m.cpp

@@ -419,20 +419,22 @@ double dipole_correction_energy(Utils::Vector3d const &box_dipole) {
 }
 } // namespace
 
-/** @details Calculate the long range electrostatics part of the stress
+/** @details Calculate the long range electrostatics part of the pressure
  *  tensor. This is part \f$\Pi_{\textrm{dir}, \alpha, \beta}\f$ eq. (2.6)
  *  in @cite essmann95a. The part \f$\Pi_{\textrm{corr}, \alpha, \beta}\f$
  *  eq. (2.8) is not present here since M is the empty set in our simulations.
  */
-Utils::Vector9d p3m_calc_kspace_stress() {
-  Utils::Vector9d node_k_space_stress{};
+Utils::Vector9d p3m_calc_kspace_pressure_tensor() {
+  Utils::Vector9d node_k_space_pressure_tensor{};
 
   if (p3m.sum_q2 > 0) {
     p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim);
     fft_perform_forw(p3m.rs_mesh.data(), p3m.fft, comm_cart);
 
+    double diagonal = 0;
     int ind = 0;
     int j[3];
+    auto const half_alpha_inv_sq = Utils::sqr(1.0 / 2.0 / p3m.params.alpha);
     for (j[0] = 0; j[0] < p3m.fft.plan[3].new_mesh[RX]; j[0]++) {
       for (j[1] = 0; j[1] < p3m.fft.plan[3].new_mesh[RY]; j[1]++) {
         for (j[2] = 0; j[2] < p3m.fft.plan[3].new_mesh[RZ]; j[2]++) {
@@ -455,40 +457,29 @@ Utils::Vector9d p3m_calc_kspace_stress() {
           ind++;
 
           auto const vterm =
-              (sqk == 0)
-                  ? 0.
-                  : -2.0 * (1 / sqk + Utils::sqr(1.0 / 2.0 / p3m.params.alpha));
-
-          node_k_space_stress[0] +=
-              node_k_space_energy *
-              (1.0 + vterm * Utils::sqr(kx)); /* sigma_xx */
-          node_k_space_stress[1] +=
-              node_k_space_energy * (vterm * kx * ky); /* sigma_xy */
-          node_k_space_stress[2] +=
-              node_k_space_energy * (vterm * kx * kz); /* sigma_xz */
-
-          node_k_space_stress[3] +=
-              node_k_space_energy * (vterm * kx * ky); /* sigma_yx */
-          node_k_space_stress[4] +=
-              node_k_space_energy *
-              (1.0 + vterm * Utils::sqr(ky)); /* sigma_yy */
-          node_k_space_stress[5] +=
-              node_k_space_energy * (vterm * ky * kz); /* sigma_yz */
-
-          node_k_space_stress[6] +=
-              node_k_space_energy * (vterm * kx * kz); /* sigma_zx */
-          node_k_space_stress[7] +=
-              node_k_space_energy * (vterm * ky * kz); /* sigma_zy */
-          node_k_space_stress[8] +=
-              node_k_space_energy *
-              (1.0 + vterm * Utils::sqr(kz)); /* sigma_zz */
+              (sqk == 0) ? 0. : -2.0 * (1 / sqk + half_alpha_inv_sq);
+
+          diagonal += node_k_space_energy;
+          auto const prefactor = node_k_space_energy * vterm;
+          node_k_space_pressure_tensor[0] += prefactor * kx * kx; /* sigma_xx */
+          node_k_space_pressure_tensor[1] += prefactor * kx * ky; /* sigma_xy */
+          node_k_space_pressure_tensor[2] += prefactor * kx * kz; /* sigma_xz */
+          node_k_space_pressure_tensor[3] += prefactor * ky * kx; /* sigma_yx */
+          node_k_space_pressure_tensor[4] += prefactor * ky * ky; /* sigma_yy */
+          node_k_space_pressure_tensor[5] += prefactor * ky * kz; /* sigma_yz */
+          node_k_space_pressure_tensor[6] += prefactor * kz * kx; /* sigma_zx */
+          node_k_space_pressure_tensor[7] += prefactor * kz * ky; /* sigma_zy */
+          node_k_space_pressure_tensor[8] += prefactor * kz * kz; /* sigma_zz */
         }
       }
     }
+    node_k_space_pressure_tensor[0] += diagonal;
+    node_k_space_pressure_tensor[4] += diagonal;
+    node_k_space_pressure_tensor[8] += diagonal;
   }
 
   auto const force_prefac = coulomb.prefactor / (2.0 * box_geo.volume());
-  return force_prefac * node_k_space_stress;
+  return force_prefac * node_k_space_pressure_tensor;
 }
 
 double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,