From d0ceeb1ffeaefe4066a2ab947940c0113445d286 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 6 Aug 2021 14:56:29 +0200 Subject: [PATCH 01/11] core: Make short-range loop kernels const --- src/core/energy.cpp | 2 +- src/core/pressure.cpp | 4 ++-- src/core/pressure_inline.hpp | 3 ++- 3 files changed, 5 insertions(+), 4 deletions(-) diff --git a/src/core/energy.cpp b/src/core/energy.cpp index 0c4e20d495a..aacc8939a21 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -70,7 +70,7 @@ void energy_calc(const double time) { } short_range_loop( - [](Particle &p1, int bond_id, Utils::Span partners) { + [](Particle const &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) { diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index cf6822297b6..9380dbd24e1 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -85,7 +85,7 @@ void pressure_calc() { } short_range_loop( - [](Particle &p1, int bond_id, Utils::Span partners) { + [](Particle const &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) { @@ -100,7 +100,7 @@ void pressure_calc() { } return true; }, - [](Particle &p1, Particle &p2, Distance const &d) { + [](Particle const &p1, Particle const &p2, Distance const &d) { add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), obs_pressure); }, diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index 54349d0cc79..eef2ff4718f 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -138,7 +138,8 @@ calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams, } inline boost::optional> -calc_bonded_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle &p1, +calc_bonded_pressure_tensor(Bonded_IA_Parameters const &iaparams, + Particle const &p1, Utils::Span partners) { switch (number_of_partners(iaparams)) { case 1: From 04f0de31a8529770206d1b32340c9a05e3a8084b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 6 Aug 2021 16:02:48 +0200 Subject: [PATCH 02/11] core: Remove intermediate Observable_stat functions --- src/core/energy.cpp | 27 ++++++++++++--------------- src/core/energy.hpp | 4 ---- src/core/pressure.cpp | 29 +++++++++++++---------------- 3 files changed, 25 insertions(+), 35 deletions(-) diff --git a/src/core/energy.cpp b/src/core/energy.cpp index aacc8939a21..d76d1add83e 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -65,7 +65,9 @@ void energy_calc(const double time) { on_observable_calc(); - for (auto const &p : cell_structure.local_particles()) { + auto const local_parts = cell_structure.local_particles(); + + for (auto const &p : local_parts) { obs_energy.kinetic[0] += calc_kinetic_energy(p); } @@ -85,9 +87,16 @@ void energy_calc(const double time) { }, maximal_cutoff(), maximal_cutoff_bonded()); - calc_long_range_energies(cell_structure.local_particles()); +#ifdef ELECTROSTATICS + /* calculate k-space part of electrostatic interaction. */ + obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(local_parts); +#endif + +#ifdef DIPOLES + /* calculate k-space part of magnetostatic interaction. */ + obs_energy.dipolar[1] = Dipole::calc_energy_long_range(local_parts); +#endif - auto local_parts = cell_structure.local_particles(); Constraints::constraints.add_energy(local_parts, time, obs_energy); #ifdef CUDA @@ -111,18 +120,6 @@ REGISTER_CALLBACK(update_energy_local) void update_energy() { mpi_call_all(update_energy_local, -1, -1); } -void calc_long_range_energies(const ParticleRange &particles) { -#ifdef ELECTROSTATICS - /* calculate k-space part of electrostatic interaction. */ - obs_energy.coulomb[1] = Coulomb::calc_energy_long_range(particles); -#endif - -#ifdef DIPOLES - /* calculate k-space part of magnetostatic interaction. */ - obs_energy.dipolar[1] = Dipole::calc_energy_long_range(particles); -#endif -} - double calculate_current_potential_energy_of_system() { update_energy(); return obs_energy.accumulate(-obs_energy.kinetic[0]); diff --git a/src/core/energy.hpp b/src/core/energy.hpp index 00c439118ac..6219e15e831 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -28,7 +28,6 @@ #define _ENERGY_H #include "Observable_stat.hpp" -#include "ParticleRange.hpp" #include "actor/ActorList.hpp" extern ActorList energyActors; @@ -42,9 +41,6 @@ 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); - /** Calculate the total energy of the system. */ double calculate_current_potential_energy_of_system(); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 9380dbd24e1..b023580d042 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -57,21 +57,7 @@ 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 - /* calculate k-space part of electrostatic interaction. */ - auto const coulomb_pressure = Coulomb::calc_pressure_long_range(particles); - boost::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9); -#endif -#ifdef DIPOLES - /* calculate k-space part of magnetostatic interaction. */ - Dipole::calc_pressure_long_range(); -#endif -} - void pressure_calc() { - auto const volume = box_geo.volume(); if (!interactions_sanity_checks()) return; @@ -80,7 +66,10 @@ void pressure_calc() { on_observable_calc(); - for (auto const &p : cell_structure.local_particles()) { + auto const volume = box_geo.volume(); + auto const local_parts = cell_structure.local_particles(); + + for (auto const &p : local_parts) { add_kinetic_virials(p, obs_pressure); } @@ -106,7 +95,15 @@ void pressure_calc() { }, maximal_cutoff(), maximal_cutoff_bonded()); - calc_long_range_virials(cell_structure.local_particles()); +#ifdef ELECTROSTATICS + /* calculate k-space part of electrostatic interaction. */ + auto const coulomb_pressure = Coulomb::calc_pressure_long_range(local_parts); + boost::copy(coulomb_pressure, obs_pressure.coulomb.begin() + 9); +#endif +#ifdef DIPOLES + /* calculate k-space part of magnetostatic interaction. */ + Dipole::calc_pressure_long_range(); +#endif #ifdef VIRTUAL_SITES if (!obs_pressure.virtual_sites.empty()) { From 7b7aaebc94873dfe751234b9d5e28ffd96b7509a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Mon, 9 Aug 2021 16:18:35 +0200 Subject: [PATCH 03/11] core: Simplify Observable_stat framework Reduce code complexity by centralizing the custom reduction code into the Observable_stat class. Remove hazardous `extern int` statements inside class methods (risk of linking errors if the symbol name changes, prefer including the relevant header file). --- src/core/CMakeLists.txt | 1 - src/core/Observable_stat.cpp | 43 ++++++++++++++---- src/core/Observable_stat.hpp | 37 +++++---------- src/core/constraints/ShapeBasedConstraint.cpp | 1 + src/core/energy.cpp | 7 +-- src/core/pressure.cpp | 7 +-- src/core/reduce_observable_stat.cpp | 45 ------------------- src/core/reduce_observable_stat.hpp | 33 -------------- 8 files changed, 48 insertions(+), 126 deletions(-) delete mode 100644 src/core/reduce_observable_stat.cpp delete mode 100644 src/core/reduce_observable_stat.hpp diff --git a/src/core/CMakeLists.txt b/src/core/CMakeLists.txt index aaedbdc1eb0..8a24cfe27f0 100644 --- a/src/core/CMakeLists.txt +++ b/src/core/CMakeLists.txt @@ -41,7 +41,6 @@ set(EspressoCore_SRC CellStructure.cpp PartCfg.cpp AtomDecomposition.cpp - reduce_observable_stat.cpp EspressoSystemStandAlone.cpp DomainDecomposition.cpp) diff --git a/src/core/Observable_stat.cpp b/src/core/Observable_stat.cpp index 6822174b292..70a2c33f2b1 100644 --- a/src/core/Observable_stat.cpp +++ b/src/core/Observable_stat.cpp @@ -24,32 +24,42 @@ #include "config.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" +#include "communication.hpp" +#include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include +#include + +#include #include #include +#include #include +inline std::size_t max_non_bonded_pairs() { + auto const n = static_cast(max_seen_particle_type); + return (n * (n + 1)) / 2; +} + Observable_stat::Observable_stat(std::size_t chunk_size) : m_chunk_size(chunk_size) { // number of chunks for different interaction types - auto constexpr n_coulomb = 2; - auto constexpr n_dipolar = 2; + constexpr std::size_t n_coulomb = 2; + constexpr std::size_t n_dipolar = 2; #ifdef VIRTUAL_SITES - auto constexpr n_vs = 1; + constexpr std::size_t n_vs = 1; #else - auto constexpr n_vs = 0; + constexpr std::size_t n_vs = 0; #endif auto const n_bonded = bonded_ia_params.size(); auto const n_non_bonded = max_non_bonded_pairs(); constexpr std::size_t n_ext_fields = 1; // reduction over all fields constexpr std::size_t n_kinetic = 1; // linear+angular kinetic contributions - // resize vector - auto const total = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb + - n_dipolar + n_vs + n_ext_fields; - m_data = std::vector(m_chunk_size * total); + auto const n_elements = n_kinetic + n_bonded + 2 * n_non_bonded + n_coulomb + + n_dipolar + n_vs + n_ext_fields; + m_data = std::vector(m_chunk_size * n_elements); // spans for the different contributions kinetic = Utils::Span(m_data.data(), m_chunk_size); @@ -65,3 +75,20 @@ Observable_stat::Observable_stat(std::size_t chunk_size) Utils::Span(non_bonded_intra.end(), n_non_bonded * m_chunk_size); assert(non_bonded_inter.end() == (m_data.data() + m_data.size())); } + +Utils::Span +Observable_stat::non_bonded_contribution(Utils::Span base_pointer, + int type1, int type2) const { + auto const offset = static_cast(Utils::upper_triangular( + std::min(type1, type2), std::max(type1, type2), max_seen_particle_type)); + return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size}; +} + +void Observable_stat::mpi_reduce() { + if (comm_cart.rank() == 0) { + std::vector temp(m_data); + boost::mpi::reduce(comm_cart, temp, m_data, std::plus<>{}, 0); + } else { + boost::mpi::reduce(comm_cart, m_data, std::plus<>{}, 0); + } +} diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index a0d44e5cc15..531fb440bb4 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -23,7 +23,6 @@ #include #include -#include #include #include @@ -38,32 +37,14 @@ class Observable_stat { /** Number of doubles per data item */ std::size_t m_chunk_size; - /** Calculate the maximal number of non-bonded interaction pairs in the - * system. - */ - static std::size_t max_non_bonded_pairs() { - extern int max_seen_particle_type; - return static_cast( - (max_seen_particle_type * (max_seen_particle_type + 1)) / 2); - } - /** Get contribution from a non-bonded interaction */ Utils::Span non_bonded_contribution(Utils::Span base_pointer, - int type1, int type2) const { - extern int max_seen_particle_type; - int offset = Utils::upper_triangular( - std::min(type1, type2), std::max(type1, type2), max_seen_particle_type); - return {base_pointer.begin() + offset * m_chunk_size, m_chunk_size}; - } + int type1, int type2) const; public: explicit Observable_stat(std::size_t chunk_size); auto chunk_size() const { return m_chunk_size; } - Utils::Span data_() { return {m_data.data(), m_data.size()}; } - Utils::Span data_() const { - return {m_data.data(), m_data.size()}; - } /** Accumulate values. * @param acc Initial value for the accumulator. @@ -107,17 +88,16 @@ class Observable_stat { Utils::Span non_bonded_inter; /** 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); + Utils::Span bonded_contribution(int bond_id) const { + auto const offset = m_chunk_size * static_cast(bond_id); + return Utils::Span(bonded.data() + offset, m_chunk_size); } void add_non_bonded_contribution(int type1, int type2, Utils::Span data) { - auto const dest = - (type1 == type2) - ? non_bonded_contribution(non_bonded_intra, type1, type2) - : non_bonded_contribution(non_bonded_inter, type1, type2); + assert(data.size() == m_chunk_size); + auto const source = (type1 == type2) ? non_bonded_intra : non_bonded_inter; + auto const dest = non_bonded_contribution(source, type1, type2); boost::transform(dest, data, dest.begin(), std::plus<>{}); } @@ -137,6 +117,9 @@ class Observable_stat { int type2) const { return non_bonded_contribution(non_bonded_inter, type1, type2); } + + /** MPI reduction. */ + void mpi_reduce(); }; #endif // ESPRESSO_OBSERVABLE_STAT_HPP diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 22b1db2f30a..f94cd656b08 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -20,6 +20,7 @@ #include "ShapeBasedConstraint.hpp" #include "BoxGeometry.hpp" +#include "Observable_stat.hpp" #include "communication.hpp" #include "config.hpp" #include "dpd.hpp" diff --git a/src/core/energy.cpp b/src/core/energy.cpp index d76d1add83e..a7cb08b360a 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -33,7 +33,6 @@ #include "forces.hpp" #include "integrate.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" -#include "reduce_observable_stat.hpp" #include "short_range_loop.hpp" @@ -107,11 +106,7 @@ void energy_calc(const double time) { obs_energy.dipolar[1] += energy_host.dipolar; #endif - /* gather data */ - auto obs_energy_res = reduce(comm_cart, obs_energy); - if (obs_energy_res) { - std::swap(obs_energy, *obs_energy_res); - } + obs_energy.mpi_reduce(); } void update_energy_local(int, int) { energy_calc(get_sim_time()); } diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index b023580d042..63f8fc829ca 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -34,7 +34,6 @@ #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "pressure_inline.hpp" -#include "reduce_observable_stat.hpp" #include "virtual_sites.hpp" #include "short_range_loop.hpp" @@ -114,11 +113,7 @@ void pressure_calc() { obs_pressure.rescale(volume); - /* gather data */ - auto obs_pressure_res = reduce(comm_cart, obs_pressure); - if (obs_pressure_res) { - std::swap(obs_pressure, *obs_pressure_res); - } + obs_pressure.mpi_reduce(); } void update_pressure_local(int, int) { pressure_calc(); } diff --git a/src/core/reduce_observable_stat.cpp b/src/core/reduce_observable_stat.cpp deleted file mode 100644 index aa43c4dbff4..00000000000 --- a/src/core/reduce_observable_stat.cpp +++ /dev/null @@ -1,45 +0,0 @@ -/* - * Copyright (C) 2010-2020 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#include "reduce_observable_stat.hpp" - -#include -#include - -#include - -/** Reduce contributions from all MPI ranks. */ -boost::optional reduce(boost::mpi::communicator const &comm, - Observable_stat const &obs) { - if (comm.rank() == 0) { - Observable_stat res{obs.chunk_size()}; - - boost::mpi::reduce(comm, obs.data_().data(), obs.data_().size(), - res.data_().data(), std::plus<>{}, 0); - - return res; - } - - boost::mpi::reduce(comm, obs.data_().data(), obs.data_().size(), - std::plus<>{}, 0); - - return {}; -} diff --git a/src/core/reduce_observable_stat.hpp b/src/core/reduce_observable_stat.hpp deleted file mode 100644 index e54c02f3b32..00000000000 --- a/src/core/reduce_observable_stat.hpp +++ /dev/null @@ -1,33 +0,0 @@ -/* - * Copyright (C) 2010-2020 The ESPResSo project - * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010 - * Max-Planck-Institute for Polymer Research, Theory Group - * - * This file is part of ESPResSo. - * - * ESPResSo is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * ESPResSo is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with this program. If not, see . - */ - -#ifndef ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP -#define ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP - -#include "Observable_stat.hpp" - -#include -#include - -/** Reduce contributions from all MPI ranks. */ -boost::optional reduce(boost::mpi::communicator const &comm, - Observable_stat const &obs); -#endif // ESPRESSO_REDUCE_OBSERVABLE_STAT_HPP From a8a84737df4f174bc11c5ba2113b165d061080ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Tue, 10 Aug 2021 09:20:18 +0200 Subject: [PATCH 04/11] core: Remove Observable_stat global variables Generate these observables on-demand. Move the private function _Observable_stat_to_dict() from the .pyx file to the .pxd file. Add the missing handle_error() guards in the Cython interface. --- src/core/Observable_stat.hpp | 4 +- src/core/energy.cpp | 36 +++-- src/core/energy.hpp | 10 +- src/core/pressure.cpp | 31 ++-- src/core/pressure.hpp | 10 +- .../EspressoSystemStandAlone_test.cpp | 23 ++- src/python/espressomd/analyze.pxd | 144 ++++++++++++++++- src/python/espressomd/analyze.pyx | 149 ++---------------- 8 files changed, 202 insertions(+), 205 deletions(-) diff --git a/src/core/Observable_stat.hpp b/src/core/Observable_stat.hpp index 531fb440bb4..72109bcf892 100644 --- a/src/core/Observable_stat.hpp +++ b/src/core/Observable_stat.hpp @@ -44,7 +44,7 @@ class Observable_stat { public: explicit Observable_stat(std::size_t chunk_size); - auto chunk_size() const { return m_chunk_size; } + auto get_chunk_size() const { return m_chunk_size; } /** Accumulate values. * @param acc Initial value for the accumulator. @@ -78,7 +78,7 @@ class Observable_stat { Utils::Span coulomb; /** Contribution(s) from dipolar interactions. */ Utils::Span dipolar; - /** Contribution(s) from virtual sites (accumulated). */ + /** Contribution from virtual sites (accumulated). */ Utils::Span virtual_sites; /** Contribution from external fields (accumulated). */ Utils::Span external_fields; diff --git a/src/core/energy.cpp b/src/core/energy.cpp index a7cb08b360a..e5b67539bb8 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -39,18 +39,18 @@ #include "electrostatics_magnetostatics/coulomb.hpp" #include "electrostatics_magnetostatics/dipole.hpp" +#include + ActorList energyActors; -/** Energy of the system */ -Observable_stat obs_energy{1}; +static std::shared_ptr calculate_energy_local() { -Observable_stat const &get_obs_energy() { return obs_energy; } + auto obs_energy_ptr = std::make_shared(1); -void energy_calc(const double time) { if (!interactions_sanity_checks()) - return; + return obs_energy_ptr; - obs_energy = Observable_stat{1}; + auto &obs_energy = *obs_energy_ptr; #ifdef CUDA clear_energy_on_GPU(); @@ -71,7 +71,8 @@ void energy_calc(const double time) { } short_range_loop( - [](Particle const &p1, int bond_id, Utils::Span partners) { + [&obs_energy](Particle const &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) { @@ -80,7 +81,7 @@ void energy_calc(const double time) { } return true; }, - [](Particle const &p1, Particle const &p2, Distance const &d) { + [&obs_energy](Particle const &p1, Particle const &p2, Distance const &d) { add_non_bonded_pair_energy(p1, p2, d.vec21, sqrt(d.dist2), d.dist2, obs_energy); }, @@ -96,7 +97,7 @@ void energy_calc(const double time) { obs_energy.dipolar[1] = Dipole::calc_energy_long_range(local_parts); #endif - Constraints::constraints.add_energy(local_parts, time, obs_energy); + Constraints::constraints.add_energy(local_parts, get_sim_time(), obs_energy); #ifdef CUDA auto const energy_host = copy_energy_from_GPU(); @@ -107,20 +108,21 @@ void energy_calc(const double time) { #endif obs_energy.mpi_reduce(); + return obs_energy_ptr; } -void update_energy_local(int, int) { energy_calc(get_sim_time()); } - -REGISTER_CALLBACK(update_energy_local) +REGISTER_CALLBACK_MASTER_RANK(calculate_energy_local) -void update_energy() { mpi_call_all(update_energy_local, -1, -1); } +std::shared_ptr calculate_energy() { + return mpi_call(Communication::Result::master_rank, calculate_energy_local); +} double calculate_current_potential_energy_of_system() { - update_energy(); - return obs_energy.accumulate(-obs_energy.kinetic[0]); + auto const obs_energy = calculate_energy(); + return obs_energy->accumulate(-obs_energy->kinetic[0]); } double observable_compute_energy() { - update_energy(); - return obs_energy.accumulate(0); + auto const obs_energy = calculate_energy(); + return obs_energy->accumulate(0); } diff --git a/src/core/energy.hpp b/src/core/energy.hpp index 6219e15e831..d9ad254c7a2 100644 --- a/src/core/energy.hpp +++ b/src/core/energy.hpp @@ -30,16 +30,12 @@ #include "Observable_stat.hpp" #include "actor/ActorList.hpp" +#include + extern ActorList energyActors; /** Parallel energy calculation. */ -void energy_calc(double time); - -/** Run @ref energy_calc in parallel. */ -void update_energy(); - -/** Return the energy observable. */ -Observable_stat const &get_obs_energy(); +std::shared_ptr calculate_energy(); /** Calculate the total energy of the system. */ double calculate_current_potential_energy_of_system(); diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 63f8fc829ca..0b7680a1683 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -49,19 +49,17 @@ #include #include #include +#include #include -/** Pressure tensor of the system */ -Observable_stat obs_pressure{9}; +static std::shared_ptr calculate_pressure_local() { -Observable_stat const &get_obs_pressure() { return obs_pressure; } - -void pressure_calc() { + auto obs_pressure_ptr = std::make_shared(9); if (!interactions_sanity_checks()) - return; + return obs_pressure_ptr; - obs_pressure = Observable_stat{9}; + auto &obs_pressure = *obs_pressure_ptr; on_observable_calc(); @@ -73,7 +71,8 @@ void pressure_calc() { } short_range_loop( - [](Particle const &p1, int bond_id, Utils::Span partners) { + [&obs_pressure](Particle const &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) { @@ -88,7 +87,8 @@ void pressure_calc() { } return true; }, - [](Particle const &p1, Particle const &p2, Distance const &d) { + [&obs_pressure](Particle const &p1, Particle const &p2, + Distance const &d) { add_non_bonded_pair_virials(p1, p2, d.vec21, sqrt(d.dist2), obs_pressure); }, @@ -114,19 +114,20 @@ void pressure_calc() { obs_pressure.rescale(volume); obs_pressure.mpi_reduce(); + return obs_pressure_ptr; } -void update_pressure_local(int, int) { pressure_calc(); } - -REGISTER_CALLBACK(update_pressure_local) +REGISTER_CALLBACK_MASTER_RANK(calculate_pressure_local) -void update_pressure() { mpi_call_all(update_pressure_local, -1, -1); } +std::shared_ptr calculate_pressure() { + return mpi_call(Communication::Result::master_rank, calculate_pressure_local); +} Utils::Vector9d observable_compute_pressure_tensor() { - update_pressure(); + auto const obs_pressure = calculate_pressure(); Utils::Vector9d pressure_tensor{}; for (std::size_t j = 0; j < 9; j++) { - pressure_tensor[j] = obs_pressure.accumulate(0, j); + pressure_tensor[j] = obs_pressure->accumulate(0, j); } return pressure_tensor; } diff --git a/src/core/pressure.hpp b/src/core/pressure.hpp index 9b08b22b340..426381ceb9e 100644 --- a/src/core/pressure.hpp +++ b/src/core/pressure.hpp @@ -29,14 +29,10 @@ #include -/** Parallel pressure calculation from a virial expansion. */ -void pressure_calc(); - -/** Run @ref pressure_calc in parallel. */ -void update_pressure(); +#include -/** Return the pressure observable. */ -Observable_stat const &get_obs_pressure(); +/** Parallel pressure calculation from a virial expansion. */ +std::shared_ptr calculate_pressure(); /** Helper function for @ref Observables::PressureTensor. */ Utils::Vector9d observable_compute_pressure_tensor(); diff --git a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp index 2ec3a6725a5..ddd295170b5 100644 --- a/src/core/unit_tests/EspressoSystemStandAlone_test.cpp +++ b/src/core/unit_tests/EspressoSystemStandAlone_test.cpp @@ -79,7 +79,6 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const box_l = 8.; auto const box_center = box_l / 2.; espresso::system->set_box_l(Utils::Vector3d::broadcast(box_l)); - auto const &obs_energy = get_obs_energy(); // particle properties auto const pid1 = 9; @@ -148,8 +147,8 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, set_particle_v(pid2, {static_cast(i), 0., 0.}); auto const &p = get_particle_data(pid2); auto const kinetic_energy = 0.5 * p.p.mass * p.m.v.norm2(); - update_energy(); - BOOST_CHECK_CLOSE(obs_energy.kinetic[0], kinetic_energy, tol); + auto const obs_energy = calculate_energy(); + BOOST_CHECK_CLOSE(obs_energy->kinetic[0], kinetic_energy, tol); BOOST_CHECK_CLOSE(observable_compute_energy(), kinetic_energy, tol); } } @@ -179,12 +178,12 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const lj_energy = 4.0 * eps * (Utils::sqr(frac6) - frac6 + shift); // measure energies - update_energy(); + auto const obs_energy = calculate_energy(); for (int i = 0; i < n_pairs; ++i) { auto const ref_inter = (i == lj_pair_ab) ? lj_energy : 0.; auto const ref_intra = (i == lj_pair_bb) ? lj_energy : 0.; - BOOST_CHECK_CLOSE(obs_energy.non_bonded_inter[i], ref_inter, 500. * tol); - BOOST_CHECK_CLOSE(obs_energy.non_bonded_intra[i], ref_intra, 500. * tol); + BOOST_CHECK_CLOSE(obs_energy->non_bonded_inter[i], ref_inter, 500. * tol); + BOOST_CHECK_CLOSE(obs_energy->non_bonded_intra[i], ref_intra, 500. * tol); } } #endif // LENNARD_JONES @@ -207,7 +206,7 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, add_particle_bond(pid2, std::vector{fene_bond_id, pid3}); // measure energies - update_energy(); + auto const obs_energy = calculate_energy(); auto const &p1 = get_particle_data(pid1); auto const &p2 = get_particle_data(pid2); auto const &p3 = get_particle_data(pid3); @@ -216,9 +215,9 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, auto const fene_energy = -0.5 * fene_bond.k * Utils::sqr(fene_bond.drmax) * std::log(1.0 - Utils::sqr((dist - fene_bond.r0) / fene_bond.drmax)); - BOOST_CHECK_CLOSE(obs_energy.bonded[none_bond_id], none_energy, 0.0); - BOOST_CHECK_CLOSE(obs_energy.bonded[harm_bond_id], harm_energy, 40. * tol); - BOOST_CHECK_CLOSE(obs_energy.bonded[fene_bond_id], fene_energy, 40. * tol); + BOOST_CHECK_CLOSE(obs_energy->bonded[none_bond_id], none_energy, 0.0); + BOOST_CHECK_CLOSE(obs_energy->bonded[harm_bond_id], harm_energy, 40. * tol); + BOOST_CHECK_CLOSE(obs_energy->bonded[fene_bond_id], fene_energy, 40. * tol); } // check electrostatics @@ -246,11 +245,11 @@ BOOST_FIXTURE_TEST_CASE(espresso_system_stand_alone, ParticleFactory, place_particle(pid2, pos2); auto const r = (pos2 - pos1).norm(); // check P3M energy - update_energy(); + auto const obs_energy = calculate_energy(); // at very short distances, the real-space contribution to // the energy is much larger than the k-space contribution auto const energy_ref = -prefactor / r; - auto const energy_p3m = obs_energy.coulomb[0] + obs_energy.coulomb[1]; + auto const energy_p3m = obs_energy->coulomb[0] + obs_energy->coulomb[1]; BOOST_CHECK_CLOSE(energy_p3m, energy_ref, 0.01); } } diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 2ab4f40e24c..221ae6df68a 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -25,6 +25,11 @@ from .utils cimport Vector3i, Vector3d, Vector9d, Span from .utils cimport create_nparray_from_double_span from libcpp.vector cimport vector # import std::vector as vector from libcpp cimport bool as cbool +from libcpp.memory cimport shared_ptr +from .interactions cimport bonded_ia_params_is_type +from .interactions cimport bonded_ia_params_size +from .interactions cimport CoreNoneBond +include "myconfig.pxi" cdef extern from "" namespace "std" nogil: cdef cppclass array4 "std::array": @@ -58,7 +63,7 @@ cdef extern from "Observable_stat.hpp": Span[double] bonded_contribution(int bond_id) Span[double] non_bonded_intra_contribution(int type1, int type2) Span[double] non_bonded_inter_contribution(int type1, int type2) - size_t chunk_size() + size_t get_chunk_size() cdef extern from "statistics.hpp": cdef void calc_structurefactor(PartCfg & , const vector[int] & p_types, int order, vector[double] & wavevectors, vector[double] & intensities) except + @@ -82,12 +87,10 @@ cdef extern from "statistics_chain.hpp": array2 calc_rh(int, int, int) cdef extern from "pressure.hpp": - cdef void update_pressure() - cdef const Observable_stat & get_obs_pressure() + cdef shared_ptr[Observable_stat] calculate_pressure() cdef extern from "energy.hpp": - cdef void update_energy() - cdef const Observable_stat & get_obs_energy() + cdef shared_ptr[Observable_stat] calculate_energy() double calculate_current_potential_energy_of_system() cdef extern from "dpd.hpp": @@ -142,3 +145,134 @@ cdef inline get_obs_contrib(Span[double] contribution, int size, if value.shape[0] == 1: return value[0] return value + +cdef inline observable_stat_matrix(size_t size, cbool calc_sp): + if size == 9 and not calc_sp: + return np.zeros((3, 3), dtype=float) + else: + return 0.0 + +cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp): + """Transform an ``Observable_stat`` object to a python dict. + + Parameters + ---------- + obs : + Core observable. + calc_sp : :obj:`bool` + Whether to calculate a scalar pressure (only relevant when + ``obs`` is a pressure tensor observable). + + Returns + ------- + :obj:`dict` + A dictionary with the following keys: + + * ``"total"``: total contribution + * ``"kinetic"``: kinetic contribution + * ``"bonded"``: total bonded contribution + * ``"bonded", ``: bonded contribution which arises from the given bond_type + * ``"nonbonded"``: total nonbonded contribution + * ``"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 + * ``"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``) + * ``"virtual_sites"``: virtual sites contribution + * ``"virtual_sites", ``: contribution from virtual site i + + """ + + cdef int i + cdef int j + + size = obs.get_chunk_size() + + # Dict to store the results + p = {} + + # Total contribution + if size == 1: + p["total"] = obs.accumulate() + else: + total = np.zeros((9,), dtype=float) + for i in range(9): + total[i] = obs.accumulate(0.0, i) + total = total.reshape((3, 3)) + if calc_sp: + p["total"] = np.einsum('ii', total) / 3 + else: + p["total"] = total + + # Kinetic + p["kinetic"] = get_obs_contrib(obs.kinetic, size, calc_sp) + + # External + p["external_fields"] = get_obs_contrib(obs.external_fields, size, calc_sp) + + # Bonded + total_bonded = observable_stat_matrix(size, calc_sp) + for i in range(bonded_ia_params_size()): + if not bonded_ia_params_is_type[CoreNoneBond](i): + val = get_obs_contrib(obs.bonded_contribution(i), size, calc_sp) + p["bonded", i] = val + total_bonded += val + p["bonded"] = total_bonded + + # Non-Bonded interactions, total as well as intra and inter molecular + total_intra = observable_stat_matrix(size, calc_sp) + total_inter = observable_stat_matrix(size, calc_sp) + for i in range(max_seen_particle_type): + for j in range(i, max_seen_particle_type): + intra = get_obs_contrib( + obs.non_bonded_intra_contribution( + i, j), size, calc_sp) + total_intra += intra + p["non_bonded_intra", i, j] = intra + inter = get_obs_contrib( + obs.non_bonded_inter_contribution( + i, j), size, calc_sp) + total_inter += inter + p["non_bonded_inter", i, j] = inter + p["non_bonded", i, j] = intra + inter + + p["non_bonded_intra"] = total_intra + p["non_bonded_inter"] = total_inter + p["non_bonded"] = total_intra + total_inter + + # Electrostatics + IF ELECTROSTATICS == 1: + cdef np.ndarray coulomb + coulomb = get_obs_contribs(obs.coulomb, size, calc_sp) + for i in range(coulomb.shape[0]): + p["coulomb", i] = coulomb[i] + p["coulomb"] = np.sum(coulomb, axis=0) + + # Dipoles + IF DIPOLES == 1: + cdef np.ndarray dipolar + dipolar = get_obs_contribs(obs.dipolar, size, calc_sp) + for i in range(dipolar.shape[0]): + p["dipolar", i] = dipolar[i] + p["dipolar"] = np.sum(dipolar, axis=0) + + # virtual sites + IF VIRTUAL_SITES == 1: + cdef np.ndarray virtual_sites + virtual_sites = get_obs_contribs(obs.virtual_sites, size, calc_sp) + for i in range(virtual_sites.shape[0]): + p["virtual_sites", i] = virtual_sites[i] + p["virtual_sites"] = np.sum(virtual_sites, axis=0) + + return p + +cdef inline get_scalar_pressure(): + return Observable_stat_to_dict(calculate_pressure().get(), True) + +cdef inline get_pressure_tensor(): + return Observable_stat_to_dict(calculate_pressure().get(), False) + +cdef inline get_energy(): + return Observable_stat_to_dict(calculate_energy().get(), False) diff --git a/src/python/espressomd/analyze.pyx b/src/python/espressomd/analyze.pyx index 312ed782a2c..50df235268a 100644 --- a/src/python/espressomd/analyze.pyx +++ b/src/python/espressomd/analyze.pyx @@ -20,16 +20,11 @@ include "myconfig.pxi" from . cimport analyze from libcpp.vector cimport vector # import std::vector as vector -from libcpp cimport bool as cbool -from .interactions cimport bonded_ia_params_is_type -from .interactions cimport bonded_ia_params_size -from .interactions cimport CoreNoneBond import numpy as np cimport numpy as np import scipy.signal from .grid cimport box_geo -from collections import OrderedDict from .system import System from .utils import array_locked, is_valid_type, handle_errors from .utils cimport Vector3i, Vector3d, Vector9d @@ -93,129 +88,6 @@ def autocorrelation(time_series): f"are supported, got shape {time_series.shape}") -cdef _Observable_stat_to_dict(Observable_stat obs, cbool calc_sp): - """Transform an Observable_stat struct to a python dict. - - Parameters - ---------- - obs : - Core observable. - calc_sp : :obj:`bool` - Whether to calculate a scalar pressure (only relevant when - ``obs`` is a pressure tensor observable). - - Returns - ------- - :obj:`dict` - A dictionary with the following keys: - - * ``"total"``: total contribution - * ``"kinetic"``: kinetic contribution - * ``"bonded"``: total bonded contribution - * ``"bonded", ``: bonded contribution which arises from the given bond_type - * ``"nonbonded"``: total nonbonded contribution - * ``"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 - * ``"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``) - * ``"virtual_sites"``: virtual sites contribution - * ``"virtual_sites", ``: contribution from virtual site i - - """ - - size = obs.chunk_size() - - def set_initial(): - if size == 9 and not calc_sp: - return np.zeros((3, 3), dtype=float) - else: - return 0.0 - - cdef int i - cdef int j - - # Dict to store the results - p = OrderedDict() - - # Total contribution - if size == 1: - p["total"] = obs.accumulate() - else: - total = np.zeros((9,), dtype=float) - for i in range(9): - total[i] = obs.accumulate(0.0, i) - total = total.reshape((3, 3)) - if calc_sp: - p["total"] = np.einsum('ii', total) / 3 - else: - p["total"] = total - - # Kinetic - p["kinetic"] = get_obs_contrib(obs.kinetic, size, calc_sp) - - # External - p["external_fields"] = get_obs_contrib(obs.external_fields, size, calc_sp) - - # Bonded - total_bonded = set_initial() - for i in range(bonded_ia_params_size()): - if not bonded_ia_params_is_type[CoreNoneBond](i): - val = get_obs_contrib(obs.bonded_contribution(i), size, calc_sp) - p["bonded", i] = val - total_bonded += val - p["bonded"] = total_bonded - - # Non-Bonded interactions, total as well as intra and inter molecular - total_intra = set_initial() - total_inter = set_initial() - for i in range(analyze.max_seen_particle_type): - for j in range(i, analyze.max_seen_particle_type): - intra = get_obs_contrib( - obs.non_bonded_intra_contribution( - i, j), size, calc_sp) - total_intra += intra - p["non_bonded_intra", i, j] = intra - inter = get_obs_contrib( - obs.non_bonded_inter_contribution( - i, j), size, calc_sp) - total_inter += inter - p["non_bonded_inter", i, j] = inter - p["non_bonded", i, j] = intra + inter - - p["non_bonded_intra"] = total_intra - p["non_bonded_inter"] = total_inter - p["non_bonded"] = total_intra + total_inter - - # Electrostatics - IF ELECTROSTATICS == 1: - cdef np.ndarray coulomb - coulomb = get_obs_contribs(obs.coulomb, size, calc_sp) - for i in range(coulomb.shape[0]): - p["coulomb", i] = coulomb[i] - p["coulomb"] = np.sum(coulomb, axis=0) - - # Dipoles - IF DIPOLES == 1: - cdef np.ndarray dipolar - dipolar = get_obs_contribs(obs.dipolar, size, calc_sp) - for i in range(dipolar.shape[0]): - p["dipolar", i] = dipolar[i] - p["dipolar"] = np.sum(dipolar, axis=0) - - # virtual sites - IF VIRTUAL_SITES == 1: - cdef np.ndarray virtual_sites - virtual_sites = get_obs_contribs(obs.virtual_sites, size, calc_sp) - for i in range(virtual_sites.shape[0]): - p["virtual_sites", i] = virtual_sites[i] - p["virtual_sites"] = np.sum(virtual_sites, axis=0) - - return p - - class Analysis: def __init__(self, system): @@ -385,10 +257,9 @@ class Analysis: """ - # Update in ESPResSo core - analyze.update_pressure() - - return _Observable_stat_to_dict(analyze.get_obs_pressure(), True) + obs = analyze.get_scalar_pressure() + handle_errors("calculate_pressure() failed") + return obs def pressure_tensor(self): """Calculate the instantaneous pressure_tensor (in parallel). This is @@ -419,10 +290,9 @@ class Analysis: """ - # Update in ESPResSo core - analyze.update_pressure() - - return _Observable_stat_to_dict(analyze.get_obs_pressure(), False) + obs = analyze.get_pressure_tensor() + handle_errors("calculate_pressure() failed") + return obs IF DPD == 1: def dpd_stress(self): @@ -470,10 +340,9 @@ class Analysis: """ - analyze.update_energy() - handle_errors("calc_long_range_energies failed") - - return _Observable_stat_to_dict(analyze.get_obs_energy(), False) + obs = analyze.get_energy() + handle_errors("calculate_energy() failed") + return obs def calc_re(self, chain_start=None, number_of_chains=None, chain_length=None): From 79255edc793121b65d6178cc6800e0e4ff93967b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Wed, 11 Aug 2021 08:25:00 +0200 Subject: [PATCH 05/11] python: Remove code duplication --- src/python/espressomd/analyze.pxd | 65 +++++++++++++------------------ 1 file changed, 28 insertions(+), 37 deletions(-) diff --git a/src/python/espressomd/analyze.pxd b/src/python/espressomd/analyze.pxd index 221ae6df68a..c42750b5fdc 100644 --- a/src/python/espressomd/analyze.pxd +++ b/src/python/espressomd/analyze.pxd @@ -96,7 +96,7 @@ 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, cbool calc_scalar_pressure): """ Convert an Observable_stat range of contributions into a correctly @@ -115,7 +115,7 @@ cdef inline get_obs_contribs(Span[double] contributions, int size, """ cdef np.ndarray value value = create_nparray_from_double_span(contributions) - if size == 9: + if contributions.size() == 9: if calc_scalar_pressure: return np.einsum('...ii', value.reshape((-1, 3, 3))) / 3 else: @@ -123,7 +123,7 @@ cdef inline get_obs_contribs(Span[double] contributions, int size, else: return value -cdef inline get_obs_contrib(Span[double] contribution, int size, +cdef inline get_obs_contrib(Span[double] contribution, cbool calc_scalar_pressure): """ Convert an Observable_stat contribution into a correctly @@ -133,15 +133,13 @@ cdef inline get_obs_contrib(Span[double] contribution, int size, ---------- contributions : (N,) array_like of :obj:`float` Flattened array of energy/pressure contributions from an observable. - size : :obj:`int`, \{1, 9\} - Dimensionality of the data. calc_scalar_pressure : :obj:`bool` Whether to calculate a scalar pressure (only relevant when ``contributions`` is a pressure tensor). """ cdef np.ndarray value - value = get_obs_contribs(contribution, size, calc_scalar_pressure) + value = get_obs_contribs(contribution, calc_scalar_pressure) if value.shape[0] == 1: return value[0] return value @@ -185,55 +183,48 @@ cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp): """ - cdef int i - cdef int j - - size = obs.get_chunk_size() + cdef size_t i + cdef size_t j + cdef size_t obs_dim = obs.get_chunk_size() + cdef size_t n_bonded = bonded_ia_params_size() + cdef size_t n_nonbonded = max_seen_particle_type + cdef double[9] total # Dict to store the results p = {} # Total contribution - if size == 1: - p["total"] = obs.accumulate() - else: - total = np.zeros((9,), dtype=float) - for i in range(9): - total[i] = obs.accumulate(0.0, i) - total = total.reshape((3, 3)) - if calc_sp: - p["total"] = np.einsum('ii', total) / 3 - else: - p["total"] = total + + for i in range(obs_dim): + total[i] = obs.accumulate(0.0, i) + p["total"] = get_obs_contrib(Span[double](total, obs_dim), calc_sp) # Kinetic - p["kinetic"] = get_obs_contrib(obs.kinetic, size, calc_sp) + p["kinetic"] = get_obs_contrib(obs.kinetic, calc_sp) # External - p["external_fields"] = get_obs_contrib(obs.external_fields, size, calc_sp) + p["external_fields"] = get_obs_contrib(obs.external_fields, calc_sp) # Bonded - total_bonded = observable_stat_matrix(size, calc_sp) - for i in range(bonded_ia_params_size()): + total_bonded = observable_stat_matrix(obs_dim, calc_sp) + for i in range(n_bonded): if not bonded_ia_params_is_type[CoreNoneBond](i): - val = get_obs_contrib(obs.bonded_contribution(i), size, calc_sp) + val = get_obs_contrib(obs.bonded_contribution(i), calc_sp) p["bonded", i] = val total_bonded += val p["bonded"] = total_bonded # Non-Bonded interactions, total as well as intra and inter molecular - total_intra = observable_stat_matrix(size, calc_sp) - total_inter = observable_stat_matrix(size, calc_sp) - for i in range(max_seen_particle_type): - for j in range(i, max_seen_particle_type): + total_intra = observable_stat_matrix(obs_dim, calc_sp) + total_inter = observable_stat_matrix(obs_dim, calc_sp) + for i in range(n_nonbonded): + for j in range(i, n_nonbonded): intra = get_obs_contrib( - obs.non_bonded_intra_contribution( - i, j), size, calc_sp) + obs.non_bonded_intra_contribution(i, j), calc_sp) total_intra += intra p["non_bonded_intra", i, j] = intra inter = get_obs_contrib( - obs.non_bonded_inter_contribution( - i, j), size, calc_sp) + obs.non_bonded_inter_contribution(i, j), calc_sp) total_inter += inter p["non_bonded_inter", i, j] = inter p["non_bonded", i, j] = intra + inter @@ -245,7 +236,7 @@ cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp): # Electrostatics IF ELECTROSTATICS == 1: cdef np.ndarray coulomb - coulomb = get_obs_contribs(obs.coulomb, size, calc_sp) + coulomb = get_obs_contribs(obs.coulomb, calc_sp) for i in range(coulomb.shape[0]): p["coulomb", i] = coulomb[i] p["coulomb"] = np.sum(coulomb, axis=0) @@ -253,7 +244,7 @@ cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp): # Dipoles IF DIPOLES == 1: cdef np.ndarray dipolar - dipolar = get_obs_contribs(obs.dipolar, size, calc_sp) + dipolar = get_obs_contribs(obs.dipolar, calc_sp) for i in range(dipolar.shape[0]): p["dipolar", i] = dipolar[i] p["dipolar"] = np.sum(dipolar, axis=0) @@ -261,7 +252,7 @@ cdef inline Observable_stat_to_dict(Observable_stat * obs, cbool calc_sp): # virtual sites IF VIRTUAL_SITES == 1: cdef np.ndarray virtual_sites - virtual_sites = get_obs_contribs(obs.virtual_sites, size, calc_sp) + virtual_sites = get_obs_contribs(obs.virtual_sites, calc_sp) for i in range(virtual_sites.shape[0]): p["virtual_sites", i] = virtual_sites[i] p["virtual_sites"] = np.sum(virtual_sites, axis=0) From fdbdb669525ff94b564c44911c8bb14507739e59 Mon Sep 17 00:00:00 2001 From: David Beyer Date: Wed, 11 Aug 2021 17:57:56 +0200 Subject: [PATCH 06/11] Set hidden particle type in reaction method sample scripts manually. --- samples/grand_canonical.py | 2 ++ samples/reaction_ensemble.py | 4 ++++ samples/reaction_ensemble_complex_reaction.py | 3 +++ samples/widom_insertion.py | 5 +++++ 4 files changed, 14 insertions(+) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 35ca8290f32..d95106ad1a6 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -97,6 +97,8 @@ print(RE.get_status()) system.setup_type_map([0, 1, 2]) +# Set the hidden particle type to the lowest possible number to speed up the simulation +RE.set_non_interacting_type(len(types)) RE.reaction(10000) diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index cd0002c1933..49ccfc5fddd 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -92,6 +92,10 @@ print(RE.get_status()) system.setup_type_map([0, 1, 2]) + +# Set the hidden particle type to the lowest possible number to speed up the simulation +RE.set_non_interacting_type(3) + for i in range(10000): RE.reaction() if i % 100 == 0: diff --git a/samples/reaction_ensemble_complex_reaction.py b/samples/reaction_ensemble_complex_reaction.py index f8f78c9f70f..090fbd2409a 100644 --- a/samples/reaction_ensemble_complex_reaction.py +++ b/samples/reaction_ensemble_complex_reaction.py @@ -95,6 +95,9 @@ numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []} +# Set the hidden particle type to the lowest possible number to speed up the simulation +RE.set_non_interacting_type(len(types)) + # warmup RE.reaction(200) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 86b42e272ea..93184e296f9 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -113,6 +113,10 @@ print(widom.get_status()) system.setup_type_map([0, 1, 2]) + +# Set the hidden particle type to the lowest possible number to speed up the simulation +widom.set_non_interacting_type(len(types)) + n_iterations = 100 for i in range(n_iterations): for j in range(50): @@ -125,5 +129,6 @@ f"A- {system.number_of_particles(type=1)}", f"H+ {system.number_of_particles(type=2)}") + print("excess chemical potential for an ion pair ", widom.measure_excess_chemical_potential(insertion_reaction_id)) From 6918faa5004bd115b26a060ff9bdf83632a4fd30 Mon Sep 17 00:00:00 2001 From: David Beyer Date: Wed, 11 Aug 2021 18:09:26 +0200 Subject: [PATCH 07/11] Set hidden particle type manually in python tests for reaction methods. --- testsuite/python/constant_pH.py | 3 +++ testsuite/python/constant_pH_stats.py | 4 ++++ testsuite/python/reaction_ensemble.py | 3 +++ testsuite/python/widom_insertion.py | 3 +++ 4 files changed, 13 insertions(+) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index accca5e22f7..a059034acaa 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -55,6 +55,9 @@ def test_ideal_alpha(self): default_charges={type_HA: 0, type_A: -1, type_H: +1}) RE.constant_pH = pH + # Set the hidden particle type to the lowest possible number to speed up the simulation + RE.set_non_interacting_type(3) + # equilibration RE.reaction(800) diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index 34e24bd6a79..293a6b01c70 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -71,6 +71,10 @@ def test_ideal_titration_curve(self): type_HA = ReactionEnsembleTest.type_HA system = ReactionEnsembleTest.system RE = ReactionEnsembleTest.RE + + # Set the hidden particle type to the lowest possible number to speed up the simulation + RE.set_non_interacting_type(3) + # chemical warmup - get close to chemical equilibrium before we start # sampling RE.reaction(40 * N0) diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 00471545038..41b4a173a64 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -90,6 +90,9 @@ def test_ideal_titration_curve(self): RE = ReactionEnsembleTest.RE target_alpha = ReactionEnsembleTest.target_alpha + # Set the hidden particle type to the lowest possible number to speed up the simulation + RE.set_non_interacting_type(3) + # chemical warmup - get close to chemical equilibrium before we start # sampling RE.reaction(5 * N0) diff --git a/testsuite/python/widom_insertion.py b/testsuite/python/widom_insertion.py index 4ee45129d58..21511975f57 100644 --- a/testsuite/python/widom_insertion.py +++ b/testsuite/python/widom_insertion.py @@ -75,6 +75,9 @@ class WidomInsertionTest(ut.TestCase): Widom = espressomd.reaction_ensemble.WidomInsertion( kT=TEMPERATURE, seed=1) + # Set the hidden particle type to the lowest possible number to speed up the simulation + Widom.set_non_interacting_type(1) + def setUp(self): self.system.part.add(pos=0.5 * self.system.box_l, type=self.TYPE_HA) From a941569ec96f0c94976feda5ec31350772f03a84 Mon Sep 17 00:00:00 2001 From: David Beyer Date: Thu, 12 Aug 2021 10:16:34 +0200 Subject: [PATCH 08/11] Introduce type lists/dicts to avoid numbers. --- samples/grand_canonical.py | 5 +-- samples/reaction_ensemble.py | 36 +++++++++++-------- samples/reaction_ensemble_complex_reaction.py | 5 +-- samples/widom_insertion.py | 5 +-- 4 files changed, 31 insertions(+), 20 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index d95106ad1a6..32ce1427ddd 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -97,8 +97,9 @@ print(RE.get_status()) system.setup_type_map([0, 1, 2]) -# Set the hidden particle type to the lowest possible number to speed up the simulation -RE.set_non_interacting_type(len(types)) +# Set the hidden particle type to the lowest possible number to speed +# up the simulation +RE.set_non_interacting_type(max(types)+1) RE.reaction(10000) diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index 49ccfc5fddd..4f99c2e36ee 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -57,9 +57,16 @@ # Particle setup ############################################################# -# type 0 = HA -# type 1 = A- -# type 2 = H+ +types = { + "HA": 0, + "A-": 1, + "H+": 2, + } +charge_dict = { + 0: 0, + 1: -1, + 2: +1, + } N0 = 50 # number of titratable units K_diss = 0.0088 @@ -75,32 +82,33 @@ kT=1, exclusion_radius=1, seed=77) - RE.add_reaction(gamma=K_diss, reactant_types=[0], reactant_coefficients=[1], - product_types=[1, 2], product_coefficients=[1, 1], - default_charges={0: 0, 1: -1, 2: +1}) + RE.add_reaction(gamma=K_diss, reactant_types=[types["HA"]], reactant_coefficients=[1], + product_types=[types["A-"], types["H+"]], product_coefficients=[1, 1], + default_charges=charge_dict) elif args.mode == "constant_pH_ensemble": RE = espressomd.reaction_ensemble.ConstantpHEnsemble( kT=1, exclusion_radius=1, seed=77) RE.constant_pH = 2 - RE.add_reaction(gamma=K_diss, reactant_types=[0], - product_types=[1, 2], - default_charges={0: 0, 1: -1, 2: +1}) + RE.add_reaction(gamma=K_diss, reactant_types=[types["HA"]], + product_types=[types["A-"], types["H+"]], + default_charges=charge_dict) else: raise RuntimeError( "Please provide either --reaction_ensemble or --constant_pH_ensemble as argument ") print(RE.get_status()) -system.setup_type_map([0, 1, 2]) +system.setup_type_map(list(types.values())) -# Set the hidden particle type to the lowest possible number to speed up the simulation -RE.set_non_interacting_type(3) +# Set the hidden particle type to the lowest possible number_of_particles +# to speed up the simulation +RE.set_non_interacting_type(max(types.values())+1) for i in range(10000): RE.reaction() if i % 100 == 0: - print("HA", system.number_of_particles(type=0), "A-", - system.number_of_particles(type=1), "H+", system.number_of_particles(type=2)) + print("HA", system.number_of_particles(type=types["HA"]), "A-", + system.number_of_particles(type=types["A-"]), "H+", system.number_of_particles(type=types["H+"])) print("reaction 0 has acceptance rate: ", RE.get_acceptance_rate_reaction(0)) print("reaction 1 has acceptance rate: ", RE.get_acceptance_rate_reaction(1)) diff --git a/samples/reaction_ensemble_complex_reaction.py b/samples/reaction_ensemble_complex_reaction.py index 090fbd2409a..b666fcca4b5 100644 --- a/samples/reaction_ensemble_complex_reaction.py +++ b/samples/reaction_ensemble_complex_reaction.py @@ -95,8 +95,9 @@ numbers = {type_A: [], type_B: [], type_C: [], type_D: [], type_E: []} -# Set the hidden particle type to the lowest possible number to speed up the simulation -RE.set_non_interacting_type(len(types)) +# Set the hidden particle type to the lowest possible number to speed +# up the simulation +RE.set_non_interacting_type(max(types)+1) # warmup RE.reaction(200) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 93184e296f9..304229617a2 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -114,8 +114,9 @@ system.setup_type_map([0, 1, 2]) -# Set the hidden particle type to the lowest possible number to speed up the simulation -widom.set_non_interacting_type(len(types)) +# Set the hidden particle type to the lowest possible number to speed +# up the simulation +widom.set_non_interacting_type(max(types)+1) n_iterations = 100 for i in range(n_iterations): From 9183526f25bd61bb22aa02b7b252c0040660dc93 Mon Sep 17 00:00:00 2001 From: David Beyer Date: Thu, 12 Aug 2021 10:40:06 +0200 Subject: [PATCH 09/11] Introduce type dicts to avoid manually setting numbers. --- testsuite/python/constant_pH.py | 32 +++++++++++++++--------- testsuite/python/constant_pH_stats.py | 36 ++++++++++++++++----------- testsuite/python/reaction_ensemble.py | 36 ++++++++++++++++----------- testsuite/python/widom_insertion.py | 3 ++- 4 files changed, 64 insertions(+), 43 deletions(-) diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index a059034acaa..f1768fca239 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -32,9 +32,16 @@ class ConstantpHTest(ut.TestCase): def test_ideal_alpha(self): N0 = 40 c0 = 0.00028 - type_HA = 0 - type_A = 1 - type_H = 5 + types = { + "HA": 0, + "A-": 1, + "H+": 5, + } + charges_dict = { + 0: 0, + 1: -1, + 5: +1, + } pH = 2.0 pKa = 2.5 box_l = (N0 / c0)**(1.0 / 3.0) @@ -42,7 +49,7 @@ def test_ideal_alpha(self): system.cell_system.skin = 0.4 system.time_step = 0.01 system.part.add(pos=np.random.random((2 * N0, 3)) * system.box_l, - type=N0 * [type_A, type_H]) + type=N0 * [types["A-"], types["H+"]]) RE = espressomd.reaction_ensemble.ConstantpHEnsemble( kT=1.0, @@ -50,13 +57,14 @@ def test_ideal_alpha(self): seed=44) RE.add_reaction( gamma=10**(-pKa), - reactant_types=[type_HA], - product_types=[type_A, type_H], - default_charges={type_HA: 0, type_A: -1, type_H: +1}) + reactant_types=[types["HA"]], + product_types=[types["A-"], types["H+"]], + default_charges=charges_dict) RE.constant_pH = pH - # Set the hidden particle type to the lowest possible number to speed up the simulation - RE.set_non_interacting_type(3) + # Set the hidden particle type to the lowest possible number to speed + # up the simulation + RE.set_non_interacting_type(max(types.values())+1) # equilibration RE.reaction(800) @@ -65,9 +73,9 @@ def test_ideal_alpha(self): alphas = [] for _ in range(80): RE.reaction(15) - num_H = system.number_of_particles(type=type_H) - num_HA = system.number_of_particles(type=type_HA) - num_A = system.number_of_particles(type=type_A) + num_H = system.number_of_particles(type=types["H+"]) + num_HA = system.number_of_particles(type=types["HA"]) + num_A = system.number_of_particles(type=types["A-"]) self.assertEqual(num_H, num_A) self.assertEqual(num_A + num_HA, N0) alphas.append(num_A / N0) diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index 293a6b01c70..23c4813d911 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -29,9 +29,16 @@ class ReactionEnsembleTest(ut.TestCase): N0 = 40 c0 = 0.00028 - type_HA = 0 - type_A = 1 - type_H = 5 + types = { + "HA": 0, + "A-": 1, + "H+": 5, + } + charges_dict = { + 0: 0, + 1: -1, + 5: +1, + } temperature = 1.0 # choose target alpha not too far from 0.5 to get good statistics in a # small number of steps @@ -51,13 +58,13 @@ class ReactionEnsembleTest(ut.TestCase): def setUpClass(cls): cls.system.part.add( pos=np.random.random((2 * cls.N0, 3)) * cls.system.box_l, - type=cls.N0 * [cls.type_A, cls.type_H]) + type=cls.N0 * [cls.types["A-"], cls.types["H+"]]) cls.RE.add_reaction( gamma=cls.Ka, - reactant_types=[cls.type_HA], - product_types=[cls.type_A, cls.type_H], - default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}) + reactant_types=[cls.types["HA"]], + product_types=[cls.types["A-"], cls.types["H+"]], + default_charges=cls.charges_dict) cls.RE.constant_pH = cls.pH @classmethod @@ -66,14 +73,13 @@ def ideal_alpha(cls, pH): def test_ideal_titration_curve(self): N0 = ReactionEnsembleTest.N0 - type_A = ReactionEnsembleTest.type_A - type_H = ReactionEnsembleTest.type_H - type_HA = ReactionEnsembleTest.type_HA + types = ReactionEnsembleTest.types system = ReactionEnsembleTest.system RE = ReactionEnsembleTest.RE - # Set the hidden particle type to the lowest possible number to speed up the simulation - RE.set_non_interacting_type(3) + # Set the hidden particle type to the lowest possible number to speed + # up the simulation + RE.set_non_interacting_type(max(types.values())+1) # chemical warmup - get close to chemical equilibrium before we start # sampling @@ -85,9 +91,9 @@ def test_ideal_titration_curve(self): num_samples = 1000 for _ in range(num_samples): RE.reaction(10) - average_NH += system.number_of_particles(type=type_H) - average_NHA += system.number_of_particles(type=type_HA) - average_NA += system.number_of_particles(type=type_A) + average_NH += system.number_of_particles(type=types["H+"]) + average_NHA += system.number_of_particles(type=types["HA"]) + average_NA += system.number_of_particles(type=types["A-"]) average_NH /= float(num_samples) average_NA /= float(num_samples) average_NHA /= float(num_samples) diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 41b4a173a64..325fd450243 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -36,9 +36,16 @@ class ReactionEnsembleTest(ut.TestCase): # ensembles are equivalent only in the thermodynamic limit N \to \infty) N0 = 40 c0 = 0.00028 - type_HA = 0 - type_A = 1 - type_H = 2 + types = { + "HA": 0, + "A-": 1, + "H+": 2, + } + charge_dict = { + 0: 0, + 1: -1, + 2: +1, + } target_alpha = 0.6 # We get best statistics at alpha=0.5 Then the test is less sensitive to # the exact sequence of random numbers and does not require hard-coded @@ -46,9 +53,9 @@ class ReactionEnsembleTest(ut.TestCase): temperature = 1.0 exclusion_radius = 1.0 # could be in this test for example anywhere in the range 0.000001 ... 9, - reactant_types = [type_HA] + reactant_types = [types["HA"]] reactant_coefficients = [1] - product_types = [type_A, type_H] + product_types = [types["A-"], types["H+"]] product_coefficients = [1, 1] nubar = 1 system = espressomd.System(box_l=np.ones(3) * (N0 / c0)**(1.0 / 3.0)) @@ -68,7 +75,7 @@ class ReactionEnsembleTest(ut.TestCase): def setUpClass(cls): cls.system.part.add( pos=np.random.random((2 * cls.N0, 3)) * cls.system.box_l, - type=cls.N0 * [cls.type_A, cls.type_H]) + type=cls.N0 * [cls.types["A-"], cls.types["H+"]]) cls.RE.add_reaction( gamma=cls.gamma, @@ -76,22 +83,21 @@ def setUpClass(cls): reactant_coefficients=cls.reactant_coefficients, product_types=cls.product_types, product_coefficients=cls.product_coefficients, - default_charges={cls.type_HA: 0, cls.type_A: -1, cls.type_H: +1}, + default_charges=cls.charge_dict, check_for_electroneutrality=True) def test_ideal_titration_curve(self): N0 = ReactionEnsembleTest.N0 - type_A = ReactionEnsembleTest.type_A - type_H = ReactionEnsembleTest.type_H - type_HA = ReactionEnsembleTest.type_HA + types = ReactionEnsembleTest.types system = ReactionEnsembleTest.system gamma = ReactionEnsembleTest.gamma RE = ReactionEnsembleTest.RE target_alpha = ReactionEnsembleTest.target_alpha - # Set the hidden particle type to the lowest possible number to speed up the simulation - RE.set_non_interacting_type(3) + # Set the hidden particle type to the lowest possible number to speed + # up the simulation + RE.set_non_interacting_type(max(types.values())+1) # chemical warmup - get close to chemical equilibrium before we start # sampling @@ -103,9 +109,9 @@ def test_ideal_titration_curve(self): num_samples = 300 for _ in range(num_samples): RE.reaction(10) - average_NH += system.number_of_particles(type=type_H) - average_NHA += system.number_of_particles(type=type_HA) - average_NA += system.number_of_particles(type=type_A) + average_NH += system.number_of_particles(type=types["H+"]) + average_NHA += system.number_of_particles(type=types["HA"]) + average_NA += system.number_of_particles(type=types["A-"]) average_NH /= num_samples average_NA /= num_samples average_NHA /= num_samples diff --git a/testsuite/python/widom_insertion.py b/testsuite/python/widom_insertion.py index 21511975f57..3046a4c82eb 100644 --- a/testsuite/python/widom_insertion.py +++ b/testsuite/python/widom_insertion.py @@ -75,7 +75,8 @@ class WidomInsertionTest(ut.TestCase): Widom = espressomd.reaction_ensemble.WidomInsertion( kT=TEMPERATURE, seed=1) - # Set the hidden particle type to the lowest possible number to speed up the simulation + # Set the hidden particle type to the lowest possible number to speed + # up the simulation Widom.set_non_interacting_type(1) def setUp(self): From c6ccb50a2f6938825a9dacbf665d8b47ee83e55f Mon Sep 17 00:00:00 2001 From: David Beyer Date: Thu, 12 Aug 2021 11:00:15 +0200 Subject: [PATCH 10/11] Fix formatting. --- samples/grand_canonical.py | 2 +- samples/reaction_ensemble.py | 18 ++++++++--------- samples/reaction_ensemble_complex_reaction.py | 2 +- samples/widom_insertion.py | 2 +- testsuite/python/constant_pH.py | 18 ++++++++--------- testsuite/python/constant_pH_stats.py | 18 ++++++++--------- testsuite/python/reaction_ensemble.py | 20 +++++++++---------- 7 files changed, 40 insertions(+), 40 deletions(-) diff --git a/samples/grand_canonical.py b/samples/grand_canonical.py index 32ce1427ddd..79549132bbd 100644 --- a/samples/grand_canonical.py +++ b/samples/grand_canonical.py @@ -99,7 +99,7 @@ # Set the hidden particle type to the lowest possible number to speed # up the simulation -RE.set_non_interacting_type(max(types)+1) +RE.set_non_interacting_type(max(types) + 1) RE.reaction(10000) diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index 4f99c2e36ee..729379a9f25 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -58,15 +58,15 @@ # Particle setup ############################################################# types = { - "HA": 0, - "A-": 1, - "H+": 2, - } + "HA": 0, + "A-": 1, + "H+": 2, +} charge_dict = { - 0: 0, - 1: -1, - 2: +1, - } + 0: 0, + 1: -1, + 2: +1, +} N0 = 50 # number of titratable units K_diss = 0.0088 @@ -102,7 +102,7 @@ # Set the hidden particle type to the lowest possible number_of_particles # to speed up the simulation -RE.set_non_interacting_type(max(types.values())+1) +RE.set_non_interacting_type(max(types.values()) + 1) for i in range(10000): RE.reaction() diff --git a/samples/reaction_ensemble_complex_reaction.py b/samples/reaction_ensemble_complex_reaction.py index b666fcca4b5..17a02e88e9f 100644 --- a/samples/reaction_ensemble_complex_reaction.py +++ b/samples/reaction_ensemble_complex_reaction.py @@ -97,7 +97,7 @@ # Set the hidden particle type to the lowest possible number to speed # up the simulation -RE.set_non_interacting_type(max(types)+1) +RE.set_non_interacting_type(max(types) + 1) # warmup RE.reaction(200) diff --git a/samples/widom_insertion.py b/samples/widom_insertion.py index 304229617a2..d3339877da8 100644 --- a/samples/widom_insertion.py +++ b/samples/widom_insertion.py @@ -116,7 +116,7 @@ # Set the hidden particle type to the lowest possible number to speed # up the simulation -widom.set_non_interacting_type(max(types)+1) +widom.set_non_interacting_type(max(types) + 1) n_iterations = 100 for i in range(n_iterations): diff --git a/testsuite/python/constant_pH.py b/testsuite/python/constant_pH.py index f1768fca239..b0e1d0f14d4 100644 --- a/testsuite/python/constant_pH.py +++ b/testsuite/python/constant_pH.py @@ -33,15 +33,15 @@ def test_ideal_alpha(self): N0 = 40 c0 = 0.00028 types = { - "HA": 0, - "A-": 1, - "H+": 5, - } + "HA": 0, + "A-": 1, + "H+": 5, + } charges_dict = { - 0: 0, - 1: -1, - 5: +1, - } + 0: 0, + 1: -1, + 5: +1, + } pH = 2.0 pKa = 2.5 box_l = (N0 / c0)**(1.0 / 3.0) @@ -64,7 +64,7 @@ def test_ideal_alpha(self): # Set the hidden particle type to the lowest possible number to speed # up the simulation - RE.set_non_interacting_type(max(types.values())+1) + RE.set_non_interacting_type(max(types.values()) + 1) # equilibration RE.reaction(800) diff --git a/testsuite/python/constant_pH_stats.py b/testsuite/python/constant_pH_stats.py index 23c4813d911..f14ae0a4779 100644 --- a/testsuite/python/constant_pH_stats.py +++ b/testsuite/python/constant_pH_stats.py @@ -30,15 +30,15 @@ class ReactionEnsembleTest(ut.TestCase): N0 = 40 c0 = 0.00028 types = { - "HA": 0, - "A-": 1, - "H+": 5, - } + "HA": 0, + "A-": 1, + "H+": 5, + } charges_dict = { - 0: 0, - 1: -1, - 5: +1, - } + 0: 0, + 1: -1, + 5: +1, + } temperature = 1.0 # choose target alpha not too far from 0.5 to get good statistics in a # small number of steps @@ -79,7 +79,7 @@ def test_ideal_titration_curve(self): # Set the hidden particle type to the lowest possible number to speed # up the simulation - RE.set_non_interacting_type(max(types.values())+1) + RE.set_non_interacting_type(max(types.values()) + 1) # chemical warmup - get close to chemical equilibrium before we start # sampling diff --git a/testsuite/python/reaction_ensemble.py b/testsuite/python/reaction_ensemble.py index 325fd450243..bc3c23ea378 100644 --- a/testsuite/python/reaction_ensemble.py +++ b/testsuite/python/reaction_ensemble.py @@ -37,15 +37,15 @@ class ReactionEnsembleTest(ut.TestCase): N0 = 40 c0 = 0.00028 types = { - "HA": 0, - "A-": 1, - "H+": 2, - } + "HA": 0, + "A-": 1, + "H+": 2, + } charge_dict = { - 0: 0, - 1: -1, - 2: +1, - } + 0: 0, + 1: -1, + 2: +1, + } target_alpha = 0.6 # We get best statistics at alpha=0.5 Then the test is less sensitive to # the exact sequence of random numbers and does not require hard-coded @@ -88,7 +88,7 @@ def setUpClass(cls): def test_ideal_titration_curve(self): N0 = ReactionEnsembleTest.N0 - types = ReactionEnsembleTest.types + types = ReactionEnsembleTest.types system = ReactionEnsembleTest.system gamma = ReactionEnsembleTest.gamma @@ -97,7 +97,7 @@ def test_ideal_titration_curve(self): # Set the hidden particle type to the lowest possible number to speed # up the simulation - RE.set_non_interacting_type(max(types.values())+1) + RE.set_non_interacting_type(max(types.values()) + 1) # chemical warmup - get close to chemical equilibrium before we start # sampling From ae1268c96e59bf45e43a588a9170ede8be597ed0 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Thu, 12 Aug 2021 12:20:11 +0200 Subject: [PATCH 11/11] Formatting --- samples/reaction_ensemble.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/samples/reaction_ensemble.py b/samples/reaction_ensemble.py index 729379a9f25..523f03c4766 100644 --- a/samples/reaction_ensemble.py +++ b/samples/reaction_ensemble.py @@ -82,8 +82,11 @@ kT=1, exclusion_radius=1, seed=77) - RE.add_reaction(gamma=K_diss, reactant_types=[types["HA"]], reactant_coefficients=[1], - product_types=[types["A-"], types["H+"]], product_coefficients=[1, 1], + RE.add_reaction(gamma=K_diss, + reactant_types=[types["HA"]], + reactant_coefficients=[1], + product_types=[types["A-"], types["H+"]], + product_coefficients=[1, 1], default_charges=charge_dict) elif args.mode == "constant_pH_ensemble": RE = espressomd.reaction_ensemble.ConstantpHEnsemble(