From 85b78a72f4eb77483be373c4a8651bc074fc0e94 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 10:44:50 +0200 Subject: [PATCH 01/11] python: Increase code reuse for Observable_stat --- src/python/espressomd/analyze.pxd | 25 +++++ src/python/espressomd/analyze.pyx | 171 ++++++++---------------------- 2 files changed, 72 insertions(+), 124 deletions(-) diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 447bc57977d..a68e41ff8c0 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -19,7 +19,9 @@ # For C-extern Analysis +cimport numpy as np from .utils cimport Vector3i, Vector3d, Vector9d, Span +from .utils cimport create_nparray_from_double_span from libcpp.vector cimport vector # import std::vector as vector from libcpp cimport bool as cbool @@ -94,3 +96,26 @@ cdef extern from "energy.hpp": cdef extern from "dpd.hpp": Vector9d dpd_stress() + +cdef inline get_obs_contribs(Span[double] contributions, int size): + """ + Convert an Observable_stat range of contributions into a correctly + shaped numpy array. + """ + cdef np.ndarray value + value = create_nparray_from_double_span(contributions) + if size == 9: + return value.reshape((-1, 3, 3)) + else: + return value + +cdef inline get_obs_contrib(Span[double] contribution, int size): + """ + Convert an Observable_stat contribution into a correctly + shaped numpy array. If the size is 1, decay to a float. + """ + cdef np.ndarray value + value = get_obs_contribs(contribution, size) + if value.shape[0] == 1: + return value[0] + return value diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 3dcd59670eb..c6b09c9a25b 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -32,12 +32,11 @@ from .utils import array_locked, is_valid_type from .utils cimport Vector3i, Vector3d, Vector9d from .utils cimport handle_errors, check_type_or_throw_except from .utils cimport create_nparray_from_double_array -from .utils cimport create_nparray_from_double_span from .particle_data cimport get_n_part -cdef _Observable_stat_to_dict(Observable_stat obs): - """Transform a scalar Observable_stat struct to a python dict. +cdef _Observable_stat_to_dict(Observable_stat obs, int size): + """Transform an Observable_stat struct to a python dict. Returns ------- @@ -52,56 +51,64 @@ cdef _Observable_stat_to_dict(Observable_stat obs): * ``"nonbonded", , ``: nonbonded contribution which arises from the interactions between type_i and type_j * ``"nonbonded_intra", , ``: nonbonded contribution between short ranged forces between type i and j and with the same mol_id * ``"nonbonded_inter", , ``: nonbonded contribution between short ranged forces between type i and j and different mol_ids - * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method. It is equivalent to 1/3 of the trace of the Coulomb contribution tensor. - For how the contribution tensor is calculated, see :ref:`Contribution Tensor`. The averaged value in an isotropic NVT simulation is equivalent to the average of - :math:`E^{coulomb}/(3V)`, see :cite:`brown95a`. + * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method * ``"coulomb", ``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) - * ``"dipolar"``: not implemented - * ``"virtual_sites"``: Contribution contribution due to virtual sites + * ``"dipolar"``: dipolar contribution, how it is calculated depends on the method + * ``"dipolar", ``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"virtual_sites"``: virtual sites contribution + * ``"virtual_sites", ``: contribution from virtual site i """ + def set_initial(): + if size == 9: + return np.zeros((3, 3), dtype=float) + else: + return 0.0 + + cdef int i + cdef int j + # Dict to store the results p = OrderedDict() - # Total pressure - p["total"] = obs.accumulate() + # Total contribution + if size == 1: + p["total"] = obs.accumulate() + else: + total = np.zeros((9,), dtype=float) + for i in range(9): + total[i] = obs.accumulate(0.0, i) + p["total"] = total.reshape((3, 3)) # Kinetic - p["kinetic"] = 0 - for kinetic_i in range(obs.kinetic.size()): - p["kinetic"] += obs.kinetic[kinetic_i] + p["kinetic"] = get_obs_contrib(obs.kinetic, size) # External - p["external_fields"] = obs.external_fields[0] + p["external_fields"] = get_obs_contrib(obs.external_fields, size) # Bonded - cdef int i - cdef double val - cdef double total_bonded - total_bonded = 0 + total_bonded = set_initial() for i in range(bonded_ia_params.size()): if bonded_ia_params[i].type != BONDED_IA_NONE: - val = obs.bonded_contribution(i)[0] + val = get_obs_contrib(obs.bonded_contribution(i), size) p["bonded", i] = val total_bonded += val p["bonded"] = total_bonded # Non-Bonded interactions, total as well as intra and inter molecular - cdef int j - cdef double total_intra - cdef double total_inter - total_inter = 0 - total_intra = 0 - + total_intra = set_initial() + total_inter = set_initial() for i in range(analyze.max_seen_particle_type): for j in range(i, analyze.max_seen_particle_type): - intra = obs.non_bonded_intra_contribution(i, j)[ - 0] + intra = get_obs_contrib( + obs.non_bonded_intra_contribution( + i, j), size) total_intra += intra p["non_bonded_intra", i, j] = intra - inter = obs.non_bonded_inter_contribution(i, j)[ - 0] + inter = get_obs_contrib( + obs.non_bonded_inter_contribution( + i, j), size) total_inter += inter p["non_bonded_inter", i, j] = inter p["non_bonded", i, j] = intra + inter @@ -113,8 +120,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs): # Electrostatics IF ELECTROSTATICS == 1: cdef np.ndarray coulomb - coulomb = create_nparray_from_double_span( - obs.coulomb) + coulomb = get_obs_contribs(obs.coulomb, size) for i in range(coulomb.shape[0]): p["coulomb", i] = coulomb[i] p["coulomb"] = np.sum(coulomb, axis=0) @@ -122,8 +128,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs): # Dipoles IF DIPOLES == 1: cdef np.ndarray dipolar - dipolar = create_nparray_from_double_span( - obs.dipolar) + dipolar = get_obs_contribs(obs.dipolar, size) for i in range(dipolar.shape[0]): p["dipolar", i] = dipolar[i] p["dipolar"] = np.sum(dipolar, axis=0) @@ -131,12 +136,10 @@ cdef _Observable_stat_to_dict(Observable_stat obs): # virtual sites IF VIRTUAL_SITES == 1: cdef np.ndarray virtual_sites - if obs.virtual_sites.size(): - virtual_sites = create_nparray_from_double_span( - obs.virtual_sites) - for i in range(virtual_sites.shape[0]): - p["virtual_sites", i] = virtual_sites[i] - p["virtual_sites"] = np.sum(virtual_sites, axis=0) + virtual_sites = get_obs_contribs(obs.virtual_sites, size) + for i in range(virtual_sites.shape[0]): + p["virtual_sites", i] = virtual_sites[i] + p["virtual_sites"] = np.sum(virtual_sites, axis=0) return p @@ -306,14 +309,14 @@ class Analysis: :math:`E^{coulomb}/(3V)`, see :cite:`brown95a`. * ``"coulomb", ``: Coulomb pressure from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) * ``"dipolar"``: not implemented - * ``"virtual_sites"``: Pressure contribution due to virtual sites + * ``"virtual_sites"``: Pressure contribution from virtual sites """ # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_scalar_pressure) + return _Observable_stat_to_dict(analyze.obs_scalar_pressure, 1) def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -340,94 +343,14 @@ class Analysis: * ``"coulomb"``: Maxwell pressure tensor, how it is calculated depends on the method * ``"coulomb", ``: Maxwell pressure tensor from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) * ``"dipolar"``: not implemented - * ``"virtual_sites"``: pressure tensor contribution for virtual sites + * ``"virtual_sites"``: pressure tensor contribution from virtual sites """ - # Dict to store the results - p = OrderedDict() - # Update in ESPResSo core analyze.update_pressure() - # Total pressure - cdef int i - total = np.zeros(9) - for i in range(9): - total[i] = analyze.obs_pressure_tensor.accumulate(0.0, i) - - p["total"] = total.reshape((3, 3)) - - # Individual components of the pressure - - # Kinetic - p["kinetic"] = create_nparray_from_double_span( - analyze.obs_pressure_tensor.kinetic).reshape((3, 3)) - - # Bonded - total_bonded = np.zeros((3, 3)) - for i in range(bonded_ia_params.size()): - if bonded_ia_params[i].type != BONDED_IA_NONE: - p["bonded", i] = np.reshape(create_nparray_from_double_span( - analyze.obs_pressure_tensor.bonded_contribution(i)), (3, 3)) - total_bonded += p["bonded", i] - p["bonded"] = total_bonded - - # Non-Bonded interactions, total as well as intra and inter molecular - cdef int j - total_non_bonded_intra = np.zeros((3, 3)) - total_non_bonded_inter = np.zeros((3, 3)) - - for i in range(analyze.max_seen_particle_type): - for j in range(i, analyze.max_seen_particle_type): - p["non_bonded_intra", i, j] = np.reshape( - create_nparray_from_double_span( - analyze.obs_pressure_tensor.non_bonded_intra_contribution(i, j)), (3, 3)) - total_non_bonded_intra += p["non_bonded_intra", i, j] - - p["non_bonded_inter", i, j] = np.reshape( - create_nparray_from_double_span( - analyze.obs_pressure_tensor.non_bonded_inter_contribution(i, j)), (3, 3)) - total_non_bonded_inter += p["non_bonded_inter", i, j] - p["non_bonded", i, j] = p["non_bonded_intra", i, j] + \ - p["non_bonded_inter", i, j] - - p["non_bonded_intra"] = total_non_bonded_intra - p["non_bonded_inter"] = total_non_bonded_inter - p["non_bonded"] = total_non_bonded_intra + total_non_bonded_inter - - # Electrostatics - IF ELECTROSTATICS == 1: - cdef np.ndarray coulomb - coulomb = np.reshape( - create_nparray_from_double_span( - analyze.obs_pressure_tensor.coulomb), (-1, 3, 3)) - for i in range(coulomb.shape[0]): - p["coulomb", i] = coulomb[i] - p["coulomb"] = np.sum(coulomb, axis=0) - - # Dipoles - IF DIPOLES == 1: - cdef np.ndarray dipolar - dipolar = np.reshape( - create_nparray_from_double_span( - analyze.obs_pressure_tensor.dipolar), (-1, 3, 3)) - for i in range(dipolar.shape[0]): - p["dipolar", i] = dipolar[i] - p["dipolar"] = np.sum(dipolar, axis=0) - - # virtual sites - IF VIRTUAL_SITES_RELATIVE == 1: - cdef np.ndarray virtual_sites - if analyze.obs_pressure_tensor.virtual_sites.size(): - virtual_sites = np.reshape( - create_nparray_from_double_span( - analyze.obs_pressure_tensor.virtual_sites), (-1, 3, 3)) - for i in range(virtual_sites.shape[0]): - p["virtual_sites", i] = virtual_sites[i] - p["virtual_sites"] = np.sum(virtual_sites, axis=0) - - return p + return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9) IF DPD == 1: def dpd_stress(self): @@ -478,7 +401,7 @@ class Analysis: analyze.update_energy() handle_errors("calc_long_range_energies failed") - return _Observable_stat_to_dict(analyze.obs_energy) + return _Observable_stat_to_dict(analyze.obs_energy, 1) def calc_re(self, chain_start=None, number_of_chains=None, chain_length=None): From 97ee461c81617906f261c2c6ea6d78a7a52a7c6d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 11:26:53 +0200 Subject: [PATCH 02/11] core: Pull out obs_scalar_pressure logic The scalar pressure can be calculated from the pressure tensor. --- src/core/Observable_stat.hpp | 2 +- .../electrostatics_magnetostatics/coulomb.cpp | 4 +-- .../electrostatics_magnetostatics/coulomb.hpp | 3 +- src/core/pressure.cpp | 28 +++++++++---------- src/core/pressure_inline.hpp | 12 ++------ 5 files changed, 19 insertions(+), 30 deletions(-) diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index 5a3c968f47b..39935ce9820 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -68,7 +68,7 @@ class Observable_stat { * @param column Which column to sum up (only relevant for multi-dimensional * observables). */ - double accumulate(double acc = 0.0, size_t column = 0) { + double accumulate(double acc = 0.0, size_t column = 0) const { assert(column < m_chunk_size); if (m_chunk_size == 1) return boost::accumulate(m_data, acc); diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 3990632f2fe..983dea0ebde 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -42,8 +42,7 @@ Coulomb_parameters coulomb; namespace Coulomb { -void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor, +void calc_pressure_long_range(Observable_stat &p_tensor, const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M @@ -60,7 +59,6 @@ void calc_pressure_long_range(Observable_stat &virials, p3m_charge_assign(particles); auto const p3m_p_tensor = p3m_calc_kspace_pressure_tensor(); std::copy_n(p3m_p_tensor.data(), 9, p_tensor.coulomb.begin() + 9); - virials.coulomb[1] = p3m_p_tensor[0] + p3m_p_tensor[4] + p3m_p_tensor[8]; break; } diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index a4a323cd0cf..3f6d8f677e4 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -60,8 +60,7 @@ struct Coulomb_parameters { extern Coulomb_parameters coulomb; namespace Coulomb { -void calc_pressure_long_range(Observable_stat &virials, - Observable_stat &p_tensor, +void calc_pressure_long_range(Observable_stat &p_tensor, const ParticleRange &particles); void sanity_checks(int &state); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 90bba9527d8..84de55ffa65 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -56,12 +56,22 @@ nptiso_struct nptiso = {0.0, false, 0}; +/** Calculate the scalar pressure from the pressure tensor. */ +Observable_stat get_scalar_pressure(Observable_stat const &pressure_tensor) { + Observable_stat scalar_pressure{1}; + auto it_scalar = scalar_pressure.data_().begin(); + auto it_tensor = pressure_tensor.data_().begin(); + auto it_tensor_end = pressure_tensor.data_().end(); + for (; it_tensor < it_tensor_end; it_tensor += 9) + *(it_scalar++) += (it_tensor[0] + it_tensor[4] + it_tensor[8]) / 3.; + return scalar_pressure; +} + /** Calculate long-range virials (P3M, ...). */ void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_pressure_long_range(obs_scalar_pressure, obs_pressure_tensor, - particles); + Coulomb::calc_pressure_long_range(obs_pressure_tensor, particles); #endif #ifdef DIPOLES /* calculate k-space part of magnetostatic interaction. */ @@ -79,7 +89,6 @@ void pressure_calc() { if (!interactions_sanity_checks()) return; - obs_scalar_pressure = Observable_stat{1}; obs_pressure_tensor = Observable_stat{9}; on_observable_calc(); @@ -97,30 +106,21 @@ void pressure_calc() { calc_long_range_virials(cell_structure.local_particles()); #ifdef VIRTUAL_SITES - if (!obs_scalar_pressure.virtual_sites.empty()) { + if (!obs_pressure_tensor.virtual_sites.empty()) { auto const vs_pressure_tensor = virtual_sites()->pressure_tensor(); - - obs_scalar_pressure.virtual_sites[0] += trace(vs_pressure_tensor); boost::copy(flatten(vs_pressure_tensor), obs_pressure_tensor.virtual_sites.begin()); } #endif - /* rescale kinetic energy (=ideal contribution) */ - obs_scalar_pressure.rescale(3.0 * volume); - obs_pressure_tensor.rescale(volume); /* gather data */ - auto obs_scalar_pressure_res = reduce(comm_cart, obs_scalar_pressure); - if (obs_scalar_pressure_res) { - std::swap(obs_scalar_pressure, *obs_scalar_pressure_res); - } - auto obs_pressure_tensor_res = reduce(comm_cart, obs_pressure_tensor); if (obs_pressure_tensor_res) { std::swap(obs_pressure_tensor, *obs_pressure_tensor_res); } + obs_scalar_pressure = get_scalar_pressure(obs_pressure_tensor); } void update_pressure() { mpi_gather_stats(GatherStats::pressure); } diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 4838023f341..e464faa2de5 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -52,14 +52,12 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, auto const type1 = p1.p.mol_id; auto const type2 = p2.p.mol_id; - obs_scalar_pressure.add_non_bonded_contribution(type1, type2, - trace(stress)); obs_pressure_tensor.add_non_bonded_contribution(type1, type2, flatten(stress)); } #ifdef ELECTROSTATICS - if (!obs_scalar_pressure.coulomb.empty()) { + if (!obs_pressure_tensor.coulomb.empty()) { /* real space Coulomb */ auto const p_coulomb = Coulomb::pair_pressure(p1, p2, d, dist); @@ -68,8 +66,6 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, obs_pressure_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; } } - - obs_scalar_pressure.coulomb[0] += trace(p_coulomb); } #endif /*ifdef ELECTROSTATICS */ @@ -150,8 +146,6 @@ inline bool add_bonded_pressure_tensor(Particle &p1, int bond_id, if (result) { auto const &tensor = result.get(); - obs_scalar_pressure.bonded_contribution(bond_id)[0] += trace(tensor); - /* pressure tensor part */ for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) @@ -171,9 +165,7 @@ inline void add_kinetic_virials(Particle const &p1) { if (p1.p.is_virtual) return; - /* kinetic energy */ - obs_scalar_pressure.kinetic[0] += p1.m.v.norm2() * p1.p.mass; - + /* kinetic pressure */ for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) obs_pressure_tensor.kinetic[k * 3 + l] += From c4a3103eea6cc30b139ac154397dc240ba6590b8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 15:32:40 +0200 Subject: [PATCH 03/11] core: Remove obs_scalar_pressure global --- src/core/pressure.cpp | 14 ---------- src/core/pressure_inline.hpp | 1 - src/python/espressomd/analyze.pxd | 37 ++++++++++++++++++++++---- src/python/espressomd/analyze.pyx | 44 +++++++++++++++++++++---------- 4 files changed, 62 insertions(+), 34 deletions(-) diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 84de55ffa65..518560d344a 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -37,8 +37,6 @@ #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" -/** Scalar pressure of the system */ -Observable_stat obs_scalar_pressure{1}; /** Pressure tensor of the system */ Observable_stat obs_pressure_tensor{9}; @@ -56,17 +54,6 @@ nptiso_struct nptiso = {0.0, false, 0}; -/** Calculate the scalar pressure from the pressure tensor. */ -Observable_stat get_scalar_pressure(Observable_stat const &pressure_tensor) { - Observable_stat scalar_pressure{1}; - auto it_scalar = scalar_pressure.data_().begin(); - auto it_tensor = pressure_tensor.data_().begin(); - auto it_tensor_end = pressure_tensor.data_().end(); - for (; it_tensor < it_tensor_end; it_tensor += 9) - *(it_scalar++) += (it_tensor[0] + it_tensor[4] + it_tensor[8]) / 3.; - return scalar_pressure; -} - /** Calculate long-range virials (P3M, ...). */ void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS @@ -120,7 +107,6 @@ void pressure_calc() { if (obs_pressure_tensor_res) { std::swap(obs_pressure_tensor, *obs_pressure_tensor_res); } - obs_scalar_pressure = get_scalar_pressure(obs_pressure_tensor); } void update_pressure() { mpi_gather_stats(GatherStats::pressure); } diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index e464faa2de5..39ea96bce25 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -32,7 +32,6 @@ #include -extern Observable_stat obs_scalar_pressure; extern Observable_stat obs_pressure_tensor; /** Calculate non bonded energies between a pair of particles. diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index a68e41ff8c0..bc8c4821dfa 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -20,6 +20,7 @@ # For C-extern Analysis cimport numpy as np +import numpy as np from .utils cimport Vector3i, Vector3d, Vector9d, Span from .utils cimport create_nparray_from_double_span from libcpp.vector cimport vector # import std::vector as vector @@ -81,7 +82,6 @@ cdef extern from "statistics_chain.hpp": array2 calc_rh(int, int, int) cdef extern from "pressure_inline.hpp": - cdef Observable_stat obs_scalar_pressure cdef Observable_stat obs_pressure_tensor cdef extern from "pressure.hpp": @@ -97,25 +97,52 @@ cdef extern from "energy.hpp": cdef extern from "dpd.hpp": Vector9d dpd_stress() -cdef inline get_obs_contribs(Span[double] contributions, int size): +cdef inline get_obs_contribs(Span[double] contributions, int size, + cbool calc_scalar_pressure): """ Convert an Observable_stat range of contributions into a correctly shaped numpy array. + + Parameters + ---------- + contributions : (N,) array_like of :obj:`float` + Flattened array of energy/pressure contributions from an observable. + size : :obj:`int`, \{1, 9\} + Dimensionality of the data. + calc_scalar_pressure : :obj:`bool` + Whether to calculate a scalar pressure (only relevant when + ``contributions`` is a pressure tensor). + """ cdef np.ndarray value value = create_nparray_from_double_span(contributions) if size == 9: - return value.reshape((-1, 3, 3)) + if calc_scalar_pressure: + return np.einsum('...ii', value.reshape((-1, 3, 3))) / 3 + else: + return value.reshape((-1, 3, 3)) else: return value -cdef inline get_obs_contrib(Span[double] contribution, int size): +cdef inline get_obs_contrib(Span[double] contribution, int size, + cbool calc_scalar_pressure): """ Convert an Observable_stat contribution into a correctly shaped numpy array. If the size is 1, decay to a float. + + Parameters + ---------- + contributions : (N,) array_like of :obj:`float` + Flattened array of energy/pressure contributions from an observable. + size : :obj:`int`, \{1, 9\} + Dimensionality of the data. + calc_scalar_pressure : :obj:`bool` + Whether to calculate a scalar pressure (only relevant when + ``contributions`` is a pressure tensor). + """ cdef np.ndarray value - value = get_obs_contribs(contribution, size) + value = get_obs_contribs(contribution, size, calc_scalar_pressure) if value.shape[0] == 1: return value[0] return value diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index c6b09c9a25b..fdf1b2a31a2 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -20,6 +20,7 @@ include "myconfig.pxi" from . cimport analyze from libcpp.vector cimport vector # import std::vector as vector +from libcpp cimport bool as cbool from .interactions cimport BONDED_IA_NONE from .interactions cimport bonded_ia_params import numpy as np @@ -35,9 +36,20 @@ from .utils cimport create_nparray_from_double_array from .particle_data cimport get_n_part -cdef _Observable_stat_to_dict(Observable_stat obs, int size): +cdef _Observable_stat_to_dict(Observable_stat obs, int size, + cbool calc_sp): """Transform an Observable_stat struct to a python dict. + Parameters + ---------- + obs : + Core observable. + size : :obj:`int`, \{1, 9\} + Dimensionality of the data. + calc_sp : :obj:`bool` + Whether to calculate a scalar pressure (only relevant when + ``obs`` is a pressure tensor observable). + Returns ------- :obj:`dict` @@ -61,7 +73,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): """ def set_initial(): - if size == 9: + if size == 9 and not calc_sp: return np.zeros((3, 3), dtype=float) else: return 0.0 @@ -79,19 +91,23 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): total = np.zeros((9,), dtype=float) for i in range(9): total[i] = obs.accumulate(0.0, i) - p["total"] = total.reshape((3, 3)) + total = total.reshape((3, 3)) + if calc_sp: + p["total"] = np.einsum('ii', total) / 3 + else: + p["total"] = total # Kinetic - p["kinetic"] = get_obs_contrib(obs.kinetic, size) + p["kinetic"] = get_obs_contrib(obs.kinetic, size, calc_sp) # External - p["external_fields"] = get_obs_contrib(obs.external_fields, size) + p["external_fields"] = get_obs_contrib(obs.external_fields, size, calc_sp) # Bonded total_bonded = set_initial() for i in range(bonded_ia_params.size()): if bonded_ia_params[i].type != BONDED_IA_NONE: - val = get_obs_contrib(obs.bonded_contribution(i), size) + val = get_obs_contrib(obs.bonded_contribution(i), size, calc_sp) p["bonded", i] = val total_bonded += val p["bonded"] = total_bonded @@ -103,12 +119,12 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): for j in range(i, analyze.max_seen_particle_type): intra = get_obs_contrib( obs.non_bonded_intra_contribution( - i, j), size) + i, j), size, calc_sp) total_intra += intra p["non_bonded_intra", i, j] = intra inter = get_obs_contrib( obs.non_bonded_inter_contribution( - i, j), size) + i, j), size, calc_sp) total_inter += inter p["non_bonded_inter", i, j] = inter p["non_bonded", i, j] = intra + inter @@ -120,7 +136,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): # Electrostatics IF ELECTROSTATICS == 1: cdef np.ndarray coulomb - coulomb = get_obs_contribs(obs.coulomb, size) + coulomb = get_obs_contribs(obs.coulomb, size, calc_sp) for i in range(coulomb.shape[0]): p["coulomb", i] = coulomb[i] p["coulomb"] = np.sum(coulomb, axis=0) @@ -128,7 +144,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): # Dipoles IF DIPOLES == 1: cdef np.ndarray dipolar - dipolar = get_obs_contribs(obs.dipolar, size) + dipolar = get_obs_contribs(obs.dipolar, size, calc_sp) for i in range(dipolar.shape[0]): p["dipolar", i] = dipolar[i] p["dipolar"] = np.sum(dipolar, axis=0) @@ -136,7 +152,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size): # virtual sites IF VIRTUAL_SITES == 1: cdef np.ndarray virtual_sites - virtual_sites = get_obs_contribs(obs.virtual_sites, size) + virtual_sites = get_obs_contribs(obs.virtual_sites, size, calc_sp) for i in range(virtual_sites.shape[0]): p["virtual_sites", i] = virtual_sites[i] p["virtual_sites"] = np.sum(virtual_sites, axis=0) @@ -316,7 +332,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_scalar_pressure, 1) + return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9, True) def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -350,7 +366,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9) + return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9, False) IF DPD == 1: def dpd_stress(self): @@ -401,7 +417,7 @@ class Analysis: analyze.update_energy() handle_errors("calc_long_range_energies failed") - return _Observable_stat_to_dict(analyze.obs_energy, 1) + return _Observable_stat_to_dict(analyze.obs_energy, 1, False) def calc_re(self, chain_start=None, number_of_chains=None, chain_length=None): From e7353221ec845578f6723b7686bc03b9b34ecc0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 16:03:09 +0200 Subject: [PATCH 04/11] core: Rename obs_pressure_tensor to obs_pressure --- src/core/Observable_stat.hpp | 2 +- src/core/pressure.cpp | 23 +++++++++++------------ src/core/pressure_inline.hpp | 15 ++++++--------- src/python/espressomd/analyze.pxd | 2 +- src/python/espressomd/analyze.pyx | 4 ++-- 5 files changed, 21 insertions(+), 25 deletions(-) diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index 39935ce9820..8b9228d43a8 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -29,7 +29,7 @@ #include #include -/** Observable for the scalar pressure, pressure tensor and energy. */ +/** Observable for the pressure and energy. */ class Observable_stat { /** Array for observables on each node. */ std::vector m_data; diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 518560d344a..72c0c95de83 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -38,7 +38,7 @@ #include "electrostatics_magnetostatics/dipole.hpp" /** Pressure tensor of the system */ -Observable_stat obs_pressure_tensor{9}; +Observable_stat obs_pressure{9}; nptiso_struct nptiso = {0.0, 0.0, @@ -58,7 +58,7 @@ nptiso_struct nptiso = {0.0, void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_pressure_long_range(obs_pressure_tensor, particles); + Coulomb::calc_pressure_long_range(obs_pressure, particles); #endif #ifdef DIPOLES /* calculate k-space part of magnetostatic interaction. */ @@ -76,7 +76,7 @@ void pressure_calc() { if (!interactions_sanity_checks()) return; - obs_pressure_tensor = Observable_stat{9}; + obs_pressure = Observable_stat{9}; on_observable_calc(); @@ -93,19 +93,18 @@ void pressure_calc() { calc_long_range_virials(cell_structure.local_particles()); #ifdef VIRTUAL_SITES - if (!obs_pressure_tensor.virtual_sites.empty()) { - auto const vs_pressure_tensor = virtual_sites()->pressure_tensor(); - boost::copy(flatten(vs_pressure_tensor), - obs_pressure_tensor.virtual_sites.begin()); + if (!obs_pressure.virtual_sites.empty()) { + auto const vs_pressure = virtual_sites()->pressure_tensor(); + boost::copy(flatten(vs_pressure), obs_pressure.virtual_sites.begin()); } #endif - obs_pressure_tensor.rescale(volume); + obs_pressure.rescale(volume); /* gather data */ - auto obs_pressure_tensor_res = reduce(comm_cart, obs_pressure_tensor); - if (obs_pressure_tensor_res) { - std::swap(obs_pressure_tensor, *obs_pressure_tensor_res); + auto obs_pressure_res = reduce(comm_cart, obs_pressure); + if (obs_pressure_res) { + std::swap(obs_pressure, *obs_pressure_res); } } @@ -115,7 +114,7 @@ Utils::Vector9d observable_compute_pressure_tensor() { update_pressure(); Utils::Vector9d pressure_tensor{}; for (size_t j = 0; j < 9; j++) { - pressure_tensor[j] = obs_pressure_tensor.accumulate(0, j); + pressure_tensor[j] = obs_pressure.accumulate(0, j); } return pressure_tensor; } diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 39ea96bce25..bb7db4e6817 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -32,7 +32,7 @@ #include -extern Observable_stat obs_pressure_tensor; +extern Observable_stat obs_pressure; /** Calculate non bonded energies between a pair of particles. * @param p1 pointer to particle 1. @@ -51,18 +51,17 @@ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, auto const type1 = p1.p.mol_id; auto const type2 = p2.p.mol_id; - obs_pressure_tensor.add_non_bonded_contribution(type1, type2, - flatten(stress)); + obs_pressure.add_non_bonded_contribution(type1, type2, flatten(stress)); } #ifdef ELECTROSTATICS - if (!obs_pressure_tensor.coulomb.empty()) { + if (!obs_pressure.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++) { - obs_pressure_tensor.coulomb[i * 3 + j] += p_coulomb[i][j]; + obs_pressure.coulomb[i * 3 + j] += p_coulomb[i][j]; } } } @@ -148,8 +147,7 @@ inline bool add_bonded_pressure_tensor(Particle &p1, int bond_id, /* pressure tensor part */ for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) - obs_pressure_tensor.bonded_contribution(bond_id)[k * 3 + l] += - tensor[k][l]; + obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l]; return false; } @@ -167,8 +165,7 @@ inline void add_kinetic_virials(Particle const &p1) { /* kinetic pressure */ for (int k = 0; k < 3; k++) for (int l = 0; l < 3; l++) - obs_pressure_tensor.kinetic[k * 3 + l] += - (p1.m.v[k]) * (p1.m.v[l]) * p1.p.mass; + obs_pressure.kinetic[k * 3 + l] += p1.m.v[k] * p1.m.v[l] * p1.p.mass; } #endif diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index bc8c4821dfa..f7b79a9c3d6 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -82,7 +82,7 @@ cdef extern from "statistics_chain.hpp": array2 calc_rh(int, int, int) cdef extern from "pressure_inline.hpp": - cdef Observable_stat obs_pressure_tensor + cdef Observable_stat obs_pressure cdef extern from "pressure.hpp": cdef void update_pressure() diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index fdf1b2a31a2..3f9581be06e 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -332,7 +332,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9, True) + return _Observable_stat_to_dict(analyze.obs_pressure, 9, True) def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -366,7 +366,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_pressure_tensor, 9, False) + return _Observable_stat_to_dict(analyze.obs_pressure, 9, False) IF DPD == 1: def dpd_stress(self): From 06a7421ad228c577bca752ecca00dd8bac1664d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 17:45:21 +0200 Subject: [PATCH 05/11] core: Make Observable_stat globals local Pass local variable obs_energy and obs_pressure to energy/pressure kernels as mutable references to remove external linkage. --- .../bonded_interactions.dox | 17 +++--- src/core/constraints/ShapeBasedConstraint.cpp | 3 +- src/core/cuda_common_cuda.cu | 3 +- src/core/cuda_interface.hpp | 3 +- src/core/energy.cpp | 11 ++-- src/core/energy.hpp | 10 ++-- src/core/energy_inline.hpp | 49 +++++++---------- src/core/pressure.cpp | 20 ++++--- src/core/pressure.hpp | 14 +++-- src/core/pressure_inline.hpp | 54 ++++++++++--------- src/python/espressomd/analyze.pxd | 8 +-- src/python/espressomd/analyze.pyx | 6 +-- 12 files changed, 98 insertions(+), 100 deletions(-) diff --git a/src/core/bonded_interactions/bonded_interactions.dox b/src/core/bonded_interactions/bonded_interactions.dox index 9bf2a84c2cd..c16c673d649 100644 --- a/src/core/bonded_interactions/bonded_interactions.dox +++ b/src/core/bonded_interactions/bonded_interactions.dox @@ -126,25 +126,20 @@ * result = fene_pair_force(iaparams, dx); * break; * @endcode - * - in function @ref add_bonded_force(): add the new interaction to the - * switch statement if it is not a pairwise bond. For the harmonic angle, - * the code looks like this: + * - in functions called by @ref add_bonded_force(): add the new interaction + * to the switch statements if it is not a pairwise bond. For the harmonic + * angle, the code looks like this: * @code{.cpp} * case BONDED_IA_ANGLE_HARMONIC: - * std::tie(force1, force2, force3) = - * angle_harmonic_force(p1->r.p, p2->r.p, p3->r.p, iaparams); - * break; + * return angle_harmonic_force(p1.r.p, p2.r.p, p3.r.p, iaparams); * @endcode * * In energy_inline.hpp: * - include the header file of the new interaction - * - in function @ref add_bonded_energy(): add the new interaction to the + * - in function @ref calc_bonded_energy(): add the new interaction to the * switch statement. For the FENE bond, e.g., the code looks like this: * @code{.cpp} - * boost::optional retval; - * // ... * case BONDED_IA_FENE: - * retval = fene_pair_energy(iaparams, dx); - * break; + * return fene_pair_energy(iaparams, dx); * @endcode * * Pressure tensor and virial calculation: if your bonded * interaction is not a pair bond or modifies the particles involved, diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 1374d83c07d..bac04384a36 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -113,7 +113,8 @@ ParticleForce ShapeBasedConstraint::force(Particle const &p, void ShapeBasedConstraint::add_energy(const Particle &p, const Utils::Vector3d &folded_pos, - double t, Observable_stat &energy) const { + double t, + Observable_stat &obs_energy) const { double dist; double nonbonded_en = 0.0; diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index eece331536f..3847dc55661 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -221,8 +221,7 @@ void clear_energy_on_GPU() { cuda_safe_mem(cudaMemset(energy_device, 0, sizeof(CUDA_energy))); } -void copy_energy_from_GPU() { - extern Observable_stat obs_energy; +void copy_energy_from_GPU(Observable_stat &obs_energy) { if (!global_part_vars_host.communication_enabled) return; cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy), diff --git a/src/core/cuda_interface.hpp b/src/core/cuda_interface.hpp index 842f1f95abc..9ba41bb5d10 100644 --- a/src/core/cuda_interface.hpp +++ b/src/core/cuda_interface.hpp @@ -24,6 +24,7 @@ #ifdef CUDA #include "CudaHostAllocator.hpp" +#include "Observable_stat.hpp" #include "ParticleRange.hpp" #include @@ -107,7 +108,7 @@ typedef struct { } CUDA_global_part_vars; void copy_forces_from_GPU(ParticleRange &particles); -void copy_energy_from_GPU(); +void copy_energy_from_GPU(Observable_stat &obs_energy); void clear_energy_on_GPU(); CUDA_global_part_vars *gpu_get_global_particle_vars_pointer_host(); diff --git a/src/core/energy.cpp b/src/core/energy.cpp index 580884f57c0..bacecfb1c11 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -43,6 +43,8 @@ ActorList energyActors; /** Energy of the system */ Observable_stat obs_energy{1}; +Observable_stat const &get_obs_energy() { return obs_energy; } + /** Reduce the system energy from all MPI ranks. */ void master_energy_calc() { mpi_gather_stats(GatherStats::energy); } @@ -65,13 +67,14 @@ void energy_calc(const double time) { on_observable_calc(); for (auto const &p : cell_structure.local_particles()) { - add_kinetic_energy(p); + add_kinetic_energy(p, obs_energy); } short_range_loop( - [](Particle &p) { add_single_particle_energy(p); }, + [](Particle &p) { add_bonded_energy(p, obs_energy); }, [](Particle const &p1, Particle const &p2, Distance const &d) { - add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2); + add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, + obs_energy); }); calc_long_range_energies(cell_structure.local_particles()); @@ -80,7 +83,7 @@ void energy_calc(const double time) { Constraints::constraints.add_energy(local_parts, time, obs_energy); #ifdef CUDA - copy_energy_from_GPU(); + copy_energy_from_GPU(obs_energy); #endif /* gather data */ diff --git a/src/core/energy.hpp b/src/core/energy.hpp index 8daa5d745ab..00c439118ac 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -27,17 +27,21 @@ #ifndef _ENERGY_H #define _ENERGY_H +#include "Observable_stat.hpp" #include "ParticleRange.hpp" #include "actor/ActorList.hpp" extern ActorList energyActors; -/** Calculate energies. */ -void update_energy(); - /** Parallel energy calculation. */ void energy_calc(double time); +/** Run @ref energy_calc in parallel. */ +void update_energy(); + +/** Return the energy observable. */ +Observable_stat const &get_obs_energy(); + /** Calculate long-range energies (P3M, ...). */ void calc_long_range_energies(const ParticleRange &particles); diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index 07f7c66ab98..e494542cd23 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -71,8 +71,6 @@ #include "electrostatics_magnetostatics/dipole_inline.hpp" #endif -extern Observable_stat obs_energy; - /** Calculate non-bonded energies between a pair of particles. * @param p1 particle 1. * @param p2 particle 2. @@ -179,10 +177,12 @@ inline double calc_non_bonded_pair_energy(Particle const &p1, * @param d vector between p1 and p2. * @param dist distance between p1 and p2. * @param dist2 distance squared between p1 and p2. + * @param[in,out] obs_energy energy observable. */ inline void add_non_bonded_pair_energy(Particle const &p1, Particle const &p2, Utils::Vector3d const &d, - double const dist, double const dist2) { + double const dist, double const dist2, + Observable_stat &obs_energy) { IA_parameters const &ia_params = *get_ia_param(p1.p.type, p2.p.type); #ifdef EXCLUSIONS @@ -279,32 +279,12 @@ calc_bonded_energy(Bonded_ia_parameters const &iaparams, Particle const &p1, throw BondInvalidSizeError(n_partners); } -/** Add bonded energies for one particle to the energy observable. - * @param[in] p1 particle for which to calculate energies - * @param[in] bond_id Numeric id of the bond - * @param[in] partners Bond partners of particle. - * - * @return True if bond was broken, false otherwise. - */ -inline bool add_bonded_energy(Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - auto const result = calc_bonded_energy(iaparams, p1, partners); - - if (result) { - obs_energy.bonded_contribution(bond_id)[0] += result.get(); - - return false; - } - - return true; -} - /** Add kinetic energies for one particle to the energy observable. * @param[in] p1 particle for which to calculate energies + * @param[out] obs_energy energy observable */ -inline void add_kinetic_energy(Particle const &p1) { +inline void add_kinetic_energy(Particle const &p1, + Observable_stat &obs_energy) { if (p1.p.is_virtual) return; @@ -323,11 +303,22 @@ inline void add_kinetic_energy(Particle const &p1) { #endif } -/** Add kinetic and bonded energies for one particle to the energy observable. +/** Add bonded energies for one particle to the energy observable. * @param[in] p particle for which to calculate energies + * @param[out] obs_energy energy observable */ -inline void add_single_particle_energy(Particle &p) { - cell_structure.execute_bond_handler(p, add_bonded_energy); +inline void add_bonded_energy(Particle &p, Observable_stat &obs_energy) { + cell_structure.execute_bond_handler( + p, [&obs_energy](Particle &p1, int bond_id, + Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + auto const result = calc_bonded_energy(iaparams, p1, partners); + if (result) { + obs_energy.bonded_contribution(bond_id)[0] += result.get(); + return false; + } + return true; + }); } #endif // ENERGY_INLINE_HPP diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 72c0c95de83..96d19c435c2 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -37,9 +37,6 @@ #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" -/** Pressure tensor of the system */ -Observable_stat obs_pressure{9}; - nptiso_struct nptiso = {0.0, 0.0, 0.0, @@ -54,6 +51,11 @@ nptiso_struct nptiso = {0.0, false, 0}; +/** Pressure tensor of the system */ +Observable_stat obs_pressure{9}; + +Observable_stat const &get_obs_pressure() { return obs_pressure; } + /** Calculate long-range virials (P3M, ...). */ void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS @@ -66,10 +68,6 @@ void calc_long_range_virials(const ParticleRange &particles) { #endif } -static void add_single_particle_virials(Particle &p) { - cell_structure.execute_bond_handler(p, add_bonded_pressure_tensor); -} - void pressure_calc() { auto const volume = box_geo.volume(); @@ -81,13 +79,13 @@ void pressure_calc() { on_observable_calc(); for (auto const &p : cell_structure.local_particles()) { - add_kinetic_virials(p); + add_kinetic_virials(p, obs_pressure); } - short_range_loop([](Particle &p) { add_single_particle_virials(p); }, + short_range_loop([](Particle &p) { add_bonded_virials(p, obs_pressure); }, [](Particle &p1, Particle &p2, Distance const &d) { - add_non_bonded_pair_virials(p1, p2, d.vec21, - sqrt(d.dist2)); + add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), + obs_pressure); }); calc_long_range_virials(cell_structure.local_particles()); diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp index 165bf5d3ee2..9b08b22b340 100644 --- a/src/core/pressure.hpp +++ b/src/core/pressure.hpp @@ -25,16 +25,20 @@ #ifndef CORE_PRESSURE_HPP #define CORE_PRESSURE_HPP +#include "Observable_stat.hpp" + #include -/** Parallel pressure calculation from a virial expansion. - */ +/** Parallel pressure calculation from a virial expansion. */ void pressure_calc(); +/** Run @ref pressure_calc in parallel. */ +void update_pressure(); + +/** Return the pressure observable. */ +Observable_stat const &get_obs_pressure(); + /** Helper function for @ref Observables::PressureTensor. */ Utils::Vector9d observable_compute_pressure_tensor(); -/** Calculate the scalar pressure and pressure tensor. */ -void update_pressure(); - #endif diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index bb7db4e6817..3b4f7050c01 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -32,16 +32,16 @@ #include -extern Observable_stat obs_pressure; - /** Calculate non bonded energies between a pair of particles. * @param p1 pointer to particle 1. * @param p2 pointer to particle 2. * @param d vector between p1 and p2. * @param dist distance between p1 and p2. + * @param[in,out] obs_pressure pressure observable. */ inline void add_non_bonded_pair_virials(Particle const &p1, Particle const &p2, - Utils::Vector3d const &d, double dist) { + Utils::Vector3d const &d, double dist, + Observable_stat &obs_pressure) { #ifdef EXCLUSIONS if (do_nonbonded(p1, p2)) #endif @@ -136,29 +136,12 @@ calc_bonded_pressure_tensor(Bonded_ia_parameters const &iaparams, Particle &p1, } } -inline bool add_bonded_pressure_tensor(Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = bonded_ia_params[bond_id]; - - auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners); - if (result) { - auto const &tensor = result.get(); - - /* pressure tensor part */ - for (int k = 0; k < 3; k++) - for (int l = 0; l < 3; l++) - obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l]; - - return false; - } - - return true; -} - /** Calculate kinetic pressure (aka energy) for one particle. - * @param p1 particle for which to calculate pressure + * @param[in] p1 particle for which to calculate pressure + * @param[out] obs_pressure pressure observable */ -inline void add_kinetic_virials(Particle const &p1) { +inline void add_kinetic_virials(Particle const &p1, + Observable_stat &obs_pressure) { if (p1.p.is_virtual) return; @@ -168,4 +151,27 @@ inline void add_kinetic_virials(Particle const &p1) { obs_pressure.kinetic[k * 3 + l] += p1.m.v[k] * p1.m.v[l] * p1.p.mass; } +/** Add bonded energies for one particle to the energy observable. + * @param[in] p particle for which to calculate pressure + * @param[out] obs_pressure pressure observable + */ +inline void add_bonded_virials(Particle &p, Observable_stat &obs_pressure) { + cell_structure.execute_bond_handler(p, [&obs_pressure]( + Particle &p1, int bond_id, + Utils::Span partners) { + auto const &iaparams = bonded_ia_params[bond_id]; + auto const result = calc_bonded_pressure_tensor(iaparams, p1, partners); + if (result) { + auto const &tensor = result.get(); + /* pressure tensor part */ + for (int k = 0; k < 3; k++) + for (int l = 0; l < 3; l++) + obs_pressure.bonded_contribution(bond_id)[k * 3 + l] += tensor[k][l]; + + return false; + } + return true; + }); +} + #endif diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index f7b79a9c3d6..b2d6f967bf5 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -81,17 +81,13 @@ cdef extern from "statistics_chain.hpp": array4 calc_rg(int, int, int) except + array2 calc_rh(int, int, int) -cdef extern from "pressure_inline.hpp": - cdef Observable_stat obs_pressure - cdef extern from "pressure.hpp": cdef void update_pressure() - -cdef extern from "energy_inline.hpp": - cdef Observable_stat obs_energy + cdef const Observable_stat & get_obs_pressure() cdef extern from "energy.hpp": cdef void update_energy() + cdef const Observable_stat & get_obs_energy() double calculate_current_potential_energy_of_system() cdef extern from "dpd.hpp": diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 3f9581be06e..edab9cd37a3 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -332,7 +332,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_pressure, 9, True) + return _Observable_stat_to_dict(analyze.get_obs_pressure(), 9, True) def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -366,7 +366,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.obs_pressure, 9, False) + return _Observable_stat_to_dict(analyze.get_obs_pressure(), 9, False) IF DPD == 1: def dpd_stress(self): @@ -417,7 +417,7 @@ class Analysis: analyze.update_energy() handle_errors("calc_long_range_energies failed") - return _Observable_stat_to_dict(analyze.obs_energy, 1, False) + return _Observable_stat_to_dict(analyze.get_obs_energy(), 1, False) def calc_re(self, chain_start=None, number_of_chains=None, chain_length=None): From acdfc60650bf2ba6721462f126129756d11d4951 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 21:40:55 +0200 Subject: [PATCH 06/11] core: Add ELC correction to the P3M energy --- src/core/Observable_stat.cpp | 2 +- src/core/electrostatics_magnetostatics/coulomb.cpp | 2 +- src/python/espressomd/analyze.pyx | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 881864f21ee..4bb3a3dc0aa 100644 --- a/src/core/Observable_stat.cpp +++ b/src/core/Observable_stat.cpp @@ -27,7 +27,7 @@ Observable_stat::Observable_stat(size_t chunk_size) : m_chunk_size(chunk_size) { // number of chunks for different interaction types - auto constexpr n_coulomb = 3; + auto constexpr n_coulomb = 2; auto constexpr n_dipolar = 3; #ifdef VIRTUAL_SITES auto constexpr n_vs = 1; diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 983dea0ebde..9594f11b8b2 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -335,7 +335,7 @@ void calc_energy_long_range(Observable_stat &energy, // restore modified sums ELC_P3M_restore_p3m_sums(particles); } - energy.coulomb[2] = ELC_energy(particles); + energy.coulomb[1] += ELC_energy(particles); break; #endif #ifdef SCAFACOS diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index edab9cd37a3..7e4bf02f577 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -64,7 +64,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size, * ``"nonbonded_intra", , ``: nonbonded contribution between short ranged forces between type i and j and with the same mol_id * ``"nonbonded_inter", , ``: nonbonded contribution between short ranged forces between type i and j and different mol_ids * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method - * ``"coulomb", ``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"coulomb", ``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: dipolar contribution, how it is calculated depends on the method * ``"dipolar", ``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) and their corrections (``i=2``) * ``"virtual_sites"``: virtual sites contribution @@ -323,7 +323,7 @@ class Analysis: * ``"coulomb"``: Coulomb pressure, how it is calculated depends on the method. It is equivalent to 1/3 of the trace of the Coulomb pressure tensor. For how the pressure tensor is calculated, see :ref:`Pressure Tensor`. The averaged value in an isotropic NVT simulation is equivalent to the average of :math:`E^{coulomb}/(3V)`, see :cite:`brown95a`. - * ``"coulomb", ``: Coulomb pressure from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"coulomb", ``: Coulomb pressure from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: not implemented * ``"virtual_sites"``: Pressure contribution from virtual sites @@ -357,7 +357,7 @@ class Analysis: * ``"nonbonded_intra", , ``: nonbonded pressure tensor between short ranged forces between type i and j and with the same mol_id * ``"nonbonded_inter", , ``: nonbonded pressure tensor between short ranged forces between type i and j and different mol_ids * ``"coulomb"``: Maxwell pressure tensor, how it is calculated depends on the method - * ``"coulomb", ``: Maxwell pressure tensor from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"coulomb", ``: Maxwell pressure tensor from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: not implemented * ``"virtual_sites"``: pressure tensor contribution from virtual sites @@ -398,7 +398,7 @@ class Analysis: * ``"nonbonded_intra", , ``: nonbonded energy between short ranged forces between type i and j and with the same mol_id * ``"nonbonded_inter", , ``: nonbonded energy between short ranged forces between type i and j and different mol_ids * ``"coulomb"``: Coulomb energy, how it is calculated depends on the method - * ``"coulomb", ``: Coulomb energy from particle pairs (``i=0``), electrostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"coulomb", ``: Coulomb energy from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: dipolar energy * ``"dipolar", ``: dipolar energy from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) and their corrections (``i=2``) From 6d1edb3f230a068b5d7d8b304473c21b78df6e83 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 23 Jun 2020 22:01:36 +0200 Subject: [PATCH 07/11] core: Add DLC correction to the DP3M/DDS energies --- src/core/Observable_stat.cpp | 2 +- src/core/electrostatics_magnetostatics/dipole.cpp | 4 ++-- src/python/espressomd/analyze.pyx | 4 ++-- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 4bb3a3dc0aa..a40377367b9 100644 --- a/src/core/Observable_stat.cpp +++ b/src/core/Observable_stat.cpp @@ -28,7 +28,7 @@ Observable_stat::Observable_stat(size_t chunk_size) : m_chunk_size(chunk_size) { // number of chunks for different interaction types auto constexpr n_coulomb = 2; - auto constexpr n_dipolar = 3; + auto constexpr n_dipolar = 2; #ifdef VIRTUAL_SITES auto constexpr n_vs = 1; #else diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index 42a5858c9db..5e88476da95 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -216,7 +216,7 @@ void calc_energy_long_range(Observable_stat &energy, case DIPOLAR_MDLC_P3M: dp3m_dipole_assign(particles); energy.dipolar[1] = dp3m_calc_kspace_forces(false, true, particles); - energy.dipolar[2] = add_mdlc_energy_corrections(particles); + energy.dipolar[1] += add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA: @@ -226,7 +226,7 @@ void calc_energy_long_range(Observable_stat &energy, case DIPOLAR_MDLC_DS: energy.dipolar[1] = magnetic_dipolar_direct_sum_calculations(false, true, particles); - energy.dipolar[2] = add_mdlc_energy_corrections(particles); + energy.dipolar[1] += add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_DS: diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 7e4bf02f577..f37e222cb07 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -66,7 +66,7 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size, * ``"coulomb"``: Coulomb contribution, how it is calculated depends on the method * ``"coulomb", ``: Coulomb contribution from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: dipolar contribution, how it is calculated depends on the method - * ``"dipolar", ``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"dipolar", ``: dipolar contribution from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) * ``"virtual_sites"``: virtual sites contribution * ``"virtual_sites", ``: contribution from virtual site i @@ -400,7 +400,7 @@ class Analysis: * ``"coulomb"``: Coulomb energy, how it is calculated depends on the method * ``"coulomb", ``: Coulomb energy from particle pairs (``i=0``), electrostatics solvers (``i=1``) * ``"dipolar"``: dipolar energy - * ``"dipolar", ``: dipolar energy from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) and their corrections (``i=2``) + * ``"dipolar", ``: dipolar energy from particle pairs and magnetic field constraints (``i=0``), magnetostatics solvers (``i=1``) Examples From 01b4003796ed1c14f1f44439d1ce9799a473c984 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Jun 2020 11:01:09 +0200 Subject: [PATCH 08/11] core: Store P3M GPU in the same slot as P3M CPU Slot 0: particle electrostatics. Slot 1: electrostatics solver. --- src/core/cuda_common_cuda.cu | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index 3847dc55661..a0ebe72a90e 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -227,7 +227,7 @@ void copy_energy_from_GPU(Observable_stat &obs_energy) { cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy), cudaMemcpyDeviceToHost)); if (!obs_energy.coulomb.empty()) - obs_energy.coulomb[0] += energy_host.coulomb; + obs_energy.coulomb[1] += energy_host.coulomb; if (obs_energy.dipolar.size() >= 2) obs_energy.dipolar[1] += energy_host.dipolar; } From f0de4eea2ef2c0c39727106ddb7942e882711b9b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Jun 2020 12:12:50 +0200 Subject: [PATCH 09/11] core: Narrow inclusion of Observable_stat --- src/core/cuda_common_cuda.cu | 10 +++----- src/core/cuda_interface.hpp | 3 +-- .../electrostatics_magnetostatics/coulomb.cpp | 23 ++++++++++--------- .../electrostatics_magnetostatics/coulomb.hpp | 3 +-- .../electrostatics_magnetostatics/dipole.cpp | 23 +++++++++---------- .../electrostatics_magnetostatics/dipole.hpp | 3 +-- src/core/energy.cpp | 17 +++++++++----- src/core/energy_inline.hpp | 14 +++++------ 8 files changed, 46 insertions(+), 50 deletions(-) diff --git a/src/core/cuda_common_cuda.cu b/src/core/cuda_common_cuda.cu index a0ebe72a90e..02b36554453 100644 --- a/src/core/cuda_common_cuda.cu +++ b/src/core/cuda_common_cuda.cu @@ -22,7 +22,6 @@ #include "debug.hpp" -#include "Observable_stat.hpp" #include "ParticleRange.hpp" #include "cuda_init.hpp" #include "cuda_interface.hpp" @@ -221,15 +220,12 @@ void clear_energy_on_GPU() { cuda_safe_mem(cudaMemset(energy_device, 0, sizeof(CUDA_energy))); } -void copy_energy_from_GPU(Observable_stat &obs_energy) { +CUDA_energy copy_energy_from_GPU() { if (!global_part_vars_host.communication_enabled) - return; + return {}; cuda_safe_mem(cudaMemcpy(&energy_host, energy_device, sizeof(CUDA_energy), cudaMemcpyDeviceToHost)); - if (!obs_energy.coulomb.empty()) - obs_energy.coulomb[1] += energy_host.coulomb; - if (obs_energy.dipolar.size() >= 2) - obs_energy.dipolar[1] += energy_host.dipolar; + return energy_host; } void _cuda_safe_mem(cudaError_t CU_err, const char *file, unsigned int line) { diff --git a/src/core/cuda_interface.hpp b/src/core/cuda_interface.hpp index 9ba41bb5d10..6bb0d710a5d 100644 --- a/src/core/cuda_interface.hpp +++ b/src/core/cuda_interface.hpp @@ -24,7 +24,6 @@ #ifdef CUDA #include "CudaHostAllocator.hpp" -#include "Observable_stat.hpp" #include "ParticleRange.hpp" #include @@ -108,7 +107,7 @@ typedef struct { } CUDA_global_part_vars; void copy_forces_from_GPU(ParticleRange &particles); -void copy_energy_from_GPU(Observable_stat &obs_energy); +CUDA_energy copy_energy_from_GPU(); 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 9594f11b8b2..314851b87e8 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -297,8 +297,8 @@ void calc_long_range_force(const ParticleRange &particles) { #endif } -void calc_energy_long_range(Observable_stat &energy, - const ParticleRange &particles) { +double calc_energy_long_range(const ParticleRange &particles) { + double energy = 0.0; switch (coulomb.method) { #ifdef P3M case COULOMB_P3M_GPU: @@ -307,46 +307,47 @@ void calc_energy_long_range(Observable_stat &energy, break; case COULOMB_P3M: p3m_charge_assign(particles); - energy.coulomb[1] = p3m_calc_kspace_forces(false, true, particles); + energy = p3m_calc_kspace_forces(false, true, particles); break; case COULOMB_ELC_P3M: // assign the original charges first // they may not have been assigned yet p3m_charge_assign(particles); if (!elc_params.dielectric_contrast_on) - energy.coulomb[1] = p3m_calc_kspace_forces(false, true, particles); + energy = p3m_calc_kspace_forces(false, true, particles); else { - energy.coulomb[1] = 0.5 * p3m_calc_kspace_forces(false, true, particles); - energy.coulomb[1] += 0.5 * coulomb.prefactor * - ELC_P3M_dielectric_layers_energy_self(particles); + energy = 0.5 * p3m_calc_kspace_forces(false, true, particles); + energy += 0.5 * coulomb.prefactor * + ELC_P3M_dielectric_layers_energy_self(particles); // assign both original and image charges now ELC_p3m_charge_assign_both(particles); ELC_P3M_modify_p3m_sums_both(particles); - energy.coulomb[1] += 0.5 * p3m_calc_kspace_forces(false, true, particles); + energy += 0.5 * p3m_calc_kspace_forces(false, true, particles); // assign only the image charges now ELC_p3m_charge_assign_image(particles); ELC_P3M_modify_p3m_sums_image(particles); - energy.coulomb[1] -= 0.5 * p3m_calc_kspace_forces(false, true, particles); + energy -= 0.5 * p3m_calc_kspace_forces(false, true, particles); // restore modified sums ELC_P3M_restore_p3m_sums(particles); } - energy.coulomb[1] += ELC_energy(particles); + energy += ELC_energy(particles); break; #endif #ifdef SCAFACOS case COULOMB_SCAFACOS: assert(!Scafacos::dipolar()); - energy.coulomb[1] += Scafacos::long_range_energy(); + energy += Scafacos::long_range_energy(); break; #endif default: break; } + return energy; } int iccp3m_sanity_check() { diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index 3f6d8f677e4..b169ade04b5 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -74,8 +74,7 @@ void init(); void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy, - const ParticleRange &particles); +double calc_energy_long_range(const ParticleRange &particles); int iccp3m_sanity_check(); diff --git a/src/core/electrostatics_magnetostatics/dipole.cpp b/src/core/electrostatics_magnetostatics/dipole.cpp index 5e88476da95..d52d70251b4 100644 --- a/src/core/electrostatics_magnetostatics/dipole.cpp +++ b/src/core/electrostatics_magnetostatics/dipole.cpp @@ -205,33 +205,31 @@ void calc_long_range_force(const ParticleRange &particles) { } } -void calc_energy_long_range(Observable_stat &energy, - const ParticleRange &particles) { +double calc_energy_long_range(const ParticleRange &particles) { + double energy = 0.; switch (dipole.method) { #ifdef DP3M case DIPOLAR_P3M: dp3m_dipole_assign(particles); - energy.dipolar[1] = dp3m_calc_kspace_forces(false, true, particles); + energy = dp3m_calc_kspace_forces(false, true, particles); break; case DIPOLAR_MDLC_P3M: dp3m_dipole_assign(particles); - energy.dipolar[1] = dp3m_calc_kspace_forces(false, true, particles); - energy.dipolar[1] += add_mdlc_energy_corrections(particles); + energy = dp3m_calc_kspace_forces(false, true, particles); + energy += add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_ALL_WITH_ALL_AND_NO_REPLICA: - energy.dipolar[1] = dawaanr_calculations(false, true, particles); + energy = dawaanr_calculations(false, true, particles); break; #ifdef DP3M case DIPOLAR_MDLC_DS: - energy.dipolar[1] = - magnetic_dipolar_direct_sum_calculations(false, true, particles); - energy.dipolar[1] += add_mdlc_energy_corrections(particles); + energy = magnetic_dipolar_direct_sum_calculations(false, true, particles); + energy += add_mdlc_energy_corrections(particles); break; #endif case DIPOLAR_DS: - energy.dipolar[1] = - magnetic_dipolar_direct_sum_calculations(false, true, particles); + energy = magnetic_dipolar_direct_sum_calculations(false, true, particles); break; case DIPOLAR_DS_GPU: // NOLINT(bugprone-branch-clone) // do nothing: it's an actor @@ -244,7 +242,7 @@ void calc_energy_long_range(Observable_stat &energy, #ifdef SCAFACOS_DIPOLES case DIPOLAR_SCAFACOS: assert(Scafacos::dipolar()); - energy.dipolar[1] = Scafacos::long_range_energy(); + energy = Scafacos::long_range_energy(); #endif case DIPOLAR_NONE: break; @@ -253,6 +251,7 @@ void calc_energy_long_range(Observable_stat &energy, << "energy calculation not implemented for dipolar method."; break; } + return energy; } int set_mesh() { diff --git a/src/core/electrostatics_magnetostatics/dipole.hpp b/src/core/electrostatics_magnetostatics/dipole.hpp index d306fc4c389..47b2f4d8ddc 100644 --- a/src/core/electrostatics_magnetostatics/dipole.hpp +++ b/src/core/electrostatics_magnetostatics/dipole.hpp @@ -80,8 +80,7 @@ void init(); void calc_long_range_force(const ParticleRange &particles); -void calc_energy_long_range(Observable_stat &energy, - const ParticleRange &particles); +double calc_energy_long_range(const ParticleRange &particles); int set_mesh(); void bcast_params(const boost::mpi::communicator &comm); diff --git a/src/core/energy.cpp b/src/core/energy.cpp index bacecfb1c11..af3e3485e15 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -67,7 +67,7 @@ void energy_calc(const double time) { on_observable_calc(); for (auto const &p : cell_structure.local_particles()) { - add_kinetic_energy(p, obs_energy); + obs_energy.kinetic[0] += calc_kinetic_energy(p); } short_range_loop( @@ -83,7 +83,11 @@ void energy_calc(const double time) { Constraints::constraints.add_energy(local_parts, time, obs_energy); #ifdef CUDA - copy_energy_from_GPU(obs_energy); + auto const energy_host = copy_energy_from_GPU(); + if (!obs_energy.coulomb.empty()) + obs_energy.coulomb[1] += energy_host.coulomb; + if (!obs_energy.dipolar.empty()) + obs_energy.dipolar[1] += energy_host.dipolar; #endif /* gather data */ @@ -98,12 +102,13 @@ void update_energy() { 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(obs_energy, particles); -#endif /* ifdef ELECTROSTATICS */ + obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(particles); +#endif #ifdef DIPOLES - Dipole::calc_energy_long_range(obs_energy, particles); -#endif /* ifdef DIPOLES */ + /* calculate k-space part of magnetostatic interaction. */ + obs_energy.dipolar[1] = Dipole::calc_energy_long_range(particles); +#endif } double calculate_current_potential_energy_of_system() { diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index e494542cd23..26242f87633 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -279,17 +279,15 @@ calc_bonded_energy(Bonded_ia_parameters const &iaparams, Particle const &p1, throw BondInvalidSizeError(n_partners); } -/** Add kinetic energies for one particle to the energy observable. +/** Calculate kinetic energies for one particle. * @param[in] p1 particle for which to calculate energies - * @param[out] obs_energy energy observable */ -inline void add_kinetic_energy(Particle const &p1, - Observable_stat &obs_energy) { +inline double calc_kinetic_energy(Particle const &p1) { if (p1.p.is_virtual) - return; + return 0.0; /* kinetic energy */ - obs_energy.kinetic[0] += 0.5 * p1.p.mass * p1.m.v.norm2(); + auto res = 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 @@ -297,10 +295,10 @@ inline void add_kinetic_energy(Particle const &p1, if (p1.p.rotation) { /* the rotational part is added to the total kinetic energy; * Here we use the rotational inertia */ - obs_energy.kinetic[0] += - 0.5 * (hadamard_product(p1.m.omega, p1.m.omega) * p1.p.rinertia); + res += 0.5 * (hadamard_product(p1.m.omega, p1.m.omega) * p1.p.rinertia); } #endif + return res; } /** Add bonded energies for one particle to the energy observable. From 6d2647c14a4d6317501142193c53dbc0774ec0c6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Jun 2020 12:45:04 +0200 Subject: [PATCH 10/11] python: Simplify function argument list --- src/python/espressomd/analyze.pxd | 1 + src/python/espressomd/analyze.pyx | 13 ++++++------- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index b2d6f967bf5..80966c4c022 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -58,6 +58,7 @@ cdef extern from "Observable_stat.hpp": Span[double] bonded_contribution(int bond_id) Span[double] non_bonded_intra_contribution(int type1, int type2) Span[double] non_bonded_inter_contribution(int type1, int type2) + size_t chunk_size() cdef extern from "statistics.hpp": cdef vector[double] calc_structurefactor(PartCfg & , const vector[int] & p_types, int order) diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index f37e222cb07..d1f80e9c9ce 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -36,16 +36,13 @@ from .utils cimport create_nparray_from_double_array from .particle_data cimport get_n_part -cdef _Observable_stat_to_dict(Observable_stat obs, int size, - cbool calc_sp): +cdef _Observable_stat_to_dict(Observable_stat obs, cbool calc_sp): """Transform an Observable_stat struct to a python dict. Parameters ---------- obs : Core observable. - size : :obj:`int`, \{1, 9\} - Dimensionality of the data. calc_sp : :obj:`bool` Whether to calculate a scalar pressure (only relevant when ``obs`` is a pressure tensor observable). @@ -72,6 +69,8 @@ cdef _Observable_stat_to_dict(Observable_stat obs, int size, """ + size = obs.chunk_size() + def set_initial(): if size == 9 and not calc_sp: return np.zeros((3, 3), dtype=float) @@ -332,7 +331,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.get_obs_pressure(), 9, True) + return _Observable_stat_to_dict(analyze.get_obs_pressure(), True) def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -366,7 +365,7 @@ class Analysis: # Update in ESPResSo core analyze.update_pressure() - return _Observable_stat_to_dict(analyze.get_obs_pressure(), 9, False) + return _Observable_stat_to_dict(analyze.get_obs_pressure(), False) IF DPD == 1: def dpd_stress(self): @@ -417,7 +416,7 @@ class Analysis: analyze.update_energy() handle_errors("calc_long_range_energies failed") - return _Observable_stat_to_dict(analyze.get_obs_energy(), 1, False) + return _Observable_stat_to_dict(analyze.get_obs_energy(), False) def calc_re(self, chain_start=None, number_of_chains=None, chain_length=None): From 88b63bd0cd969dffa729233d5c8c9802107811d0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 24 Jun 2020 16:56:47 +0200 Subject: [PATCH 11/11] core: Narrow inclusion of Observable_stat --- src/core/electrostatics_magnetostatics/coulomb.cpp | 9 +++------ src/core/electrostatics_magnetostatics/coulomb.hpp | 4 +--- src/core/electrostatics_magnetostatics/dipole.hpp | 2 -- src/core/pressure.cpp | 3 ++- 4 files changed, 6 insertions(+), 12 deletions(-) diff --git a/src/core/electrostatics_magnetostatics/coulomb.cpp b/src/core/electrostatics_magnetostatics/coulomb.cpp index 314851b87e8..41f5f3122dd 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.cpp +++ b/src/core/electrostatics_magnetostatics/coulomb.cpp @@ -42,8 +42,7 @@ Coulomb_parameters coulomb; namespace Coulomb { -void calc_pressure_long_range(Observable_stat &p_tensor, - const ParticleRange &particles) { +Utils::Vector9d calc_pressure_long_range(const ParticleRange &particles) { switch (coulomb.method) { #ifdef P3M case COULOMB_ELC_P3M: @@ -57,10 +56,7 @@ void calc_pressure_long_range(Observable_stat &p_tensor, break; case COULOMB_P3M: { p3m_charge_assign(particles); - auto const p3m_p_tensor = p3m_calc_kspace_pressure_tensor(); - std::copy_n(p3m_p_tensor.data(), 9, p_tensor.coulomb.begin() + 9); - - break; + return p3m_calc_kspace_pressure_tensor(); } #endif case COULOMB_MMM1D: @@ -72,6 +68,7 @@ void calc_pressure_long_range(Observable_stat &p_tensor, default: break; } + return {}; } void sanity_checks(int &state) { diff --git a/src/core/electrostatics_magnetostatics/coulomb.hpp b/src/core/electrostatics_magnetostatics/coulomb.hpp index b169ade04b5..0b83db00409 100644 --- a/src/core/electrostatics_magnetostatics/coulomb.hpp +++ b/src/core/electrostatics_magnetostatics/coulomb.hpp @@ -24,7 +24,6 @@ #include #ifdef ELECTROSTATICS -#include "Observable_stat.hpp" #include #include @@ -60,8 +59,7 @@ struct Coulomb_parameters { extern Coulomb_parameters coulomb; namespace Coulomb { -void calc_pressure_long_range(Observable_stat &p_tensor, - const ParticleRange &particles); +Utils::Vector9d calc_pressure_long_range(const ParticleRange &particles); void sanity_checks(int &state); double cutoff(const Utils::Vector3d &box_l); diff --git a/src/core/electrostatics_magnetostatics/dipole.hpp b/src/core/electrostatics_magnetostatics/dipole.hpp index 47b2f4d8ddc..7f557057bbc 100644 --- a/src/core/electrostatics_magnetostatics/dipole.hpp +++ b/src/core/electrostatics_magnetostatics/dipole.hpp @@ -24,8 +24,6 @@ #include #ifdef DIPOLES -#include "Observable_stat.hpp" - #include #include diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 96d19c435c2..db90822ff55 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -60,7 +60,8 @@ Observable_stat const &get_obs_pressure() { return obs_pressure; } void calc_long_range_virials(const ParticleRange &particles) { #ifdef ELECTROSTATICS /* calculate k-space part of electrostatic interaction. */ - Coulomb::calc_pressure_long_range(obs_pressure, particles); + auto const coulomb_pressure = Coulomb::calc_pressure_long_range(particles); + boost::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9); #endif #ifdef DIPOLES /* calculate k-space part of magnetostatic interaction. */