Refactor Observable_stat framework

src/core/AtomDecomposition.hpp

@@ -38,7 +38,7 @@
  * by just evenly distributing the particles over
  * all part-taking nodes. Pairs are found by just
  * considering all pairs independent of logical or
- * physical location, it has therefor quadratic time
+ * physical location, it has therefore quadratic time
  * complexity in the number of particles.
  *
  * For a more detailed discussion please see @cite plimpton1995a.

src/core/CMakeLists.txt

@@ -31,6 +31,7 @@ set(EspressoCore_SRC
     reaction_ensemble.cpp
     rotate_system.cpp
     rotation.cpp
+    Observable_stat.cpp
     RuntimeErrorCollector.cpp
     RuntimeError.cpp
     RuntimeErrorStream.cpp

src/core/Observable_stat.cpp

@@ -0,0 +1,94 @@
+/*
+ * Copyright (C) 2010-2019 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 "Observable_stat.hpp"
+
+#include "config.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 =
+      m_pressure_obs ? Coulomb::pressure_n() : Coulomb::energy_n();
+  auto const n_dipolar =
+      m_pressure_obs ? Dipole::pressure_n() : Dipole::energy_n();
+  size_t n_vs = 0;
+#ifdef VIRTUAL_SITES
+  n_vs = 1;
+#endif
+  auto const n_bonded = bonded_ia_params.size();
+  auto const n_non_bonded = max_non_bonded_pairs();
+  constexpr size_t n_ext_fields = 1; // energies from all fields: accumulated
+  constexpr size_t n_kinetic = 1; // linear+angular kinetic energy: accumulated
+
+  // resize vector
+  auto const total = n_kinetic + n_bonded + n_non_bonded + n_coulomb +
+                     n_dipolar + n_vs + n_ext_fields;
+  data.resize(m_chunk_size * total);
+
+  // spans for the different contributions
+  kinetic = Utils::Span<double>(data.data(), m_chunk_size);
+  bonded = Utils::Span<double>(kinetic.end(), n_bonded * m_chunk_size);
+  non_bonded = Utils::Span<double>(bonded.end(), n_non_bonded * m_chunk_size);
+  coulomb = Utils::Span<double>(non_bonded.end(), n_coulomb * m_chunk_size);
+  dipolar = Utils::Span<double>(coulomb.end(), n_dipolar * m_chunk_size);
+  virtual_sites = Utils::Span<double>(dipolar.end(), n_vs * m_chunk_size);
+  external_fields =
+      Utils::Span<double>(virtual_sites.end(), n_ext_fields * m_chunk_size);
+}
+
+void Observable_stat_non_bonded::resize() {
+  // number of chunks for different interaction types
+  auto const n_non_bonded = max_non_bonded_pairs();
+  // resize vector
+  data.resize(m_chunk_size * 2 * n_non_bonded);
+  // spans for the different contributions
+  auto const span_size = n_non_bonded * m_chunk_size;
+  non_bonded_intra = Utils::Span<double>(data.data(), span_size);
+  non_bonded_inter = Utils::Span<double>(non_bonded_intra.end(), span_size);
+}
+
+void Observable_stat_base::reduce(double *out) const {
+  MPI_Reduce(data.data(), out, data.size(), MPI_DOUBLE, MPI_SUM, 0, comm_cart);
+}

src/core/Observable_stat.hpp

@@ -19,59 +19,214 @@
 #ifndef ESPRESSO_OBSERVABLE_STAT_HPP
 #define ESPRESSO_OBSERVABLE_STAT_HPP
 
-#include <vector>
+#include <boost/range/numeric.hpp>
 
-struct Observable_stat {
-  /** Status flag for observable calculation. For 'analyze energy': 0
-   *  re-initialize observable struct, else everything is fine,
-   *  calculation can start. For 'analyze pressure' and 'analyze p_inst':
-   *  0 or !(1+v_comp) re-initialize, else all OK.
-   */
-  int init_status;
+#include <utils/Span.hpp>
+#include <utils/index.hpp>
+
+#include <algorithm>
+#include <utility>
+#include <vector>
 
+/** Cache for system statistics.
+ *  Store unidimensional (energy, scalar pressure) and multi-dimensional
+ *  (stress tensors) properties of the system and provide accumulation and
+ *  reduction functionality.
+ */
+class Observable_stat_base {
+protected:
   /** Array for observables on each node. */
   std::vector<double> data;
+  /** Number of doubles per data item */
+  size_t m_chunk_size;
+
+public:
+  /** @param chunk_size Dimensionality of the data
+   */
+  explicit Observable_stat_base(size_t chunk_size) : m_chunk_size(chunk_size) {}
+
+  /** Resize the observable */
+  virtual void resize() = 0;
+
+  /** Resize and zero out the observable */
+  void resize_and_clear() {
+    resize();
+    std::fill(data.begin(), data.end(), 0);
+  }
+
+  /** Gather the contributions from the current MPI rank.
+   *  @param[out] out Destination of the reduction.
+   */
+  void reduce(double *out) const;
 
-  /** number of Coulomb interactions */
-  int n_coulomb;
-  /** number of dipolar interactions */
-  int n_dipolar;
-  /** Number of virtual sites relative (rigid body) contributions */
-  int n_virtual_sites;
-  /** Number of external field contributions */
-  const static int n_external_field = 1;
-
-  /** start of bonded interactions. Right after the special ones */
-  double *bonded;
-  /** start of observables for non-bonded interactions. */
-  double *non_bonded;
-  /** start of observables for Coulomb interaction. */
-  double *coulomb;
-  /** start of observables for Coulomb interaction. */
-  double *dipolar;
-  /** Start of observables for virtual sites relative (rigid bodies) */
-  double *virtual_sites;
-  /** Start of observables for external fields */
-  double *external_fields;
-
-  /** number of doubles per data item */
-  int chunk_size;
+  /** Accumulate values.
+   *  @param acc    Initial value for the accumulator.
+   *  @param column Which column to sum up (only relevant for multi-dimensional
+   *                observables).
+   */
+  double accumulate(double acc = 0.0, size_t column = 0) {
+    assert(column < m_chunk_size);
+    if (m_chunk_size == 1)
+      return boost::accumulate(data, acc);
+
+    for (auto it = data.begin() + column; it < data.end(); it += m_chunk_size)
+      acc += *it;
+    return acc;
+  }
+
+  /** Rescale values */
+  void rescale(double volume) {
+    auto const factor = 1. / volume;
+    for (auto &e : data) {
+      e *= factor;
+    }
+  }
+
+protected:
+  /** Get contribution from a non-bonded interaction */
+  Utils::Span<double> non_bonded_contribution(double *base_pointer, int type1,
+                                              int type2) const {
+    extern int max_seen_particle_type;
+    if (type1 > type2) {
+      using std::swap;
+      swap(type1, type2);
+    }
+    int offset = Utils::upper_triangular(type1, type2, max_seen_particle_type);
+    return Utils::Span<double>(base_pointer + offset * m_chunk_size,
+                               m_chunk_size);
+  }
+
+  /** Calculate the maximal number of non-bonded interaction pairs in the
+   *  system.
+   */
+  size_t max_non_bonded_pairs() {
+    extern int max_seen_particle_type;
+    return static_cast<size_t>(
+        (max_seen_particle_type * (max_seen_particle_type + 1)) / 2);
+  }
+};
+
+/** Structure used to cache the results of the scalar pressure, stress tensor
+ *  and energy calculations.
+ */
+class Observable_stat : public Observable_stat_base {
+private:
+  /** Whether this observable is a pressure or energy observable */
+  bool m_pressure_obs;
+
+public:
+  explicit Observable_stat(size_t chunk_size, bool pressure_obs = true)
+      : Observable_stat_base(chunk_size), m_pressure_obs(pressure_obs) {
+    resize_and_clear();
+  }
+
+  /** Contribution from linear and angular kinetic energy (accumulated). */
+  Utils::Span<double> kinetic;
+  /** Contribution(s) from bonded interactions. */
+  Utils::Span<double> bonded;
+  /** Contribution(s) from non-bonded interactions. */
+  Utils::Span<double> non_bonded;
+  /** Contribution(s) from Coulomb interactions. */
+  Utils::Span<double> coulomb;
+  /** Contribution(s) from dipolar interactions. */
+  Utils::Span<double> dipolar;
+  /** Contribution(s) from virtual sites (accumulated). */
+  Utils::Span<double> virtual_sites;
+  /** Contribution from external fields (accumulated). */
+  Utils::Span<double> external_fields;
+
+  /** Resize the observable */
+  void resize() final;
+
+  /** 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);
+  }
+
+  /** Get contribution from a non-bonded interaction */
+  Utils::Span<double> non_bonded_contribution(int type1, int type2) const {
+    return Observable_stat_base::non_bonded_contribution(non_bonded.data(),
+                                                         type1, type2);
+  }
 };
 
 /** Structure used only in the pressure and stress tensor calculation to
-   distinguish
-    non-bonded intra- and inter- molecular contributions. */
-struct Observable_stat_non_bonded {
-  /** Array for observables on each node. */
-  std::vector<double> data_nb;
+ *  distinguish non-bonded intra- and inter- molecular contributions.
+ */
+class Observable_stat_non_bonded : public Observable_stat_base {
+public:
+  explicit Observable_stat_non_bonded(size_t chunk_size)
+      : Observable_stat_base(chunk_size) {
+    resize_and_clear();
+  }
+
+  /** Contribution(s) from non-bonded intramolecular interactions. */
+  Utils::Span<double> non_bonded_intra;
+  /** Contribution(s) from non-bonded intermolecular interactions. */
+  Utils::Span<double> non_bonded_inter;
+
+  /** Resize the observable */
+  void resize() final;
+
+  /** Get contribution from a non-bonded intramolecular interaction */
+  Utils::Span<double> non_bonded_intra_contribution(int type1,
+                                                    int type2) const {
+    return non_bonded_contribution(non_bonded_intra.data(), type1, type2);
+  }
+
+  /** Get contribution from a non-bonded intermolecular interaction */
+  Utils::Span<double> non_bonded_inter_contribution(int type1,
+                                                    int type2) const {
+    return non_bonded_contribution(non_bonded_inter.data(), type1, type2);
+  }
+};
+
+/** 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.
+   */
+  bool v_comp;
+
+  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();
+  }
+
+  /** 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 {
+public:
+  /** Observed statistic for the current MPI rank. */
+  Observable_stat_non_bonded local;
 
-  /** start of observables for non-bonded intramolecular interactions. */
-  double *non_bonded_intra;
-  /** start of observables for non-bonded intermolecular interactions. */
-  double *non_bonded_inter;
+  explicit Observable_stat_non_bonded_wrapper(size_t chunk_size)
+      : Observable_stat_non_bonded{chunk_size}, local{chunk_size} {}
 
-  /** number of doubles per data item */
-  int chunk_size_nb;
+  /** Gather the contributions from all MPI ranks. */
+  void reduce() { local.reduce(data.data()); }
 };
 
 #endif // ESPRESSO_OBSERVABLE_STAT_HPP

src/core/actor/ActorList.hpp

@@ -28,4 +28,4 @@ class ActorList : public std::vector<Actor *> {
   void remove(Actor *actor);
 };
 
-#endif /* ACTORLIST_HPP_ */
+#endif /* _ACTOR_ACTORLIST_HPP */