diff --git a/src/core/AtomDecomposition.hpp b/src/core/AtomDecomposition.hpp
index 85ead10c6e9..795d57bb3bc 100644
--- a/src/core/AtomDecomposition.hpp
+++ b/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.
diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt
index 8a9ac0d1d66..186b0a5f28e 100644
--- a/src/core/CMakeLists.txt
+++ b/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
diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp
new file mode 100644
index 00000000000..ce4629fd0f6
--- /dev/null
+++ b/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 .
+ */
+
+#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
+
+extern boost::mpi::communicator comm_cart;
+
+/** Tracker of observables */
+std::vector ®istered_observables() {
+ static std::vector 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(data.data(), m_chunk_size);
+ bonded = Utils::Span(kinetic.end(), n_bonded * m_chunk_size);
+ non_bonded = Utils::Span(bonded.end(), n_non_bonded * m_chunk_size);
+ coulomb = Utils::Span(non_bonded.end(), n_coulomb * m_chunk_size);
+ dipolar = Utils::Span(coulomb.end(), n_dipolar * m_chunk_size);
+ virtual_sites = Utils::Span(dipolar.end(), n_vs * m_chunk_size);
+ external_fields =
+ Utils::Span(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(data.data(), span_size);
+ non_bonded_inter = Utils::Span(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);
+}
diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp
index 54c8de2f230..8bf031ced69 100644
--- a/src/core/Observable_stat.hpp
+++ b/src/core/Observable_stat.hpp
@@ -19,59 +19,214 @@
#ifndef ESPRESSO_OBSERVABLE_STAT_HPP
#define ESPRESSO_OBSERVABLE_STAT_HPP
-#include
+#include
-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
+#include
+
+#include
+#include
+#include
+/** 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 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 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(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(
+ (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 kinetic;
+ /** Contribution(s) from bonded interactions. */
+ Utils::Span bonded;
+ /** Contribution(s) from non-bonded interactions. */
+ Utils::Span non_bonded;
+ /** Contribution(s) from Coulomb interactions. */
+ Utils::Span coulomb;
+ /** Contribution(s) from dipolar interactions. */
+ Utils::Span dipolar;
+ /** Contribution(s) from virtual sites (accumulated). */
+ Utils::Span virtual_sites;
+ /** Contribution from external fields (accumulated). */
+ Utils::Span external_fields;
+
+ /** Resize the observable */
+ void resize() final;
+
+ /** Get contribution from a bonded interaction */
+ Utils::Span bonded_contribution(int bond_id) {
+ return Utils::Span(bonded.data() + m_chunk_size * bond_id,
+ m_chunk_size);
+ }
+
+ /** Get contribution from a non-bonded interaction */
+ Utils::Span 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 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 non_bonded_intra;
+ /** Contribution(s) from non-bonded intermolecular interactions. */
+ Utils::Span non_bonded_inter;
+
+ /** Resize the observable */
+ void resize() final;
+
+ /** Get contribution from a non-bonded intramolecular interaction */
+ Utils::Span 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 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
diff --git a/src/core/actor/ActorList.hpp b/src/core/actor/ActorList.hpp
index e8abcd2c948..3de187d9491 100644
--- a/src/core/actor/ActorList.hpp
+++ b/src/core/actor/ActorList.hpp
@@ -28,4 +28,4 @@ class ActorList : public std::vector {
void remove(Actor *actor);
};
-#endif /* ACTORLIST_HPP_ */
+#endif /* _ACTOR_ACTORLIST_HPP */
diff --git a/src/core/communication.cpp b/src/core/communication.cpp
index 4998ed283fc..df92012f2d4 100644
--- a/src/core/communication.cpp
+++ b/src/core/communication.cpp
@@ -57,7 +57,6 @@
#include "particle_data.hpp"
#include "pressure.hpp"
#include "rotation.hpp"
-#include "statistics.hpp"
#include "virtual_sites.hpp"
#include "electrostatics_magnetostatics/coulomb.hpp"
@@ -365,69 +364,52 @@ void mpi_bcast_max_seen_particle_type(int ns) {
}
/*************** GATHER ************/
-void mpi_gather_stats(int job, void *result, void *result_t, void *result_nb,
- void *result_t_nb) {
+void mpi_gather_stats(GatherStats job, double *result) {
+ auto job_slave = static_cast(job);
switch (job) {
- case 1:
- mpi_call(mpi_gather_stats_slave, -1, 1);
- energy_calc((double *)result, sim_time);
+ case GatherStats::energy:
+ mpi_call(mpi_gather_stats_slave, -1, job_slave);
+ energy_calc(sim_time);
break;
- case 2:
- /* calculate and reduce (sum up) virials for 'analyze pressure' or
- 'analyze stress_tensor' */
- mpi_call(mpi_gather_stats_slave, -1, 2);
- pressure_calc((double *)result, (double *)result_t, (double *)result_nb,
- (double *)result_t_nb, 0);
+ case GatherStats::pressure:
+ case GatherStats::pressure_v_comp:
+ mpi_call(mpi_gather_stats_slave, -1, job_slave);
+ pressure_calc(job == GatherStats::pressure_v_comp);
break;
- case 3:
- mpi_call(mpi_gather_stats_slave, -1, 3);
- pressure_calc((double *)result, (double *)result_t, (double *)result_nb,
- (double *)result_t_nb, 1);
- break;
- case 6:
- mpi_call(mpi_gather_stats_slave, -1, 6);
- lb_calc_fluid_momentum((double *)result, lbpar, lbfields, lblattice);
- break;
- case 7:
+ case GatherStats::lb_fluid_momentum:
+ mpi_call(mpi_gather_stats_slave, -1, job_slave);
+ lb_calc_fluid_momentum(result, lbpar, lbfields, lblattice);
break;
#ifdef LB_BOUNDARIES
- case 8:
- mpi_call(mpi_gather_stats_slave, -1, 8);
- lb_collect_boundary_forces((double *)result);
+ case GatherStats::lb_boundary_forces:
+ mpi_call(mpi_gather_stats_slave, -1, job_slave);
+ lb_collect_boundary_forces(result);
break;
#endif
default:
fprintf(
stderr,
"%d: INTERNAL ERROR: illegal request %d for mpi_gather_stats_slave\n",
- this_node, job);
+ this_node, job_slave);
errexit();
}
}
-void mpi_gather_stats_slave(int, int job) {
+void mpi_gather_stats_slave(int, int job_slave) {
+ auto job = static_cast(job_slave);
switch (job) {
- case 1:
- /* calculate and reduce (sum up) energies */
- energy_calc(nullptr, sim_time);
+ case GatherStats::energy:
+ energy_calc(sim_time);
break;
- case 2:
- /* calculate and reduce (sum up) virials for 'analyze pressure' or 'analyze
- * stress_tensor'*/
- pressure_calc(nullptr, nullptr, nullptr, nullptr, 0);
+ case GatherStats::pressure:
+ case GatherStats::pressure_v_comp:
+ pressure_calc(job == GatherStats::pressure_v_comp);
break;
- case 3:
- /* calculate and reduce (sum up) virials, revert velocities half a timestep
- * for 'analyze p_inst' */
- pressure_calc(nullptr, nullptr, nullptr, nullptr, 1);
- break;
- case 6:
+ case GatherStats::lb_fluid_momentum:
lb_calc_fluid_momentum(nullptr, lbpar, lbfields, lblattice);
break;
- case 7:
- break;
#ifdef LB_BOUNDARIES
- case 8:
+ case GatherStats::lb_boundary_forces:
lb_collect_boundary_forces(nullptr);
break;
#endif
@@ -435,7 +417,7 @@ void mpi_gather_stats_slave(int, int job) {
fprintf(
stderr,
"%d: INTERNAL ERROR: illegal request %d for mpi_gather_stats_slave\n",
- this_node, job);
+ this_node, job_slave);
errexit();
}
}
diff --git a/src/core/communication.hpp b/src/core/communication.hpp
index 79c15390888..0e8a93ce31f 100644
--- a/src/core/communication.hpp
+++ b/src/core/communication.hpp
@@ -65,6 +65,14 @@ extern int this_node;
extern int n_nodes;
/** The communicator */
extern boost::mpi::communicator comm_cart;
+/** Statistics to calculate */
+enum class GatherStats : int {
+ energy,
+ pressure,
+ pressure_v_comp,
+ lb_fluid_momentum,
+ lb_boundary_forces
+};
/**
* Default MPI tag used by callbacks.
@@ -208,47 +216,22 @@ void mpi_bcast_ia_params(int i, int j);
void mpi_bcast_max_seen_particle_type(int s);
/** Gather data for analysis.
- * \todo update parameter descriptions
- * \param job what to do:
- * \arg \c 1 calculate and reduce (sum up) energies,
- * using \ref energy_calc.
- * \arg \c 2 calculate and reduce (sum up) pressure, stress tensor,
- * using \ref pressure_calc.
- * \arg \c 3 calculate and reduce (sum up) instantaneous pressure,
- * using \ref pressure_calc.
- * \arg \c 6 use \ref lb_calc_fluid_momentum
- * \arg \c 8 use \ref lb_collect_boundary_forces
- * \param result where to store the gathered value(s):
- * \arg for \c job=1 unused (the results are stored in a global
- * energy array of type \ref Observable_stat)
- * \arg for \c job=2 unused (the results are stored in a global
- * virials array of type \ref Observable_stat)
- * \arg for \c job=3 unused (the results are stored in a global
- * virials array of type \ref Observable_stat)
- * \param result_t where to store the gathered value(s):
- * \arg for \c job=1 unused (the results are stored in a global
- * energy array of type \ref Observable_stat)
- * \arg for \c job=2 unused (the results are stored in a global
- * p_tensor tensor of type \ref Observable_stat)
- * \arg for \c job=3 unused (the results are stored in a global
- * p_tensor tensor of type \ref Observable_stat)
- * \param result_nb where to store the gathered value(s):
- * \arg for \c job=1 unused (the results are stored in a global
- * energy array of type \ref Observable_stat_non_bonded)
- * \arg for \c job=2 unused (the results are stored in a global
- * virials_non_bonded array of type \ref Observable_stat_non_bonded)
- * \arg for \c job=3 unused (the results are stored in a global
- * virials_non_bonded array of type \ref Observable_stat_non_bonded)
- * \param result_t_nb where to store the gathered value(s):
- * \arg for \c job=1 unused (the results are stored in a global
- * energy array of type \ref Observable_stat_non_bonded)
- * \arg for \c job=2 unused (the results are stored in a global
- * p_tensor_non_bonded tensor of type \ref Observable_stat_non_bonded)
- * \arg for \c job=3 unused (the results are stored in a global
- * p_tensor_non_bonded tensor of type \ref Observable_stat_non_bonded)
+ * \param[in] job what to do:
+ * \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.
+ * \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
+ * \ref lb_calc_fluid_momentum.
+ * \arg for \ref GatherStats::lb_boundary_forces, use
+ * \ref lb_collect_boundary_forces.
+ * \param[out] result where to store values gathered by
+ * \ref GatherStats::lb_fluid_momentum,
+ * \ref GatherStats::lb_boundary_forces
*/
-void mpi_gather_stats(int job, void *result, void *result_t, void *result_nb,
- void *result_t_nb);
+void mpi_gather_stats(GatherStats job, double *result = nullptr);
/** Send new \ref time_step and rescale the velocities accordingly. */
void mpi_set_time_step(double time_step);
diff --git a/src/core/constraints/Constraint.hpp b/src/core/constraints/Constraint.hpp
index e363bb3fa2b..678bae90a2e 100644
--- a/src/core/constraints/Constraint.hpp
+++ b/src/core/constraints/Constraint.hpp
@@ -19,6 +19,7 @@
#ifndef CONSTRAINTS_CONSTRAINT_HPP
#define CONSTRAINTS_CONSTRAINT_HPP
+#include "Observable_stat.hpp"
#include "Particle.hpp"
#include "energy.hpp"
@@ -26,14 +27,12 @@ namespace Constraints {
class Constraint {
public:
/**
- * @brief Add energy contribution of this constraints to energy.
- *
- * Add constraint energy for particle to observable.
+ * @brief Add constraint energy for particle to the energy observable.
*
* @param[in] p The particle to add the energy for.
* @param[in] folded_pos Folded position of the particle.
* @param[in] t The time at which the energy should be calculated.
- * @param[out] energy to add the energy to.
+ * @param[out] energy Energy observable to add this constraint's energy to.
*/
virtual void add_energy(const Particle &p, const Utils::Vector3d &folded_pos,
double t, Observable_stat &energy) const = 0;
@@ -41,8 +40,6 @@ class Constraint {
/**
* @brief Calculate the force of the constraint on a particle.
*
- * Add constraint energy for particle to observable.
- *
* @param[in] p The particle to calculate the force for.
* @param[in] folded_pos Folded position of the particle.
* @param[in] t The time at which the force should be calculated.
diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp
index a330094b611..72fe5bac3e7 100644
--- a/src/core/constraints/ShapeBasedConstraint.cpp
+++ b/src/core/constraints/ShapeBasedConstraint.cpp
@@ -136,6 +136,7 @@ void ShapeBasedConstraint::add_energy(const Particle &p,
}
}
if (part_rep.p.type >= 0)
- *obsstat_nonbonded(&energy, p.p.type, part_rep.p.type) += nonbonded_en;
+ energy.non_bonded_contribution(p.p.type, part_rep.p.type)[0] +=
+ nonbonded_en;
}
} // namespace Constraints
diff --git a/src/core/constraints/ShapeBasedConstraint.hpp b/src/core/constraints/ShapeBasedConstraint.hpp
index 2517e19daaf..9d9c246316a 100644
--- a/src/core/constraints/ShapeBasedConstraint.hpp
+++ b/src/core/constraints/ShapeBasedConstraint.hpp
@@ -23,6 +23,7 @@
#include "Constraint.hpp"
#include "Particle.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include
#include
diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu
index 770b288681f..f27d938a90b 100644
--- a/src/core/cuda_common_cuda.cu
+++ b/src/core/cuda_common_cuda.cu
@@ -22,6 +22,7 @@
#include "debug.hpp"
+#include "Observable_stat.hpp"
#include "ParticleRange.hpp"
#include "cuda_init.hpp"
#include "cuda_interface.hpp"
@@ -221,15 +222,21 @@ void clear_energy_on_GPU() {
}
void copy_energy_from_GPU() {
+ extern Observable_stat_wrapper obs_energy;
if (!global_part_vars_host.communication_enabled)
return;
cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy),
cudaMemcpyDeviceToHost));
- copy_CUDA_energy_to_energy(energy_host);
+ if (!obs_energy.local.bonded.empty())
+ obs_energy.local.bonded[0] += energy_host.bonded;
+ if (!obs_energy.local.non_bonded.empty())
+ obs_energy.local.non_bonded[0] += energy_host.non_bonded;
+ if (!obs_energy.local.coulomb.empty())
+ obs_energy.local.coulomb[0] += energy_host.coulomb;
+ if (obs_energy.local.dipolar.size() >= 2)
+ obs_energy.local.dipolar[1] += energy_host.dipolar;
}
-/** @name Generic copy functions from and to device */
-
void _cuda_safe_mem(cudaError_t CU_err, const char *file, unsigned int line) {
if (cudaSuccess != CU_err) {
fprintf(stderr, "Cuda Memory error at %s:%u.\n", file, line);
@@ -254,5 +261,3 @@ void _cuda_safe_mem(cudaError_t CU_err, const char *file, unsigned int line) {
}
}
}
-
-/*@}*/
diff --git a/src/core/cuda_interface.cpp b/src/core/cuda_interface.cpp
index d068684a918..909a3da4c0a 100644
--- a/src/core/cuda_interface.cpp
+++ b/src/core/cuda_interface.cpp
@@ -23,7 +23,6 @@
#include "communication.hpp"
#include "config.hpp"
-#include "energy.hpp"
#include "grid.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include "serialization/CUDA_particle_data.hpp"
@@ -161,19 +160,6 @@ void cuda_mpi_send_forces(const ParticleRange &particles,
}
}
-/** Takes a CUDA_energy struct and adds it to the core energy struct.
-This cannot be done from inside cuda_common_cuda.cu:copy_energy_from_GPU()
-because energy.hpp indirectly includes on mpi.h while .cu files may not depend
-on mpi.h. */
-void copy_CUDA_energy_to_energy(CUDA_energy energy_host) {
- energy.bonded[0] += energy_host.bonded;
- energy.non_bonded[0] += energy_host.non_bonded;
- if (energy.n_coulomb >= 1)
- energy.coulomb[0] += energy_host.coulomb;
- if (energy.n_dipolar >= 2)
- energy.dipolar[1] += energy_host.dipolar;
-}
-
void cuda_bcast_global_part_params() { mpi_bcast_cuda_global_part_vars(); }
#endif /* ifdef CUDA */
diff --git a/src/core/cuda_interface.hpp b/src/core/cuda_interface.hpp
index 76887922869..994e7fc0e3d 100644
--- a/src/core/cuda_interface.hpp
+++ b/src/core/cuda_interface.hpp
@@ -108,7 +108,6 @@ typedef struct {
void copy_forces_from_GPU(ParticleRange &particles);
void copy_energy_from_GPU();
-void copy_CUDA_energy_to_energy(CUDA_energy energy_host);
void clear_energy_on_GPU();
CUDA_global_part_vars *gpu_get_global_particle_vars_pointer_host();
diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp
index 8f2e7210741..f7b916bc7aa 100644
--- a/src/core/electrostatics_magnetostatics/coulomb.cpp
+++ b/src/core/electrostatics_magnetostatics/coulomb.cpp
@@ -19,9 +19,6 @@
#include "electrostatics_magnetostatics/coulomb.hpp"
-// Real space cutoff
-double coulomb_cutoff;
-
#ifdef ELECTROSTATICS
#include "cells.hpp"
#include "communication.hpp"
@@ -45,21 +42,6 @@ double coulomb_cutoff;
Coulomb_parameters coulomb;
namespace Coulomb {
-
-void pressure_n(int &n_coulomb) {
- switch (coulomb.method) {
- case COULOMB_NONE:
- n_coulomb = 0;
- break;
- case COULOMB_P3M_GPU:
- case COULOMB_P3M:
- n_coulomb = 2;
- break;
- default:
- n_coulomb = 1;
- }
-}
-
void calc_pressure_long_range(Observable_stat &virials,
Observable_stat &p_tensor,
const ParticleRange &particles) {
@@ -77,7 +59,7 @@ void calc_pressure_long_range(Observable_stat &virials,
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 + 9);
+ std::copy_n(p3m_stress.data(), 9, p_tensor.coulomb.begin() + 9);
virials.coulomb[1] = p3m_stress[0] + p3m_stress[4] + p3m_stress[8];
break;
@@ -340,7 +322,7 @@ void calc_energy_long_range(Observable_stat &energy,
energy.coulomb[1] +=
0.5 * ELC_P3M_dielectric_layers_energy_self(particles);
- // assign both original and image charges now
+ // assign both original and image charges now
ELC_p3m_charge_assign_both(particles);
ELC_P3M_modify_p3m_sums_both(particles);
@@ -369,21 +351,6 @@ void calc_energy_long_range(Observable_stat &energy,
}
}
-int energy_n() {
- switch (coulomb.method) {
- case COULOMB_NONE:
- return 0;
- case COULOMB_ELC_P3M:
- return 3;
- case COULOMB_P3M_GPU:
- case COULOMB_P3M:
- case COULOMB_SCAFACOS:
- return 2;
- default:
- return 1;
- }
-}
-
int iccp3m_sanity_check() {
switch (coulomb.method) {
#ifdef P3M
@@ -426,10 +393,7 @@ int elc_sanity_check() {
return ES_ERROR;
}
case COULOMB_ELC_P3M:
-
case COULOMB_P3M:
- p3m.params.epsilon = P3M_EPSILON_METALLIC;
- coulomb.method = COULOMB_ELC_P3M;
return ES_OK;
default:
break;
@@ -480,8 +444,6 @@ int set_prefactor(double prefactor) {
return ES_OK;
}
-/** @brief Deactivates the current Coulomb method
- */
void deactivate_method() {
coulomb.prefactor = 0;
diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp
index c99d7bf0c3e..71884d0a96a 100644
--- a/src/core/electrostatics_magnetostatics/coulomb.hpp
+++ b/src/core/electrostatics_magnetostatics/coulomb.hpp
@@ -21,21 +21,19 @@
#include "config.hpp"
+#include
+
#ifdef ELECTROSTATICS
#include "Observable_stat.hpp"
#include
#include
-/** \name Type codes for the type of Coulomb interaction
- Enumeration of implemented methods for the electrostatic
- interaction.
-*/
-/************************************************************/
-/*@{*/
-
+/** Type codes for the type of %Coulomb interaction.
+ * Enumeration of implemented methods for the electrostatic interaction.
+ */
enum CoulombMethod {
- COULOMB_NONE, ///< %Coulomb interaction switched off (NONE)
+ COULOMB_NONE, ///< %Coulomb interaction switched off
COULOMB_DH, ///< %Coulomb method is Debye-Hueckel
COULOMB_P3M, ///< %Coulomb method is P3M
COULOMB_P3M_GPU, ///< %Coulomb method is P3M with GPU-based long-range part
@@ -43,31 +41,54 @@ enum CoulombMethod {
COULOMB_MMM1D, ///< %Coulomb method is one-dimensional MMM
COULOMB_RF, ///< %Coulomb method is Reaction-Field
COULOMB_MMM1D_GPU, ///< %Coulomb method is one-dimensional MMM running on GPU
- COULOMB_SCAFACOS, ///< %Coulomb method is scafacos
+ COULOMB_SCAFACOS, ///< %Coulomb method is ScaFaCoS
};
-/*@}*/
-
-/** \name Compounds for Coulomb interactions */
-/*@{*/
-/** field containing the interaction parameters for
- * the Coulomb interaction. */
+/** Interaction parameters for the %Coulomb interaction. */
struct Coulomb_parameters {
- /** bjerrum length times temperature. */
+ /** Bjerrum length times temperature. */
double prefactor = 0.;
double field_induced = 0.;
double field_applied = 0.;
- /** Method to treat Coulomb interaction. */
+ /** Method to treat %Coulomb interaction. */
CoulombMethod method = COULOMB_NONE;
};
-/** Structure containing the Coulomb parameters. */
+/** Structure containing the %Coulomb parameters. */
extern Coulomb_parameters coulomb;
namespace Coulomb {
-void pressure_n(int &n_coulomb);
+/** Number of electrostatic contributions to the system pressure calculation. */
+inline size_t pressure_n() {
+ switch (coulomb.method) {
+ case COULOMB_NONE:
+ return 0;
+ case COULOMB_P3M_GPU:
+ case COULOMB_P3M:
+ return 2;
+ default:
+ return 1;
+ }
+}
+
+/** Number of electrostatic contributions to the system energy calculation. */
+inline size_t energy_n() {
+ switch (coulomb.method) {
+ case COULOMB_NONE:
+ return 0;
+ case COULOMB_ELC_P3M:
+ return 3;
+ case COULOMB_P3M_GPU:
+ case COULOMB_P3M:
+ case COULOMB_SCAFACOS:
+ return 2;
+ default:
+ return 1;
+ }
+}
+
void calc_pressure_long_range(Observable_stat &virials,
Observable_stat &p_tensor,
const ParticleRange &particles);
@@ -85,7 +106,6 @@ void calc_long_range_force(const ParticleRange &particles);
void calc_energy_long_range(Observable_stat &energy,
const ParticleRange &particles);
-int energy_n();
int iccp3m_sanity_check();
@@ -96,15 +116,18 @@ void bcast_coulomb_params();
/** @brief Set the electrostatics prefactor */
int set_prefactor(double prefactor);
-/** @brief Deactivates the current Coulomb method
- This was part of coulomb_set_bjerrum()
-*/
+/** @brief Deactivates the current %Coulomb method. */
void deactivate_method();
-/** @brief Update particles with properties depending on other particles
- * e.g., charges of ICC particles
+/** @brief Update particles with properties depending on other particles,
+ * namely ICC charges.
*/
void update_dependent_particles();
} // namespace Coulomb
+#else // ELECTROSTATICS
+namespace Coulomb {
+constexpr size_t pressure_n() { return 0; }
+constexpr size_t energy_n() { return 0; }
+} // namespace Coulomb
#endif // ELECTROSTATICS
#endif // ESPRESSO_COULOMB_HPP
diff --git a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
index 4747685661a..ac3d579a37c 100644
--- a/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
+++ b/src/core/electrostatics_magnetostatics/coulomb_inline.hpp
@@ -102,10 +102,10 @@ pair_force(Particle const &p1, Particle const &p2, Utils::Vector3d const &d,
* If supported by the method, this returns the virial
* contribution to the pressure tensor for this pair.
*
- * @param p1 Particle
- * @param p2 Particle
- * @param d Distance
- * @param dist |d|
+ * @param p1 particle
+ * @param p2 particle
+ * @param d distance vector
+ * @param dist distance norm
* @return Contribution to the pressure tensor.
*/
inline Utils::Vector pair_pressure(Particle const &p1,
diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp
index 828db98d2ab..42a5858c9db 100644
--- a/src/core/electrostatics_magnetostatics/dipole.cpp
+++ b/src/core/electrostatics_magnetostatics/dipole.cpp
@@ -19,9 +19,6 @@
#include "electrostatics_magnetostatics/dipole.hpp"
-// Real space cutoff of long range methods
-double dipolar_cutoff;
-
#ifdef DIPOLES
#include "actor/DipolarBarnesHut.hpp"
@@ -48,8 +45,6 @@ Dipole_parameters dipole = {
};
namespace Dipole {
-int pressure_n() { return 0; }
-
void calc_pressure_long_range() {
switch (dipole.method) {
case DIPOLAR_NONE:
@@ -260,30 +255,6 @@ void calc_energy_long_range(Observable_stat &energy,
}
}
-void energy_n(int &n_dipolar) {
- switch (dipole.method) {
- case DIPOLAR_NONE:
- n_dipolar = 1; // because there may be an external magnetic field
- break;
- case DIPOLAR_P3M:
- case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA:
- case DIPOLAR_DS:
- case DIPOLAR_DS_GPU:
-#ifdef DIPOLAR_BARNES_HUT
- case DIPOLAR_BH_GPU:
-#endif
- case DIPOLAR_SCAFACOS:
- n_dipolar = 2;
- break;
- case DIPOLAR_MDLC_P3M:
- case DIPOLAR_MDLC_DS:
- n_dipolar = 3;
- break;
- default:
- break;
- }
-}
-
int set_mesh() {
switch (dipole.method) {
#ifdef DP3M
@@ -340,7 +311,7 @@ void set_method_local(DipolarInteraction method) {
if ((dipole.method == DIPOLAR_BH_GPU) && (method != DIPOLAR_BH_GPU)) {
deactivate_dipolar_barnes_hut();
}
-#endif // BARNES_HUT
+#endif // DIPOLAR_ BARNES_HUT
dipole.method = method;
}
diff --git a/src/core/electrostatics_magnetostatics/dipole.hpp b/src/core/electrostatics_magnetostatics/dipole.hpp
index 15db5234b89..db477d24c65 100644
--- a/src/core/electrostatics_magnetostatics/dipole.hpp
+++ b/src/core/electrostatics_magnetostatics/dipole.hpp
@@ -21,7 +21,7 @@
#include "config.hpp"
-extern double dipolar_cutoff;
+#include
#ifdef DIPOLES
#include "Observable_stat.hpp"
@@ -31,52 +31,68 @@ extern double dipolar_cutoff;
#include
#include
-/** \name Compounds for Dipole interactions */
-/*@{*/
-
-/** \name Type codes for the type of dipolar interaction
- Enumeration of implemented methods for the magnetostatic
- interaction.
+/** Type codes for the type of dipolar interaction.
+ * Enumeration of implemented methods for the magnetostatic interaction.
*/
-/************************************************************/
-/*@{*/
enum DipolarInteraction {
- /** dipolar interaction switched off (NONE). */
+ /** Dipolar interaction switched off. */
DIPOLAR_NONE = 0,
- /** dipolar method is P3M. */
+ /** Dipolar method is P3M. */
DIPOLAR_P3M,
/** Dipolar method is P3M plus DLC. */
DIPOLAR_MDLC_P3M,
- /** Dipolar method is all with all and no replicas */
+ /** Dipolar method is all with all and no replicas. */
DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA,
- /** Dipolar method is magnetic dipolar direct sum */
+ /** Dipolar method is magnetic dipolar direct summation. */
DIPOLAR_DS,
- /** Dipolar method is direct sum plus DLC. */
+ /** Dipolar method is direct summation plus DLC. */
DIPOLAR_MDLC_DS,
- /** Direct summation on gpu */
+ /** Dipolar method is direct summation on GPU. */
DIPOLAR_DS_GPU,
#ifdef DIPOLAR_BARNES_HUT
- /** Direct summation on gpu by Barnes-Hut algorithm */
+ /** Dipolar method is direct summation on GPU by Barnes-Hut algorithm. */
DIPOLAR_BH_GPU,
#endif
- /** Scafacos library */
+ /** Dipolar method is ScaFaCoS. */
DIPOLAR_SCAFACOS
};
-/** field containing the interaction parameters for
- * the Dipole interaction. */
+/** Interaction parameters for the %dipole interaction. */
struct Dipole_parameters {
double prefactor;
DipolarInteraction method;
};
-/*@}*/
-/** Structure containing the Dipole parameters. */
+/** Structure containing the %dipole parameters. */
extern Dipole_parameters dipole;
namespace Dipole {
-int pressure_n();
+/** 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. */
+inline size_t energy_n() {
+ switch (dipole.method) {
+ case DIPOLAR_NONE:
+ return 1; // because there may be an external magnetic field
+ case DIPOLAR_P3M:
+ case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA:
+ case DIPOLAR_DS:
+ case DIPOLAR_DS_GPU:
+#ifdef DIPOLAR_BARNES_HUT
+ case DIPOLAR_BH_GPU:
+#endif
+ case DIPOLAR_SCAFACOS:
+ return 2;
+ case DIPOLAR_MDLC_P3M:
+ case DIPOLAR_MDLC_DS:
+ return 3;
+ default:
+ return 0;
+ }
+}
+
void calc_pressure_long_range();
void nonbonded_sanity_check(int &state);
@@ -91,7 +107,6 @@ void calc_long_range_force(const ParticleRange &particles);
void calc_energy_long_range(Observable_stat &energy,
const ParticleRange &particles);
-void energy_n(int &n_dipolar);
int set_mesh();
void bcast_params(const boost::mpi::communicator &comm);
@@ -100,8 +115,11 @@ void bcast_params(const boost::mpi::communicator &comm);
int set_Dprefactor(double prefactor);
void set_method_local(DipolarInteraction method);
-
} // namespace Dipole
-
+#else // DIPOLES
+namespace Dipole {
+constexpr size_t pressure_n() { return 0; }
+constexpr size_t energy_n() { return 0; }
+} // namespace Dipole
#endif // DIPOLES
#endif // ESPRESSO_DIPOLE_HPP
diff --git a/src/core/electrostatics_magnetostatics/dipole_inline.hpp b/src/core/electrostatics_magnetostatics/dipole_inline.hpp
index 7941c66d9aa..344bd739910 100644
--- a/src/core/electrostatics_magnetostatics/dipole_inline.hpp
+++ b/src/core/electrostatics_magnetostatics/dipole_inline.hpp
@@ -56,7 +56,6 @@ inline double pair_energy(Particle const &p1, Particle const &p2,
Utils::Vector3d const &d, double dist, double dist2) {
double ret = 0;
if (dipole.method != DIPOLAR_NONE) {
- // ret=0;
switch (dipole.method) {
#ifdef DP3M
case DIPOLAR_MDLC_P3M:
diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp
index 3f1c326ce37..c3451d335c7 100644
--- a/src/core/electrostatics_magnetostatics/elc.cpp
+++ b/src/core/electrostatics_magnetostatics/elc.cpp
@@ -1229,6 +1229,9 @@ int ELC_set_params(double maxPWerror, double gap_size, double far_cut,
Coulomb::elc_sanity_check();
+ p3m.params.epsilon = P3M_EPSILON_METALLIC;
+ coulomb.method = COULOMB_ELC_P3M;
+
elc_params.far_cut = far_cut;
if (far_cut != -1) {
elc_params.far_cut2 = Utils::sqr(far_cut);
diff --git a/src/core/electrostatics_magnetostatics/icc.cpp b/src/core/electrostatics_magnetostatics/icc.cpp
index 894165019da..ec77dbbd577 100644
--- a/src/core/electrostatics_magnetostatics/icc.cpp
+++ b/src/core/electrostatics_magnetostatics/icc.cpp
@@ -65,7 +65,7 @@ void init_forces_iccp3m(const ParticleRange &particles,
void force_calc_iccp3m(const ParticleRange &particles,
const ParticleRange &ghost_particles);
-/** Variant of @ref add_non_bonded_pair_force where only Coulomb
+/** Variant of @ref add_non_bonded_pair_force where only %Coulomb
* contributions are calculated
*/
inline void add_non_bonded_pair_force_iccp3m(Particle &p1, Particle &p2,
@@ -97,7 +97,7 @@ int iccp3m_iteration(const ParticleRange &particles,
Coulomb::iccp3m_sanity_check();
- if ((iccp3m_cfg.eout <= 0)) {
+ if (iccp3m_cfg.eout <= 0) {
runtimeErrorMsg()
<< "ICCP3M: nonpositive dielectric constant is not allowed.";
}
@@ -160,7 +160,7 @@ int iccp3m_iteration(const ParticleRange &particles,
/* check if the charge now is more than 1e6, to determine if ICC still
* leads to reasonable results */
/* this is kind of an arbitrary measure but does a good job spotting
- * divergence !*/
+ * divergence! */
if (std::abs(p.p.q) > 1e6) {
runtimeErrorMsg()
<< "too big charge assignment in iccp3m! q >1e6 , assigned "
@@ -199,7 +199,7 @@ void force_calc_iccp3m(const ParticleRange &particles,
short_range_loop(Utils::NoOp{}, [](Particle &p1, Particle &p2,
Distance const &d) {
- /* calc non bonded interactions */
+ /* calc non-bonded interactions */
add_non_bonded_pair_force_iccp3m(p1, p2, d.vec21, sqrt(d.dist2), d.dist2);
});
diff --git a/src/core/energy.cpp b/src/core/energy.cpp
index 14a80957c6d..335d433e2eb 100644
--- a/src/core/energy.cpp
+++ b/src/core/energy.cpp
@@ -23,6 +23,7 @@
*/
#include "EspressoSystemInterface.hpp"
+#include "Observable_stat.hpp"
#include "communication.hpp"
#include "constraints.hpp"
#include "cuda_interface.hpp"
@@ -30,7 +31,6 @@
#include "energy_inline.hpp"
#include "event.hpp"
#include "forces.hpp"
-#include
#include "short_range_loop.hpp"
@@ -39,46 +39,20 @@
ActorList energyActors;
-Observable_stat energy{};
-Observable_stat total_energy{};
-
-/************************************************************/
-
-void init_energies(Observable_stat *stat) {
- int n_pre, n_non_bonded, n_dipolar(0);
-
- n_pre = 1;
- n_non_bonded = (max_seen_particle_type * (max_seen_particle_type + 1)) / 2;
-
-#ifdef ELECTROSTATICS
- auto const n_coulomb = Coulomb::energy_n();
-#else
- int const n_coulomb = 0;
-#endif
-#ifdef DIPOLES
- Dipole::energy_n(n_dipolar);
-#endif
-
- obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded,
- n_coulomb, n_dipolar, 0, 1);
- stat->init_status = 0;
-}
-
-/************************************************************/
+/** Energy of the system */
+Observable_stat_wrapper obs_energy{1, false};
+/** Reduce the system energy from all MPI ranks. */
void master_energy_calc() {
- mpi_gather_stats(1, total_energy.data.data(), nullptr, nullptr, nullptr);
-
- total_energy.init_status = 1;
+ mpi_gather_stats(GatherStats::energy);
+ obs_energy.is_initialized = true;
}
-/************************************************************/
-
-void energy_calc(double *result, const double time) {
+void energy_calc(const double time) {
if (!interactions_sanity_checks())
return;
- init_energies(&energy);
+ obs_energy.local.resize_and_clear();
#ifdef CUDA
clear_energy_on_GPU();
@@ -105,36 +79,35 @@ void energy_calc(double *result, const double time) {
calc_long_range_energies(cell_structure.local_particles());
auto local_parts = cell_structure.local_particles();
- Constraints::constraints.add_energy(local_parts, time, energy);
+ Constraints::constraints.add_energy(local_parts, time, obs_energy.local);
#ifdef CUDA
copy_energy_from_GPU();
#endif
/* gather data */
- MPI_Reduce(energy.data.data(), result, energy.data.size(), MPI_DOUBLE,
- MPI_SUM, 0, comm_cart);
+ obs_energy.reduce();
}
-/************************************************************/
+void update_energy() {
+ if (!obs_energy.is_initialized) {
+ obs_energy.resize();
+ master_energy_calc();
+ }
+}
void calc_long_range_energies(const ParticleRange &particles) {
#ifdef ELECTROSTATICS
/* calculate k-space part of electrostatic interaction. */
- Coulomb::calc_energy_long_range(energy, particles);
+ Coulomb::calc_energy_long_range(obs_energy.local, particles);
#endif /* ifdef ELECTROSTATICS */
#ifdef DIPOLES
- Dipole::calc_energy_long_range(energy, particles);
+ Dipole::calc_energy_long_range(obs_energy.local, particles);
#endif /* ifdef DIPOLES */
}
double calculate_current_potential_energy_of_system() {
- // calculate potential energy
- if (total_energy.init_status == 0) {
- init_energies(&total_energy);
- master_energy_calc();
- }
-
- return boost::accumulate(total_energy.data, -total_energy.data[0]);
+ update_energy();
+ return obs_energy.accumulate(-obs_energy.kinetic[0]);
}
diff --git a/src/core/energy.hpp b/src/core/energy.hpp
index d11a0458fa7..73fcb102df9 100644
--- a/src/core/energy.hpp
+++ b/src/core/energy.hpp
@@ -27,41 +27,21 @@
#ifndef _ENERGY_H
#define _ENERGY_H
-/* include the energy files */
#include "ParticleRange.hpp"
#include "actor/ActorList.hpp"
-#include "statistics.hpp"
-
-/** \name Exported Variables */
-/************************************************************/
-/*@{*/
-///
-extern Observable_stat energy, total_energy;
extern ActorList energyActors;
-/*@}*/
-
-/** \name Exported Functions */
-/************************************************************/
-/*@{*/
-/** allocate energy arrays and initialize with zero */
-void init_energies(Observable_stat *stat);
+/** Recalculate energies (only if necessary). */
+void update_energy();
-/** on the master node: calc energies only if necessary */
-void master_energy_calc();
+/** Parallel energy calculation. */
+void energy_calc(double time);
-/** parallel energy calculation.
- @param result non-zero only on master node; will contain the cumulative over
- all nodes. */
-void energy_calc(double *result, double time);
-
-/** Calculate long range energies (P3M, ...). */
+/** Calculate long-range energies (P3M, ...). */
void calc_long_range_energies(const ParticleRange &particles);
-/** Calculate the total energy */
+/** Calculate the total energy of the system. */
double calculate_current_potential_energy_of_system();
-/*@}*/
-
#endif
diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp
index 6e82b18a422..15d9cd3d6c9 100644
--- a/src/core/energy_inline.hpp
+++ b/src/core/energy_inline.hpp
@@ -27,6 +27,7 @@
#include "config.hpp"
#include
+#include "Observable_stat.hpp"
#include "bonded_interactions/angle_cosine.hpp"
#include "bonded_interactions/angle_cossquare.hpp"
#include "bonded_interactions/angle_harmonic.hpp"
@@ -63,7 +64,6 @@
#endif
#include "cells.hpp"
#include "exclusions.hpp"
-#include "statistics.hpp"
#include "energy.hpp"
@@ -71,6 +71,8 @@
#include "electrostatics_magnetostatics/dipole_inline.hpp"
#endif
+extern Observable_stat_wrapper obs_energy;
+
/** Calculate non-bonded energies between a pair of particles.
* @param p1 particle 1.
* @param p2 particle 2.
@@ -175,7 +177,7 @@ inline double calc_non_bonded_pair_energy(Particle const &p1,
}
/** Add non-bonded and short-range Coulomb energies between a pair of particles
- * to the @ref energy observable.
+ * to the energy observable.
* @param p1 particle 1.
* @param p2 particle 2.
* @param d vector between p1 and p2.
@@ -190,16 +192,18 @@ inline void add_non_bonded_pair_energy(Particle const &p1, Particle const &p2,
#ifdef EXCLUSIONS
if (do_nonbonded(p1, p2))
#endif
- *obsstat_nonbonded(&energy, p1.p.type, p2.p.type) +=
+ obs_energy.local.non_bonded_contribution(p1.p.type, p2.p.type)[0] +=
calc_non_bonded_pair_energy(p1, p2, ia_params, d, dist);
#ifdef ELECTROSTATICS
- energy.coulomb[0] +=
- Coulomb::pair_energy(p1, p2, p1.p.q * p2.p.q, d, dist, dist2);
+ if (!obs_energy.local.coulomb.empty())
+ obs_energy.local.coulomb[0] +=
+ Coulomb::pair_energy(p1, p2, p1.p.q * p2.p.q, d, dist, dist2);
#endif
#ifdef DIPOLES
- energy.dipolar[0] += Dipole::pair_energy(p1, p2, d, dist, dist2);
+ if (!obs_energy.local.dipolar.empty())
+ obs_energy.local.dipolar[0] += Dipole::pair_energy(p1, p2, d, dist, dist2);
#endif
}
@@ -278,7 +282,7 @@ calc_bonded_energy(Bonded_ia_parameters const &iaparams, Particle const &p1,
throw BondInvalidSizeError(n_partners);
}
-/** Add bonded energies for one particle to the @ref energy observable.
+/** 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.
@@ -292,7 +296,7 @@ inline bool add_bonded_energy(Particle &p1, int bond_id,
auto const result = calc_bonded_energy(iaparams, p1, partners);
if (result) {
- *obsstat_bonded(&energy, bond_id) += result.get();
+ obs_energy.local.bonded_contribution(bond_id)[0] += result.get();
return false;
}
@@ -300,7 +304,7 @@ inline bool add_bonded_energy(Particle &p1, int bond_id,
return true;
}
-/** Add kinetic energies for one particle to the @ref energy observable.
+/** Add kinetic energies for one particle to the energy observable.
* @param[in] p1 particle for which to calculate energies
*/
inline void add_kinetic_energy(Particle const &p1) {
@@ -308,24 +312,23 @@ inline void add_kinetic_energy(Particle const &p1) {
return;
/* kinetic energy */
- if (not p1.p.is_virtual)
- energy.data[0] += 0.5 * p1.p.mass * p1.m.v.norm2();
+ obs_energy.local.kinetic[0] += 0.5 * p1.p.mass * p1.m.v.norm2();
- // Note that rotational degrees of virtual sites are integrated
- // and therefore can contribute to kinetic energy
+ // Note that rotational degrees of virtual sites are integrated
+ // and therefore can contribute to kinetic energy
#ifdef ROTATION
if (p1.p.rotation) {
/* the rotational part is added to the total kinetic energy;
- Here we use the rotational inertia */
- energy.data[0] += 0.5 * (Utils::sqr(p1.m.omega[0]) * p1.p.rinertia[0] +
- Utils::sqr(p1.m.omega[1]) * p1.p.rinertia[1] +
- Utils::sqr(p1.m.omega[2]) * p1.p.rinertia[2]);
+ * Here we use the rotational inertia */
+ obs_energy.local.kinetic[0] +=
+ 0.5 * (Utils::sqr(p1.m.omega[0]) * p1.p.rinertia[0] +
+ Utils::sqr(p1.m.omega[1]) * p1.p.rinertia[1] +
+ Utils::sqr(p1.m.omega[2]) * p1.p.rinertia[2]);
}
#endif
}
-/** Add kinetic and bonded energies for one particle to the @ref energy
- * observable.
+/** Add kinetic and bonded energies for one particle to the energy observable.
* @param[in] p particle for which to calculate energies
*/
inline void add_single_particle_energy(Particle &p) {
diff --git a/src/core/event.cpp b/src/core/event.cpp
index 898e239eee2..76684d89325 100644
--- a/src/core/event.cpp
+++ b/src/core/event.cpp
@@ -25,6 +25,7 @@
*/
#include "event.hpp"
+#include "Observable_stat.hpp"
#include "Particle.hpp"
#include "bonded_interactions/thermalized_bond.hpp"
#include "cells.hpp"
@@ -42,7 +43,6 @@
#include "npt.hpp"
#include "partCfg_global.hpp"
#include "particle_data.hpp"
-#include "statistics.hpp"
#include "thermostat.hpp"
#include "virtual_sites.hpp"
@@ -50,6 +50,7 @@
#include "electrostatics_magnetostatics/coulomb.hpp"
#include "electrostatics_magnetostatics/dipole.hpp"
+#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#ifdef SCAFACOS
#include "electrostatics_magnetostatics/scafacos.hpp"
diff --git a/src/core/event.hpp b/src/core/event.hpp
index ed099bfeb05..b49a9c61fcd 100644
--- a/src/core/event.hpp
+++ b/src/core/event.hpp
@@ -107,8 +107,8 @@ unsigned global_ghost_flags();
/** called every time the walls for the lb fluid are changed */
void on_lbboundary_change();
-/** @brief Update particles with properties depending on other particles
- * e.g., virtual sites, ICC charges
+/** @brief Update particles with properties depending on other particles,
+ * namely virtual sites and ICC charges.
*/
void update_dependent_particles();
diff --git a/src/core/grid_based_algorithms/lb_boundaries.cpp b/src/core/grid_based_algorithms/lb_boundaries.cpp
index 3b5feec4c75..c404dbddbce 100644
--- a/src/core/grid_based_algorithms/lb_boundaries.cpp
+++ b/src/core/grid_based_algorithms/lb_boundaries.cpp
@@ -271,7 +271,7 @@ Utils::Vector3d lbboundary_get_force(LBBoundary const *lbb) {
#endif
} else if (lattice_switch == ActiveLB::CPU) {
#if defined(LB_BOUNDARIES)
- mpi_gather_stats(8, forces.data(), nullptr, nullptr, nullptr);
+ mpi_gather_stats(GatherStats::lb_boundary_forces, forces.data());
#endif
}
auto const container_index = std::distance(lbboundaries.begin(), it);
diff --git a/src/core/grid_based_algorithms/lb_interface.cpp b/src/core/grid_based_algorithms/lb_interface.cpp
index 5ad9b0c48d7..d0184e5e6e9 100644
--- a/src/core/grid_based_algorithms/lb_interface.cpp
+++ b/src/core/grid_based_algorithms/lb_interface.cpp
@@ -1223,7 +1223,7 @@ Utils::Vector3d lb_lbfluid_calc_fluid_momentum() {
lb_calc_fluid_momentum_GPU(fluid_momentum.data());
#endif
} else if (lattice_switch == ActiveLB::CPU) {
- mpi_gather_stats(6, fluid_momentum.data(), nullptr, nullptr, nullptr);
+ mpi_gather_stats(GatherStats::lb_fluid_momentum, fluid_momentum.data());
}
return fluid_momentum;
}
diff --git a/src/core/grid_based_algorithms/lbgpu.cpp b/src/core/grid_based_algorithms/lbgpu.cpp
index 4da999a5eed..f7d477bfca8 100644
--- a/src/core/grid_based_algorithms/lbgpu.cpp
+++ b/src/core/grid_based_algorithms/lbgpu.cpp
@@ -36,7 +36,6 @@
#include "grid_based_algorithms/lbgpu.hpp"
#include "integrate.hpp"
#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "statistics.hpp"
#include
diff --git a/src/core/io/writer/h5md_core.cpp b/src/core/io/writer/h5md_core.cpp
index 4b0f2aaeb8d..ab64a602d33 100644
--- a/src/core/io/writer/h5md_core.cpp
+++ b/src/core/io/writer/h5md_core.cpp
@@ -35,7 +35,7 @@ static void backup_file(const std::string &from, const std::string &to) {
if (this_node == 0) {
/*
* If the file itself *and* a backup file exists, something must
- * have went wrong.
+ * have gone wrong.
*/
boost::filesystem::path pfrom(from), pto(to);
try {
diff --git a/src/core/observables/StressTensor.hpp b/src/core/observables/StressTensor.hpp
index 4dad362baf4..9ff584bf364 100644
--- a/src/core/observables/StressTensor.hpp
+++ b/src/core/observables/StressTensor.hpp
@@ -30,7 +30,7 @@ class StressTensor : public Observable {
std::vector shape() const override { return {3, 3}; }
std::vector operator()() const override {
std::vector res(n_values());
- observable_compute_stress_tensor(1, res.data());
+ observable_compute_stress_tensor(res.data());
return res;
}
};
diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp
index f90fd1a9d71..827dde7e407 100644
--- a/src/core/pressure.cpp
+++ b/src/core/pressure.cpp
@@ -22,6 +22,7 @@
* Implementation of pressure.hpp.
*/
+#include "Observable_stat.hpp"
#include "cells.hpp"
#include "communication.hpp"
#include "event.hpp"
@@ -36,17 +37,18 @@
#include "electrostatics_magnetostatics/coulomb.hpp"
#include "electrostatics_magnetostatics/dipole.hpp"
-Observable_stat virials{};
-Observable_stat total_pressure{};
-Observable_stat p_tensor{};
-Observable_stat total_p_tensor{};
-
-/* Observables used in the calculation of intra- and inter- molecular
- non-bonded contributions to pressure and to stress tensor */
-Observable_stat_non_bonded virials_non_bonded{};
-Observable_stat_non_bonded total_pressure_non_bonded{};
-Observable_stat_non_bonded p_tensor_non_bonded{};
-Observable_stat_non_bonded total_p_tensor_non_bonded{};
+/** Scalar pressure of the system */
+Observable_stat_wrapper obs_scalar_pressure{1};
+/** Stress tensor of the system */
+Observable_stat_wrapper obs_stress_tensor{9};
+/** Contribution from the intra- and inter-molecular non-bonded interactions
+ * to the scalar pressure of the system.
+ */
+Observable_stat_non_bonded_wrapper obs_scalar_pressure_non_bonded{1};
+/** Contribution from the intra- and inter-molecular non-bonded interactions
+ * to the stress tensor of the system.
+ */
+Observable_stat_non_bonded_wrapper obs_stress_tensor_non_bonded{9};
nptiso_struct nptiso = {0.0,
0.0,
@@ -63,54 +65,33 @@ nptiso_struct nptiso = {0.0,
false,
0};
-/************************************************************/
-/* local prototypes */
-/************************************************************/
-
-/** Calculate long range virials (P3M, ...). */
-void calc_long_range_virials(const ParticleRange &particles);
-
-/** Initializes a virials Observable stat. */
-void init_virials(Observable_stat *stat);
-
-/** Initializes a virials Observable stat. */
-void init_virials_non_bonded(Observable_stat_non_bonded *stat_nb);
-
-/** on the master node: calc energies only if necessary
- * @param v_comp flag which enables (1) compensation of the velocities required
- * for deriving a pressure reflecting \ref nptiso_struct::p_inst
- * (hence it only works with domain decomposition); naturally it
- * therefore doesn't make sense to use it without NpT.
- */
-void master_pressure_calc(int v_comp);
-
-/** Initializes stat to be used by \ref pressure_calc. */
-void init_p_tensor(Observable_stat *stat);
-
-/** Initializes stat_nb to be used by \ref pressure_calc. */
-void init_p_tensor_non_bonded(Observable_stat_non_bonded *stat_nb);
+/** 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.local,
+ obs_stress_tensor.local, particles);
+#endif
+#ifdef DIPOLES
+ /* calculate k-space part of magnetostatic interaction. */
+ Dipole::calc_pressure_long_range();
+#endif
+}
-/*********************************/
-/* Scalar and Tensorial Pressure */
-/*********************************/
static void add_single_particle_virials(Particle &p) {
cell_structure.execute_bond_handler(p, add_bonded_stress);
}
-void pressure_calc(double *result, double *result_t, double *result_nb,
- double *result_t_nb, int v_comp) {
+void pressure_calc(bool v_comp) {
auto const volume = box_geo.volume();
if (!interactions_sanity_checks())
return;
- init_virials(&virials);
-
- init_p_tensor(&p_tensor);
-
- init_virials_non_bonded(&virials_non_bonded);
-
- init_p_tensor_non_bonded(&p_tensor_non_bonded);
+ obs_scalar_pressure.local.resize_and_clear();
+ obs_scalar_pressure_non_bonded.local.resize_and_clear();
+ obs_stress_tensor.local.resize_and_clear();
+ obs_stress_tensor_non_bonded.local.resize_and_clear();
on_observable_calc();
@@ -124,183 +105,89 @@ void pressure_calc(double *result, double *result_t, double *result_nb,
sqrt(d.dist2));
});
- /* rescale kinetic energy (=ideal contribution) */
- virials.data[0] /= (3.0 * volume * time_step * time_step);
-
calc_long_range_virials(cell_structure.local_particles());
#ifdef VIRTUAL_SITES
- {
+ if (!obs_scalar_pressure.local.virtual_sites.empty()) {
auto const vs_stress = virtual_sites()->stress_tensor();
- *virials.virtual_sites += trace(vs_stress);
- boost::copy(flatten(vs_stress), p_tensor.virtual_sites);
+ obs_scalar_pressure.local.virtual_sites[0] += trace(vs_stress);
+ boost::copy(flatten(vs_stress),
+ obs_stress_tensor.local.virtual_sites.begin());
}
#endif
- for (size_t n = 1; n < virials.data.size(); n++)
- virials.data[n] /= 3.0 * volume;
-
- for (int i = 0; i < 9; i++)
- p_tensor.data[i] /= (volume * time_step * time_step);
+ /* rescale kinetic energy (=ideal contribution) */
+ obs_scalar_pressure.local.rescale(3.0 * volume);
- for (size_t i = 9; i < p_tensor.data.size(); i++)
- p_tensor.data[i] /= volume;
+ obs_stress_tensor.local.rescale(volume);
/* Intra- and Inter- part of nonbonded interaction */
- for (auto &value : virials_non_bonded.data_nb) {
- value /= 3.0 * volume;
- }
+ obs_scalar_pressure_non_bonded.local.rescale(3.0 * volume);
- for (auto &value : p_tensor_non_bonded.data_nb) {
- value /= volume;
- }
+ obs_stress_tensor_non_bonded.local.rescale(volume);
/* gather data */
- MPI_Reduce(virials.data.data(), result, virials.data.size(), MPI_DOUBLE,
- MPI_SUM, 0, comm_cart);
- MPI_Reduce(p_tensor.data.data(), result_t, p_tensor.data.size(), MPI_DOUBLE,
- MPI_SUM, 0, comm_cart);
-
- MPI_Reduce(virials_non_bonded.data_nb.data(), result_nb,
- virials_non_bonded.data_nb.size(), MPI_DOUBLE, MPI_SUM, 0,
- comm_cart);
- MPI_Reduce(p_tensor_non_bonded.data_nb.data(), result_t_nb,
- p_tensor_non_bonded.data_nb.size(), MPI_DOUBLE, MPI_SUM, 0,
- comm_cart);
-}
-
-/************************************************************/
-
-void calc_long_range_virials(const ParticleRange &particles) {
-#ifdef ELECTROSTATICS
- /* calculate k-space part of electrostatic interaction. */
- Coulomb::calc_pressure_long_range(virials, p_tensor, particles);
-#endif /*ifdef ELECTROSTATICS */
-
-#ifdef DIPOLES
- /* calculate k-space part of magnetostatic interaction. */
- Dipole::calc_pressure_long_range();
-#endif /*ifdef DIPOLES */
-}
-
-/* Initialize the virials used in the calculation of the scalar pressure */
-/************************************************************/
-void init_virials(Observable_stat *stat) {
- // Determine number of contribution for different interaction types
- // bonded, nonbonded, Coulomb, dipolar, rigid bodies
- int n_pre, n_non_bonded, n_coulomb(0), n_dipolar(0), n_vs(0);
-
- n_pre = 1;
- n_non_bonded = (max_seen_particle_type * (max_seen_particle_type + 1)) / 2;
-
-#ifdef ELECTROSTATICS
- Coulomb::pressure_n(n_coulomb);
-#endif
-#ifdef DIPOLES
- Dipole::pressure_n();
-#endif
-#ifdef VIRTUAL_SITES
- n_vs = 1;
-#endif
-
- // Allocate memory for the data
- obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded,
- n_coulomb, n_dipolar, n_vs, 1);
- stat->init_status = 0;
-}
-
-/************************************************************/
-void init_virials_non_bonded(Observable_stat_non_bonded *stat_nb) {
- int n_non_bonded;
-
- n_non_bonded = (max_seen_particle_type * (max_seen_particle_type + 1)) / 2;
-
- obsstat_realloc_and_clear_non_bonded(stat_nb, n_non_bonded, 1);
+ obs_scalar_pressure.reduce();
+ obs_stress_tensor.reduce();
+ obs_scalar_pressure_non_bonded.reduce();
+ obs_stress_tensor_non_bonded.reduce();
}
-/* Initialize the p_tensor */
-/***************************/
-void init_p_tensor(Observable_stat *stat) {
- // Determine number of contribution for different interaction types
- // bonded, nonbonded, Coulomb, dipolar, rigid bodies
- int n_pre, n_non_bonded, n_coulomb(0), n_vs(0);
-
- n_pre = 1;
- n_non_bonded = (max_seen_particle_type * (max_seen_particle_type + 1)) / 2;
-
-#ifdef ELECTROSTATICS
- Coulomb::pressure_n(n_coulomb);
-#endif
-#ifdef DIPOLES
- auto const n_dipolar = Dipole::pressure_n();
-#else
- auto const n_dipolar = 0;
-#endif
-#ifdef VIRTUAL_SITES
- n_vs = 1;
-#endif
-
- obsstat_realloc_and_clear(stat, n_pre, bonded_ia_params.size(), n_non_bonded,
- n_coulomb, n_dipolar, n_vs, 9);
- stat->init_status = 0;
-}
-
-/***************************/
-void init_p_tensor_non_bonded(Observable_stat_non_bonded *stat_nb) {
- int n_nonbonded;
- n_nonbonded = (max_seen_particle_type * (max_seen_particle_type + 1)) / 2;
-
- obsstat_realloc_and_clear_non_bonded(stat_nb, n_nonbonded, 9);
+/** Reduce the system scalar pressure and stress tensor from all MPI ranks.
+ * @param v_comp flag which enables compensation of the velocities required
+ * for deriving a pressure reflecting \ref nptiso_struct::p_inst
+ * (hence it only works with domain decomposition); naturally it
+ * therefore doesn't make sense to use it without NpT.
+ */
+void master_pressure_calc(bool v_comp) {
+ mpi_gather_stats(v_comp ? GatherStats::pressure_v_comp
+ : GatherStats::pressure);
+
+ obs_scalar_pressure.is_initialized = true;
+ obs_stress_tensor.is_initialized = true;
+ obs_scalar_pressure.v_comp = v_comp;
+ obs_stress_tensor.v_comp = v_comp;
}
-/************************************************************/
-void master_pressure_calc(int v_comp) {
- if (v_comp)
- mpi_gather_stats(3, total_pressure.data.data(), total_p_tensor.data.data(),
- total_pressure_non_bonded.data_nb.data(),
- total_p_tensor_non_bonded.data_nb.data());
- else
- mpi_gather_stats(2, total_pressure.data.data(), total_p_tensor.data.data(),
- total_pressure_non_bonded.data_nb.data(),
- total_p_tensor_non_bonded.data_nb.data());
-
- total_pressure.init_status = 1 + v_comp;
- total_p_tensor.init_status = 1 + v_comp;
+/** Calculate the sum of the ideal gas components of the instantaneous
+ * pressure, rescaled by the box dimensions.
+ */
+double calculate_npt_p_vel_rescaled() {
+ Utils::Vector3d p_vel;
+ MPI_Reduce(nptiso.p_vel.data(), p_vel.data(), 3, MPI_DOUBLE, MPI_SUM, 0,
+ MPI_COMM_WORLD);
+ double acc = 0.0;
+ for (int i = 0; i < 3; i++)
+ if (nptiso.geometry & nptiso.nptgeom_dir[i])
+ acc += p_vel[i];
+ return acc / (nptiso.dimension * nptiso.volume);
}
-/************************************************************/
-int observable_compute_stress_tensor(int v_comp, double *A) {
- /* if desired (v_comp==1) replace ideal component with instantaneous one */
- if (total_pressure.init_status != 1 + v_comp) {
- init_virials(&total_pressure);
- init_p_tensor(&total_p_tensor);
-
- init_virials_non_bonded(&total_pressure_non_bonded);
- init_p_tensor_non_bonded(&total_p_tensor_non_bonded);
+void update_pressure(bool v_comp) {
+ if (obs_scalar_pressure.is_initialized &&
+ obs_scalar_pressure.v_comp == v_comp)
+ return;
- if (v_comp && (integ_switch == INTEG_METHOD_NPT_ISO) &&
- !(nptiso.invalidate_p_vel)) {
- if (total_pressure.init_status == 0)
- master_pressure_calc(0);
- p_tensor.data[0] = 0.0;
- Utils::Vector3d p_vel;
- MPI_Reduce(nptiso.p_vel.data(), p_vel.data(), 3, MPI_DOUBLE, MPI_SUM, 0,
- MPI_COMM_WORLD);
- for (int i = 0; i < 3; i++)
- if (nptiso.geometry & nptiso.nptgeom_dir[i])
- p_tensor.data[0] += p_vel[i];
- p_tensor.data[0] /= (nptiso.dimension * nptiso.volume);
- total_pressure.init_status = 1 + v_comp;
- } else
- master_pressure_calc(v_comp);
+ obs_scalar_pressure.resize();
+ obs_scalar_pressure_non_bonded.resize();
+ obs_stress_tensor.resize();
+ obs_stress_tensor_non_bonded.resize();
+
+ if (v_comp && (integ_switch == INTEG_METHOD_NPT_ISO) &&
+ !(nptiso.invalidate_p_vel)) {
+ if (!obs_scalar_pressure.is_initialized)
+ master_pressure_calc(false);
+ obs_scalar_pressure.kinetic[0] = calculate_npt_p_vel_rescaled();
+ obs_scalar_pressure.v_comp = true;
+ } else {
+ master_pressure_calc(v_comp);
}
+}
- for (int j = 0; j < 9; j++) {
- double value = total_p_tensor.data[j];
- for (size_t i = 1; i < total_p_tensor.data.size() / 9; i++)
- value += total_p_tensor.data[9 * i + j];
- A[j] = value;
+void observable_compute_stress_tensor(double *output) {
+ update_pressure(true);
+ for (size_t j = 0; j < 9; j++) {
+ output[j] = obs_stress_tensor.accumulate(0, j);
}
- return 0;
}
diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp
index c70632378ff..10111fdaefa 100644
--- a/src/core/pressure.hpp
+++ b/src/core/pressure.hpp
@@ -25,48 +25,27 @@
#ifndef CORE_PRESSURE_HPP
#define CORE_PRESSURE_HPP
-#include "statistics.hpp"
-
-/** \name Exported Variables */
-/************************************************************/
-/*@{*/
-///
-extern Observable_stat virials, total_pressure, p_tensor, total_p_tensor;
-///
-extern Observable_stat_non_bonded virials_non_bonded, total_pressure_non_bonded,
- p_tensor_non_bonded, total_p_tensor_non_bonded;
-/*@}*/
-
-/** \name Exported Functions */
-/************************************************************/
-/*@{*/
-void init_virials(Observable_stat *stat);
-void init_virials_non_bonded(Observable_stat_non_bonded *stat_nb);
-void init_p_tensor_non_bonded(Observable_stat_non_bonded *stat_nb);
-void init_p_tensor(Observable_stat *stat);
-void master_pressure_calc(int v_comp);
-
-/** Calculates the pressure in the system from a virial expansion.
- * @param[out] result Calculated scalar pressure
- * @param[out] result_t Calculated stress tensor
- * @param[out] result_nb Calculated intra- and inter-molecular nonbonded
- * contributions to the scalar pressure
- * @param[out] result_t_nb Calculated intra- and inter-molecular nonbonded
- * contributions to the stress tensor
- * @param[in] v_comp flag which enables (1) compensation of the velocities
+/** Parallel pressure calculation from a virial expansion.
+ * @param[in] v_comp flag which enables compensation of the velocities
* required for deriving a pressure reflecting
* \ref nptiso_struct::p_inst (hence it only works with
* domain decomposition); naturally it therefore doesn't
* make sense to use it without NpT.
*/
-void pressure_calc(double *result, double *result_t, double *result_nb,
- double *result_t_nb, int v_comp);
+void pressure_calc(bool v_comp);
-/** Function to calculate stress tensor for the observables */
-int observable_compute_stress_tensor(int v_comp, double *A);
-
-void update_pressure(int v_comp);
+/** Helper function for @ref Observables::StressTensor.
+ * @param[out] output Destination array.
+ */
+void observable_compute_stress_tensor(double *output);
-/*@}*/
+/** Recalculate the virials, pressure and stress tensors (only if necessary).
+ * @param[in] v_comp flag which enables compensation of the velocities
+ * required for deriving a pressure reflecting
+ * \ref nptiso_struct::p_inst (hence it only works with
+ * domain decomposition); naturally it therefore doesn't
+ * make sense to use it without NpT.
+ */
+void update_pressure(bool v_comp);
#endif
diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp
index 82b285deb6a..7ea538f9bb6 100644
--- a/src/core/pressure_inline.hpp
+++ b/src/core/pressure_inline.hpp
@@ -25,6 +25,7 @@
#ifndef CORE_PRESSURE_INLINE_HPP
#define CORE_PRESSURE_INLINE_HPP
+#include "Observable_stat.hpp"
#include "exclusions.hpp"
#include "forces_inline.hpp"
#include "integrate.hpp"
@@ -33,6 +34,11 @@
#include
+extern Observable_stat_wrapper obs_scalar_pressure;
+extern Observable_stat_wrapper obs_stress_tensor;
+extern Observable_stat_non_bonded_wrapper obs_scalar_pressure_non_bonded;
+extern Observable_stat_non_bonded_wrapper obs_stress_tensor_non_bonded;
+
/** Calculate non bonded energies between a pair of particles.
* @param p1 pointer to particle 1.
* @param p2 pointer to particle 2.
@@ -46,47 +52,48 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2,
#endif
{
auto const force = calc_non_bonded_pair_force(p1, p2, d, dist);
- *obsstat_nonbonded(&virials, p1.p.type, p2.p.type) += d * force;
+ obs_scalar_pressure.local.non_bonded_contribution(
+ p1.p.type, p2.p.type)[0] += d * force;
/* stress tensor part */
for (int k = 0; k < 3; k++)
for (int l = 0; l < 3; l++)
- obsstat_nonbonded(&p_tensor, p1.p.type, p2.p.type)[k * 3 + l] +=
- force[k] * d[l];
+ obs_stress_tensor.local.non_bonded_contribution(
+ p1.p.type, p2.p.type)[k * 3 + l] += force[k] * d[l];
auto const p1molid = p1.p.mol_id;
auto const p2molid = p2.p.mol_id;
if (p1molid == p2molid) {
- *obsstat_nonbonded_intra(&virials_non_bonded, p1.p.type, p2.p.type) +=
- d * force;
+ obs_scalar_pressure_non_bonded.local.non_bonded_intra_contribution(
+ p1.p.type, p2.p.type)[0] += d * force;
for (int k = 0; k < 3; k++)
for (int l = 0; l < 3; l++)
- obsstat_nonbonded_intra(&p_tensor_non_bonded, p1.p.type,
- p2.p.type)[k * 3 + l] += force[k] * d[l];
+ obs_stress_tensor_non_bonded.local.non_bonded_intra_contribution(
+ p1.p.type, p2.p.type)[k * 3 + l] += force[k] * d[l];
} else {
- *obsstat_nonbonded_inter(&virials_non_bonded, p1.p.type, p2.p.type) +=
- d * force;
+ obs_scalar_pressure_non_bonded.local.non_bonded_inter_contribution(
+ p1.p.type, p2.p.type)[0] += d * force;
for (int k = 0; k < 3; k++)
for (int l = 0; l < 3; l++)
- obsstat_nonbonded_inter(&p_tensor_non_bonded, p1.p.type,
- p2.p.type)[k * 3 + l] += force[k] * d[l];
+ obs_stress_tensor_non_bonded.local.non_bonded_inter_contribution(
+ p1.p.type, p2.p.type)[k * 3 + l] += force[k] * d[l];
}
}
#ifdef ELECTROSTATICS
- {
+ if (!obs_scalar_pressure.local.coulomb.empty()) {
/* real space Coulomb */
auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist);
for (int i = 0; i < 3; i++) {
for (int j = 0; j < 3; j++) {
- p_tensor.coulomb[i * 3 + j] += p_coulomb[i][j];
+ obs_stress_tensor.local.coulomb[i * 3 + j] += p_coulomb[i][j];
}
}
- virials.coulomb[0] += trace(p_coulomb);
+ obs_scalar_pressure.local.coulomb[0] += trace(p_coulomb);
}
#endif /*ifdef ELECTROSTATICS */
@@ -167,12 +174,13 @@ inline bool add_bonded_stress(Particle &p1, int bond_id,
if (result) {
auto const &stress = result.get();
- *obsstat_bonded(&virials, bond_id) += trace(stress);
+ obs_scalar_pressure.local.bonded_contribution(bond_id)[0] += trace(stress);
/* stress tensor part */
for (int k = 0; k < 3; k++)
for (int l = 0; l < 3; l++)
- obsstat_bonded(&p_tensor, bond_id)[k * 3 + l] += stress[k][l];
+ obs_stress_tensor.local.bonded_contribution(bond_id)[k * 3 + l] +=
+ stress[k][l];
return false;
}
@@ -183,27 +191,26 @@ inline bool add_bonded_stress(Particle &p1, int bond_id,
/** Calculate kinetic pressure (aka energy) for one particle.
* @param p1 particle for which to calculate pressure
* @param v_comp flag which enables compensation of the velocities required
- * for deriving a pressure reflecting \ref
- * nptiso_struct::p_inst (hence it only works with domain decomposition);
- * naturally it therefore doesn't make sense to use it without NpT.
+ * for deriving a pressure reflecting \ref nptiso_struct::p_inst
+ * (hence it only works with domain decomposition); naturally it
+ * therefore doesn't make sense to use it without NpT.
*/
-inline void add_kinetic_virials(Particle const &p1, int v_comp) {
- if (not p1.p.is_virtual) {
- /* kinetic energy */
- if (v_comp) {
- virials.data[0] += ((p1.m.v * time_step) -
- (p1.f.f * (0.5 * Utils::sqr(time_step) / p1.p.mass)))
- .norm2() *
- p1.p.mass;
- } else {
- virials.data[0] += Utils::sqr(time_step) * p1.m.v.norm2() * p1.p.mass;
- }
-
- for (int k = 0; k < 3; k++)
- for (int l = 0; l < 3; l++)
- p_tensor.data[k * 3 + l] +=
- (p1.m.v[k] * time_step) * (p1.m.v[l] * time_step) * p1.p.mass;
+inline void add_kinetic_virials(Particle const &p1, bool v_comp) {
+ if (p1.p.is_virtual)
+ return;
+
+ /* kinetic energy */
+ if (v_comp) {
+ obs_scalar_pressure.local.kinetic[0] +=
+ (p1.m.v - (p1.f.f * (0.5 * time_step / p1.p.mass))).norm2() * p1.p.mass;
+ } else {
+ obs_scalar_pressure.local.kinetic[0] += p1.m.v.norm2() * p1.p.mass;
}
+
+ for (int k = 0; k < 3; k++)
+ for (int l = 0; l < 3; l++)
+ obs_stress_tensor.local.kinetic[k * 3 + l] +=
+ (p1.m.v[k]) * (p1.m.v[l]) * p1.p.mass;
}
#endif
diff --git a/src/core/reaction_ensemble.cpp b/src/core/reaction_ensemble.cpp
index fdf46fe5871..282a85b3084 100644
--- a/src/core/reaction_ensemble.cpp
+++ b/src/core/reaction_ensemble.cpp
@@ -25,6 +25,7 @@
#include "grid.hpp"
#include "integrate.hpp"
#include "partCfg_global.hpp"
+#include "statistics.hpp"
#include
#include
diff --git a/src/core/short_range_loop.hpp b/src/core/short_range_loop.hpp
index b531886ba15..b826e822449 100644
--- a/src/core/short_range_loop.hpp
+++ b/src/core/short_range_loop.hpp
@@ -55,9 +55,9 @@ struct EuclidianDistance {
};
/**
- * @brief Decided which distance function to use depending on the
- cell system, and call the pair code.
-*/
+ * @brief Decide which distance function to use depending on the
+ * cell system, and call the pair code.
+ */
template
void decide_distance(CellIterator first, CellIterator last,
@@ -82,8 +82,7 @@ void decide_distance(CellIterator first, CellIterator last,
}
/**
- * @brief Functor that returns true for
- * any arguments.
+ * @brief Functor that returns true for any argument.
*/
struct True {
template bool operator()(T...) const { return true; }
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 963ec2bf748..763ed8dcce4 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -27,23 +27,17 @@
#include "statistics.hpp"
#include "Particle.hpp"
-#include "bonded_interactions/bonded_interaction_data.hpp"
+#include "cells.hpp"
#include "communication.hpp"
-#include "energy.hpp"
#include "errorhandling.hpp"
#include "grid.hpp"
#include "grid_based_algorithms/lb_interface.hpp"
-#include "integrate.hpp"
-#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
-#include "npt.hpp"
#include "partCfg_global.hpp"
-#include "pressure.hpp"
-#include "short_range_loop.hpp"
-#include
#include
#include
#include
+#include
#include
#include
@@ -513,86 +507,3 @@ void analyze_append(PartCfg &partCfg) {
}
configs.emplace_back(config);
}
-
-/****************************************************************************************
- * Observables handling
- ****************************************************************************************/
-
-void obsstat_realloc_and_clear(Observable_stat *stat, int n_pre, int n_bonded,
- int n_non_bonded, int n_coulomb, int n_dipolar,
- int n_vs, int c_size) {
-
- // Number of doubles to store pressure in
- const int total =
- c_size *
- (n_pre + static_cast(bonded_ia_params.size()) + n_non_bonded +
- n_coulomb + n_dipolar + n_vs + Observable_stat::n_external_field);
-
- // Allocate mem for the double list
- stat->data.resize(total);
-
- // Number of doubles per interaction (pressure=1, stress tensor=9,...)
- stat->chunk_size = c_size;
-
- // Number of chunks for different interaction types
- stat->n_coulomb = n_coulomb;
- stat->n_dipolar = n_dipolar;
- stat->n_virtual_sites = n_vs;
- // Pointers to the start of different contributions
- stat->bonded = stat->data.data() + c_size * n_pre;
- stat->non_bonded = stat->bonded + c_size * bonded_ia_params.size();
- stat->coulomb = stat->non_bonded + c_size * n_non_bonded;
- stat->dipolar = stat->coulomb + c_size * n_coulomb;
- stat->virtual_sites = stat->dipolar + c_size * n_dipolar;
- stat->external_fields = stat->virtual_sites + c_size * n_vs;
-
- // Set all observables to zero
- for (int i = 0; i < total; i++)
- stat->data[i] = 0.0;
-}
-
-void obsstat_realloc_and_clear_non_bonded(Observable_stat_non_bonded *stat_nb,
- int n_nonbonded, int c_size) {
- auto const total = c_size * (n_nonbonded + n_nonbonded);
-
- stat_nb->data_nb.resize(total);
- stat_nb->chunk_size_nb = c_size;
- stat_nb->non_bonded_intra = stat_nb->data_nb.data();
- stat_nb->non_bonded_inter = stat_nb->non_bonded_intra + c_size * n_nonbonded;
-
- for (int i = 0; i < total; i++)
- stat_nb->data_nb[i] = 0.0;
-}
-
-void invalidate_obs() {
- total_energy.init_status = 0;
- total_pressure.init_status = 0;
- total_p_tensor.init_status = 0;
-}
-
-void update_pressure(int v_comp) {
- Utils::Vector3d p_vel;
- /* if desired (v_comp==1) replace ideal component with instantaneous one */
- if (total_pressure.init_status != 1 + v_comp) {
- init_virials(&total_pressure);
- init_p_tensor(&total_p_tensor);
-
- init_virials_non_bonded(&total_pressure_non_bonded);
- init_p_tensor_non_bonded(&total_p_tensor_non_bonded);
-
- if (v_comp && (integ_switch == INTEG_METHOD_NPT_ISO) &&
- !(nptiso.invalidate_p_vel)) {
- if (total_pressure.init_status == 0)
- master_pressure_calc(0);
- total_pressure.data[0] = 0.0;
- MPI_Reduce(nptiso.p_vel.data(), p_vel.data(), 3, MPI_DOUBLE, MPI_SUM, 0,
- MPI_COMM_WORLD);
- for (int i = 0; i < 3; i++)
- if (nptiso.geometry & nptiso.nptgeom_dir[i])
- total_pressure.data[0] += p_vel[i];
- total_pressure.data[0] /= (nptiso.dimension * nptiso.volume);
- total_pressure.init_status = 1 + v_comp;
- } else
- master_pressure_calc(v_comp);
- }
-}
diff --git a/src/core/statistics.hpp b/src/core/statistics.hpp
index 52486cb9ba7..c350de5db43 100644
--- a/src/core/statistics.hpp
+++ b/src/core/statistics.hpp
@@ -26,10 +26,8 @@
* Implementation in statistics.cpp.
*/
-#include "Observable_stat.hpp"
#include "PartCfg.hpp"
#include "Particle.hpp"
-#include "nonbonded_interactions/nonbonded_interaction_data.hpp"
#include