diff --git a/src/core/BoxGeometry.hpp b/src/core/BoxGeometry.hpp index 7224ec4a2b7..de8b32e21d2 100644 --- a/src/core/BoxGeometry.hpp +++ b/src/core/BoxGeometry.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef ESPRESSO_SRC_CORE_BOX_GEOMETRY_HPP -#define ESPRESSO_SRC_CORE_BOX_GEOMETRY_HPP + +#pragma once #include "algorithm/periodic_fold.hpp" #include "lees_edwards/LeesEdwardsBC.hpp" @@ -294,5 +294,3 @@ inline Utils::Vector3d folded_position(const Utils::Vector3d &p, return p_folded; } - -#endif diff --git a/src/core/EspressoSystemStandAlone.cpp b/src/core/EspressoSystemStandAlone.cpp index 2e9062dea94..d997a478756 100644 --- a/src/core/EspressoSystemStandAlone.cpp +++ b/src/core/EspressoSystemStandAlone.cpp @@ -19,6 +19,7 @@ #include "config/config.hpp" +#include "BoxGeometry.hpp" #include "EspressoSystemStandAlone.hpp" #include "cell_system/CellStructure.hpp" #include "cell_system/CellStructureType.hpp" @@ -58,7 +59,8 @@ EspressoSystemStandAlone::EspressoSystemStandAlone(int argc, char **argv) { } void EspressoSystemStandAlone::set_box_l(Utils::Vector3d const &box_l) const { - set_box_length(box_l); + System::get_system().box_geo->set_length(box_l); + on_boxl_change(); } void EspressoSystemStandAlone::set_node_grid( diff --git a/src/core/LocalBox.hpp b/src/core/LocalBox.hpp index 538a8529038..02154a0de57 100644 --- a/src/core/LocalBox.hpp +++ b/src/core/LocalBox.hpp @@ -16,8 +16,8 @@ * You should have received a copy of the GNU General Public License * along with this program. If not, see . */ -#ifndef ESPRESSO_SRC_CORE_LOCALBOX_HPP -#define ESPRESSO_SRC_CORE_LOCALBOX_HPP + +#pragma once #include "cell_system/CellStructureType.hpp" @@ -66,5 +66,3 @@ class LocalBox { m_cell_structure_type = cell_structure_type; } }; - -#endif diff --git a/src/core/PartCfg.cpp b/src/core/PartCfg.cpp index 55227795edb..dcec7738093 100644 --- a/src/core/PartCfg.cpp +++ b/src/core/PartCfg.cpp @@ -19,8 +19,10 @@ #include "PartCfg.hpp" +#include "BoxGeometry.hpp" #include "grid.hpp" #include "particle_node.hpp" +#include "system/System.hpp" #include @@ -31,6 +33,8 @@ void PartCfg::update() { if (m_valid) return; + auto const &box_geo = *System::get_system().box_geo; + m_parts.clear(); auto const ids = get_particle_ids(); diff --git a/src/core/analysis/statistics.cpp b/src/core/analysis/statistics.cpp index eadc654d986..96930c13587 100644 --- a/src/core/analysis/statistics.cpp +++ b/src/core/analysis/statistics.cpp @@ -26,6 +26,7 @@ #include "analysis/statistics.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" #include "errorhandling.hpp" @@ -48,6 +49,7 @@ double mindist(PartCfg &partCfg, std::vector const &set1, std::vector const &set2) { using Utils::contains; + auto const &box_geo = *System::get_system().box_geo; auto mindist_sq = std::numeric_limits::infinity(); for (auto jt = partCfg.begin(); jt != partCfg.end(); ++jt) { @@ -142,6 +144,7 @@ std::vector nbhood(PartCfg &partCfg, Utils::Vector3d const &pos, double dist) { std::vector ids; auto const dist_sq = dist * dist; + auto const &box_geo = *System::get_system().box_geo; for (auto const &p : partCfg) { auto const r_sq = box_geo.get_mi_vector(pos, p.pos()).norm2(); @@ -158,6 +161,7 @@ calc_part_distribution(PartCfg &partCfg, std::vector const &p1_types, std::vector const &p2_types, double r_min, double r_max, int r_bins, bool log_flag, bool int_flag) { + auto const &box_geo = *System::get_system().box_geo; auto const r_max2 = Utils::sqr(r_max); auto const r_min2 = Utils::sqr(r_min); auto const start_dist2 = Utils::sqr(r_max + 1.); @@ -240,6 +244,7 @@ structure_factor(PartCfg &partCfg, std::vector const &p_types, int order) { if (order < 1) throw std::domain_error("order has to be a strictly positive number"); + auto const &box_geo = *System::get_system().box_geo; auto const order_sq = Utils::sqr(static_cast(order)); auto const twoPI_L = 2. * Utils::pi() * box_geo.length_inv()[0]; std::vector ff(2 * order_sq + 1); diff --git a/src/core/analysis/statistics_chain.cpp b/src/core/analysis/statistics_chain.cpp index 8787674c8f1..0f0f34ba49a 100644 --- a/src/core/analysis/statistics_chain.cpp +++ b/src/core/analysis/statistics_chain.cpp @@ -24,9 +24,11 @@ #include "analysis/statistics_chain.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "grid.hpp" #include "particle_node.hpp" +#include "system/System.hpp" #include #include @@ -36,6 +38,7 @@ #include std::array calc_re(int chain_start, int n_chains, int chain_length) { + auto const &box_geo = *System::get_system().box_geo; double dist = 0.0, dist2 = 0.0, dist4 = 0.0; std::array re; @@ -61,6 +64,7 @@ std::array calc_re(int chain_start, int n_chains, int chain_length) { } std::array calc_rg(int chain_start, int n_chains, int chain_length) { + auto const &box_geo = *System::get_system().box_geo; double r_G = 0.0, r_G2 = 0.0, r_G4 = 0.0; std::array rg; @@ -101,6 +105,7 @@ std::array calc_rg(int chain_start, int n_chains, int chain_length) { } std::array calc_rh(int chain_start, int n_chains, int chain_length) { + auto const &box_geo = *System::get_system().box_geo; double r_H = 0.0, r_H2 = 0.0; std::array rh; diff --git a/src/core/bonded_interactions/angle_common.hpp b/src/core/bonded_interactions/angle_common.hpp index 3da4b2e7d45..9b188bfaab5 100644 --- a/src/core/bonded_interactions/angle_common.hpp +++ b/src/core/bonded_interactions/angle_common.hpp @@ -24,6 +24,7 @@ * Common code for functions calculating angle forces. */ +#include "BoxGeometry.hpp" #include "config/config.hpp" #include "grid.hpp" @@ -37,6 +38,7 @@ * @f$ \vec{r_{ij}} @f$ (from particle @f$ j @f$ to particle @f$ i @f$) * and @f$ \vec{r_{kj}} @f$, and their normalization constants. * + * @param[in] box_geo Box geometry. * @param[in] r_mid Position of second/middle particle. * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. @@ -47,7 +49,8 @@ * @f$ \cos(\theta_{ijk}) @f$ */ inline std::tuple -calc_vectors_and_cosine(Utils::Vector3d const &r_mid, +calc_vectors_and_cosine(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right, bool sanitize_cosine = false) { @@ -75,6 +78,7 @@ calc_vectors_and_cosine(Utils::Vector3d const &r_mid, * See the details in @ref bondedIA_angle_force. The @f$ K(\theta_{ijk}) @f$ * term is provided as a lambda function in @p forceFactor. * + * @param[in] box_geo Box geometry. * @param[in] r_mid Position of second/middle particle. * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. @@ -86,11 +90,12 @@ calc_vectors_and_cosine(Utils::Vector3d const &r_mid, */ template std::tuple -angle_generic_force(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, +angle_generic_force(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right, ForceFactor forceFactor, bool sanitize_cosine) { auto const [vec1, vec2, d1i, d2i, cosine] = - calc_vectors_and_cosine(r_mid, r_left, r_right, sanitize_cosine); + calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, sanitize_cosine); /* force factor */ auto const fac = forceFactor(cosine); /* distribute forces */ diff --git a/src/core/bonded_interactions/angle_cosine.hpp b/src/core/bonded_interactions/angle_cosine.hpp index b9918fd102f..11beb86991a 100644 --- a/src/core/bonded_interactions/angle_cosine.hpp +++ b/src/core/bonded_interactions/angle_cosine.hpp @@ -26,6 +26,7 @@ * @ref bondedIA_angle_cosine. */ +#include "BoxGeometry.hpp" #include "angle_common.hpp" #include @@ -52,9 +53,10 @@ struct AngleCosineBond { AngleCosineBond(double bend, double phi0); std::tuple - forces(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; - double energy(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; + double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; private: @@ -75,7 +77,8 @@ struct AngleCosineBond { * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleCosineBond::forces(Utils::Vector3d const &r_mid, +AngleCosineBond::forces(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { @@ -86,7 +89,8 @@ AngleCosineBond::forces(Utils::Vector3d const &r_mid, return -bend * (sin_phi * cos_phi0 - cos_phi * sin_phi0) / sin_phi; }; - return angle_generic_force(r_mid, r_left, r_right, forceFactor, false); + return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, + false); } /** Computes the three-body angle interaction energy. @@ -94,10 +98,12 @@ AngleCosineBond::forces(Utils::Vector3d const &r_mid, * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. */ -inline double AngleCosineBond::energy(Utils::Vector3d const &r_mid, +inline double AngleCosineBond::energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { - auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true); + auto const vectors = + calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); auto const cos_phi = std::get<4>(vectors); auto const sin_phi = sqrt(1 - Utils::sqr(cos_phi)); // potential: U(phi) = k * [1 - cos(phi - phi0)] diff --git a/src/core/bonded_interactions/angle_cossquare.hpp b/src/core/bonded_interactions/angle_cossquare.hpp index 473ea7817d2..1db675825e9 100644 --- a/src/core/bonded_interactions/angle_cossquare.hpp +++ b/src/core/bonded_interactions/angle_cossquare.hpp @@ -49,9 +49,10 @@ struct AngleCossquareBond { AngleCossquareBond(double bend, double phi0); std::tuple - forces(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; - double energy(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; + double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; private: @@ -71,7 +72,8 @@ struct AngleCossquareBond { * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleCossquareBond::forces(Utils::Vector3d const &r_mid, +AngleCossquareBond::forces(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { @@ -79,7 +81,8 @@ AngleCossquareBond::forces(Utils::Vector3d const &r_mid, return bend * (cos_phi - cos_phi0); }; - return angle_generic_force(r_mid, r_left, r_right, forceFactor, false); + return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, + false); } /** Computes the three-body angle interaction energy. @@ -87,10 +90,12 @@ AngleCossquareBond::forces(Utils::Vector3d const &r_mid, * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. */ -inline double AngleCossquareBond::energy(Utils::Vector3d const &r_mid, +inline double AngleCossquareBond::energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { - auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true); + auto const vectors = + calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); auto const cos_phi = std::get<4>(vectors); return 0.5 * bend * Utils::sqr(cos_phi - cos_phi0); } diff --git a/src/core/bonded_interactions/angle_harmonic.hpp b/src/core/bonded_interactions/angle_harmonic.hpp index a1f53d08da2..d1c5ce03dd3 100644 --- a/src/core/bonded_interactions/angle_harmonic.hpp +++ b/src/core/bonded_interactions/angle_harmonic.hpp @@ -51,9 +51,10 @@ struct AngleHarmonicBond { } std::tuple - forces(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; - double energy(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; + double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; private: @@ -72,7 +73,8 @@ struct AngleHarmonicBond { * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -AngleHarmonicBond::forces(Utils::Vector3d const &r_mid, +AngleHarmonicBond::forces(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { @@ -82,7 +84,8 @@ AngleHarmonicBond::forces(Utils::Vector3d const &r_mid, return -bend * (phi - phi0) / sin_phi; }; - return angle_generic_force(r_mid, r_left, r_right, forceFactor, true); + return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, + true); } /** Compute the three-body angle interaction energy. @@ -90,10 +93,12 @@ AngleHarmonicBond::forces(Utils::Vector3d const &r_mid, * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. */ -inline double AngleHarmonicBond::energy(Utils::Vector3d const &r_mid, +inline double AngleHarmonicBond::energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { - auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true); + auto const vectors = + calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); auto const cos_phi = std::get<4>(vectors); auto const phi = acos(cos_phi); return 0.5 * bend * Utils::sqr(phi - phi0); diff --git a/src/core/bonded_interactions/bonded_tab.hpp b/src/core/bonded_interactions/bonded_tab.hpp index 0725aad8fe1..515ae986340 100644 --- a/src/core/bonded_interactions/bonded_tab.hpp +++ b/src/core/bonded_interactions/bonded_tab.hpp @@ -30,6 +30,7 @@ #include "config/config.hpp" +#include "BoxGeometry.hpp" #include "TabulatedPotential.hpp" #include "angle_common.hpp" #include "bonded_interactions/dihedral.hpp" @@ -92,9 +93,10 @@ struct TabulatedAngleBond : public TabulatedBond { TabulatedAngleBond(double min, double max, std::vector const &energy, std::vector const &force); std::tuple - forces(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, - Utils::Vector3d const &r_right) const; - double energy(Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; + double energy(BoxGeometry const &box_geo, Utils::Vector3d const &r_mid, + Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const; }; @@ -109,9 +111,11 @@ struct TabulatedDihedralBond : public TabulatedBond { std::vector const &force); boost::optional> - forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, Utils::Vector3d const &r4) const; - boost::optional energy(Utils::Vector3d const &r1, + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, + Utils::Vector3d const &r4) const; + boost::optional energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r1, Utils::Vector3d const &r2, Utils::Vector3d const &r3, Utils::Vector3d const &r4) const; @@ -155,13 +159,15 @@ TabulatedDistanceBond::energy(Utils::Vector3d const &dx) const { } /** Compute the three-body angle interaction force. + * @param[in] box_geo Box geometry. * @param[in] r_mid Position of second/middle particle. * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. * @return Forces on the second, first and third particles, in that order. */ inline std::tuple -TabulatedAngleBond::forces(Utils::Vector3d const &r_mid, +TabulatedAngleBond::forces(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { @@ -177,21 +183,25 @@ TabulatedAngleBond::forces(Utils::Vector3d const &r_mid, return -gradient / sin_phi; }; - return angle_generic_force(r_mid, r_left, r_right, forceFactor, true); + return angle_generic_force(box_geo, r_mid, r_left, r_right, forceFactor, + true); } /** Compute the three-body angle interaction energy. * It is assumed that the potential is tabulated * for all angles between 0 and Pi. * + * @param[in] box_geo Box geometry. * @param[in] r_mid Position of second/middle particle. * @param[in] r_left Position of first/left particle. * @param[in] r_right Position of third/right particle. */ -inline double TabulatedAngleBond::energy(Utils::Vector3d const &r_mid, +inline double TabulatedAngleBond::energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r_mid, Utils::Vector3d const &r_left, Utils::Vector3d const &r_right) const { - auto const vectors = calc_vectors_and_cosine(r_mid, r_left, r_right, true); + auto const vectors = + calc_vectors_and_cosine(box_geo, r_mid, r_left, r_right, true); auto const cos_phi = std::get<4>(vectors); /* calculate phi */ #ifdef TABANGLEMINUS @@ -206,6 +216,7 @@ inline double TabulatedAngleBond::energy(Utils::Vector3d const &r_mid, * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$ * (see @cite swope92a page 592). * + * @param[in] box_geo Box geometry. * @param[in] r1 Position of the first particle. * @param[in] r2 Position of the second particle. * @param[in] r3 Position of the third particle. @@ -214,7 +225,8 @@ inline double TabulatedAngleBond::energy(Utils::Vector3d const &r_mid, */ inline boost::optional> -TabulatedDihedralBond::forces(Utils::Vector3d const &r1, +TabulatedDihedralBond::forces(BoxGeometry const &box_geo, + Utils::Vector3d const &r1, Utils::Vector3d const &r2, Utils::Vector3d const &r3, Utils::Vector3d const &r4) const { @@ -226,8 +238,8 @@ TabulatedDihedralBond::forces(Utils::Vector3d const &r1, /* dihedral angle */ auto const angle_is_undefined = - calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23, - v23Xv34, l_v23Xv34, cos_phi, phi); + calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, + l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - force zero */ if (angle_is_undefined) { return {}; @@ -255,22 +267,24 @@ TabulatedDihedralBond::forces(Utils::Vector3d const &r1, /** Compute the four-body dihedral interaction energy. * The energy doesn't have any singularity if the angle phi is well-defined. * + * @param[in] box_geo Box geometry. * @param[in] r1 Position of the first particle. * @param[in] r2 Position of the second particle. * @param[in] r3 Position of the third particle. * @param[in] r4 Position of the fourth particle. */ inline boost::optional TabulatedDihedralBond::energy( - Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, Utils::Vector3d const &r4) const { + BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, + Utils::Vector3d const &r4) const { /* vectors for dihedral calculations. */ Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; double l_v12Xv23, l_v23Xv34; /* dihedral angle, cosine of the dihedral angle */ double phi, cos_phi; auto const angle_is_undefined = - calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23, - v23Xv34, l_v23Xv34, cos_phi, phi); + calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, + l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - energy zero */ if (angle_is_undefined) { return {}; diff --git a/src/core/bonded_interactions/dihedral.hpp b/src/core/bonded_interactions/dihedral.hpp index 63c28fef1cd..54a266cb0e2 100644 --- a/src/core/bonded_interactions/dihedral.hpp +++ b/src/core/bonded_interactions/dihedral.hpp @@ -57,10 +57,12 @@ struct DihedralBond { boost::optional> - forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, Utils::Vector3d const &r4) const; + forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, + Utils::Vector3d const &r4) const; - inline boost::optional energy(Utils::Vector3d const &r1, + inline boost::optional energy(BoxGeometry const &box_geo, + Utils::Vector3d const &r1, Utils::Vector3d const &r2, Utils::Vector3d const &r3, Utils::Vector3d const &r4) const; @@ -85,6 +87,7 @@ struct DihedralBond { * If the a,b or b,c are parallel the dihedral angle is not defined in which * case the function returns true. Calling functions should check for that. * + * @param[in] box_geo Box geometry. * @param[in] r1 , r2 , r3 , r4 Positions of the particles forming the dihedral * @param[out] a Vector from @p p1 to @p p2 * @param[out] b Vector from @p p2 to @p p3 @@ -98,9 +101,10 @@ struct DihedralBond { * @return Whether the angle is undefined. */ inline bool -calc_dihedral_angle(Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, Utils::Vector3d const &r4, - Utils::Vector3d &a, Utils::Vector3d &b, Utils::Vector3d &c, +calc_dihedral_angle(BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, + Utils::Vector3d const &r4, Utils::Vector3d &a, + Utils::Vector3d &b, Utils::Vector3d &c, Utils::Vector3d &aXb, double &l_aXb, Utils::Vector3d &bXc, double &l_bXc, double &cosphi, double &phi) { a = box_geo.get_mi_vector(r2, r1); @@ -141,6 +145,7 @@ calc_dihedral_angle(Utils::Vector3d const &r1, Utils::Vector3d const &r2, * The forces have a singularity at @f$ \phi = 0 @f$ and @f$ \phi = \pi @f$ * (see @cite swope92a page 592). * + * @param[in] box_geo Box geometry. * @param[in] r1 Position of the first particle. * @param[in] r2 Position of the second particle. * @param[in] r3 Position of the third particle. @@ -149,8 +154,8 @@ calc_dihedral_angle(Utils::Vector3d const &r1, Utils::Vector3d const &r2, */ inline boost::optional> -DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, +DihedralBond::forces(BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, Utils::Vector3d const &r4) const { /* vectors for dihedral angle calculation */ Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; @@ -160,8 +165,8 @@ DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2, /* dihedral angle */ auto const angle_is_undefined = - calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23, - v23Xv34, l_v23Xv34, cos_phi, phi); + calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, + l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - force zero */ if (angle_is_undefined) { return {}; @@ -199,14 +204,15 @@ DihedralBond::forces(Utils::Vector3d const &r1, Utils::Vector3d const &r2, /** Compute the four-body dihedral interaction energy. * The energy doesn't have any singularity if the angle phi is well-defined. * + * @param[in] box_geo Box geometry. * @param[in] r1 Position of the first particle. * @param[in] r2 Position of the second particle. * @param[in] r3 Position of the third particle. * @param[in] r4 Position of the fourth particle. */ inline boost::optional -DihedralBond::energy(Utils::Vector3d const &r1, Utils::Vector3d const &r2, - Utils::Vector3d const &r3, +DihedralBond::energy(BoxGeometry const &box_geo, Utils::Vector3d const &r1, + Utils::Vector3d const &r2, Utils::Vector3d const &r3, Utils::Vector3d const &r4) const { /* vectors for dihedral calculations. */ Utils::Vector3d v12, v23, v34, v12Xv23, v23Xv34; @@ -215,8 +221,8 @@ DihedralBond::energy(Utils::Vector3d const &r1, Utils::Vector3d const &r2, double phi, cos_phi; auto const angle_is_undefined = - calc_dihedral_angle(r1, r2, r3, r4, v12, v23, v34, v12Xv23, l_v12Xv23, - v23Xv34, l_v23Xv34, cos_phi, phi); + calc_dihedral_angle(box_geo, r1, r2, r3, r4, v12, v23, v34, v12Xv23, + l_v12Xv23, v23Xv34, l_v23Xv34, cos_phi, phi); /* dihedral angle not defined - energy zero */ if (angle_is_undefined) { return {}; diff --git a/src/core/cell_system/RegularDecomposition.cpp b/src/core/cell_system/RegularDecomposition.cpp index 25902264ec8..f7ef28bc426 100644 --- a/src/core/cell_system/RegularDecomposition.cpp +++ b/src/core/cell_system/RegularDecomposition.cpp @@ -27,6 +27,7 @@ #include "error_handling/RuntimeErrorStream.hpp" #include "errorhandling.hpp" #include "grid.hpp" +#include "system/System.hpp" #include #include @@ -396,6 +397,7 @@ template auto make_flat_set(Comparator &&comp) { void RegularDecomposition::init_cell_interactions() { + auto const &box_geo = *System::get_system().box_geo; auto const halo = Utils::Vector3i{1, 1, 1}; auto const cart_info = Utils::Mpi::cart_get<3>(m_comm); auto const &node_pos = cart_info.coords; diff --git a/src/core/cells.cpp b/src/core/cells.cpp index 99f4644fe63..c5cb38f2d1e 100644 --- a/src/core/cells.cpp +++ b/src/core/cells.cpp @@ -205,13 +205,19 @@ std::vector get_neighbor_pids() { void set_hybrid_decomposition(std::set n_square_types, double cutoff_regular) { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &local_geo = *system.local_geo; + auto &cell_structure = *system.cell_structure; cell_structure.set_hybrid_decomposition(comm_cart, cutoff_regular, box_geo, local_geo, n_square_types); on_cell_structure_change(); } void cells_re_init(CellStructure &cell_structure, CellStructureType new_cs) { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &local_geo = *system.local_geo; switch (new_cs) { case CellStructureType::CELL_STRUCTURE_REGULAR: cell_structure.set_regular_decomposition(comm_cart, interaction_range(), @@ -249,7 +255,9 @@ void check_resort_particles() { } void cells_update_ghosts(unsigned data_parts) { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; /* data parts that are only updated on resort */ auto constexpr resort_only_parts = diff --git a/src/core/cluster_analysis/Cluster.cpp b/src/core/cluster_analysis/Cluster.cpp index b29d0340535..bf60fba4b9a 100644 --- a/src/core/cluster_analysis/Cluster.cpp +++ b/src/core/cluster_analysis/Cluster.cpp @@ -30,6 +30,7 @@ #include "errorhandling.hpp" #include "grid.hpp" #include "particle_node.hpp" +#include "system/System.hpp" #include @@ -52,6 +53,7 @@ Utils::Vector3d Cluster::center_of_mass() { Utils::Vector3d Cluster::center_of_mass_subcluster(std::vector const &particle_ids) { sanity_checks(); + auto const &box_geo = *System::get_system().box_geo; Utils::Vector3d com{}; // The distances between the particles are "folded", such that all distances @@ -82,6 +84,7 @@ Cluster::center_of_mass_subcluster(std::vector const &particle_ids) { double Cluster::longest_distance() { sanity_checks(); + auto const &box_geo = *System::get_system().box_geo; double ld = 0.; for (auto a = particles.begin(); a != particles.end(); a++) { for (auto b = a; ++b != particles.end();) { @@ -105,6 +108,7 @@ double Cluster::radius_of_gyration() { double Cluster::radius_of_gyration_subcluster(std::vector const &particle_ids) { sanity_checks(); + auto const &box_geo = *System::get_system().box_geo; // Center of mass Utils::Vector3d com = center_of_mass_subcluster(particle_ids); double sum_sq_dist = 0.; @@ -133,6 +137,7 @@ std::vector sort_indices(const std::vector &v) { std::pair Cluster::fractal_dimension(double dr) { #ifdef GSL sanity_checks(); + auto const &box_geo = *System::get_system().box_geo; Utils::Vector3d com = center_of_mass(); // calculate Df using linear regression on the logarithms of the radii of // gyration against the number of particles in sub-clusters. Particles are @@ -182,7 +187,8 @@ std::pair Cluster::fractal_dimension(double dr) { } void Cluster::sanity_checks() const { - if (::box_geo.type() != BoxType::CUBOID) { + auto const &box_geo = *System::get_system().box_geo; + if (box_geo.type() != BoxType::CUBOID) { throw std::runtime_error( "Cluster analysis is not compatible with non-cuboid box types"); } diff --git a/src/core/cluster_analysis/ClusterStructure.cpp b/src/core/cluster_analysis/ClusterStructure.cpp index 6cd7a03cd3b..0f32c1a6519 100644 --- a/src/core/cluster_analysis/ClusterStructure.cpp +++ b/src/core/cluster_analysis/ClusterStructure.cpp @@ -22,9 +22,9 @@ #include "Cluster.hpp" #include "PartCfg.hpp" #include "errorhandling.hpp" -#include "grid.hpp" #include "partCfg_global.hpp" #include "particle_node.hpp" +#include "system/System.hpp" #include @@ -195,7 +195,7 @@ int ClusterStructure::get_next_free_cluster_id() { } void ClusterStructure::sanity_checks() const { - if (::box_geo.type() != BoxType::CUBOID) { + if (System::get_system().box_geo->type() != BoxType::CUBOID) { throw std::runtime_error( "Cluster analysis is not compatible with non-cuboid box types"); } diff --git a/src/core/collision.cpp b/src/core/collision.cpp index d03d3d8abd8..da8a0b41c90 100644 --- a/src/core/collision.cpp +++ b/src/core/collision.cpp @@ -20,6 +20,7 @@ #include "collision.hpp" #ifdef COLLISION_DETECTION +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" #include "cell_system/Cell.hpp" @@ -30,6 +31,7 @@ #include "event.hpp" #include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "system/System.hpp" #include "virtual_sites.hpp" #include @@ -235,45 +237,46 @@ void queue_collision(const int part1, const int part2) { /** @brief Calculate position of vs for GLUE_TO_SURFACE mode. * Returns id of particle to bind vs to. */ -const Particle &glue_to_surface_calc_vs_pos(const Particle &p1, - const Particle &p2, - Utils::Vector3d &pos) { +static auto const &glue_to_surface_calc_vs_pos(Particle const &p1, + Particle const &p2, + BoxGeometry const &box_geo, + Utils::Vector3d &pos) { double c; auto const vec21 = box_geo.get_mi_vector(p1.pos(), p2.pos()); - const double dist_betw_part = vec21.norm(); + auto const dist = vec21.norm(); // Find out, which is the particle to be glued. if ((p1.type() == collision_params.part_type_to_be_glued) and (p2.type() == collision_params.part_type_to_attach_vs_to)) { - c = 1 - collision_params.dist_glued_part_to_vs / dist_betw_part; + c = 1. - collision_params.dist_glued_part_to_vs / dist; } else if ((p2.type() == collision_params.part_type_to_be_glued) and (p1.type() == collision_params.part_type_to_attach_vs_to)) { - c = collision_params.dist_glued_part_to_vs / dist_betw_part; + c = collision_params.dist_glued_part_to_vs / dist; } else { throw std::runtime_error("This should never be thrown. Bug."); } - for (int i = 0; i < 3; i++) { - pos[i] = p2.pos()[i] + vec21[i] * c; - } + pos = p2.pos() + vec21 * c; if (p1.type() == collision_params.part_type_to_attach_vs_to) return p1; return p2; } -void bind_at_point_of_collision_calc_vs_pos(const Particle *const p1, - const Particle *const p2, - Utils::Vector3d &pos1, - Utils::Vector3d &pos2) { - auto const vec21 = box_geo.get_mi_vector(p1->pos(), p2->pos()); - pos1 = p1->pos() - vec21 * collision_params.vs_placement; - pos2 = p1->pos() - vec21 * (1. - collision_params.vs_placement); +static void bind_at_point_of_collision_calc_vs_pos(Particle const &p1, + Particle const &p2, + BoxGeometry const &box_geo, + Utils::Vector3d &pos1, + Utils::Vector3d &pos2) { + auto const vec21 = box_geo.get_mi_vector(p1.pos(), p2.pos()); + pos1 = p1.pos() - vec21 * collision_params.vs_placement; + pos2 = p1.pos() - vec21 * (1. - collision_params.vs_placement); } // Considers three particles for three_particle_binding and performs // the binding if the criteria are met -void coldet_do_three_particle_bond(Particle &p, Particle const &p1, - Particle const &p2) { +static void coldet_do_three_particle_bond(Particle &p, Particle const &p1, + Particle const &p2, + BoxGeometry const &box_geo) { // If p1 and p2 are not closer or equal to the cutoff distance, skip // p1: if (box_geo.get_mi_vector(p.pos(), p1.pos()).norm() > @@ -339,6 +342,7 @@ void coldet_do_three_particle_bond(Particle &p, Particle const &p1, #ifdef VIRTUAL_SITES_RELATIVE static void place_vs_and_relate_to_particle(CellStructure &cell_structure, + BoxGeometry const &box_geo, int const current_vs_pid, Utils::Vector3d const &pos, int const relate_to) { @@ -346,7 +350,7 @@ static void place_vs_and_relate_to_particle(CellStructure &cell_structure, new_part.id() = current_vs_pid; new_part.pos() = pos; auto p_vs = cell_structure.add_particle(std::move(new_part)); - vs_relate_to(*p_vs, get_part(cell_structure, relate_to)); + vs_relate_to(*p_vs, get_part(cell_structure, relate_to), box_geo); p_vs->set_virtual(true); p_vs->type() = collision_params.vs_particle_type; } @@ -407,29 +411,30 @@ std::vector gather_global_collision_queue() { } static void three_particle_binding_do_search(Cell *basecell, Particle &p1, - Particle &p2) { - auto handle_cell = [&p1, &p2](Cell *c) { - for (auto &P : c->particles()) { + Particle &p2, + BoxGeometry const &box_geo) { + auto handle_cell = [&p1, &p2, &box_geo](Cell *c) { + for (auto &p : c->particles()) { // Skip collided particles themselves - if ((P.id() == p1.id()) || (P.id() == p2.id())) { + if ((p.id() == p1.id()) or (p.id() == p2.id())) { continue; } // We need all cyclical permutations, here (bond is placed on 1st // particle, order of bond partners does not matter, so we don't need // non-cyclic permutations). - // coldet_do_three_particle_bond checks the bonding criterion and if + // @ref coldet_do_three_particle_bond checks the bonding criterion and if // the involved particles are not already bonded before it binds them. - if (!P.is_ghost()) { - coldet_do_three_particle_bond(P, p1, p2); + if (!p.is_ghost()) { + coldet_do_three_particle_bond(p, p1, p2, box_geo); } if (!p1.is_ghost()) { - coldet_do_three_particle_bond(p1, P, p2); + coldet_do_three_particle_bond(p1, p, p2, box_geo); } if (!p2.is_ghost()) { - coldet_do_three_particle_bond(p2, P, p1); + coldet_do_three_particle_bond(p2, p, p1, box_geo); } } }; @@ -447,7 +452,7 @@ static void three_particle_binding_do_search(Cell *basecell, Particle &p1, // looks for a third particle by using the domain decomposition // cell system. If found, it performs three particle binding static void three_particle_binding_domain_decomposition( - CellStructure &cell_structure, + CellStructure &cell_structure, BoxGeometry const &box_geo, std::vector const &gathered_queue) { for (auto &c : gathered_queue) { @@ -461,9 +466,9 @@ static void three_particle_binding_domain_decomposition( auto cell2 = cell_structure.find_current_cell(p2); if (cell1) - three_particle_binding_do_search(cell1, p1, p2); + three_particle_binding_do_search(cell1, p1, p2, box_geo); if (cell2 and cell1 != cell2) - three_particle_binding_do_search(cell2, p1, p2); + three_particle_binding_do_search(cell2, p1, p2, box_geo); } // If local particles exist } // Loop over total collisions @@ -471,6 +476,7 @@ static void three_particle_binding_domain_decomposition( // Handle the collisions stored in the queue void handle_collisions(CellStructure &cell_structure) { + auto const &box_geo = *System::get_system().box_geo; // Note that the glue to surface mode adds bonds between the centers // but does so later in the process. This is needed to guarantee that // a particle can only be glued once, even if queued twice in a single @@ -551,12 +557,12 @@ void handle_collisions(CellStructure &cell_structure) { p2->set_can_rotate_all_axes(); // Positions of the virtual sites - bind_at_point_of_collision_calc_vs_pos(p1, p2, pos1, pos2); + bind_at_point_of_collision_calc_vs_pos(*p1, *p2, box_geo, pos1, pos2); auto handle_particle = [&](Particle *p, Utils::Vector3d const &pos) { if (not p->is_ghost()) { - place_vs_and_relate_to_particle(cell_structure, current_vs_pid, - pos, p->id()); + place_vs_and_relate_to_particle(cell_structure, box_geo, + current_vs_pid, pos, p->id()); // Particle storage locations may have changed due to // added particle p1 = cell_structure.get_local_particle(c.pp1); @@ -593,8 +599,8 @@ void handle_collisions(CellStructure &cell_structure) { } Utils::Vector3d pos; - const Particle &attach_vs_to = - glue_to_surface_calc_vs_pos(*p1, *p2, pos); + Particle const &attach_vs_to = + glue_to_surface_calc_vs_pos(*p1, *p2, box_geo, pos); // Add a bond between the centers of the colliding particles // The bond is placed on the node that has p1 @@ -615,7 +621,8 @@ void handle_collisions(CellStructure &cell_structure) { // Vs placement happens on the node that has p1 if (!attach_vs_to.is_ghost()) { - place_vs_and_relate_to_particle(cell_structure, current_vs_pid, pos, + place_vs_and_relate_to_particle(cell_structure, box_geo, + current_vs_pid, pos, attach_vs_to.id()); // Particle storage locations may have changed due to // added particle @@ -647,7 +654,8 @@ void handle_collisions(CellStructure &cell_structure) { // three-particle-binding part if (collision_params.mode == CollisionModeType::BIND_THREE_PARTICLES) { auto gathered_queue = gather_global_collision_queue(); - three_particle_binding_domain_decomposition(cell_structure, gathered_queue); + three_particle_binding_domain_decomposition(cell_structure, box_geo, + gathered_queue); } // if TPB local_collision_queue.clear(); diff --git a/src/core/constraints/Constraints.hpp b/src/core/constraints/Constraints.hpp index ead490ea029..276574e57e4 100644 --- a/src/core/constraints/Constraints.hpp +++ b/src/core/constraints/Constraints.hpp @@ -19,8 +19,10 @@ #ifndef CORE_CONSTRAINTS_CONSTRAINTS_HPP #define CORE_CONSTRAINTS_CONSTRAINTS_HPP +#include "BoxGeometry.hpp" #include "Observable_stat.hpp" #include "grid.hpp" +#include "system/System.hpp" #include #include @@ -50,6 +52,7 @@ template class Constraints { public: void add(std::shared_ptr const &constraint) { + auto const &box_geo = *System::get_system().box_geo; if (not constraint->fits_in_box(box_geo.length())) { throw std::runtime_error("Constraint not compatible with box size."); } @@ -73,7 +76,8 @@ template class Constraints { const_iterator begin() const { return m_constraints.begin(); } const_iterator end() const { return m_constraints.end(); } - void add_forces(ParticleRange &particles, double t) const { + void add_forces(BoxGeometry const &box_geo, ParticleRange &particles, + double time) const { if (m_constraints.empty()) return; @@ -83,16 +87,16 @@ template class Constraints { auto const pos = folded_position(p.pos(), box_geo); ParticleForce force{}; for (auto const &constraint : *this) { - force += constraint->force(p, pos, t); + force += constraint->force(p, pos, time); } p.force_and_torque() += force; } } - void add_energy(const ParticleRange &particles, double time, - Observable_stat &obs_energy) const { - for (auto &p : particles) { + void add_energy(BoxGeometry const &box_geo, ParticleRange const &particles, + double time, Observable_stat &obs_energy) const { + for (auto const &p : particles) { auto const pos = folded_position(p.pos(), box_geo); for (auto const &constraint : *this) { diff --git a/src/core/constraints/ShapeBasedConstraint.cpp b/src/core/constraints/ShapeBasedConstraint.cpp index 8a3a1b85542..1e2b7f3750e 100644 --- a/src/core/constraints/ShapeBasedConstraint.cpp +++ b/src/core/constraints/ShapeBasedConstraint.cpp @@ -51,12 +51,13 @@ double ShapeBasedConstraint::total_normal_force() const { } double ShapeBasedConstraint::min_dist(const ParticleRange &particles) { + auto const &box_geo = *System::get_system().box_geo; double global_mindist = std::numeric_limits::infinity(); auto const local_mindist = std::accumulate( particles.begin(), particles.end(), std::numeric_limits::infinity(), - [this](double min, Particle const &p) { + [this, &box_geo](double min, Particle const &p) { auto const &ia_params = get_ia_param(p.type(), part_rep.type()); if (checkIfInteraction(ia_params)) { double dist; diff --git a/src/core/dpd.cpp b/src/core/dpd.cpp index ead9419beba..99bb694e525 100644 --- a/src/core/dpd.cpp +++ b/src/core/dpd.cpp @@ -27,6 +27,7 @@ #include "dpd.hpp" +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" #include "event.hpp" @@ -108,6 +109,8 @@ Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, return {}; } + auto const &box_geo = *System::get_system().box_geo; + auto const v21 = box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v()); auto const noise_vec = @@ -127,29 +130,32 @@ Utils::Vector3d dpd_pair_force(Particle const &p1, Particle const &p2, } static auto dpd_viscous_stress_local() { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; on_observable_calc(); Utils::Matrix stress{}; - cell_structure.non_bonded_loop( - [&stress](const Particle &p1, const Particle &p2, Distance const &d) { - auto const v21 = - box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v()); + cell_structure.non_bonded_loop([&stress, &box_geo](Particle const &p1, + Particle const &p2, + Distance const &d) { + auto const v21 = + box_geo.velocity_difference(p1.pos(), p2.pos(), p1.v(), p2.v()); - auto const &ia_params = get_ia_param(p1.type(), p2.type()); - auto const dist = std::sqrt(d.dist2); + auto const &ia_params = get_ia_param(p1.type(), p2.type()); + auto const dist = std::sqrt(d.dist2); - auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, {}); - auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, {}); + auto const f_r = dpd_pair_force(ia_params.dpd.radial, v21, dist, {}); + auto const f_t = dpd_pair_force(ia_params.dpd.trans, v21, dist, {}); - /* Projection operator to radial direction */ - auto const P = tensor_product(d.vec21 / d.dist2, d.vec21); - /* This is equivalent to P * f_r + (1 - P) * f_t, but with - * doing only one matrix-vector multiplication */ - auto const f = P * (f_r - f_t) + f_t; + /* Projection operator to radial direction */ + auto const P = tensor_product(d.vec21 / d.dist2, d.vec21); + /* This is equivalent to P * f_r + (1 - P) * f_t, but with + * doing only one matrix-vector multiplication */ + auto const f = P * (f_r - f_t) + f_t; - stress += tensor_product(d.vec21, f); - }); + stress += tensor_product(d.vec21, f); + }); return stress; } @@ -170,6 +176,7 @@ static auto dpd_viscous_stress_local() { * @return Stress tensor contribution. */ Utils::Vector9d dpd_stress(boost::mpi::communicator const &comm) { + auto const &box_geo = *System::get_system().box_geo; auto const local_stress = dpd_viscous_stress_local(); std::remove_const_t global_stress; diff --git a/src/core/ek/EKWalberla.cpp b/src/core/ek/EKWalberla.cpp index 2a23fd24348..a828440c5b8 100644 --- a/src/core/ek/EKWalberla.cpp +++ b/src/core/ek/EKWalberla.cpp @@ -21,10 +21,10 @@ #ifdef WALBERLA +#include "BoxGeometry.hpp" #include "ek/EKReactions.hpp" #include "ek/EKWalberla.hpp" #include "errorhandling.hpp" -#include "grid.hpp" #include "integrate.hpp" #include "lb/Implementation.hpp" #include "lb/LBWalberla.hpp" @@ -110,8 +110,9 @@ void EKWalberla::veto_time_step(double time_step) const { } void EKWalberla::sanity_checks() const { + auto const &box_geo = *System::get_system().box_geo; auto const &lattice = ek_container->get_lattice(); - auto const agrid = ::box_geo.length()[0] / lattice.get_grid_dimensions()[0]; + auto const agrid = box_geo.length()[0] / lattice.get_grid_dimensions()[0]; auto [my_left, my_right] = lattice.get_local_domain(); my_left *= agrid; my_right *= agrid; diff --git a/src/core/electrostatics/coulomb_inline.hpp b/src/core/electrostatics/coulomb_inline.hpp index 81b6cecec67..e3365c9ba33 100644 --- a/src/core/electrostatics/coulomb_inline.hpp +++ b/src/core/electrostatics/coulomb_inline.hpp @@ -26,6 +26,7 @@ #include "electrostatics/solver.hpp" #include "Particle.hpp" +#include "system/System.hpp" #include #include @@ -85,9 +86,11 @@ struct ShortRangeForceCorrectionsKernel { result_type operator()(std::shared_ptr const &ptr) const { auto const &actor = *ptr; - return kernel_type{[&actor](Particle &p1, Particle &p2, double q1q2) { - actor.add_pair_force_corrections(p1, p2, q1q2); - }}; + auto const &box_geo = *System::get_system().box_geo; + return kernel_type{ + [&actor, &box_geo](Particle &p1, Particle &p2, double q1q2) { + actor.add_pair_force_corrections(p1, p2, q1q2, box_geo); + }}; } #endif // P3M }; @@ -134,15 +137,16 @@ struct ShortRangeEnergyKernel { result_type operator()(std::shared_ptr const &ptr) const { auto const &actor = *ptr; + auto const &box_geo = *System::get_system().box_geo; auto const energy_kernel = std::visit(*this, actor.base_solver); - return kernel_type{[&actor, energy_kernel]( + return kernel_type{[&actor, &box_geo, energy_kernel]( Particle const &p1, Particle const &p2, double q1q2, Utils::Vector3d const &d, double dist) { auto energy = 0.; if (energy_kernel) { energy = (*energy_kernel)(p1, p2, q1q2, d, dist); } - return energy + actor.pair_energy_correction(q1q2, p1, p2); + return energy + actor.pair_energy_correction(p1, p2, q1q2, box_geo); }}; } #endif // P3M diff --git a/src/core/electrostatics/elc.cpp b/src/core/electrostatics/elc.cpp index 97eff73cce1..09da66d5190 100644 --- a/src/core/electrostatics/elc.cpp +++ b/src/core/electrostatics/elc.cpp @@ -30,6 +30,7 @@ #include "electrostatics/p3m.hpp" #include "electrostatics/p3m_gpu.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "cell_system/CellStructure.hpp" @@ -125,7 +126,8 @@ static std::vector calc_sc_cache(ParticleRange const &particles, } static std::pair -prepare_sc_cache(ParticleRange const &particles, double far_cut) { +prepare_sc_cache(ParticleRange const &particles, BoxGeometry const &box_geo, + double far_cut) { assert(far_cut >= 0.); auto const n_freq_x = static_cast(std::ceil(far_cut * box_geo.length()[0]) + 1.); @@ -201,6 +203,7 @@ void ElectrostaticLayerCorrection::check_gap(Particle const &p) const { void ElectrostaticLayerCorrection::add_dipole_force( ParticleRange const &particles) const { constexpr std::size_t size = 3; + auto const &box_geo = *System::get_system().box_geo; auto const pref = prefactor * 4. * Utils::pi() / box_geo.volume(); /* for non-neutral systems, this shift gives the background contribution @@ -266,6 +269,7 @@ void ElectrostaticLayerCorrection::add_dipole_force( double ElectrostaticLayerCorrection::dipole_energy( ParticleRange const &particles) const { constexpr std::size_t size = 7; + auto const &box_geo = *System::get_system().box_geo; auto const pref = prefactor * 2. * Utils::pi() / box_geo.volume(); auto const lz = box_geo.length()[2]; /* for nonneutral systems, this shift gives the background contribution @@ -363,6 +367,7 @@ struct ImageSum { double ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const { constexpr std::size_t size = 4; + auto const &box_geo = *System::get_system().box_geo; auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1]; auto const pref = prefactor * 2. * Utils::pi() * xy_area_inv; auto const delta = elc.delta_mid_top * elc.delta_mid_bot; @@ -438,6 +443,7 @@ ElectrostaticLayerCorrection::z_energy(ParticleRange const &particles) const { void ElectrostaticLayerCorrection::add_z_force( ParticleRange const &particles) const { constexpr std::size_t size = 1; + auto const &box_geo = *System::get_system().box_geo; auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1]; auto const pref = prefactor * 2. * Utils::pi() * xy_area_inv; auto const delta = elc.delta_mid_top * elc.delta_mid_bot; @@ -498,6 +504,7 @@ void setup_PoQ(elc_data const &elc, double prefactor, std::size_t index, double omega, ParticleRange const &particles) { assert(index >= 1); constexpr std::size_t size = 4; + auto const &box_geo = *System::get_system().box_geo; auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1]; auto const pref_di = prefactor * 4. * Utils::pi() * xy_area_inv; auto const pref = -pref_di / expm1(omega * box_geo.length()[2]); @@ -639,6 +646,7 @@ static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p, assert(index_p >= 1); assert(index_q >= 1); constexpr std::size_t size = 8; + auto const &box_geo = *System::get_system().box_geo; auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1]; auto const pref_di = prefactor * 8 * Utils::pi() * xy_area_inv; auto const pref = -pref_di / expm1(omega * box_geo.length()[2]); @@ -763,7 +771,8 @@ static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p, } static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega, - const ParticleRange &particles) { + ParticleRange const &particles, + BoxGeometry const &box_geo) { auto constexpr c_2pi = 2. * Utils::pi(); auto const pref_x = c_2pi * box_geo.length_inv()[0] * static_cast(index_p) / omega; @@ -823,7 +832,8 @@ static double PQ_energy(double omega, std::size_t n_part) { void ElectrostaticLayerCorrection::add_force( ParticleRange const &particles) const { auto constexpr c_2pi = 2. * Utils::pi(); - auto const n_freqs = prepare_sc_cache(particles, elc.far_cut); + auto const &box_geo = *System::get_system().box_geo; + auto const n_freqs = prepare_sc_cache(particles, box_geo, elc.far_cut); auto const n_scxcache = std::get<0>(n_freqs); auto const n_scycache = std::get<1>(n_freqs); partblk.resize(particles.size() * 8); @@ -869,7 +879,7 @@ void ElectrostaticLayerCorrection::add_force( Utils::sqr(box_geo.length_inv()[1] * static_cast(q))); setup_PQ(elc, prefactor, p, q, omega, particles); distribute(8); - add_PQ_force(p, q, omega, particles); + add_PQ_force(p, q, omega, particles, box_geo); } } } @@ -877,8 +887,9 @@ void ElectrostaticLayerCorrection::add_force( double ElectrostaticLayerCorrection::calc_energy( ParticleRange const &particles) const { auto constexpr c_2pi = 2. * Utils::pi(); + auto const &box_geo = *System::get_system().box_geo; auto energy = dipole_energy(particles) + z_energy(particles); - auto const n_freqs = prepare_sc_cache(particles, elc.far_cut); + auto const n_freqs = prepare_sc_cache(particles, box_geo, elc.far_cut); auto const n_scxcache = std::get<0>(n_freqs); auto const n_scycache = std::get<1>(n_freqs); @@ -933,6 +944,7 @@ double ElectrostaticLayerCorrection::calc_energy( double ElectrostaticLayerCorrection::tune_far_cut() const { // Largest reasonable cutoff for far formula auto constexpr maximal_far_cut = 50.; + auto const &box_geo = *System::get_system().box_geo; auto const box_l_x_inv = box_geo.length_inv()[0]; auto const box_l_y_inv = box_geo.length_inv()[1]; auto const min_inv_boxl = std::min(box_l_x_inv, box_l_y_inv); @@ -972,6 +984,7 @@ static auto calc_total_charge() { } void ElectrostaticLayerCorrection::sanity_checks_periodicity() const { + auto const &box_geo = *System::get_system().box_geo; if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) { throw std::runtime_error("ELC: requires periodicity (True, True, True)"); } @@ -1006,6 +1019,7 @@ void ElectrostaticLayerCorrection::adapt_solver() { } void ElectrostaticLayerCorrection::recalc_box_h() { + auto const &box_geo = *System::get_system().box_geo; auto const new_box_h = box_geo.length()[2] - elc.gap_size; if (new_box_h < 0.) { throw std::runtime_error("ELC gap size (" + std::to_string(elc.gap_size) + @@ -1043,9 +1057,9 @@ elc_data::elc_data(double maxPWerror, double gap_size, double far_cut, bool neutralize, double delta_top, double delta_bot, bool with_const_pot, double potential_diff) : maxPWerror{maxPWerror}, gap_size{gap_size}, - box_h{box_geo.length()[2] - gap_size}, far_cut{far_cut}, far_cut2{-1.}, - far_calculated{far_cut == -1.}, dielectric_contrast_on{delta_top != 0. or - delta_bot != 0.}, + box_h{System::get_system().box_geo->length()[2] - gap_size}, + far_cut{far_cut}, far_cut2{-1.}, far_calculated{far_cut == -1.}, + dielectric_contrast_on{delta_top != 0. or delta_bot != 0.}, const_pot{with_const_pot and dielectric_contrast_on}, neutralize{neutralize and !dielectric_contrast_on}, delta_mid_top{std::clamp(delta_top, -1., +1.)}, delta_mid_bot{std::clamp( @@ -1094,11 +1108,6 @@ ElectrostaticLayerCorrection::ElectrostaticLayerCorrection( adapt_solver(); } -Utils::Vector3d elc_data::get_mi_vector(Utils::Vector3d const &a, - Utils::Vector3d const &b) const { - return box_geo.get_mi_vector(a, b); -} - static void p3m_assign_image_charge(elc_data const &elc, CoulombP3M &p3m, double q, Utils::Vector3d const &pos) { if (pos[2] < elc.space_layer) { @@ -1181,6 +1190,7 @@ double ElectrostaticLayerCorrection::long_range_energy( auto const energy = std::visit( [this, &particles](auto const &solver_ptr) { auto &solver = *solver_ptr; + auto const &box_geo = *System::get_system().box_geo; // assign the original charges (they may not have been assigned yet) solver.charge_assign(particles); @@ -1191,7 +1201,8 @@ double ElectrostaticLayerCorrection::long_range_energy( auto energy = 0.; energy += 0.5 * solver.long_range_energy(particles); - energy += 0.5 * elc.dielectric_layers_self_energy(solver, particles); + energy += + 0.5 * elc.dielectric_layers_self_energy(solver, box_geo, particles); // assign both original and image charges charge_assign(elc, solver, particles); @@ -1218,9 +1229,10 @@ void ElectrostaticLayerCorrection::add_long_range_forces( [this, &particles](auto const &solver_ptr) { auto &solver = *solver_ptr; if (elc.dielectric_contrast_on) { + auto const &box_geo = *System::get_system().box_geo; modify_p3m_sums(elc, solver, particles); charge_assign(elc, solver, particles); - elc.dielectric_layers_self_forces(solver, particles); + elc.dielectric_layers_self_forces(solver, box_geo, particles); } else { solver.charge_assign(particles); } diff --git a/src/core/electrostatics/elc.hpp b/src/core/electrostatics/elc.hpp index 197aa5fe0f4..d0071881ea3 100644 --- a/src/core/electrostatics/elc.hpp +++ b/src/core/electrostatics/elc.hpp @@ -40,6 +40,7 @@ #include "electrostatics/p3m.hpp" #include "electrostatics/p3m_gpu.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" @@ -118,29 +119,31 @@ struct elc_data { /// pairwise contributions from the lowest and top layers template void dielectric_layers_contribution(CoulombP3M const &p3m, + BoxGeometry const &box_geo, Utils::Vector3d const &pos1, Utils::Vector3d const &pos2, double q1q2, Kernel &&kernel) const { if (pos1[2] < space_layer) { auto const q_eff = delta_mid_bot * q1q2; - auto const d = get_mi_vector(pos2, {pos1[0], pos1[1], -pos1[2]}); + auto const d = box_geo.get_mi_vector(pos2, {pos1[0], pos1[1], -pos1[2]}); kernel(q_eff, d); } if (pos1[2] > (box_h - space_layer)) { auto const q_eff = delta_mid_top * q1q2; - auto const l = 2. * box_h; - auto const d = get_mi_vector(pos2, {pos1[0], pos1[1], l - pos1[2]}); + auto const z = 2. * box_h - pos1[2]; + auto const d = box_geo.get_mi_vector(pos2, {pos1[0], pos1[1], z}); kernel(q_eff, d); } } /// self energies of top and bottom layers with their virtual images double dielectric_layers_self_energy(CoulombP3M const &p3m, + BoxGeometry const &box_geo, ParticleRange const &particles) const { auto energy = 0.; for (auto const &p : particles) { dielectric_layers_contribution( - p3m, p.pos(), p.pos(), Utils::sqr(p.q()), + p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), [&](double q1q2, Utils::Vector3d const &d) { energy += p3m.pair_energy(q1q2, d.norm()); }); @@ -150,19 +153,16 @@ struct elc_data { /// forces of particles in border layers with themselves void dielectric_layers_self_forces(CoulombP3M const &p3m, + BoxGeometry const &box_geo, ParticleRange const &particles) const { for (auto &p : particles) { dielectric_layers_contribution( - p3m, p.pos(), p.pos(), Utils::sqr(p.q()), + p3m, box_geo, p.pos(), p.pos(), Utils::sqr(p.q()), [&](double q1q2, Utils::Vector3d const &d) { p.force() += p3m.pair_force(q1q2, d, d.norm()); }); } } - -private: - Utils::Vector3d get_mi_vector(Utils::Vector3d const &a, - Utils::Vector3d const &b) const; }; struct ElectrostaticLayerCorrection @@ -252,23 +252,23 @@ struct ElectrostaticLayerCorrection } /** @brief Calculate short-range pair energy correction. */ - double pair_energy_correction(double q1q2, Particle const &p1, - Particle const &p2) const { + double pair_energy_correction(Particle const &p1, Particle const &p2, + double q1q2, BoxGeometry const &box_geo) const { double energy = 0.; if (elc.dielectric_contrast_on) { energy = std::visit( - [this, &p1, &p2, q1q2](auto &p3m_ptr) { + [this, &p1, &p2, q1q2, &box_geo](auto &p3m_ptr) { auto const &pos1 = p1.pos(); auto const &pos2 = p2.pos(); auto const &p3m = *p3m_ptr; auto energy = 0.; elc.dielectric_layers_contribution( - p3m, pos1, pos2, q1q2, + p3m, box_geo, pos1, pos2, q1q2, [&](double q_eff, Utils::Vector3d const &d) { energy += p3m.pair_energy(q_eff, d.norm()); }); elc.dielectric_layers_contribution( - p3m, pos2, pos1, q1q2, + p3m, box_geo, pos2, pos1, q1q2, [&](double q_eff, Utils::Vector3d const &d) { energy += p3m.pair_energy(q_eff, d.norm()); }); @@ -280,21 +280,21 @@ struct ElectrostaticLayerCorrection } /** @brief Add short-range pair force corrections. */ - void add_pair_force_corrections(Particle &p1, Particle &p2, - double q1q2) const { + void add_pair_force_corrections(Particle &p1, Particle &p2, double q1q2, + BoxGeometry const &box_geo) const { if (elc.dielectric_contrast_on) { std::visit( - [this, &p1, &p2, q1q2](auto &p3m_ptr) { + [this, &p1, &p2, q1q2, &box_geo](auto &p3m_ptr) { auto const &pos1 = p1.pos(); auto const &pos2 = p2.pos(); auto const &p3m = *p3m_ptr; elc.dielectric_layers_contribution( - p3m, pos1, pos2, q1q2, + p3m, box_geo, pos1, pos2, q1q2, [&](double q_eff, Utils::Vector3d const &d) { p1.force() += p3m.pair_force(q_eff, d, d.norm()); }); elc.dielectric_layers_contribution( - p3m, pos2, pos1, q1q2, + p3m, box_geo, pos2, pos1, q1q2, [&](double q_eff, Utils::Vector3d const &d) { p2.force() += p3m.pair_force(q_eff, d, d.norm()); }); diff --git a/src/core/electrostatics/mmm1d.cpp b/src/core/electrostatics/mmm1d.cpp index 8bb7aa1782f..b3ebd37355f 100644 --- a/src/core/electrostatics/mmm1d.cpp +++ b/src/core/electrostatics/mmm1d.cpp @@ -29,12 +29,15 @@ #include "electrostatics/mmm-common.hpp" #include "electrostatics/mmm-modpsi.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "Particle.hpp" #include "cell_system/CellStructureType.hpp" #include "errorhandling.hpp" #include "event.hpp" #include "grid.hpp" #include "specfunc.hpp" +#include "system/System.hpp" #include "tuning.hpp" #include @@ -56,6 +59,7 @@ #endif static double far_error(int P, double minrad) { + auto const &box_geo = *System::get_system().box_geo; auto const wavenumber = 2. * Utils::pi() * box_geo.length_inv()[2]; // this uses an upper bound to all force components and the potential auto const rhores = wavenumber * minrad; @@ -67,6 +71,7 @@ static double far_error(int P, double minrad) { static double determine_minrad(double maxPWerror, int P) { // bisection to search for where the error is maxPWerror auto constexpr min_rad = 0.01; + auto const &box_geo = *System::get_system().box_geo; auto const rgranularity = min_rad * box_geo.length()[2]; auto rmin = rgranularity; auto rmax = std::min(box_geo.length()[0], box_geo.length()[1]); @@ -138,12 +143,14 @@ CoulombMMM1D::CoulombMMM1D(double prefactor, double maxPWerror, } void CoulombMMM1D::sanity_checks_periodicity() const { + auto const &box_geo = *System::get_system().box_geo; if (box_geo.periodic(0) || box_geo.periodic(1) || !box_geo.periodic(2)) { throw std::runtime_error("MMM1D requires periodicity (False, False, True)"); } } void CoulombMMM1D::sanity_checks_cell_structure() const { + auto const &local_geo = *System::get_system().local_geo; if (local_geo.cell_structure_type() != CellStructureType::CELL_STRUCTURE_NSQUARE) { throw std::runtime_error("MMM1D requires the N-square cellsystem"); @@ -151,6 +158,8 @@ void CoulombMMM1D::sanity_checks_cell_structure() const { } void CoulombMMM1D::recalc_boxl_parameters() { + auto const &box_geo = *System::get_system().box_geo; + if (far_switch_radius_sq >= Utils::sqr(box_geo.length()[2])) far_switch_radius_sq = 0.8 * Utils::sqr(box_geo.length()[2]); @@ -165,6 +174,7 @@ void CoulombMMM1D::recalc_boxl_parameters() { Utils::Vector3d CoulombMMM1D::pair_force(double q1q2, Utils::Vector3d const &d, double dist) const { auto constexpr c_2pi = 2. * Utils::pi(); + auto const &box_geo = *System::get_system().box_geo; auto const n_modPsi = static_cast(modPsi.size()) >> 1; auto const rxy2 = d[0] * d[0] + d[1] * d[1]; auto const rxy2_d = rxy2 * uz2; @@ -258,6 +268,7 @@ double CoulombMMM1D::pair_energy(double const q1q2, Utils::Vector3d const &d, return 0.; auto constexpr c_2pi = 2. * Utils::pi(); + auto const &box_geo = *System::get_system().box_geo; auto const n_modPsi = static_cast(modPsi.size()) >> 1; auto const rxy2 = d[0] * d[0] + d[1] * d[1]; auto const rxy2_d = rxy2 * uz2; @@ -321,6 +332,7 @@ void CoulombMMM1D::tune() { recalc_boxl_parameters(); if (far_switch_radius_sq < 0.) { + auto const &box_geo = *System::get_system().box_geo; auto const maxrad = box_geo.length()[2]; auto min_time = std::numeric_limits::infinity(); auto min_rad = -1.; diff --git a/src/core/electrostatics/mmm1d_gpu.cpp b/src/core/electrostatics/mmm1d_gpu.cpp index 9d73f35c67e..e16943e1f99 100644 --- a/src/core/electrostatics/mmm1d_gpu.cpp +++ b/src/core/electrostatics/mmm1d_gpu.cpp @@ -23,6 +23,8 @@ #include "electrostatics/mmm1d_gpu.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "cell_system/CellStructureType.hpp" #include "communication.hpp" #include "event.hpp" @@ -38,6 +40,7 @@ CoulombMMM1DGpu::CoulombMMM1DGpu(double prefactor, double maxPWerror, far_switch_radius_sq{-1.}, bessel_cutoff{bessel_cutoff}, m_is_tuned{ false} { set_prefactor(prefactor); + auto const &box_geo = *System::get_system().box_geo; if (maxPWerror <= 0.) { throw std::domain_error("Parameter 'maxPWerror' must be > 0"); } @@ -62,12 +65,14 @@ CoulombMMM1DGpu::CoulombMMM1DGpu(double prefactor, double maxPWerror, } void CoulombMMM1DGpu::sanity_checks_periodicity() const { + auto const &box_geo = *System::get_system().box_geo; if (box_geo.periodic(0) || box_geo.periodic(1) || !box_geo.periodic(2)) { throw std::runtime_error("MMM1D requires periodicity (False, False, True)"); } } void CoulombMMM1DGpu::sanity_checks_cell_structure() const { + auto const &local_geo = *System::get_system().local_geo; if (local_geo.cell_structure_type() != CellStructureType::CELL_STRUCTURE_NSQUARE) { throw std::runtime_error("MMM1D requires the N-square cellsystem"); diff --git a/src/core/electrostatics/p3m.cpp b/src/core/electrostatics/p3m.cpp index 5538f5d2a68..d61f5744d1e 100644 --- a/src/core/electrostatics/p3m.cpp +++ b/src/core/electrostatics/p3m.cpp @@ -38,6 +38,8 @@ #include "p3m/fft.hpp" #include "p3m/influence_function.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "actor/visitors.hpp" @@ -107,6 +109,7 @@ void CoulombP3M::count_charged_particles() { * @cite deserno98a @cite deserno98b. */ void CoulombP3M::calc_influence_function_force() { + auto const &box_geo = *System::get_system().box_geo; auto const start = Utils::Vector3i{p3m.fft.plan[3].start}; auto const size = Utils::Vector3i{p3m.fft.plan[3].new_mesh}; @@ -118,6 +121,7 @@ void CoulombP3M::calc_influence_function_force() { * self energy correction. */ void CoulombP3M::calc_influence_function_energy() { + auto const &box_geo = *System::get_system().box_geo; auto const start = Utils::Vector3i{p3m.fft.plan[3].start}; auto const size = Utils::Vector3i{p3m.fft.plan[3].new_mesh}; @@ -170,6 +174,7 @@ static void p3m_tune_aliasing_sums(int nx, int ny, int nz, */ static double p3m_real_space_error(double pref, double r_cut_iL, int n_c_part, double sum_q2, double alpha_L) { + auto const &box_geo = *System::get_system().box_geo; return (2. * pref * sum_q2 * exp(-Utils::sqr(r_cut_iL * alpha_L))) / sqrt(static_cast(n_c_part) * r_cut_iL * box_geo.length()[0] * box_geo.volume()); @@ -190,6 +195,7 @@ static double p3m_real_space_error(double pref, double r_cut_iL, int n_c_part, static double p3m_k_space_error(double pref, Utils::Vector3i const &mesh, int cao, int n_c_part, double sum_q2, double alpha_L) { + auto const &box_geo = *System::get_system().box_geo; auto const mesh_i = Utils::hadamard_division(Utils::Vector3d::broadcast(1.), mesh); auto const alpha_L_i = 1. / alpha_L; @@ -226,6 +232,7 @@ static double p3m_k_space_error(double pref, Utils::Vector3i const &mesh, static double p3mgpu_k_space_error(double prefactor, Utils::Vector3i const &mesh, int cao, int npart, double sum_q2, double alpha_L) { + auto const &box_geo = *System::get_system().box_geo; auto ks_err = 0.; if (this_node == 0) { ks_err = p3m_k_space_error_gpu(prefactor, mesh.data(), cao, npart, sum_q2, @@ -241,12 +248,16 @@ void CoulombP3M::init() { assert(p3m.params.cao >= 1 and p3m.params.cao <= 7); assert(p3m.params.alpha > 0.); + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; + p3m.params.cao3 = Utils::int_pow<3>(p3m.params.cao); p3m.params.recalc_a_ai_cao_cut(box_geo.length()); sanity_checks(); - auto const &solver = System::get_system().coulomb.impl->solver; + auto const &solver = system.coulomb.impl->solver; double elc_layer = 0.; if (auto actor = get_actor_by_type(solver)) { elc_layer = actor->elc.space_layer; @@ -398,6 +409,8 @@ auto calc_dipole_moment(boost::mpi::communicator const &comm, Utils::Vector9d CoulombP3M::p3m_calc_kspace_pressure_tensor() { using namespace detail::FFT_indexing; + auto const &box_geo = *System::get_system().box_geo; + Utils::Vector9d node_k_space_pressure_tensor{}; if (p3m.sum_q2 > 0.) { @@ -453,6 +466,8 @@ Utils::Vector9d CoulombP3M::p3m_calc_kspace_pressure_tensor() { double CoulombP3M::long_range_kernel(bool force_flag, bool energy_flag, ParticleRange const &particles) { + auto const &box_geo = *System::get_system().box_geo; + /* Gather information for FFT grid inside the nodes domain (inner local mesh) * and perform forward 3D FFT (Charge Assignment Mesh). */ p3m.sm.gather_grid(p3m.rs_mesh.data(), comm_cart, p3m.local_mesh.dim); @@ -568,8 +583,10 @@ class CoulombTuningAlgorithm : public TuningAlgorithm { void on_solver_change() const override { on_coulomb_change(); } void setup_logger(bool verbose) override { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; #ifdef CUDA - auto const &solver = System::get_system().coulomb.impl->solver; + auto const &solver = system.coulomb.impl->solver; auto const on_gpu = has_actor_of_type(solver); #else auto const on_gpu = false; @@ -629,6 +646,7 @@ class CoulombTuningAlgorithm : public TuningAlgorithm { } void determine_mesh_limits() override { + auto const &box_geo = *System::get_system().box_geo; auto const mesh_density = static_cast(p3m.params.mesh[0]) * box_geo.length_inv()[0]; @@ -661,10 +679,12 @@ class CoulombTuningAlgorithm : public TuningAlgorithm { } TuningAlgorithm::Parameters get_time() override { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &solver = system.coulomb.impl->solver; auto tuned_params = TuningAlgorithm::Parameters{}; auto time_best = time_sentinel; auto mesh_density = m_mesh_density_min; - auto const &solver = System::get_system().coulomb.impl->solver; while (mesh_density <= m_mesh_density_max) { auto trial_params = TuningAlgorithm::Parameters{}; if (m_tune_mesh) { @@ -708,6 +728,7 @@ class CoulombTuningAlgorithm : public TuningAlgorithm { }; void CoulombP3M::tune() { + auto const &box_geo = *System::get_system().box_geo; if (p3m.params.alpha_L == 0. and p3m.params.alpha != 0.) { p3m.params.alpha_L = p3m.params.alpha * box_geo.length()[0]; } @@ -740,7 +761,10 @@ void CoulombP3M::tune() { } void CoulombP3M::sanity_checks_boxl() const { - for (unsigned int i = 0; i < 3; i++) { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; + for (unsigned int i = 0u; i < 3u; i++) { /* check k-space cutoff */ if (p3m.params.cao_cut[i] >= box_geo.length_half()[i]) { std::stringstream msg; @@ -768,6 +792,7 @@ void CoulombP3M::sanity_checks_boxl() const { } void CoulombP3M::sanity_checks_periodicity() const { + auto const &box_geo = *System::get_system().box_geo; if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) { throw std::runtime_error( "CoulombP3M: requires periodicity (True, True, True)"); @@ -775,6 +800,7 @@ void CoulombP3M::sanity_checks_periodicity() const { } void CoulombP3M::sanity_checks_cell_structure() const { + auto const &local_geo = *System::get_system().local_geo; if (local_geo.cell_structure_type() != CellStructureType::CELL_STRUCTURE_REGULAR && local_geo.cell_structure_type() != @@ -799,6 +825,7 @@ void CoulombP3M::sanity_checks_node_grid() const { } void CoulombP3M::scaleby_box_l() { + auto const &box_geo = *System::get_system().box_geo; p3m.params.r_cut = p3m.params.r_cut_iL * box_geo.length()[0]; p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0]; p3m.params.recalc_a_ai_cao_cut(box_geo.length()); diff --git a/src/core/electrostatics/scafacos_impl.cpp b/src/core/electrostatics/scafacos_impl.cpp index ed96208788f..fe540bff7f8 100644 --- a/src/core/electrostatics/scafacos_impl.cpp +++ b/src/core/electrostatics/scafacos_impl.cpp @@ -26,6 +26,8 @@ #include "electrostatics/scafacos.hpp" #include "electrostatics/scafacos_impl.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" #include "event.hpp" @@ -52,7 +54,9 @@ make_coulomb_scafacos(std::string const &method, } void CoulombScafacosImpl::update_particle_data() { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; positions.clear(); charges.clear(); @@ -90,6 +94,9 @@ double CoulombScafacosImpl::time_r_cut(double r_cut) { void CoulombScafacosImpl::tune_r_cut() { auto constexpr convergence_threshold = 1e-3; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; auto const min_box_l = *boost::min_element(box_geo.length()); auto const min_local_box_l = *boost::min_element(local_geo.length()); diff --git a/src/core/energy.cpp b/src/core/energy.cpp index 69fdc93f1e2..b91fcfa7645 100644 --- a/src/core/energy.cpp +++ b/src/core/energy.cpp @@ -22,6 +22,7 @@ * Energy calculation. */ +#include "BoxGeometry.hpp" #include "Observable_stat.hpp" #include "cells.hpp" #include "communication.hpp" @@ -60,6 +61,7 @@ std::shared_ptr calculate_energy() { #endif on_observable_calc(); + auto const &box_geo = *system.box_geo; auto const local_parts = cell_structure.local_particles(); auto const n_nodes = ::communicator.size; @@ -71,11 +73,11 @@ std::shared_ptr calculate_energy() { auto const dipoles_kernel = system.dipoles.pair_energy_kernel(); short_range_loop( - [&obs_energy, coulomb_kernel_ptr = get_ptr(coulomb_kernel)]( + [&obs_energy, coulomb_kernel_ptr = get_ptr(coulomb_kernel), &box_geo]( Particle const &p1, int bond_id, Utils::Span partners) { auto const &iaparams = *bonded_ia_params.at(bond_id); - auto const result = - calc_bonded_energy(iaparams, p1, partners, coulomb_kernel_ptr); + auto const result = calc_bonded_energy(iaparams, p1, partners, box_geo, + coulomb_kernel_ptr); if (result) { obs_energy.bonded_contribution(bond_id)[0] += result.get(); return false; @@ -101,7 +103,8 @@ std::shared_ptr calculate_energy() { obs_energy.dipolar[1] = system.dipoles.calc_energy_long_range(local_parts); #endif - Constraints::constraints.add_energy(local_parts, get_sim_time(), obs_energy); + Constraints::constraints.add_energy(box_geo, local_parts, get_sim_time(), + obs_energy); #if defined(CUDA) and (defined(ELECTROSTATICS) or defined(DIPOLES)) auto const energy_host = gpu_particle_data.copy_energy_to_host(); diff --git a/src/core/energy_inline.hpp b/src/core/energy_inline.hpp index eaa24054fd3..f70bf8d578e 100644 --- a/src/core/energy_inline.hpp +++ b/src/core/energy_inline.hpp @@ -49,6 +49,7 @@ #include "nonbonded_interactions/thole.hpp" #include "nonbonded_interactions/wca.hpp" +#include "BoxGeometry.hpp" #include "Observable_stat.hpp" #include "Particle.hpp" #include "bond_error.hpp" @@ -202,7 +203,7 @@ inline void add_non_bonded_pair_energy( inline boost::optional calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, - Utils::Span partners, + Utils::Span partners, BoxGeometry const &box_geo, Coulomb::ShortRangeEnergyKernel::kernel_type const *kernel) { auto const n_partners = static_cast(partners.size()); @@ -246,16 +247,16 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, } // 1 partner if (n_partners == 2) { if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p1.pos(), p2->pos(), p3->pos()); + return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p1.pos(), p2->pos(), p3->pos()); + return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p1.pos(), p2->pos(), p3->pos()); + return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p1.pos(), p2->pos(), p3->pos()); + return iap->energy(box_geo, p1.pos(), p2->pos(), p3->pos()); } if (boost::get(&iaparams)) { runtimeWarningMsg() << "Unsupported bond type " + @@ -267,10 +268,10 @@ calc_bonded_energy(Bonded_IA_Parameters const &iaparams, Particle const &p1, } // 2 partners if (n_partners == 3) { if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p2->pos(), p1.pos(), p3->pos(), p4->pos()); + return iap->energy(box_geo, p2->pos(), p1.pos(), p3->pos(), p4->pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->energy(p2->pos(), p1.pos(), p3->pos(), p4->pos()); + return iap->energy(box_geo, p2->pos(), p1.pos(), p3->pos(), p4->pos()); } if (boost::get(&iaparams)) { runtimeWarningMsg() << "Unsupported bond type " + diff --git a/src/core/event.cpp b/src/core/event.cpp index f322a06ac01..6a8e49fe7be 100644 --- a/src/core/event.cpp +++ b/src/core/event.cpp @@ -98,7 +98,7 @@ void on_integration_start(double time_step) { } #ifdef NPT - npt_ensemble_init(box_geo); + npt_ensemble_init(*system.box_geo); #endif partCfg().invalidate(); @@ -220,11 +220,15 @@ void on_constraint_change() { recalc_forces = true; } void on_lb_boundary_conditions_change() { recalc_forces = true; } +static void update_local_geo(System::System &system) { + *system.local_geo = + regular_decomposition(*system.box_geo, ::communicator.calc_node_index(), + ::communicator.node_grid); +} + void on_boxl_change(bool skip_method_adaption) { auto &system = System::get_system(); - grid_changed_box_l(box_geo); - /* Electrostatics cutoffs mostly depend on the system size, - * therefore recalculate them. */ + update_local_geo(system); cells_re_init(*system.cell_structure); /* Now give methods a chance to react to the change in box length */ @@ -278,6 +282,7 @@ void on_periodicity_change() { #ifdef STOKESIAN_DYNAMICS if (integ_switch == INTEG_METHOD_SD) { + auto const &box_geo = *System::get_system().box_geo; if (box_geo.periodic(0) || box_geo.periodic(1) || box_geo.periodic(2)) runtimeErrorMsg() << "Stokesian Dynamics requires periodicity " << "(False, False, False)\n"; @@ -304,9 +309,10 @@ void on_timestep_change() { void on_forcecap_change() { recalc_forces = true; } void on_node_grid_change() { - grid_changed_node_grid(); auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; + ::communicator.init_comm_cart(); + update_local_geo(system); system.lb.on_node_grid_change(); system.ek.on_node_grid_change(); #ifdef ELECTROSTATICS diff --git a/src/core/forces.cpp b/src/core/forces.cpp index cc1e99fc8f3..fa6a783ec13 100644 --- a/src/core/forces.cpp +++ b/src/core/forces.cpp @@ -24,6 +24,7 @@ * The corresponding header file is forces.hpp. */ +#include "BoxGeometry.hpp" #include "bond_breakage/bond_breakage.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" @@ -153,6 +154,7 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { #endif auto &system = System::get_system(); + auto const &box_geo = *system.box_geo; auto const n_nodes = ::communicator.size; #ifdef CUDA #ifdef CALIPER @@ -199,9 +201,10 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { #endif short_range_loop( - [coulomb_kernel_ptr = get_ptr(coulomb_kernel)]( - Particle &p1, int bond_id, Utils::Span partners) { - return add_bonded_force(p1, bond_id, partners, coulomb_kernel_ptr); + [coulomb_kernel_ptr = get_ptr(coulomb_kernel), + &box_geo](Particle &p1, int bond_id, Utils::Span partners) { + return add_bonded_force(p1, bond_id, partners, box_geo, + coulomb_kernel_ptr); }, [coulomb_kernel_ptr = get_ptr(coulomb_kernel), dipoles_kernel_ptr = get_ptr(dipoles_kernel), @@ -219,27 +222,28 @@ void force_calc(CellStructure &cell_structure, double time_step, double kT) { VerletCriterion<>{skin, interaction_range(), coulomb_cutoff, dipole_cutoff, collision_detection_cutoff()}); - Constraints::constraints.add_forces(particles, get_sim_time()); + Constraints::constraints.add_forces(box_geo, particles, get_sim_time()); if (max_oif_objects) { // There are two global quantities that need to be evaluated: // object's surface and object's volume. for (int i = 0; i < max_oif_objects; i++) { auto const area_volume = boost::mpi::all_reduce( - comm_cart, calc_oif_global(i, cell_structure), std::plus<>()); + comm_cart, calc_oif_global(i, box_geo, cell_structure), + std::plus<>()); auto const oif_part_area = std::abs(area_volume[0]); auto const oif_part_vol = std::abs(area_volume[1]); if (oif_part_area < 1e-100 and oif_part_vol < 1e-100) { break; } - add_oif_global_forces(area_volume, i, cell_structure); + add_oif_global_forces(area_volume, i, box_geo, cell_structure); } } // Must be done here. Forces need to be ghost-communicated immersed_boundaries.volume_conservation(cell_structure); - if (System::get_system().lb.is_solver_set()) { + if (system.lb.is_solver_set()) { LB::couple_particles(thermo_virtual, particles, ghost_particles, time_step); } diff --git a/src/core/forces_inline.hpp b/src/core/forces_inline.hpp index e18a18c7868..387ce2fe3d3 100644 --- a/src/core/forces_inline.hpp +++ b/src/core/forces_inline.hpp @@ -28,6 +28,7 @@ #include "forces.hpp" +#include "BoxGeometry.hpp" #include "actor/visitors.hpp" #include "bond_breakage/bond_breakage.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" @@ -324,6 +325,7 @@ inline boost::optional calc_bond_pair_force( inline bool add_bonded_two_body_force( Bonded_IA_Parameters const &iaparams, Particle &p1, Particle &p2, + BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { auto const dx = box_geo.get_mi_vector(p1.pos(), p2.pos()); @@ -355,35 +357,37 @@ inline bool add_bonded_two_body_force( inline boost::optional< std::tuple> calc_bonded_three_body_force(Bonded_IA_Parameters const &iaparams, - Particle const &p1, Particle const &p2, - Particle const &p3) { + BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3) { if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p1.pos(), p2.pos(), p3.pos()); + return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p1.pos(), p2.pos(), p3.pos()); + return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); } if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p1.pos(), p2.pos(), p3.pos()); + return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); } #ifdef TABULATED if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p1.pos(), p2.pos(), p3.pos()); + return iap->forces(box_geo, p1.pos(), p2.pos(), p3.pos()); } #endif if (auto const *iap = boost::get(&iaparams)) { - return iap->calc_forces(p1, p2, p3); + return iap->calc_forces(box_geo, p1, p2, p3); } throw BondUnknownTypeError(); } inline bool add_bonded_three_body_force(Bonded_IA_Parameters const &iaparams, + BoxGeometry const &box_geo, Particle &p1, Particle &p2, Particle &p3) { if (boost::get(&iaparams)) { return false; } - auto const result = calc_bonded_three_body_force(iaparams, p1, p2, p3); + auto const result = + calc_bonded_three_body_force(iaparams, box_geo, p1, p2, p3); if (result) { auto const &forces = result.get(); @@ -399,29 +403,32 @@ inline bool add_bonded_three_body_force(Bonded_IA_Parameters const &iaparams, inline boost::optional> calc_bonded_four_body_force(Bonded_IA_Parameters const &iaparams, - Particle const &p1, Particle const &p2, - Particle const &p3, Particle const &p4) { + BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3, + Particle const &p4) { if (auto const *iap = boost::get(&iaparams)) { - return iap->calc_forces(p1, p2, p3, p4); + return iap->calc_forces(box_geo, p1, p2, p3, p4); } if (auto const *iap = boost::get(&iaparams)) { - return iap->calc_forces(p1, p2, p3, p4); + return iap->calc_forces(box_geo, p1, p2, p3, p4); } if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p2.pos(), p1.pos(), p3.pos(), p4.pos()); + return iap->forces(box_geo, p2.pos(), p1.pos(), p3.pos(), p4.pos()); } #ifdef TABULATED if (auto const *iap = boost::get(&iaparams)) { - return iap->forces(p2.pos(), p1.pos(), p3.pos(), p4.pos()); + return iap->forces(box_geo, p2.pos(), p1.pos(), p3.pos(), p4.pos()); } #endif throw BondUnknownTypeError(); } inline bool add_bonded_four_body_force(Bonded_IA_Parameters const &iaparams, - Particle &p1, Particle &p2, Particle &p3, + BoxGeometry const &box_geo, Particle &p1, + Particle &p2, Particle &p3, Particle &p4) { - auto const result = calc_bonded_four_body_force(iaparams, p1, p2, p3, p4); + auto const result = + calc_bonded_four_body_force(iaparams, box_geo, p1, p2, p3, p4); if (result) { auto const &forces = result.get(); @@ -438,6 +445,7 @@ inline bool add_bonded_four_body_force(Bonded_IA_Parameters const &iaparams, inline bool add_bonded_force(Particle &p1, int bond_id, Utils::Span partners, + BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { // Consider for bond breakage @@ -463,13 +471,14 @@ add_bonded_force(Particle &p1, int bond_id, Utils::Span partners, case 0: return false; case 1: - return add_bonded_two_body_force(iaparams, p1, *partners[0], kernel); + return add_bonded_two_body_force(iaparams, p1, *partners[0], box_geo, + kernel); case 2: - return add_bonded_three_body_force(iaparams, p1, *partners[0], + return add_bonded_three_body_force(iaparams, box_geo, p1, *partners[0], *partners[1]); case 3: - return add_bonded_four_body_force(iaparams, p1, *partners[0], *partners[1], - *partners[2]); + return add_bonded_four_body_force(iaparams, box_geo, p1, *partners[0], + *partners[1], *partners[2]); default: throw BondInvalidSizeError{number_of_partners(iaparams)}; } diff --git a/src/core/galilei/Galilei.cpp b/src/core/galilei/Galilei.cpp index b446107ec23..c5c88c7e047 100644 --- a/src/core/galilei/Galilei.cpp +++ b/src/core/galilei/Galilei.cpp @@ -21,6 +21,7 @@ #include "galilei/Galilei.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "cell_system/CellStructure.hpp" @@ -69,14 +70,15 @@ void Galilei::kill_particle_forces(bool torque) const { } Utils::Vector3d Galilei::calc_system_cms_position() const { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_l = system.box_geo->length(); + auto &cell_structure = *system.cell_structure; auto total_mass = 0.; auto cms_pos = Utils::Vector3d{}; for (auto const &p : cell_structure.local_particles()) { if (not p.is_virtual()) { total_mass += p.mass(); - cms_pos += p.mass() * - unfolded_position(p.pos(), p.image_box(), box_geo.length()); + cms_pos += p.mass() * unfolded_position(p.pos(), p.image_box(), box_l); } } total_mass = boost::mpi::all_reduce(comm_cart, total_mass, std::plus<>()); diff --git a/src/core/ghosts.cpp b/src/core/ghosts.cpp index 6a22941130b..28a4d1b52ff 100644 --- a/src/core/ghosts.cpp +++ b/src/core/ghosts.cpp @@ -30,8 +30,10 @@ */ #include "ghosts.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "grid.hpp" +#include "system/System.hpp" #include #include @@ -122,6 +124,7 @@ template static void serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts, ReductionPolicy policy, SerializationDirection direction, + BoxGeometry const &box_geo, Utils::Vector3d const *ghost_shift) { if (data_parts & GHOSTTRANS_PROPRTS) { ar &p.id() & p.mol_id() & p.type(); @@ -171,7 +174,7 @@ serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts, /* ok, this is not nice, but perhaps fast */ auto pos = p.pos() + *ghost_shift; auto img = p.image_box(); - fold_position(pos, img, ::box_geo); + fold_position(pos, img, box_geo); ar &pos; ar &img; } else { @@ -225,16 +228,18 @@ serialize_and_reduce(Archive &ar, Particle &p, unsigned int data_parts, #endif } -static std::size_t calc_transmit_size(unsigned data_parts) { +static auto calc_transmit_size(BoxGeometry const &box_geo, + unsigned data_parts) { SerializationSizeCalculator sizeof_archive; Particle p{}; serialize_and_reduce(sizeof_archive, p, data_parts, ReductionPolicy::MOVE, - SerializationDirection::SAVE, nullptr); + SerializationDirection::SAVE, box_geo, nullptr); return sizeof_archive.size(); } -static std::size_t calc_transmit_size(const GhostCommunication &ghost_comm, - unsigned int data_parts) { +static auto calc_transmit_size(GhostCommunication const &ghost_comm, + BoxGeometry const &box_geo, + unsigned int data_parts) { if (data_parts & GHOSTTRANS_PARTNUM) return sizeof(unsigned int) * ghost_comm.part_lists.size(); @@ -242,14 +247,16 @@ static std::size_t calc_transmit_size(const GhostCommunication &ghost_comm, ghost_comm.part_lists, std::size_t{0}, [](std::size_t sum, auto part_list) { return sum + part_list->size(); }); - return n_part * calc_transmit_size(data_parts); + return n_part * calc_transmit_size(box_geo, data_parts); } static void prepare_send_buffer(CommBuf &send_buffer, const GhostCommunication &ghost_comm, unsigned int data_parts) { + auto const &box_geo = *System::get_system().box_geo; + /* reallocate send buffer */ - send_buffer.resize(calc_transmit_size(ghost_comm, data_parts)); + send_buffer.resize(calc_transmit_size(ghost_comm, box_geo, data_parts)); send_buffer.bonds().clear(); auto archiver = Utils::MemcpyOArchive{Utils::make_span(send_buffer)}; @@ -269,7 +276,8 @@ static void prepare_send_buffer(CommBuf &send_buffer, } else { for (auto &p : *part_list) { serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE, - SerializationDirection::SAVE, &ghost_comm.shift); + SerializationDirection::SAVE, box_geo, + &ghost_comm.shift); if (data_parts & GHOSTTRANS_BONDS) { bond_archiver << p.bonds(); } @@ -293,8 +301,9 @@ static void prepare_ghost_cell(ParticleList *cell, std::size_t size) { static void prepare_recv_buffer(CommBuf &recv_buffer, const GhostCommunication &ghost_comm, unsigned int data_parts) { + auto const &box_geo = *System::get_system().box_geo; /* reallocate recv buffer */ - recv_buffer.resize(calc_transmit_size(ghost_comm, data_parts)); + recv_buffer.resize(calc_transmit_size(ghost_comm, box_geo, data_parts)); /* clear bond buffer */ recv_buffer.bonds().clear(); } @@ -304,6 +313,7 @@ static void put_recv_buffer(CommBuf &recv_buffer, unsigned int data_parts) { /* put back data */ auto archiver = Utils::MemcpyIArchive{Utils::make_span(recv_buffer)}; + auto const &box_geo = *System::get_system().box_geo; if (data_parts & GHOSTTRANS_PARTNUM) { for (auto part_list : ghost_comm.part_lists) { @@ -315,7 +325,7 @@ static void put_recv_buffer(CommBuf &recv_buffer, for (auto part_list : ghost_comm.part_lists) { for (auto &p : *part_list) { serialize_and_reduce(archiver, p, data_parts, ReductionPolicy::MOVE, - SerializationDirection::LOAD, nullptr); + SerializationDirection::LOAD, box_geo, nullptr); } } if (data_parts & GHOSTTRANS_BONDS) { @@ -368,9 +378,10 @@ static void add_forces_from_recv_buffer(CommBuf &recv_buffer, static void cell_cell_transfer(const GhostCommunication &ghost_comm, unsigned int data_parts) { + auto const &box_geo = *System::get_system().box_geo; CommBuf buffer; if (!(data_parts & GHOSTTRANS_PARTNUM)) { - buffer.resize(calc_transmit_size(data_parts)); + buffer.resize(calc_transmit_size(box_geo, data_parts)); } /* transfer data */ auto const offset = ghost_comm.part_lists.size() / 2; @@ -391,9 +402,10 @@ static void cell_cell_transfer(const GhostCommunication &ghost_comm, auto &p1 = src_part.begin()[i]; auto &p2 = dst_part.begin()[i]; serialize_and_reduce(ar_out, p1, data_parts, ReductionPolicy::UPDATE, - SerializationDirection::SAVE, &ghost_comm.shift); + SerializationDirection::SAVE, box_geo, + &ghost_comm.shift); serialize_and_reduce(ar_in, p2, data_parts, ReductionPolicy::UPDATE, - SerializationDirection::LOAD, nullptr); + SerializationDirection::LOAD, box_geo, nullptr); if (data_parts & GHOSTTRANS_BONDS) { p2.bonds() = p1.bonds(); } diff --git a/src/core/grid.cpp b/src/core/grid.cpp index 69c840de124..20cdb8849cf 100644 --- a/src/core/grid.cpp +++ b/src/core/grid.cpp @@ -26,17 +26,13 @@ #include "grid.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "communication.hpp" #include "event.hpp" #include -#include -#include - -BoxGeometry box_geo; -LocalBox local_geo; - LocalBox regular_decomposition(BoxGeometry const &box, Utils::Vector3i const &node_index, Utils::Vector3i const &node_grid) { @@ -56,25 +52,7 @@ LocalBox regular_decomposition(BoxGeometry const &box, CellStructureType::CELL_STRUCTURE_REGULAR}; } -void grid_changed_box_l(const BoxGeometry &box) { - local_geo = regular_decomposition(box, ::communicator.calc_node_index(), - ::communicator.node_grid); -} - -void grid_changed_node_grid(bool update_box_geo) { - ::communicator.init_comm_cart(); - - if (update_box_geo) { - grid_changed_box_l(box_geo); - } -} - void set_node_grid(Utils::Vector3i const &value) { ::communicator.node_grid = value; on_node_grid_change(); } - -void set_box_length(Utils::Vector3d const &value) { - ::box_geo.set_length(value); - on_boxl_change(); -} diff --git a/src/core/grid.hpp b/src/core/grid.hpp index e255d69f77d..d65bf7b4816 100644 --- a/src/core/grid.hpp +++ b/src/core/grid.hpp @@ -32,18 +32,10 @@ */ #include "BoxGeometry.hpp" - #include "LocalBox.hpp" #include -extern BoxGeometry box_geo; -extern LocalBox local_geo; - -void grid_changed_node_grid(bool update_box_geo = true); - -void grid_changed_box_l(const BoxGeometry &box); - /** @brief Calculate image box shift vector. * @param image_box image box offset * @param box box parameters (side lengths) @@ -79,4 +71,3 @@ LocalBox regular_decomposition(BoxGeometry const &box, Utils::Vector3i const &node_grid); void set_node_grid(Utils::Vector3i const &value); -void set_box_length(Utils::Vector3d const &value); diff --git a/src/core/immersed_boundary/ImmersedBoundaries.cpp b/src/core/immersed_boundary/ImmersedBoundaries.cpp index 1a5a09c196d..b2317f0fe6a 100644 --- a/src/core/immersed_boundary/ImmersedBoundaries.cpp +++ b/src/core/immersed_boundary/ImmersedBoundaries.cpp @@ -25,6 +25,7 @@ #include "communication.hpp" #include "grid.hpp" #include "ibm_volcons.hpp" +#include "system/System.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" @@ -99,12 +100,14 @@ void ImmersedBoundaries::calc_volumes(CellStructure &cs) { if (!BoundariesFound) return; + auto const &box_geo = *System::get_system().box_geo; + // Partial volumes for each soft particle, to be summed up std::vector tempVol(VolumesCurrent.size()); // Loop over all particles on local node - cs.bond_loop([&tempVol](Particle &p1, int bond_id, - Utils::Span partners) { + cs.bond_loop([&tempVol, &box_geo](Particle &p1, int bond_id, + Utils::Span partners) { auto const vol_cons_params = vol_cons_parameters(p1); if (vol_cons_params && @@ -156,8 +159,10 @@ void ImmersedBoundaries::calc_volume_force(CellStructure &cs) { if (!BoundariesFound) return; - cs.bond_loop([this](Particle &p1, int bond_id, - Utils::Span partners) { + auto const &box_geo = *System::get_system().box_geo; + + cs.bond_loop([this, &box_geo](Particle &p1, int bond_id, + Utils::Span partners) { if (boost::get(bonded_ia_params.at(bond_id).get()) != nullptr) { // Check if particle has an IBM Triel bonded interaction and an // IBM VolCons bonded interaction. Basically this loops over all diff --git a/src/core/immersed_boundary/ibm_tribend.cpp b/src/core/immersed_boundary/ibm_tribend.cpp index 4d627ed9373..9bb6ee67241 100644 --- a/src/core/immersed_boundary/ibm_tribend.cpp +++ b/src/core/immersed_boundary/ibm_tribend.cpp @@ -22,6 +22,7 @@ #include "BoxGeometry.hpp" #include "grid.hpp" #include "ibm_common.hpp" +#include "system/System.hpp" #include @@ -30,8 +31,9 @@ #include std::tuple -IBMTribend::calc_forces(Particle const &p1, Particle const &p2, - Particle const &p3, Particle const &p4) const { +IBMTribend::calc_forces(BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3, + Particle const &p4) const { // Get vectors making up the two triangles auto const dx1 = box_geo.get_mi_vector(p1.pos(), p3.pos()); @@ -94,6 +96,8 @@ IBMTribend::calc_forces(Particle const &p1, Particle const &p2, IBMTribend::IBMTribend(const int ind1, const int ind2, const int ind3, const int ind4, const double kb, const bool flat) { + auto const &box_geo = *System::get_system().box_geo; + // Compute theta0 if (flat) { theta0 = 0; diff --git a/src/core/immersed_boundary/ibm_tribend.hpp b/src/core/immersed_boundary/ibm_tribend.hpp index 1445f84702b..3a013d040db 100644 --- a/src/core/immersed_boundary/ibm_tribend.hpp +++ b/src/core/immersed_boundary/ibm_tribend.hpp @@ -22,6 +22,7 @@ #include "config/config.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include @@ -52,8 +53,8 @@ struct IBMTribend { * @return forces on @p p1, @p p2, @p p3, @p p4 */ std::tuple - calc_forces(Particle const &p1, Particle const &p2, Particle const &p3, - Particle const &p4) const; + calc_forces(BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3, Particle const &p4) const; private: friend boost::serialization::access; diff --git a/src/core/immersed_boundary/ibm_triel.cpp b/src/core/immersed_boundary/ibm_triel.cpp index 63426c01875..da9edc58c9b 100644 --- a/src/core/immersed_boundary/ibm_triel.cpp +++ b/src/core/immersed_boundary/ibm_triel.cpp @@ -23,6 +23,7 @@ #include "Particle.hpp" #include "grid.hpp" #include "ibm_common.hpp" +#include "system/System.hpp" #include #include @@ -73,8 +74,8 @@ RotateForces(Utils::Vector2d const &f1_rot, Utils::Vector2d const &f2_rot, } // namespace boost::optional> -IBMTriel::calc_forces(Particle const &p1, Particle const &p2, - Particle const &p3) const { +IBMTriel::calc_forces(BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3) const { // Calculate the current shape of the triangle (l,lp,cos(phi),sin(phi)); // l = length between 1 and 3 @@ -199,6 +200,8 @@ IBMTriel::IBMTriel(const int ind1, const int ind2, const int ind3, const double maxDist, const tElasticLaw elasticLaw, const double k1, const double k2) { + auto const &box_geo = *System::get_system().box_geo; + // collect particles from nodes auto const pos1 = get_ibm_particle_position(ind1); auto const pos2 = get_ibm_particle_position(ind2); diff --git a/src/core/immersed_boundary/ibm_triel.hpp b/src/core/immersed_boundary/ibm_triel.hpp index 9039f094890..78c52027a67 100644 --- a/src/core/immersed_boundary/ibm_triel.hpp +++ b/src/core/immersed_boundary/ibm_triel.hpp @@ -20,6 +20,7 @@ #ifndef IBM_TRIEL_H #define IBM_TRIEL_H +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "config/config.hpp" @@ -69,7 +70,8 @@ struct IBMTriel { * @return the forces on @p p1, @p p2, @p p3 */ boost::optional> - calc_forces(Particle const &p1, Particle const &p2, Particle const &p3) const; + calc_forces(BoxGeometry const &box_geo, Particle const &p1, + Particle const &p2, Particle const &p3) const; private: friend boost::serialization::access; diff --git a/src/core/integrate.cpp b/src/core/integrate.cpp index 287e2caa009..0d3df5f54f0 100644 --- a/src/core/integrate.cpp +++ b/src/core/integrate.cpp @@ -33,6 +33,7 @@ #include "integrators/velocity_verlet_inline.hpp" #include "integrators/velocity_verlet_npt.hpp" +#include "BoxGeometry.hpp" #include "ParticleRange.hpp" #include "accumulators.hpp" #include "bond_breakage/bond_breakage.hpp" @@ -121,7 +122,7 @@ std::weak_ptr get_protocol() { return protocol; } * @brief Update the Lees-Edwards parameters of the box geometry * for the current simulation time. */ -static void update_box_params() { +static void update_box_params(BoxGeometry &box_geo) { if (box_geo.type() == BoxType::LEES_EDWARDS) { assert(protocol != nullptr); box_geo.lees_edwards_update(get_pos_offset(sim_time, *protocol), @@ -132,9 +133,10 @@ static void update_box_params() { void set_protocol(std::shared_ptr new_protocol) { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; + auto &box_geo = *system.box_geo; box_geo.set_type(BoxType::LEES_EDWARDS); protocol = std::move(new_protocol); - LeesEdwards::update_box_params(); + LeesEdwards::update_box_params(box_geo); ::recalc_forces = true; cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } @@ -142,13 +144,14 @@ void set_protocol(std::shared_ptr new_protocol) { void unset_protocol() { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; + auto &box_geo = *system.box_geo; protocol = nullptr; box_geo.set_type(BoxType::CUBOID); ::recalc_forces = true; cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } -template void run_kernel() { +template void run_kernel(BoxGeometry const &box_geo) { if (box_geo.type() == BoxType::LEES_EDWARDS) { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; @@ -179,7 +182,7 @@ void integrator_sanity_checks() { case INTEG_METHOD_NPT_ISO: if (thermo_switch != THERMO_OFF and thermo_switch != THERMO_NPT_ISO) runtimeErrorMsg() << "The NpT integrator requires the NpT thermostat"; - if (box_geo.type() == BoxType::LEES_EDWARDS) + if (System::get_system().box_geo->type() == BoxType::LEES_EDWARDS) runtimeErrorMsg() << "The NpT integrator cannot use Lees-Edwards"; break; #endif @@ -228,8 +231,9 @@ void walberla_agrid_sanity_checks(std::string method, Utils::Vector3d const &lattice_right, double agrid) { // waLBerla and ESPResSo must agree on domain decomposition - auto const &my_left = ::local_geo.my_left(); - auto const &my_right = ::local_geo.my_right(); + auto const &system = System::get_system(); + auto const &my_left = system.local_geo->my_left(); + auto const &my_right = system.local_geo->my_right(); auto const tol = agrid / 1E6; if ((lattice_left - my_left).norm2() > tol or (lattice_right - my_right).norm2() > tol) { @@ -253,7 +257,7 @@ static void resort_particles_if_needed(ParticleRange const &particles) { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; auto const offset = LeesEdwards::verlet_list_offset( - box_geo, cell_structure.get_le_pos_offset_at_last_resort()); + *system.box_geo, cell_structure.get_le_pos_offset_at_last_resort()); if (cell_structure.check_resort_required(particles, skin, offset)) { cell_structure.set_resort_particles(Cells::RESORT_LOCAL); } @@ -327,6 +331,7 @@ int integrate(int n_steps, int reuse_forces) { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; + auto &box_geo = *system.box_geo; // Prepare particle structure and run sanity checks of all active algorithms on_integration_start(time_step); @@ -404,12 +409,12 @@ int integrate(int n_steps, int reuse_forces) { save_old_position(particles, cell_structure.ghost_particles()); #endif - LeesEdwards::update_box_params(); + LeesEdwards::update_box_params(box_geo); bool early_exit = integrator_step_1(particles); if (early_exit) break; - LeesEdwards::run_kernel(); + LeesEdwards::run_kernel(box_geo); #ifdef NPT if (integ_switch != INTEG_METHOD_NPT_ISO) @@ -424,7 +429,7 @@ int integrate(int n_steps, int reuse_forces) { #ifdef BOND_CONSTRAINT // Correct particle positions that participate in a rigid/constrained bond if (n_rigidbonds) { - correct_position_shake(cell_structure); + correct_position_shake(cell_structure, box_geo); } #endif @@ -446,11 +451,11 @@ int integrate(int n_steps, int reuse_forces) { virtual_sites()->after_force_calc(time_step); #endif integrator_step_2(particles, temperature); - LeesEdwards::run_kernel(); + LeesEdwards::run_kernel(box_geo); #ifdef BOND_CONSTRAINT // SHAKE velocity updates if (n_rigidbonds) { - correct_velocity_shake(cell_structure); + correct_velocity_shake(cell_structure, box_geo); } #endif @@ -522,7 +527,7 @@ int integrate(int n_steps, int reuse_forces) { } } // for-loop over integration steps - LeesEdwards::update_box_params(); + LeesEdwards::update_box_params(box_geo); #ifdef CALIPER CALI_CXX_MARK_LOOP_END(integration_loop); #endif @@ -654,7 +659,9 @@ void set_skin(double value) { void set_time(double value) { ::sim_time = value; ::recalc_forces = true; - LeesEdwards::update_box_params(); + auto &system = System::get_system(); + auto &box_geo = *system.box_geo; + LeesEdwards::update_box_params(box_geo); } void set_integ_switch(int value) { diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp index 4b0890ceeba..21b3f4756d4 100644 --- a/src/core/integrators/velocity_verlet_npt.cpp +++ b/src/core/integrators/velocity_verlet_npt.cpp @@ -22,6 +22,7 @@ #ifdef NPT #include "velocity_verlet_npt.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "cell_system/CellStructure.hpp" @@ -91,6 +92,10 @@ void velocity_verlet_npt_finalize_p_inst(double time_step) { void velocity_verlet_npt_propagate_pos(const ParticleRange &particles, double time_step) { + + auto const &system = System::get_system(); + auto &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; Utils::Vector3d scal{}; double L_new = 0.0; @@ -138,7 +143,6 @@ void velocity_verlet_npt_propagate_pos(const ParticleRange &particles, } } - auto &cell_structure = *System::get_system().cell_structure; cell_structure.set_resort_particles(Cells::RESORT_LOCAL); /* Apply new volume to the box-length, communicate it, and account for diff --git a/src/core/lb/Solver.cpp b/src/core/lb/Solver.cpp index 5e7e7c9fce6..7d481fb9f4a 100644 --- a/src/core/lb/Solver.cpp +++ b/src/core/lb/Solver.cpp @@ -153,7 +153,8 @@ Solver::get_interpolated_velocity(Utils::Vector3d const &pos) const { return std::visit( [&](auto &ptr) { auto const agrid = ptr->get_agrid(); - auto const lb_pos = folded_position(pos, ::box_geo) / agrid; + auto const &box_geo = *System::get_system().box_geo; + auto const lb_pos = folded_position(pos, box_geo) / agrid; return ptr->get_velocity_at_pos(lb_pos, false); }, *impl->solver); @@ -164,7 +165,8 @@ Solver::get_interpolated_density(Utils::Vector3d const &pos) const { return std::visit( [&](auto &ptr) { auto const agrid = ptr->get_agrid(); - auto const lb_pos = folded_position(pos, ::box_geo) / agrid; + auto const &box_geo = *System::get_system().box_geo; + auto const lb_pos = folded_position(pos, box_geo) / agrid; return ptr->get_density_at_pos(lb_pos, false); }, *impl->solver); @@ -213,7 +215,8 @@ void Solver::set(std::shared_ptr lb_fluid, assert(not impl->solver.has_value()); auto lb_instance = std::make_shared(lb_fluid, lb_params); lb_instance->sanity_checks(); - auto const &lebc = ::box_geo.lees_edwards_bc(); + auto const &box_geo = *System::get_system().box_geo; + auto const &lebc = box_geo.lees_edwards_bc(); lb_fluid->check_lebc(lebc.shear_direction, lebc.shear_plane_normal); impl->solver = lb_instance; } diff --git a/src/core/lb/particle_coupling.cpp b/src/core/lb/particle_coupling.cpp index edad9319ee7..4f36c0980dc 100644 --- a/src/core/lb/particle_coupling.cpp +++ b/src/core/lb/particle_coupling.cpp @@ -19,6 +19,7 @@ #include "lb/particle_coupling.hpp" #include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" @@ -151,6 +152,7 @@ Utils::Vector3d lb_drag_force(LB::Solver const &lb, Particle const &p, */ static bool in_local_domain(Utils::Vector3d const &pos, double halo = 0.) { auto const halo_vec = Utils::Vector3d::broadcast(halo); + auto const &local_geo = *System::get_system().local_geo; auto const lower_corner = local_geo.my_left() - halo_vec; auto const upper_corner = local_geo.my_right() + halo_vec; @@ -179,6 +181,9 @@ bool in_local_halo(Utils::Vector3d const &pos) { std::vector positions_in_halo(Utils::Vector3d const &pos, BoxGeometry const &box, double agrid) { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; auto const halo = 0.5 * agrid; auto const halo_vec = Utils::Vector3d::broadcast(halo); auto const fully_inside_lower = local_geo.my_left() + 2. * halo_vec; @@ -241,6 +246,7 @@ void ParticleCoupling::kernel(Particle &p) { return; auto const agrid = m_lb.get_agrid(); + auto const &box_geo = *System::get_system().box_geo; // Calculate coupling force Utils::Vector3d force_on_particle = {}; diff --git a/src/core/magnetostatics/dipolar_direct_sum.cpp b/src/core/magnetostatics/dipolar_direct_sum.cpp index ab3f6ba37cd..5384ff1e080 100644 --- a/src/core/magnetostatics/dipolar_direct_sum.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum.cpp @@ -25,10 +25,12 @@ #include "magnetostatics/dipolar_direct_sum.hpp" +#include "BoxGeometry.hpp" #include "cells.hpp" #include "communication.hpp" #include "errorhandling.hpp" #include "grid.hpp" +#include "system/System.hpp" #include #include @@ -50,7 +52,6 @@ #include #include -namespace { /** * @brief Pair force of two interacting dipoles. * @@ -60,8 +61,8 @@ namespace { * * @return Resulting force. */ -auto pair_force(Utils::Vector3d const &d, Utils::Vector3d const &m1, - Utils::Vector3d const &m2) { +static auto pair_force(Utils::Vector3d const &d, Utils::Vector3d const &m1, + Utils::Vector3d const &m2) { auto const pe2 = m1 * d; auto const pe3 = m2 * d; @@ -90,8 +91,8 @@ auto pair_force(Utils::Vector3d const &d, Utils::Vector3d const &m1, * * @return Interaction energy. */ -auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1, - Utils::Vector3d const &m2) { +static auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1, + Utils::Vector3d const &m2) { auto const r2 = d * d; auto const r = sqrt(r2); auto const r3 = r2 * r; @@ -113,7 +114,7 @@ auto pair_potential(Utils::Vector3d const &d, Utils::Vector3d const &m1, * * @return Utils::Vector3d containing dipole field components. */ -auto dipole_field(Utils::Vector3d const &d, Utils::Vector3d const &m1) { +static auto dipole_field(Utils::Vector3d const &d, Utils::Vector3d const &m1) { auto const r2 = d * d; auto const r = sqrt(r2); auto const r3 = r2 * r; @@ -185,7 +186,7 @@ struct PosMom { * at all. If false, distances are calculated as Euclidean * distances, and not using minimum image convention. * @param ncut Number of replicas in each direction. - * @param box_l Box dimensions. + * @param box_geo Box geometry. * @param init Initial value of the sum. * @param f Binary operation mapping distance and moment of the * interaction partner to the value to be summed up for this pair. @@ -195,13 +196,14 @@ struct PosMom { template T image_sum(InputIterator begin, InputIterator end, InputIterator it, bool with_replicas, Utils::Vector3i const &ncut, - Utils::Vector3d const &box_l, T init, F f) { + BoxGeometry const &box_geo, T init, F f) { + auto const &box_l = box_geo.length(); for (auto jt = begin; jt != end; ++jt) { auto const exclude_primary = (it == jt); - auto const primary_distance = - (with_replicas) ? (it->pos - jt->pos) - : ::box_geo.get_mi_vector(it->pos, jt->pos); + auto const primary_distance = (with_replicas) + ? (it->pos - jt->pos) + : box_geo.get_mi_vector(it->pos, jt->pos); for_each_image(ncut, [&](int nx, int ny, int nz) { if (!(exclude_primary && nx == 0 && ny == 0 && nz == 0)) { @@ -216,7 +218,9 @@ T image_sum(InputIterator begin, InputIterator end, InputIterator it, return init; } -auto gather_particle_data(ParticleRange const &particles, int n_replicas) { +static auto gather_particle_data(BoxGeometry const &box_geo, + ParticleRange const &particles, + int n_replicas) { auto const &comm = ::comm_cart; std::vector local_particles; std::vector local_posmom; @@ -230,7 +234,7 @@ auto gather_particle_data(ParticleRange const &particles, int n_replicas) { if (p.dipm() != 0.0) { local_particles.emplace_back(&p); local_posmom.emplace_back( - PosMom{folded_position(p.pos(), ::box_geo), p.calc_dip()}); + PosMom{folded_position(p.pos(), box_geo), p.calc_dip()}); } } @@ -255,14 +259,12 @@ auto gather_particle_data(ParticleRange const &particles, int n_replicas) { std::move(reqs), offset); } -auto get_n_cut(int n_replicas) { - return n_replicas * Utils::Vector3i{static_cast(::box_geo.periodic(0)), - static_cast(::box_geo.periodic(1)), - static_cast(::box_geo.periodic(2))}; +static auto get_n_cut(BoxGeometry const &box_geo, int n_replicas) { + return n_replicas * Utils::Vector3i{static_cast(box_geo.periodic(0)), + static_cast(box_geo.periodic(1)), + static_cast(box_geo.periodic(2))}; } -} // namespace - /** * @brief Calculate and add the interaction forces/torques to the particles. * @@ -287,12 +289,13 @@ auto get_n_cut(int n_replicas) { */ void DipolarDirectSum::add_long_range_forces( ParticleRange const &particles) const { - auto const &box_l = ::box_geo.length(); + auto const &box_geo = *System::get_system().box_geo; + auto const &box_l = box_geo.length(); auto [local_particles, all_posmom, reqs, offset] = - gather_particle_data(particles, n_replicas); + gather_particle_data(box_geo, particles, n_replicas); /* Number of image boxes considered */ - auto const ncut = get_n_cut(n_replicas); + auto const ncut = get_n_cut(box_geo, n_replicas); auto const with_replicas = (ncut.norm2() > 0); /* Range of particles we calculate the ia for on this node */ @@ -307,7 +310,7 @@ void DipolarDirectSum::add_long_range_forces( for (auto it = local_posmom_begin; it != local_posmom_end; ++it, ++p) { /* IA with own images */ auto fi = image_sum( - it, std::next(it), it, with_replicas, ncut, box_l, ParticleForce{}, + it, std::next(it), it, with_replicas, ncut, box_geo, ParticleForce{}, [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) { return pair_force(rn, it->m, mj); }); @@ -315,9 +318,8 @@ void DipolarDirectSum::add_long_range_forces( /* IA with other local particles */ auto q = std::next(p); for (auto jt = std::next(it); jt != local_posmom_end; ++jt, ++q) { - auto const d = (with_replicas) - ? (it->pos - jt->pos) - : ::box_geo.get_mi_vector(it->pos, jt->pos); + auto const d = (with_replicas) ? (it->pos - jt->pos) + : box_geo.get_mi_vector(it->pos, jt->pos); ParticleForce fij{}; ParticleForce fji{}; @@ -352,14 +354,14 @@ void DipolarDirectSum::add_long_range_forces( // red particles auto fi = image_sum(all_posmom.begin(), local_posmom_begin, it, with_replicas, - ncut, box_l, ParticleForce{}, + ncut, box_geo, ParticleForce{}, [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) { return pair_force(rn, it->m, mj); }); // black particles fi += image_sum(local_posmom_end, all_posmom.end(), it, with_replicas, ncut, - box_l, ParticleForce{}, + box_geo, ParticleForce{}, [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) { return pair_force(rn, it->m, mj); }); @@ -376,12 +378,12 @@ void DipolarDirectSum::add_long_range_forces( */ double DipolarDirectSum::long_range_energy(ParticleRange const &particles) const { - auto const &box_l = ::box_geo.length(); + auto const &box_geo = *System::get_system().box_geo; auto [local_particles, all_posmom, reqs, offset] = - gather_particle_data(particles, n_replicas); + gather_particle_data(box_geo, particles, n_replicas); /* Number of image boxes considered */ - auto const ncut = get_n_cut(n_replicas); + auto const ncut = get_n_cut(box_geo, n_replicas); auto const with_replicas = (ncut.norm2() > 0); /* Wait for the rest of the data to arrive */ @@ -394,7 +396,7 @@ DipolarDirectSum::long_range_energy(ParticleRange const &particles) const { auto u = 0.; for (auto it = local_posmom_begin; it != local_posmom_end; ++it) { - u = image_sum(it, all_posmom.end(), it, with_replicas, ncut, box_l, u, + u = image_sum(it, all_posmom.end(), it, with_replicas, ncut, box_geo, u, [it](Utils::Vector3d const &rn, Utils::Vector3d const &mj) { return pair_potential(rn, it->m, mj); }); @@ -415,12 +417,12 @@ DipolarDirectSum::long_range_energy(ParticleRange const &particles) const { #ifdef DIPOLE_FIELD_TRACKING void DipolarDirectSum::dipole_field_at_part( ParticleRange const &particles) const { - auto const &box_l = ::box_geo.length(); + auto const &box_geo = *System::get_system().box_geo; /* collect particle data */ auto [local_particles, all_posmom, reqs, offset] = - gather_particle_data(particles, n_replicas); + gather_particle_data(box_geo, particles, n_replicas); - auto const ncut = get_n_cut(n_replicas); + auto const ncut = get_n_cut(box_geo, n_replicas); auto const with_replicas = (ncut.norm2() > 0); boost::mpi::wait_all(reqs.begin(), reqs.end()); @@ -433,7 +435,7 @@ void DipolarDirectSum::dipole_field_at_part( auto p = local_particles.begin(); for (auto pi = local_posmom_begin; pi != local_posmom_end; ++pi, ++p) { auto const u = image_sum( - all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_l, + all_posmom.begin(), all_posmom.end(), pi, with_replicas, ncut, box_geo, u_init, [](Utils::Vector3d const &rn, Utils::Vector3d const &mj) { return dipole_field(rn, mj); }); diff --git a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp index d2603adc95a..7f9783bcc99 100644 --- a/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp +++ b/src/core/magnetostatics/dipolar_direct_sum_gpu.cpp @@ -24,13 +24,15 @@ #include "magnetostatics/dipolar_direct_sum_gpu.hpp" #include "magnetostatics/dipolar_direct_sum_gpu_cuda.cuh" +#include "BoxGeometry.hpp" #include "communication.hpp" #include "grid.hpp" #include "system/GpuParticleData.hpp" #include "system/System.hpp" -static void get_simulation_box(float *box, int *per) { - for (int i = 0; i < 3; i++) { +static void get_simulation_box(BoxGeometry const &box_geo, float *box, + int *per) { + for (unsigned int i = 0u; i < 3u; i++) { box[i] = static_cast(box_geo.length()[i]); per[i] = box_geo.periodic(i); } @@ -46,14 +48,15 @@ DipolarDirectSumGpu::DipolarDirectSumGpu(double prefactor) } void DipolarDirectSumGpu::add_long_range_forces() const { - auto &gpu = System::get_system().gpu; + auto &system = System::get_system(); + auto &gpu = system.gpu; gpu.update(); if (this_node != 0) { return; } float box[3]; int periodicity[3]; - get_simulation_box(box, periodicity); + get_simulation_box(*system.box_geo, box, periodicity); auto const npart = static_cast(gpu.n_particles()); auto const forces_device = gpu.get_particle_forces_device(); auto const torques_device = gpu.get_particle_torques_device(); @@ -65,14 +68,15 @@ void DipolarDirectSumGpu::add_long_range_forces() const { } void DipolarDirectSumGpu::long_range_energy() const { - auto &gpu = System::get_system().gpu; + auto &system = System::get_system(); + auto &gpu = system.gpu; gpu.update(); if (this_node != 0) { return; } float box[3]; int periodicity[3]; - get_simulation_box(box, periodicity); + get_simulation_box(*system.box_geo, box, periodicity); auto const npart = static_cast(gpu.n_particles()); auto const energy_device = &(gpu.get_energy_device()->dipolar); auto const positions_device = gpu.get_particle_positions_device(); diff --git a/src/core/magnetostatics/dlc.cpp b/src/core/magnetostatics/dlc.cpp index 807d9b183ae..e0b82a0e6c0 100644 --- a/src/core/magnetostatics/dlc.cpp +++ b/src/core/magnetostatics/dlc.cpp @@ -30,6 +30,7 @@ #include "magnetostatics/dp3m.hpp" #include "p3m/common.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" @@ -113,6 +114,7 @@ static void dipolar_force_corrections(int kcut, std::vector &fs, std::vector &ts, ParticleRange const &particles) { + auto const &box_geo = *System::get_system().box_geo; auto const facux = 2. * Utils::pi() * box_geo.length_inv()[0]; auto const facuy = 2. * Utils::pi() * box_geo.length_inv()[1]; @@ -246,6 +248,7 @@ static void dipolar_force_corrections(int kcut, */ static double dipolar_energy_correction(int kcut, ParticleRange const &particles) { + auto const &box_geo = *System::get_system().box_geo; auto const facux = 2. * Utils::pi() * box_geo.length_inv()[0]; auto const facuy = 2. * Utils::pi() * box_geo.length_inv()[1]; @@ -303,6 +306,7 @@ static double dipolar_energy_correction(int kcut, void DipolarLayerCorrection::add_force_corrections( ParticleRange const &particles) const { assert(dlc.far_cut > 0.); + auto const &box_geo = *System::get_system().box_geo; auto const volume = box_geo.volume(); auto const correc = 4. * Utils::pi() / volume; @@ -354,6 +358,7 @@ void DipolarLayerCorrection::add_force_corrections( double DipolarLayerCorrection::energy_correction( ParticleRange const &particles) const { assert(dlc.far_cut > 0.); + auto const &box_geo = *System::get_system().box_geo; auto const volume = box_geo.volume(); auto const pref = prefactor * 2. * Utils::pi() / volume; @@ -411,6 +416,7 @@ static int count_magnetic_particles() { * which there is no particles): Lz = h + gap_size */ double DipolarLayerCorrection::tune_far_cut() const { + auto const &box_geo = *System::get_system().box_geo; /* we take the maximum dipole in the system, to be sure that the errors * in the other case will be equal or less than for this one */ auto const mu_max_sq = Utils::sqr(calc_mu_max()); @@ -476,14 +482,15 @@ void DipolarLayerCorrection::adapt_solver() { } void DipolarLayerCorrection::sanity_checks_node_grid() const { + auto const &box_geo = *System::get_system().box_geo; if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) { throw std::runtime_error("DLC: requires periodicity (True, True, True)"); } } dlc_data::dlc_data(double maxPWerror, double gap_size, double far_cut) - : maxPWerror{maxPWerror}, gap_size{gap_size}, box_h{box_geo.length()[2] - - gap_size}, + : maxPWerror{maxPWerror}, gap_size{gap_size}, + box_h{System::get_system().box_geo->length()[2] - gap_size}, far_cut{far_cut}, far_calculated{far_cut == -1.} { if (far_cut <= 0. and not far_calculated) { throw std::domain_error("Parameter 'far_cut' must be > 0"); @@ -497,6 +504,7 @@ dlc_data::dlc_data(double maxPWerror, double gap_size, double far_cut) } void DipolarLayerCorrection::recalc_box_h() { + auto const &box_geo = *System::get_system().box_geo; auto const new_box_h = box_geo.length()[2] - dlc.gap_size; if (new_box_h < 0.) { throw std::runtime_error("DLC gap size (" + std::to_string(dlc.gap_size) + diff --git a/src/core/magnetostatics/dp3m.cpp b/src/core/magnetostatics/dp3m.cpp index 58db93e000c..2ba55c63e5a 100644 --- a/src/core/magnetostatics/dp3m.cpp +++ b/src/core/magnetostatics/dp3m.cpp @@ -39,6 +39,8 @@ #include "p3m/interpolation.hpp" #include "p3m/send_mesh.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "cell_system/CellStructure.hpp" @@ -107,6 +109,7 @@ double DipolarP3M::calc_average_self_energy_k_space() const { auto const start = Utils::Vector3i{dp3m.fft.plan[3].start}; auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh}; + auto const &box_geo = *System::get_system().box_geo; auto const node_phi = grid_influence_function_self_energy( dp3m.params, start, start + size, dp3m.g_energy); @@ -121,6 +124,10 @@ void DipolarP3M::init() { assert(dp3m.params.cao >= 1 and dp3m.params.cao <= 7); assert(dp3m.params.alpha > 0.); + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; + dp3m.params.cao3 = Utils::int_pow<3>(dp3m.params.cao); dp3m.params.recalc_a_ai_cao_cut(box_geo.length()); dp3m.local_mesh.calc_local_ca_mesh(dp3m.params, local_geo, skin, 0.); @@ -258,6 +265,7 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag, double k_space_energy_dip = 0.; double tmp0, tmp1; + auto const &box_geo = *System::get_system().box_geo; auto const dipole_prefac = prefactor / Utils::int_pow<3>(dp3m.params.mesh[0]); if (dp3m.sum_mu2 > 0) { @@ -506,6 +514,7 @@ double DipolarP3M::kernel(bool force_flag, bool energy_flag, double DipolarP3M::calc_surface_term(bool force_flag, bool energy_flag, ParticleRange const &particles) { + auto const &box_geo = *System::get_system().box_geo; auto const pref = prefactor * 4. * Utils::pi() / box_geo.volume() / (2. * dp3m.params.epsilon + 1.); auto const n_local_part = particles.size(); @@ -573,6 +582,7 @@ double DipolarP3M::calc_surface_term(bool force_flag, bool energy_flag, void DipolarP3M::calc_influence_function_force() { auto const start = Utils::Vector3i{dp3m.fft.plan[3].start}; auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh}; + auto const &box_geo = *System::get_system().box_geo; dp3m.g_force = grid_influence_function<3>(dp3m.params, start, start + size, box_geo.length()); @@ -581,6 +591,7 @@ void DipolarP3M::calc_influence_function_force() { void DipolarP3M::calc_influence_function_energy() { auto const start = Utils::Vector3i{dp3m.fft.plan[3].start}; auto const size = Utils::Vector3i{dp3m.fft.plan[3].new_mesh}; + auto const &box_geo = *System::get_system().box_geo; dp3m.g_energy = grid_influence_function<2>(dp3m.params, start, start + size, box_geo.length()); @@ -605,6 +616,7 @@ class DipolarTuningAlgorithm : public TuningAlgorithm { } void setup_logger(bool verbose) override { + auto const &box_geo = *System::get_system().box_geo; m_logger = std::make_unique( verbose and this_node == 0, "DipolarP3M", TuningLogger::Mode::Dipolar); m_logger->tuning_goals(dp3m.params.accuracy, m_prefactor, @@ -618,6 +630,7 @@ class DipolarTuningAlgorithm : public TuningAlgorithm { double r_cut_iL) const override { double alpha_L, rs_err, ks_err; + auto const &box_geo = *System::get_system().box_geo; /* calc maximal real space error for setting */ rs_err = dp3m_real_space_error(box_geo.length()[0], r_cut_iL, @@ -698,6 +711,7 @@ class DipolarTuningAlgorithm : public TuningAlgorithm { }; void DipolarP3M::tune() { + auto const &box_geo = *System::get_system().box_geo; if (dp3m.params.alpha_L == 0. and dp3m.params.alpha != 0.) { dp3m.params.alpha_L = dp3m.params.alpha * box_geo.length()[0]; } @@ -861,7 +875,10 @@ double dp3m_rtbisection(double box_size, double r_cut_iL, int n_c_part, } void DipolarP3M::sanity_checks_boxl() const { - for (unsigned int i = 0; i < 3; i++) { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; + for (unsigned int i = 0u; i < 3u; i++) { /* check k-space cutoff */ if (dp3m.params.cao_cut[i] >= box_geo.length_half()[i]) { std::stringstream msg; @@ -884,6 +901,7 @@ void DipolarP3M::sanity_checks_boxl() const { } void DipolarP3M::sanity_checks_periodicity() const { + auto const &box_geo = *System::get_system().box_geo; if (!box_geo.periodic(0) or !box_geo.periodic(1) or !box_geo.periodic(2)) { throw std::runtime_error( "DipolarP3M: requires periodicity (True, True, True)"); @@ -891,6 +909,7 @@ void DipolarP3M::sanity_checks_periodicity() const { } void DipolarP3M::sanity_checks_cell_structure() const { + auto const &local_geo = *System::get_system().local_geo; if (local_geo.cell_structure_type() != CellStructureType::CELL_STRUCTURE_REGULAR && local_geo.cell_structure_type() != @@ -915,6 +934,7 @@ void DipolarP3M::sanity_checks_node_grid() const { } void DipolarP3M::scaleby_box_l() { + auto const &box_geo = *System::get_system().box_geo; dp3m.params.r_cut = dp3m.params.r_cut_iL * box_geo.length()[0]; dp3m.params.alpha = dp3m.params.alpha_L * box_geo.length_inv()[0]; dp3m.params.recalc_a_ai_cao_cut(box_geo.length()); @@ -926,6 +946,7 @@ void DipolarP3M::scaleby_box_l() { } void DipolarP3M::calc_energy_correction() { + auto const &box_geo = *System::get_system().box_geo; auto const Ukp3m = calc_average_self_energy_k_space() * box_geo.volume(); auto const Ewald_volume = Utils::int_pow<3>(dp3m.params.alpha_L); auto const Eself = -2. * Ewald_volume * Utils::sqrt_pi_i() / 3.; diff --git a/src/core/magnetostatics/scafacos_impl.cpp b/src/core/magnetostatics/scafacos_impl.cpp index 22b8eeab091..0a7a8b5ad14 100644 --- a/src/core/magnetostatics/scafacos_impl.cpp +++ b/src/core/magnetostatics/scafacos_impl.cpp @@ -26,6 +26,7 @@ #include "magnetostatics/scafacos.hpp" #include "magnetostatics/scafacos_impl.hpp" +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" #include "grid.hpp" @@ -46,7 +47,9 @@ make_dipolar_scafacos(std::string const &method, } void DipolarScafacosImpl::update_particle_data() { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; positions.clear(); dipoles.clear(); diff --git a/src/core/object-in-fluid/oif_global_forces.cpp b/src/core/object-in-fluid/oif_global_forces.cpp index 6cdc35d3a13..641b305e6ef 100644 --- a/src/core/object-in-fluid/oif_global_forces.cpp +++ b/src/core/object-in-fluid/oif_global_forces.cpp @@ -33,13 +33,14 @@ int max_oif_objects = 0; -Utils::Vector2d calc_oif_global(int molType, CellStructure &cs) { +Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo, + CellStructure &cs) { // first-fold-then-the-same approach double partArea = 0.0; // z volume double VOL_partVol = 0.; - cs.bond_loop([&partArea, &VOL_partVol, + cs.bond_loop([&partArea, &VOL_partVol, &box_geo, molType](Particle &p1, int bond_id, Utils::Span partners) { if (p1.mol_id() != molType) @@ -70,13 +71,14 @@ Utils::Vector2d calc_oif_global(int molType, CellStructure &cs) { } void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, - CellStructure &cs) { + BoxGeometry const &box_geo, CellStructure &cs) { // first-fold-then-the-same approach double area = area_volume[0]; double VOL_volume = area_volume[1]; - cs.bond_loop([area, VOL_volume, molType](Particle &p1, int bond_id, - Utils::Span partners) { + cs.bond_loop([&box_geo, area, VOL_volume, + molType](Particle &p1, int bond_id, + Utils::Span partners) { if (p1.mol_id() != molType) return false; diff --git a/src/core/object-in-fluid/oif_global_forces.hpp b/src/core/object-in-fluid/oif_global_forces.hpp index a0bd8924350..908a349e136 100644 --- a/src/core/object-in-fluid/oif_global_forces.hpp +++ b/src/core/object-in-fluid/oif_global_forces.hpp @@ -23,6 +23,7 @@ * for a particle triple (triangle from mesh). See @cite dupin07a. */ +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "oif_global_forces_params.hpp" @@ -35,11 +36,12 @@ * - MPI synchronization with all reduce * - !!! loop over particles from regular_decomposition !!! */ -Utils::Vector2d calc_oif_global(int molType, CellStructure &cs); +Utils::Vector2d calc_oif_global(int molType, BoxGeometry const &box_geo, + CellStructure &cs); /** Distribute the OIF global forces to all particles in the mesh. */ void add_oif_global_forces(Utils::Vector2d const &area_volume, int molType, - CellStructure &cs); + BoxGeometry const &box_geo, CellStructure &cs); extern int max_oif_objects; diff --git a/src/core/object-in-fluid/oif_local_forces.hpp b/src/core/object-in-fluid/oif_local_forces.hpp index 3a07c67595e..52e0ebbd8a8 100644 --- a/src/core/object-in-fluid/oif_local_forces.hpp +++ b/src/core/object-in-fluid/oif_local_forces.hpp @@ -79,8 +79,8 @@ struct OifLocalForcesBond { } std::tuple - calc_forces(Particle const &p2, Particle const &p1, Particle const &p3, - Particle const &p4) const; + calc_forces(BoxGeometry const &box_geo, Particle const &p2, + Particle const &p1, Particle const &p3, Particle const &p4) const; private: friend boost::serialization::access; @@ -107,8 +107,9 @@ struct OifLocalForcesBond { */ inline std::tuple -OifLocalForcesBond::calc_forces(Particle const &p2, Particle const &p1, - Particle const &p3, Particle const &p4) const { +OifLocalForcesBond::calc_forces(BoxGeometry const &box_geo, Particle const &p2, + Particle const &p1, Particle const &p3, + Particle const &p4) const { // first-fold-then-the-same approach auto const fp2 = diff --git a/src/core/observables/BondAngles.hpp b/src/core/observables/BondAngles.hpp index 56500e1a378..981e68ff570 100644 --- a/src/core/observables/BondAngles.hpp +++ b/src/core/observables/BondAngles.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "PidObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include @@ -58,6 +58,7 @@ class BondAngles : public PidObservable { return {}; } + auto const &box_geo = *System::get_system().box_geo; std::vector res(n_values()); auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]); auto n1 = v1.norm(); diff --git a/src/core/observables/BondDihedrals.hpp b/src/core/observables/BondDihedrals.hpp index 864628cd915..0d246313509 100644 --- a/src/core/observables/BondDihedrals.hpp +++ b/src/core/observables/BondDihedrals.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "PidObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include @@ -62,6 +62,7 @@ class BondDihedrals : public PidObservable { return {}; } + auto const &box_geo = *System::get_system().box_geo; std::vector res(n_values()); auto v1 = box_geo.get_mi_vector(positions_sorted[1], positions_sorted[0]); auto v2 = box_geo.get_mi_vector(positions_sorted[2], positions_sorted[1]); diff --git a/src/core/observables/CosPersistenceAngles.hpp b/src/core/observables/CosPersistenceAngles.hpp index 10c7f892c1d..84b75d878ae 100644 --- a/src/core/observables/CosPersistenceAngles.hpp +++ b/src/core/observables/CosPersistenceAngles.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "PidObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include @@ -59,6 +59,7 @@ class CosPersistenceAngles : public PidObservable { return {}; } + auto const &box_geo = *System::get_system().box_geo; auto const no_of_angles = n_values(); auto const no_of_bonds = no_of_angles + 1; std::vector angles(no_of_angles); diff --git a/src/core/observables/CylindricalDensityProfile.hpp b/src/core/observables/CylindricalDensityProfile.hpp index 4b0c18f7871..93ad7837159 100644 --- a/src/core/observables/CylindricalDensityProfile.hpp +++ b/src/core/observables/CylindricalDensityProfile.hpp @@ -22,7 +22,7 @@ #include "CylindricalPidProfileObservable.hpp" #include "BoxGeometry.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -40,6 +40,7 @@ class CylindricalDensityProfile : public CylindricalPidProfileObservable { ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { using pos_type = Utils::Vector3d; + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; local_folded_positions.reserve(local_particles.size()); diff --git a/src/core/observables/CylindricalFluxDensityProfile.hpp b/src/core/observables/CylindricalFluxDensityProfile.hpp index 9e27a6c3394..ed36325549d 100644 --- a/src/core/observables/CylindricalFluxDensityProfile.hpp +++ b/src/core/observables/CylindricalFluxDensityProfile.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "CylindricalPidProfileObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -42,6 +42,7 @@ class CylindricalFluxDensityProfile : public CylindricalPidProfileObservable { const ParticleObservables::traits &traits) const override { using pos_type = decltype(traits.position(std::declval())); using vel_type = decltype(traits.velocity(std::declval())); + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; local_folded_positions.reserve(local_particles.size()); diff --git a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp index 9c4b06e8ed5..074ebaf3ea2 100644 --- a/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBFluxDensityProfileAtParticlePositions.cpp @@ -19,7 +19,6 @@ #include "CylindricalLBFluxDensityProfileAtParticlePositions.hpp" #include "BoxGeometry.hpp" -#include "grid.hpp" #include "system/System.hpp" #include "utils_histogram.hpp" @@ -46,7 +45,9 @@ CylindricalLBFluxDensityProfileAtParticlePositions::evaluate( local_folded_positions.reserve(local_particles.size()); local_flux_densities.reserve(local_particles.size()); - auto const &lb = System::get_system().lb; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &lb = system.lb; auto const vel_conv = lb.get_lattice_speed(); for (auto const &p : local_particles) { diff --git a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp index cce02d6f13e..943d55f758d 100644 --- a/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp +++ b/src/core/observables/CylindricalLBVelocityProfileAtParticlePositions.cpp @@ -19,7 +19,6 @@ #include "CylindricalLBVelocityProfileAtParticlePositions.hpp" #include "BoxGeometry.hpp" -#include "grid.hpp" #include "system/System.hpp" #include "utils_histogram.hpp" @@ -42,7 +41,9 @@ std::vector CylindricalLBVelocityProfileAtParticlePositions::evaluate( local_folded_positions.reserve(local_particles.size()); local_velocities.reserve(local_particles.size()); - auto const &lb = System::get_system().lb; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &lb = system.lb; auto const vel_conv = lb.get_lattice_speed(); for (auto const &p : local_particles) { diff --git a/src/core/observables/CylindricalVelocityProfile.hpp b/src/core/observables/CylindricalVelocityProfile.hpp index 394df7cb034..20a408e88fd 100644 --- a/src/core/observables/CylindricalVelocityProfile.hpp +++ b/src/core/observables/CylindricalVelocityProfile.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "CylindricalPidProfileObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -41,6 +41,7 @@ class CylindricalVelocityProfile : public CylindricalPidProfileObservable { const ParticleObservables::traits &traits) const override { using pos_type = Utils::Vector3d; using vel_type = Utils::Vector3d; + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; std::vector local_velocities{}; diff --git a/src/core/observables/DensityProfile.hpp b/src/core/observables/DensityProfile.hpp index 4ebad035dca..82d803d96fc 100644 --- a/src/core/observables/DensityProfile.hpp +++ b/src/core/observables/DensityProfile.hpp @@ -22,7 +22,7 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" #include "PidProfileObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -41,6 +41,7 @@ class DensityProfile : public PidProfileObservable { ParticleReferenceRange const &local_particles, const ParticleObservables::traits &traits) const override { using pos_type = decltype(traits.position(std::declval())); + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; local_folded_positions.reserve(local_particles.size()); diff --git a/src/core/observables/FluxDensityProfile.hpp b/src/core/observables/FluxDensityProfile.hpp index 1887c8d8c61..fa0f5eeb367 100644 --- a/src/core/observables/FluxDensityProfile.hpp +++ b/src/core/observables/FluxDensityProfile.hpp @@ -22,7 +22,7 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" #include "PidProfileObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -46,6 +46,7 @@ class FluxDensityProfile : public PidProfileObservable { const ParticleObservables::traits &traits) const override { using pos_type = decltype(traits.position(std::declval())); using vel_type = decltype(traits.velocity(std::declval())); + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; std::vector local_velocities{}; diff --git a/src/core/observables/ForceDensityProfile.hpp b/src/core/observables/ForceDensityProfile.hpp index bd3ded0f03e..1b90464de37 100644 --- a/src/core/observables/ForceDensityProfile.hpp +++ b/src/core/observables/ForceDensityProfile.hpp @@ -22,7 +22,7 @@ #include "BoxGeometry.hpp" #include "Particle.hpp" #include "PidProfileObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "utils_histogram.hpp" #include @@ -50,6 +50,7 @@ class ForceDensityProfile : public PidProfileObservable { const ParticleObservables::traits &traits) const override { using pos_type = decltype(traits.position(std::declval())); using force_type = decltype(traits.force(std::declval())); + auto const &box_geo = *System::get_system().box_geo; std::vector local_folded_positions{}; local_folded_positions.reserve(local_particles.size()); diff --git a/src/core/observables/ParticleDistances.hpp b/src/core/observables/ParticleDistances.hpp index 1db61a1a690..5329244ce09 100644 --- a/src/core/observables/ParticleDistances.hpp +++ b/src/core/observables/ParticleDistances.hpp @@ -21,7 +21,7 @@ #include "BoxGeometry.hpp" #include "PidObservable.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include #include @@ -55,6 +55,7 @@ class ParticleDistances : public PidObservable { return {}; } + auto const &box_geo = *System::get_system().box_geo; std::vector res(n_values()); for (std::size_t i = 0, end = n_values(); i < end; i++) { diff --git a/src/core/observables/ParticleTraits.hpp b/src/core/observables/ParticleTraits.hpp index fe57ed0d419..91df8817b64 100644 --- a/src/core/observables/ParticleTraits.hpp +++ b/src/core/observables/ParticleTraits.hpp @@ -24,6 +24,7 @@ #include "config/config.hpp" #include "grid.hpp" #include "rotation.hpp" +#include "system/System.hpp" namespace ParticleObservables { /** @@ -34,6 +35,7 @@ namespace ParticleObservables { template <> struct traits { auto id(Particle const &p) const { return p.id(); } auto position(Particle const &p) const { + auto const &box_geo = *System::get_system().box_geo; return unfolded_position(p.pos(), p.image_box(), box_geo.length()); } auto position_folded(Particle const &p) const { return p.pos(); } diff --git a/src/core/observables/RDF.cpp b/src/core/observables/RDF.cpp index c4975ca4e66..e631fefc783 100644 --- a/src/core/observables/RDF.cpp +++ b/src/core/observables/RDF.cpp @@ -20,7 +20,7 @@ #include "BoxGeometry.hpp" #include "fetch_particles.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include #include @@ -61,12 +61,13 @@ RDF::evaluate(boost::mpi::communicator const &comm, return {}; } + auto const &box_geo = *System::get_system().box_geo; auto const bin_width = (max_r - min_r) / static_cast(n_r_bins); auto const inv_bin_width = 1.0 / bin_width; std::vector res(n_values(), 0.0); long int cnt = 0; - auto op = [this, inv_bin_width, &cnt, &res](auto const &pos1, - auto const &pos2) { + auto op = [this, inv_bin_width, &cnt, &res, &box_geo](auto const &pos1, + auto const &pos2) { auto const dist = box_geo.get_mi_vector(pos1, pos2).norm(); if (dist > min_r && dist < max_r) { auto const ind = diff --git a/src/core/p3m/TuningAlgorithm.cpp b/src/core/p3m/TuningAlgorithm.cpp index 9fa852e36ed..21502e82d3d 100644 --- a/src/core/p3m/TuningAlgorithm.cpp +++ b/src/core/p3m/TuningAlgorithm.cpp @@ -28,9 +28,12 @@ #include "tuning.hpp" +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" #include "communication.hpp" #include "grid.hpp" #include "integrate.hpp" +#include "system/System.hpp" #include @@ -55,6 +58,9 @@ static auto constexpr P3M_TUNE_ACCURACY_TOO_LARGE = 3.; static auto constexpr P3M_RCUT_PREC = 1e-3; void TuningAlgorithm::determine_r_cut_limits() { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; auto const r_cut_iL = get_params().r_cut_iL; if (r_cut_iL == 0.) { auto const min_box_l = *boost::min_element(box_geo.length()); @@ -84,6 +90,7 @@ void TuningAlgorithm::determine_cao_limits(int initial_cao) { void TuningAlgorithm::commit(Utils::Vector3i const &mesh, int cao, double r_cut_iL, double alpha_L) { + auto const &box_geo = *System::get_system().box_geo; auto &p3m_params = get_params(); p3m_params.r_cut = r_cut_iL * box_geo.length()[0]; p3m_params.r_cut_iL = r_cut_iL; @@ -113,6 +120,9 @@ double TuningAlgorithm::get_mc_time(Utils::Vector3i const &mesh, int cao, double &tuned_r_cut_iL, double &tuned_alpha_L, double &tuned_accuracy) { + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto const &local_geo = *system.local_geo; auto const target_accuracy = get_params().accuracy; double rs_err, ks_err; double r_cut_iL_min = m_r_cut_iL_min; diff --git a/src/core/pair_criteria/DistanceCriterion.hpp b/src/core/pair_criteria/DistanceCriterion.hpp index 2cf2ee72378..2d0a935b56b 100644 --- a/src/core/pair_criteria/DistanceCriterion.hpp +++ b/src/core/pair_criteria/DistanceCriterion.hpp @@ -21,7 +21,8 @@ #include "pair_criteria/PairCriterion.hpp" -#include "grid.hpp" +#include "BoxGeometry.hpp" +#include "system/System.hpp" namespace PairCriteria { /** @@ -31,6 +32,7 @@ namespace PairCriteria { class DistanceCriterion : public PairCriterion { public: bool decide(const Particle &p1, const Particle &p2) const override { + auto const &box_geo = *System::get_system().box_geo; return box_geo.get_mi_vector(p1.pos(), p2.pos()).norm() <= m_cut_off; } double get_cut_off() { return m_cut_off; } diff --git a/src/core/pair_criteria/EnergyCriterion.hpp b/src/core/pair_criteria/EnergyCriterion.hpp index af331e4e8ed..e5a09dae382 100644 --- a/src/core/pair_criteria/EnergyCriterion.hpp +++ b/src/core/pair_criteria/EnergyCriterion.hpp @@ -21,9 +21,9 @@ #include "pair_criteria/PairCriterion.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "energy_inline.hpp" -#include "grid.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" #include "system/System.hpp" @@ -34,6 +34,8 @@ namespace PairCriteria { class EnergyCriterion : public PairCriterion { public: bool decide(Particle const &p1, Particle const &p2) const override { + auto const &box_geo = *System::get_system().box_geo; + // Distance between particles auto const d = box_geo.get_mi_vector(p1.pos(), p2.pos()); diff --git a/src/core/particle_node.cpp b/src/core/particle_node.cpp index b2a11a00976..12ae5d4add1 100644 --- a/src/core/particle_node.cpp +++ b/src/core/particle_node.cpp @@ -21,6 +21,7 @@ #include "particle_node.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" @@ -438,6 +439,7 @@ static int calculate_max_seen_id() { * @return Whether the particle was created on that node. */ static bool maybe_insert_particle(int p_id, Utils::Vector3d const &pos) { + auto const &box_geo = *System::get_system().box_geo; auto folded_pos = pos; auto image_box = Utils::Vector3i{}; fold_position(folded_pos, image_box, box_geo); @@ -458,7 +460,9 @@ static bool maybe_insert_particle(int p_id, Utils::Vector3d const &pos) { * @return Whether the particle was moved from that node. */ static bool maybe_move_particle(int p_id, Utils::Vector3d const &pos) { - auto p = get_cell_structure().get_local_particle(p_id); + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto p = system.cell_structure->get_local_particle(p_id); if (p == nullptr) { return false; } diff --git a/src/core/polymer.cpp b/src/core/polymer.cpp index d17e2f5ea5c..6adfec25dd2 100644 --- a/src/core/polymer.cpp +++ b/src/core/polymer.cpp @@ -33,8 +33,8 @@ #include "constraints.hpp" #include "constraints/Constraints.hpp" #include "constraints/ShapeBasedConstraint.hpp" -#include "grid.hpp" #include "random.hpp" +#include "system/System.hpp" #include #include @@ -48,7 +48,8 @@ #include #include -template static Utils::Vector3d random_position(RNG &rng) { +template +static Utils::Vector3d random_position(BoxGeometry const &box_geo, RNG &rng) { Utils::Vector3d v; for (int i = 0; i < 3; ++i) v[i] = box_geo.length()[i] * rng(); @@ -73,6 +74,7 @@ template static Utils::Vector3d random_unit_vector(RNG &rng) { * @return the minimal distance of a particle to coordinates @p pos */ static double distto(PartCfg &partCfg, Utils::Vector3d const &pos) { + auto const &box_geo = *System::get_system().box_geo; auto mindist_sq = std::numeric_limits::infinity(); for (auto const &p : partCfg) { @@ -97,8 +99,8 @@ static double distto(PartCfg &partCfg, Utils::Vector3d const &pos) { static bool is_valid_position(Utils::Vector3d const &pos, std::vector> const &positions, - PartCfg &partCfg, double const min_distance, - int const respect_constraints) { + PartCfg &partCfg, BoxGeometry const &box_geo, + double const min_distance, int const respect_constraints) { // check if constraint is violated if (respect_constraints) { Utils::Vector3d const folded_pos = folded_position(pos, box_geo); @@ -144,6 +146,8 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers, double const min_distance, int const max_tries, int const use_bond_angle, double const bond_angle, int const respect_constraints, int const seed) { + + auto const &box_geo = *System::get_system().box_geo; auto rng = [mt = Random::mt19937(static_cast(seed)), dist = std::uniform_real_distribution( 0.0, 1.0)]() mutable { return dist(mt); }; @@ -153,9 +157,9 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers, p.reserve(beads_per_chain); } - auto is_valid_pos = [&positions, &partCfg, min_distance, + auto is_valid_pos = [&positions, &partCfg, min_distance, &box_geo, respect_constraints](Utils::Vector3d const &v) { - return is_valid_position(v, positions, partCfg, min_distance, + return is_valid_position(v, positions, partCfg, box_geo, min_distance, respect_constraints); }; @@ -172,7 +176,7 @@ draw_polymer_positions(PartCfg &partCfg, int const n_polymers, auto draw_monomer_position = [&](int p, int m) { if (m == 0) { return (p < start_positions.size()) ? start_positions[p] - : random_position(rng); + : random_position(box_geo, rng); } if (not use_bond_angle or m < 2) { diff --git a/src/core/pressure.cpp b/src/core/pressure.cpp index 21ea649d5d3..ae698501de1 100644 --- a/src/core/pressure.cpp +++ b/src/core/pressure.cpp @@ -23,11 +23,11 @@ */ #include "pressure.hpp" +#include "BoxGeometry.hpp" #include "Observable_stat.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" -//#include "cells.hpp" #include "communication.hpp" #include "electrostatics/coulomb.hpp" #include "event.hpp" @@ -63,7 +63,7 @@ std::shared_ptr calculate_pressure() { auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; - auto const volume = box_geo.volume(); + auto const volume = system.box_geo->volume(); auto const local_parts = cell_structure.local_particles(); auto const n_nodes = ::communicator.size; @@ -76,11 +76,12 @@ std::shared_ptr calculate_pressure() { auto const coulomb_pressure_kernel = coulomb.pair_pressure_kernel(); short_range_loop( - [&obs_pressure, coulomb_force_kernel_ptr = get_ptr(coulomb_force_kernel)]( - Particle const &p1, int bond_id, Utils::Span partners) { + [&obs_pressure, coulomb_force_kernel_ptr = get_ptr(coulomb_force_kernel), + box_geo = system.box_geo](Particle const &p1, int bond_id, + Utils::Span partners) { auto const &iaparams = *bonded_ia_params.at(bond_id); auto const result = calc_bonded_pressure_tensor( - iaparams, p1, partners, coulomb_force_kernel_ptr); + iaparams, p1, partners, *box_geo, coulomb_force_kernel_ptr); if (result) { auto const &tensor = result.get(); /* pressure tensor part */ diff --git a/src/core/pressure_inline.hpp b/src/core/pressure_inline.hpp index a1f6354eb89..a761cdd7e2f 100644 --- a/src/core/pressure_inline.hpp +++ b/src/core/pressure_inline.hpp @@ -29,6 +29,7 @@ #include "magnetostatics/dipoles.hpp" #include "nonbonded_interactions/nonbonded_interaction_data.hpp" +#include "BoxGeometry.hpp" #include "Observable_stat.hpp" #include "Particle.hpp" #include "errorhandling.hpp" @@ -100,7 +101,7 @@ inline void add_non_bonded_pair_virials( boost::optional> calc_bonded_virial_pressure_tensor( Bonded_IA_Parameters const &iaparams, Particle const &p1, - Particle const &p2, + Particle const &p2, BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { auto const dx = box_geo.get_mi_vector(p1.pos(), p2.pos()); auto const result = calc_bond_pair_force(p1, p2, iaparams, dx, kernel); @@ -116,7 +117,8 @@ boost::optional> calc_bonded_virial_pressure_tensor( boost::optional> calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams, Particle const &p1, Particle const &p2, - Particle const &p3) { + Particle const &p3, + BoxGeometry const &box_geo) { if ((boost::get(&iaparams) != nullptr) || (boost::get(&iaparams) != nullptr) || #ifdef TABULATED @@ -126,7 +128,8 @@ calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams, auto const dx21 = -box_geo.get_mi_vector(p1.pos(), p2.pos()); auto const dx31 = box_geo.get_mi_vector(p3.pos(), p1.pos()); - auto const result = calc_bonded_three_body_force(iaparams, p1, p2, p3); + auto const result = + calc_bonded_three_body_force(iaparams, box_geo, p1, p2, p3); if (result) { Utils::Vector3d force2, force3; std::tie(std::ignore, force2, force3) = result.get(); @@ -146,15 +149,15 @@ calc_bonded_three_body_pressure_tensor(Bonded_IA_Parameters const &iaparams, inline boost::optional> calc_bonded_pressure_tensor( Bonded_IA_Parameters const &iaparams, Particle const &p1, - Utils::Span partners, + Utils::Span partners, BoxGeometry const &box_geo, Coulomb::ShortRangeForceKernel::kernel_type const *kernel) { switch (number_of_partners(iaparams)) { case 1: return calc_bonded_virial_pressure_tensor(iaparams, p1, *partners[0], - kernel); + box_geo, kernel); case 2: return calc_bonded_three_body_pressure_tensor(iaparams, p1, *partners[0], - *partners[1]); + *partners[1], box_geo); default: runtimeWarningMsg() << "Unsupported bond type " + std::to_string(iaparams.which()) + diff --git a/src/core/rattle.cpp b/src/core/rattle.cpp index aede9e40437..57f7c7444a0 100644 --- a/src/core/rattle.cpp +++ b/src/core/rattle.cpp @@ -23,6 +23,7 @@ #ifdef BOND_CONSTRAINT +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" #include "bonded_interactions/bonded_interaction_data.hpp" @@ -73,6 +74,7 @@ static void init_correction_vector(const ParticleRange &particles, * @return True if there was a correction. */ static bool calculate_positional_correction(RigidBond const &ia_params, + BoxGeometry const &box_geo, Particle &p1, Particle &p2) { auto const r_ij = box_geo.get_mi_vector(p1.pos(), p2.pos()); auto const r_ij2 = r_ij.norm2(); @@ -102,21 +104,24 @@ static bool calculate_positional_correction(RigidBond const &ia_params, * @return True if correction is necessary */ template -static bool compute_correction_vector(CellStructure &cs, Kernel kernel) { +static bool compute_correction_vector(CellStructure &cs, + BoxGeometry const &box_geo, + Kernel kernel) { bool correction = false; - cs.bond_loop([&correction, &kernel](Particle &p1, int bond_id, - Utils::Span partners) { - auto const &iaparams = *bonded_ia_params.at(bond_id); - - if (auto const *bond = boost::get(&iaparams)) { - auto const corrected = kernel(*bond, p1, *partners[0]); - if (corrected) - correction = true; - } - - /* Rigid bonds cannot break */ - return false; - }); + cs.bond_loop( + [&correction, &kernel, &box_geo](Particle &p1, int bond_id, + Utils::Span partners) { + auto const &iaparams = *bonded_ia_params.at(bond_id); + + if (auto const *bond = boost::get(&iaparams)) { + auto const corrected = kernel(*bond, box_geo, p1, *partners[0]); + if (corrected) + correction = true; + } + + /* Rigid bonds cannot break */ + return false; + }); return correction; } @@ -133,7 +138,7 @@ static void apply_positional_correction(const ParticleRange &particles) { }); } -void correct_position_shake(CellStructure &cs) { +void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo) { cells_update_ghosts(Cells::DATA_PART_POSITION | Cells::DATA_PART_PROPERTIES); auto particles = cs.local_particles(); @@ -143,7 +148,7 @@ void correct_position_shake(CellStructure &cs) { for (cnt = 0; cnt < SHAKE_MAX_ITERATIONS; ++cnt) { init_correction_vector(particles, ghost_particles); bool const repeat_ = - compute_correction_vector(cs, calculate_positional_correction); + compute_correction_vector(cs, box_geo, calculate_positional_correction); bool const repeat = boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or()); @@ -176,6 +181,7 @@ void correct_position_shake(CellStructure &cs) { * @return True if there was a correction. */ static bool calculate_velocity_correction(RigidBond const &ia_params, + BoxGeometry const &box_geo, Particle &p1, Particle &p2) { auto const v_ij = p1.v() - p2.v(); auto const r_ij = box_geo.get_mi_vector(p1.pos(), p2.pos()); @@ -200,12 +206,13 @@ static bool calculate_velocity_correction(RigidBond const &ia_params, * * @param particles particle range */ -static void apply_velocity_correction(const ParticleRange &particles) { +static void apply_velocity_correction(ParticleRange const &particles, + BoxGeometry const &box_geo) { boost::for_each(particles, [](Particle &p) { p.v() += p.rattle_params().correction; }); } -void correct_velocity_shake(CellStructure &cs) { +void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo) { cs.ghosts_update(Cells::DATA_PART_POSITION | Cells::DATA_PART_MOMENTUM); auto particles = cs.local_particles(); @@ -215,7 +222,7 @@ void correct_velocity_shake(CellStructure &cs) { for (cnt = 0; cnt < SHAKE_MAX_ITERATIONS; ++cnt) { init_correction_vector(particles, ghost_particles); bool const repeat_ = - compute_correction_vector(cs, calculate_velocity_correction); + compute_correction_vector(cs, box_geo, calculate_velocity_correction); bool const repeat = boost::mpi::all_reduce(comm_cart, repeat_, std::logical_or()); @@ -225,7 +232,7 @@ void correct_velocity_shake(CellStructure &cs) { cs.ghosts_reduce_rattle_correction(); - apply_velocity_correction(particles); + apply_velocity_correction(particles, box_geo); cs.ghosts_update(Cells::DATA_PART_MOMENTUM); } diff --git a/src/core/rattle.hpp b/src/core/rattle.hpp index 41ceb8cb037..5502bd62558 100644 --- a/src/core/rattle.hpp +++ b/src/core/rattle.hpp @@ -31,6 +31,7 @@ #ifdef BOND_CONSTRAINT +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" /** Transfer the current particle positions from @ref Particle::pos @@ -43,17 +44,13 @@ void save_old_position(const ParticleRange &particles, /** * @brief Propagate velocity and position while using SHAKE algorithm for bond * constraint. - * - * @param cs cell structure */ -void correct_position_shake(CellStructure &cs); +void correct_position_shake(CellStructure &cs, BoxGeometry const &box_geo); /** * @brief Correction of current velocities using RATTLE algorithm. - * - * @param cs cell structure */ -void correct_velocity_shake(CellStructure &cs); +void correct_velocity_shake(CellStructure &cs, BoxGeometry const &box_geo); #endif #endif diff --git a/src/core/reaction_methods/ReactionAlgorithm.cpp b/src/core/reaction_methods/ReactionAlgorithm.cpp index 8035a1eb11c..67d23b6783a 100644 --- a/src/core/reaction_methods/ReactionAlgorithm.cpp +++ b/src/core/reaction_methods/ReactionAlgorithm.cpp @@ -21,11 +21,11 @@ #include "reaction_methods/ReactionAlgorithm.hpp" +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" #include "energy.hpp" #include "event.hpp" -#include "grid.hpp" #include "particle_node.hpp" #include "system/System.hpp" @@ -85,6 +85,7 @@ void ReactionAlgorithm::add_reaction( */ void ReactionAlgorithm::restore_old_system_state() { auto const &old_state = get_old_system_state(); + auto const &box_geo = *System::get_system().box_geo; // restore the properties of changed and hidden particles for (auto const &state : {old_state.changed, old_state.hidden}) { for (auto const &[p_id, p_type] : state) { @@ -125,7 +126,10 @@ void ReactionAlgorithm::restore_old_system_state() { * Automatically sets the volume which is used by the reaction ensemble to the * volume of a cuboid box. */ -void ReactionAlgorithm::update_volume() { volume = box_geo.volume(); } +void ReactionAlgorithm::update_volume() { + auto const &box_geo = *System::get_system().box_geo; + volume = box_geo.volume(); +} /** * Checks whether all particles exist for the provided reaction. @@ -339,6 +343,7 @@ void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) { if (p1_ptr != nullptr) { auto &p1 = *p1_ptr; auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; auto &cell_structure = *system.cell_structure; /* Check if the inserted particle within the exclusion radius of any other @@ -357,7 +362,7 @@ void ReactionAlgorithm::check_exclusion_range(int p_id, int p_type) { exclusion_radius_per_type[p2.type()]; } - auto const d_min = ::box_geo.get_mi_vector(p2.pos(), p1.pos()).norm(); + auto const d_min = box_geo.get_mi_vector(p2.pos(), p1.pos()).norm(); if (d_min < excluded_distance) { particle_inside_exclusion_range_touched = true; @@ -410,6 +415,7 @@ void ReactionAlgorithm::delete_particle(int p_id) { void ReactionAlgorithm::set_cyl_constraint(double center_x, double center_y, double radius) { + auto const &box_geo = *System::get_system().box_geo; if (center_x < 0. or center_x > box_geo.length()[0]) throw std::domain_error("center_x is outside the box"); if (center_y < 0. or center_y > box_geo.length()[1]) @@ -424,6 +430,7 @@ void ReactionAlgorithm::set_cyl_constraint(double center_x, double center_y, void ReactionAlgorithm::set_slab_constraint(double slab_start_z, double slab_end_z) { + auto const &box_geo = *System::get_system().box_geo; if (slab_start_z < 0. or slab_start_z > box_geo.length()[2]) throw std::domain_error("slab_start_z is outside the box"); if (slab_end_z < 0. or slab_end_z > box_geo.length()[2]) @@ -439,6 +446,7 @@ void ReactionAlgorithm::set_slab_constraint(double slab_start_z, * Writes a random position inside the central box into the provided array. */ Utils::Vector3d ReactionAlgorithm::get_random_position_in_box() { + auto const &box_geo = *System::get_system().box_geo; Utils::Vector3d out_pos{}; if (m_reaction_constraint == ReactionConstraint::CYL_Z) { diff --git a/src/core/scafacos/ScafacosContext.cpp b/src/core/scafacos/ScafacosContext.cpp index ef461e947a4..fc4fa1937d5 100644 --- a/src/core/scafacos/ScafacosContext.cpp +++ b/src/core/scafacos/ScafacosContext.cpp @@ -25,9 +25,9 @@ #include "scafacos/ScafacosContext.hpp" +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "communication.hpp" -#include "grid.hpp" #include "system/System.hpp" #include @@ -40,7 +40,9 @@ namespace detail { std::tuple get_system_params() { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; auto periodicity = Utils::Vector3i{static_cast(box_geo.periodic(0)), static_cast(box_geo.periodic(1)), static_cast(box_geo.periodic(2))}; diff --git a/src/core/stokesian_dynamics/sd_interface.cpp b/src/core/stokesian_dynamics/sd_interface.cpp index aaee5dbf88a..656e5833cd8 100644 --- a/src/core/stokesian_dynamics/sd_interface.cpp +++ b/src/core/stokesian_dynamics/sd_interface.cpp @@ -24,9 +24,10 @@ #include "stokesian_dynamics/sd_cpu.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "ParticleRange.hpp" -#include "grid.hpp" +#include "system/System.hpp" #include "thermostat.hpp" #include @@ -77,7 +78,8 @@ static double sd_kT = 0.0; static std::vector v_sd{}; void register_integrator(StokesianDynamicsParameters const &obj) { - if (::box_geo.periodic(0) or ::box_geo.periodic(1) or ::box_geo.periodic(2)) { + auto const &box_geo = *System::get_system().box_geo; + if (box_geo.periodic(0) or box_geo.periodic(1) or box_geo.periodic(2)) { throw std::runtime_error( "Stokesian Dynamics requires periodicity (False, False, False)"); } diff --git a/src/core/system/GpuParticleData.cpp b/src/core/system/GpuParticleData.cpp index 03718111fd0..c6129887344 100644 --- a/src/core/system/GpuParticleData.cpp +++ b/src/core/system/GpuParticleData.cpp @@ -84,10 +84,10 @@ BOOST_SERIALIZATION_SPLIT_FREE(GpuParticleData::GpuParticle) static void pack_particles(ParticleRange const &particles, GpuParticleData::GpuParticle *buffer) { - auto const &box_l = ::box_geo; + auto const &box = *System::get_system().box_geo; unsigned long int i = 0u; for (auto const &p : particles) { - buffer[i].p = static_cast(folded_position(p.pos(), box_l)); + buffer[i].p = static_cast(folded_position(p.pos(), box)); #ifdef DIPOLES buffer[i].dip = static_cast(p.calc_dip()); #endif diff --git a/src/core/system/System.cpp b/src/core/system/System.cpp index 20a15f0c949..cbb462559e3 100644 --- a/src/core/system/System.cpp +++ b/src/core/system/System.cpp @@ -21,7 +21,6 @@ #include "System.hpp" #include "System.impl.hpp" -#include "grid.hpp" #include @@ -32,7 +31,9 @@ namespace System { static std::shared_ptr instance = std::make_shared(); System::System() { - cell_structure = std::make_shared(::box_geo); + box_geo = std::make_shared(); + local_geo = std::make_shared(); + cell_structure = std::make_shared(*box_geo); } bool is_system_set() { return instance != nullptr; } @@ -45,7 +46,7 @@ void set_system(std::shared_ptr new_instance) { System &get_system() { return *instance; } -Utils::Vector3d System::box() const { return ::box_geo.length(); } +Utils::Vector3d System::box() const { return box_geo->length(); } void System::init() { #ifdef CUDA diff --git a/src/core/system/System.hpp b/src/core/system/System.hpp index 18a37c6b2d9..afae6c3807d 100644 --- a/src/core/system/System.hpp +++ b/src/core/system/System.hpp @@ -34,6 +34,8 @@ #include +class BoxGeometry; +class LocalBox; struct CellStructure; namespace System { @@ -53,6 +55,8 @@ class System { Dipoles::Solver dipoles; LB::Solver lb; EK::Solver ek; + std::shared_ptr box_geo; + std::shared_ptr local_geo; std::shared_ptr cell_structure; }; diff --git a/src/core/system/System.impl.hpp b/src/core/system/System.impl.hpp index a053342b436..4bf2d44044f 100644 --- a/src/core/system/System.impl.hpp +++ b/src/core/system/System.impl.hpp @@ -26,3 +26,6 @@ #include "lb/Implementation.hpp" #include "cell_system/CellStructure.hpp" + +#include "BoxGeometry.hpp" +#include "LocalBox.hpp" diff --git a/src/core/tuning.cpp b/src/core/tuning.cpp index 3781e67755d..c4d8f5395a9 100644 --- a/src/core/tuning.cpp +++ b/src/core/tuning.cpp @@ -138,7 +138,7 @@ void tune_skin(double min_skin, double max_skin, double tol, int int_steps, double const max_permissible_skin = std::min(*boost::min_element(cell_structure.max_cutoff()) - maximal_cutoff(::communicator.size), - 0.5 * *boost::max_element(box_geo.length())); + 0.5 * *boost::max_element(system.box_geo->length())); if (adjust_max_skin and max_skin > max_permissible_skin) b = max_permissible_skin; diff --git a/src/core/unit_tests/lb_particle_coupling_test.cpp b/src/core/unit_tests/lb_particle_coupling_test.cpp index d4f8cde86a8..32a829a9820 100644 --- a/src/core/unit_tests/lb_particle_coupling_test.cpp +++ b/src/core/unit_tests/lb_particle_coupling_test.cpp @@ -401,6 +401,7 @@ BOOST_DATA_TEST_CASE_F(CleanupActorLB, coupling_particle_lattice_ia, // check box shifts if (rank == 0) { + auto const &box_geo = *system.box_geo; auto constexpr reference_shifts = std::array{{{{0, 0, 0}}, {{0, 0, 8}}, diff --git a/src/core/virtual_sites.cpp b/src/core/virtual_sites.cpp index 28d341f4093..4311c841bed 100644 --- a/src/core/virtual_sites.cpp +++ b/src/core/virtual_sites.cpp @@ -25,6 +25,7 @@ #include "virtual_sites.hpp" +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "communication.hpp" #include "errorhandling.hpp" @@ -55,10 +56,10 @@ void set_virtual_sites(std::shared_ptr const &v) { /** Calculate the rotation quaternion and distance between two particles */ std::tuple, double> -calculate_vs_relate_to_params(Particle const &p_vs, - Particle const &p_relate_to) { +calculate_vs_relate_to_params(Particle const &p_vs, Particle const &p_relate_to, + BoxGeometry const &box_geo) { // get the distance between the particles - auto d = ::box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos()); + auto d = box_geo.get_mi_vector(p_vs.pos(), p_relate_to.pos()); // Check if the distance between virtual and non-virtual particles is larger // than minimum global cutoff. If so, warn user. diff --git a/src/core/virtual_sites.hpp b/src/core/virtual_sites.hpp index bf540b3beef..6aec6f6b55e 100644 --- a/src/core/virtual_sites.hpp +++ b/src/core/virtual_sites.hpp @@ -36,26 +36,30 @@ void set_virtual_sites(std::shared_ptr const &v); #ifdef VIRTUAL_SITES_RELATIVE +#include "BoxGeometry.hpp" + #include #include std::tuple, double> calculate_vs_relate_to_params(Particle const &p_current, - Particle const &p_relate_to); + Particle const &p_relate_to, + BoxGeometry const &box_geo); /** * @brief Setup a virtual site to track a real particle. * @param[in,out] p_vs Virtual site. * @param[in] p_relate_to Real particle to follow. */ -inline void vs_relate_to(Particle &p_vs, Particle const &p_relate_to) { +inline void vs_relate_to(Particle &p_vs, Particle const &p_relate_to, + BoxGeometry const &box_geo) { // Set the particle id of the particle we want to relate to, the distance // and the relative orientation auto &vs_relative = p_vs.vs_relative(); vs_relative.to_particle_id = p_relate_to.id(); std::tie(vs_relative.rel_orientation, vs_relative.distance) = - calculate_vs_relate_to_params(p_vs, p_relate_to); + calculate_vs_relate_to_params(p_vs, p_relate_to, box_geo); } #endif // VIRTUAL_SITES_RELATIVE diff --git a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp index cfb0fd1a1ee..a27268b95e1 100644 --- a/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp +++ b/src/core/virtual_sites/VirtualSitesInertialessTracers.cpp @@ -22,11 +22,11 @@ #include "VirtualSitesInertialessTracers.hpp" +#include "BoxGeometry.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" #include "errorhandling.hpp" #include "forces.hpp" -#include "grid.hpp" #include "integrate.hpp" #include "lb/particle_coupling.hpp" #include "system/System.hpp" @@ -52,6 +52,7 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { DeferredActiveLBChecks check_lb_solver_set{}; auto &system = System::get_system(); auto &cell_structure = *system.cell_structure; + auto const &box_geo = *system.box_geo; auto const agrid = (check_lb_solver_set()) ? system.lb.get_agrid() : 0.; auto const to_lb_units = (check_lb_solver_set()) ? 1. / agrid : 0.; @@ -72,7 +73,7 @@ void VirtualSitesInertialessTracers::after_force_calc(double time_step) { return; } if (bookkeeping.should_be_coupled(p)) { - for (auto pos : positions_in_halo(p.pos(), ::box_geo, agrid)) { + for (auto pos : positions_in_halo(p.pos(), box_geo, agrid)) { add_md_force(system.lb, pos * to_lb_units, p.force(), time_step); } } diff --git a/src/core/virtual_sites/VirtualSitesRelative.cpp b/src/core/virtual_sites/VirtualSitesRelative.cpp index c16e5ce8a16..1dd21f49943 100644 --- a/src/core/virtual_sites/VirtualSitesRelative.cpp +++ b/src/core/virtual_sites/VirtualSitesRelative.cpp @@ -21,12 +21,12 @@ #ifdef VIRTUAL_SITES_RELATIVE +#include "BoxGeometry.hpp" #include "Particle.hpp" #include "cell_system/CellStructure.hpp" #include "cells.hpp" #include "errorhandling.hpp" #include "forces.hpp" -#include "grid.hpp" #include "integrate.hpp" #include "rotation.hpp" #include "system/System.hpp" @@ -113,7 +113,9 @@ static auto constraint_stress(Particle const &p_ref, Particle const &p_vs) { } void VirtualSitesRelative::update() const { - auto &cell_structure = *System::get_system().cell_structure; + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + auto &cell_structure = *system.cell_structure; cell_structure.ghosts_update(Cells::DATA_PART_POSITION | Cells::DATA_PART_MOMENTUM); diff --git a/src/script_interface/analysis/Analysis.cpp b/src/script_interface/analysis/Analysis.cpp index d7f9f20f530..3152e946a21 100644 --- a/src/script_interface/analysis/Analysis.cpp +++ b/src/script_interface/analysis/Analysis.cpp @@ -19,16 +19,17 @@ #include "Analysis.hpp" +#include "core/BoxGeometry.hpp" #include "core/analysis/statistics.hpp" #include "core/analysis/statistics_chain.hpp" #include "core/cells.hpp" #include "core/dpd.hpp" #include "core/energy.hpp" #include "core/event.hpp" -#include "core/grid.hpp" #include "core/nonbonded_interactions/nonbonded_interaction_data.hpp" #include "core/partCfg_global.hpp" #include "core/particle_node.hpp" +#include "core/system/System.hpp" #include "script_interface/communication.hpp" @@ -217,9 +218,9 @@ Variant Analysis::do_call_method(std::string const &name, return make_vector_of_variants(result); } if (name == "distribution") { + auto const &box_l = System::get_system().box_geo->length(); auto const r_max_limit = - 0.5 * std::min(std::min(::box_geo.length()[0], ::box_geo.length()[1]), - ::box_geo.length()[2]); + 0.5 * std::min(std::min(box_l[0], box_l[1]), box_l[2]); auto const r_min = get_value_or(parameters, "r_min", 0.); auto const r_max = get_value_or(parameters, "r_max", r_max_limit); auto const r_bins = get_value_or(parameters, "r_bins", 100); diff --git a/src/script_interface/cell_system/CellSystem.cpp b/src/script_interface/cell_system/CellSystem.cpp index b03f6f32ce6..37336516a49 100644 --- a/src/script_interface/cell_system/CellSystem.cpp +++ b/src/script_interface/cell_system/CellSystem.cpp @@ -21,6 +21,7 @@ #include "script_interface/ScriptInterface.hpp" +#include "core/BoxGeometry.hpp" #include "core/bonded_interactions/bonded_interaction_data.hpp" #include "core/cell_system/CellStructure.hpp" #include "core/cell_system/HybridDecomposition.hpp" @@ -28,7 +29,6 @@ #include "core/cells.hpp" #include "core/communication.hpp" #include "core/event.hpp" -#include "core/grid.hpp" #include "core/integrate.hpp" #include "core/nonbonded_interactions/nonbonded_interaction_data.hpp" #include "core/particle_node.hpp" @@ -245,7 +245,7 @@ Variant CellSystem::do_call_method(std::string const &name, std::vector CellSystem::mpi_resort_particles(bool global_flag) const { auto &cell_structure = get_cell_structure(); - cell_structure.resort_particles(global_flag, box_geo); + cell_structure.resort_particles(global_flag, *::System::get_system().box_geo); clear_particle_node(); auto const size = static_cast(cell_structure.local_particles().size()); std::vector n_part_per_node; diff --git a/src/script_interface/lees_edwards/LeesEdwards.hpp b/src/script_interface/lees_edwards/LeesEdwards.hpp index 92a030cf4ca..579b4493f6f 100644 --- a/src/script_interface/lees_edwards/LeesEdwards.hpp +++ b/src/script_interface/lees_edwards/LeesEdwards.hpp @@ -21,7 +21,7 @@ #include "Protocol.hpp" -#include "core/grid.hpp" +#include "core/BoxGeometry.hpp" #include "core/lees_edwards/LeesEdwardsBC.hpp" #include "core/lees_edwards/lees_edwards.hpp" #include "core/system/System.hpp" @@ -37,21 +37,23 @@ namespace LeesEdwards { class LeesEdwards : public AutoParameters { std::shared_ptr m_protocol; - LeesEdwardsBC const &m_lebc = ::box_geo.lees_edwards_bc(); + LeesEdwardsBC const &m_lebc; public: - LeesEdwards() : m_protocol{nullptr} { + LeesEdwards() + : m_protocol{nullptr}, + m_lebc{System::get_system().box_geo->lees_edwards_bc()} { add_parameters( {{"protocol", [this](Variant const &value) { if (is_none(value)) { - context()->parallel_try_catch([]() { + auto const &system = System::get_system(); + context()->parallel_try_catch([&system]() { auto constexpr invalid_dir = LeesEdwardsBC::invalid_dir; - auto const &system = System::get_system(); system.lb.lebc_sanity_checks(invalid_dir, invalid_dir); }); m_protocol = nullptr; - ::box_geo.set_lees_edwards_bc(LeesEdwardsBC{}); + system.box_geo->set_lees_edwards_bc(LeesEdwardsBC{}); ::LeesEdwards::unset_protocol(); return; } @@ -104,7 +106,7 @@ class LeesEdwards : public AutoParameters { auto const &system = System::get_system(); system.lb.lebc_sanity_checks(shear_direction, shear_plane_normal); // update box geometry and cell structure - ::box_geo.set_lees_edwards_bc( + system.box_geo->set_lees_edwards_bc( LeesEdwardsBC{0., 0., shear_direction, shear_plane_normal}); ::LeesEdwards::set_protocol(m_protocol->protocol()); }); diff --git a/src/script_interface/particle_data/ParticleHandle.cpp b/src/script_interface/particle_data/ParticleHandle.cpp index d2c179877c5..cb451eb60b1 100644 --- a/src/script_interface/particle_data/ParticleHandle.cpp +++ b/src/script_interface/particle_data/ParticleHandle.cpp @@ -24,6 +24,7 @@ #include "script_interface/Variant.hpp" #include "script_interface/get_value.hpp" +#include "core/BoxGeometry.hpp" #include "core/bonded_interactions/bonded_interaction_data.hpp" #include "core/cell_system/CellStructure.hpp" #include "core/event.hpp" @@ -197,7 +198,8 @@ ParticleHandle::ParticleHandle() { auto const p = get_particle_data(m_pid); auto const pos = p.pos(); auto const image_box = p.image_box(); - return unfolded_position(pos, image_box, ::box_geo.length()); + auto const &system = System::get_system(); + return unfolded_position(pos, image_box, system.box_geo->length()); }}, {"v", [this](Variant const &value) { @@ -394,7 +396,9 @@ ParticleHandle::ParticleHandle() { #endif // THERMOSTAT_PER_PARTICLE {"pos_folded", AutoParameter::read_only, [this]() { - return folded_position(get_particle_data(m_pid).pos(), ::box_geo); + auto const &system = System::get_system(); + auto const &box_geo = *system.box_geo; + return folded_position(get_particle_data(m_pid).pos(), box_geo); }}, {"lees_edwards_offset", @@ -611,8 +615,8 @@ Variant ParticleHandle::do_call_method(std::string const &name, */ auto const &p_current = get_particle_data(m_pid); auto const &p_relate_to = get_particle_data(other_pid); - auto const [quat, dist] = - calculate_vs_relate_to_params(p_current, p_relate_to); + auto const [quat, dist] = calculate_vs_relate_to_params( + p_current, p_relate_to, *::System::get_system().box_geo); set_parameter("vs_relative", Variant{std::vector{ {other_pid, dist, quat2vector(quat)}}}); set_parameter("virtual", true); diff --git a/src/script_interface/system/System.cpp b/src/script_interface/system/System.cpp index 31c4d80f433..b1ada840926 100644 --- a/src/script_interface/system/System.cpp +++ b/src/script_interface/system/System.cpp @@ -19,11 +19,11 @@ #include "System.hpp" +#include "core/BoxGeometry.hpp" #include "core/cell_system/CellStructure.hpp" #include "core/cell_system/CellStructureType.hpp" #include "core/cells.hpp" #include "core/event.hpp" -#include "core/grid.hpp" #include "core/nonbonded_interactions/nonbonded_interaction_data.hpp" #include "core/object-in-fluid/oif_global_forces.hpp" #include "core/particle_node.hpp" @@ -45,12 +45,14 @@ namespace System { static bool system_created = false; static Utils::Vector3b get_periodicity() { - return {::box_geo.periodic(0), ::box_geo.periodic(1), ::box_geo.periodic(2)}; + auto const &box_geo = *::System::get_system().box_geo; + return {box_geo.periodic(0), box_geo.periodic(1), box_geo.periodic(2)}; } static void set_periodicity(Utils::Vector3b const &value) { + auto &box_geo = *::System::get_system().box_geo; for (int i = 0; i < 3; ++i) { - ::box_geo.set_periodic(i, value[i]); + box_geo.set_periodic(i, value[i]); } on_periodicity_change(); } @@ -64,11 +66,12 @@ System::System() { if (not(new_value > Utils::Vector3d::broadcast(0.))) { throw std::domain_error("Attribute 'box_l' must be > 0"); } - box_geo.set_length(new_value); + puts("box_geo.set_length()"); + m_instance->box_geo->set_length(new_value); on_boxl_change(); }); }, - []() { return ::box_geo.length(); }}, + []() { return ::System::get_system().box_geo->length(); }}, {"min_global_cut", [this](Variant const &v) { context()->parallel_try_catch([&]() { @@ -104,6 +107,7 @@ void System::do_construct(VariantMap const ¶ms) { std::set const skip{{"box_l", "periodicity", "min_global_cut"}}; std::set const allowed{valid_params.begin(), valid_params.end()}; + puts("System::do_construct()"); m_instance = std::make_shared<::System::System>(); ::System::set_system(m_instance); m_instance->init(); @@ -112,6 +116,7 @@ void System::do_construct(VariantMap const ¶ms) { auto &cell_structure = *m_instance->cell_structure; cells_re_init(cell_structure, CellStructureType::CELL_STRUCTURE_NSQUARE); do_set_parameter("box_l", params.at("box_l")); + std::cout << " boxl=[" << m_instance->box() << "]\n"; cells_re_init(cell_structure, CellStructureType::CELL_STRUCTURE_REGULAR); if (params.count("periodicity")) { @@ -154,10 +159,11 @@ Variant System::do_call_method(std::string const &name, return {}; } if (name == "rescale_boxl") { + auto &box_geo = *::System::get_system().box_geo; auto const coord = get_value(parameters, "coord"); auto const length = get_value(parameters, "length"); - auto const scale = (coord == 3) ? length * ::box_geo.length_inv()[0] - : length * ::box_geo.length_inv()[coord]; + auto const scale = (coord == 3) ? length * box_geo.length_inv()[0] + : length * box_geo.length_inv()[coord]; context()->parallel_try_catch([&]() { if (length <= 0.) { throw std::domain_error("Parameter 'd_new' be > 0"); @@ -167,14 +173,15 @@ Variant System::do_call_method(std::string const &name, if (coord == 3) { new_value = Utils::Vector3d::broadcast(length); } else { - new_value = ::box_geo.length(); + new_value = box_geo.length(); new_value[coord] = length; } // when shrinking, rescale the particles first if (scale <= 1.) { rescale_particles(*m_instance->cell_structure, coord, scale); } - set_box_length(new_value); + m_instance->box_geo->set_length(new_value); + on_boxl_change(); if (scale > 1.) { rescale_particles(*m_instance->cell_structure, coord, scale); } @@ -192,16 +199,18 @@ Variant System::do_call_method(std::string const &name, return ::number_of_particles_with_type(type); } if (name == "velocity_difference") { + auto const &box_geo = *::System::get_system().box_geo; auto const pos1 = get_value(parameters, "pos1"); auto const pos2 = get_value(parameters, "pos2"); auto const v1 = get_value(parameters, "v1"); auto const v2 = get_value(parameters, "v2"); - return ::box_geo.velocity_difference(pos2, pos1, v2, v1); + return box_geo.velocity_difference(pos2, pos1, v2, v1); } if (name == "distance_vec") { + auto const &box_geo = *::System::get_system().box_geo; auto const pos1 = get_value(parameters, "pos1"); auto const pos2 = get_value(parameters, "pos2"); - return ::box_geo.get_mi_vector(pos2, pos1); + return box_geo.get_mi_vector(pos2, pos1); } if (name == "rotate_system") { rotate_system(*m_instance->cell_structure, diff --git a/src/script_interface/walberla/LBFluid.cpp b/src/script_interface/walberla/LBFluid.cpp index 7c1a5a93a44..7dfa6db2c78 100644 --- a/src/script_interface/walberla/LBFluid.cpp +++ b/src/script_interface/walberla/LBFluid.cpp @@ -26,7 +26,6 @@ #include "core/BoxGeometry.hpp" #include "core/event.hpp" -#include "core/grid.hpp" #include "core/integrate.hpp" #include "core/lb/LBWalberla.hpp" #include "core/lees_edwards/lees_edwards.hpp" @@ -83,6 +82,7 @@ Variant LBFluid::do_call_method(std::string const &name, return {}; } if (name == "add_force_at_pos") { + auto const &box_geo = *::System::get_system().box_geo; auto const pos = get_value(params, "pos"); auto const f = get_value(params, "force"); auto const folded_pos = folded_position(pos, box_geo); @@ -177,7 +177,7 @@ void LBFluid::do_construct(VariantMap const ¶ms) { throw std::runtime_error( "Lees-Edwards LB doesn't support thermalization"); } - auto const &le_bc = ::box_geo.lees_edwards_bc(); + auto const &le_bc = ::System::get_system().box_geo->lees_edwards_bc(); auto lees_edwards_object = std::make_unique( le_bc.shear_direction, le_bc.shear_plane_normal, [this, le_protocol]() { @@ -211,6 +211,7 @@ std::vector LBFluid::get_average_pressure_tensor() const { } Variant LBFluid::get_interpolated_velocity(Utils::Vector3d const &pos) const { + auto const &box_geo = *::System::get_system().box_geo; auto const lb_pos = folded_position(pos, box_geo) * m_conv_dist; auto const result = m_instance->get_velocity_at_pos(lb_pos); return Utils::Mpi::reduce_optional(context()->get_comm(), result) / diff --git a/src/script_interface/walberla/LatticeWalberla.hpp b/src/script_interface/walberla/LatticeWalberla.hpp index 2a159f9f54e..6c9c99a9de8 100644 --- a/src/script_interface/walberla/LatticeWalberla.hpp +++ b/src/script_interface/walberla/LatticeWalberla.hpp @@ -23,8 +23,9 @@ #ifdef WALBERLA +#include "core/BoxGeometry.hpp" #include "core/communication.hpp" -#include "core/grid.hpp" +#include "core/system/System.hpp" #include @@ -56,8 +57,9 @@ class LatticeWalberla : public AutoParameters { } void do_construct(VariantMap const &args) override { + auto const &box_geo = *System::get_system().box_geo; m_agrid = get_value(args, "agrid"); - m_box_l = get_value_or(args, "_box_l", ::box_geo.length()); + m_box_l = get_value_or(args, "_box_l", box_geo.length()); auto const n_ghost_layers = get_value(args, "n_ghost_layers"); context()->parallel_try_catch([&]() {