Encapsulate local_box and box_geo globals

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 <http://www.gnu.org/licenses/>.
  */
-#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

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(

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 <http://www.gnu.org/licenses/>.
  */
-#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

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 <utils/Span.hpp>
 
@@ -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();

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<int> const &set1,
                std::vector<int> const &set2) {
   using Utils::contains;
 
+  auto const &box_geo = *System::get_system().box_geo;
   auto mindist_sq = std::numeric_limits<double>::infinity();
 
   for (auto jt = partCfg.begin(); jt != partCfg.end(); ++jt) {
@@ -142,6 +144,7 @@ std::vector<int> nbhood(PartCfg &partCfg, Utils::Vector3d const &pos,
                         double dist) {
   std::vector<int> 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<int> const &p1_types,
                        std::vector<int> 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<int> 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<std::size_t>(order));
   auto const twoPI_L = 2. * Utils::pi() * box_geo.length_inv()[0];
   std::vector<double> ff(2 * order_sq + 1);

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 <utils/Vector.hpp>
 #include <utils/math/sqr.hpp>
@@ -36,6 +38,7 @@
 #include <stdexcept>
 
 std::array<double, 4> 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<double, 4> re;
 
@@ -61,6 +64,7 @@ std::array<double, 4> calc_re(int chain_start, int n_chains, int chain_length) {
 }
 
 std::array<double, 4> 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<double, 4> rg;
 
@@ -101,6 +105,7 @@ std::array<double, 4> calc_rg(int chain_start, int n_chains, int chain_length) {
 }
 
 std::array<double, 2> 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<double, 2> rh;
 

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<Utils::Vector3d, Utils::Vector3d, double, double, double>
-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 <typename ForceFactor>
 std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-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 */

src/core/bonded_interactions/angle_cosine.hpp

@@ -26,6 +26,7 @@
  *  @ref bondedIA_angle_cosine.
  */
 
+#include "BoxGeometry.hpp"
 #include "angle_common.hpp"
 
 #include <utils/Vector.hpp>
@@ -52,9 +53,10 @@ struct AngleCosineBond {
   AngleCosineBond(double bend, double phi0);
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-  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<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-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,18 +89,21 @@ 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.
  *  @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 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)]

src/core/bonded_interactions/angle_cossquare.hpp

@@ -49,9 +49,10 @@ struct AngleCossquareBond {
   AngleCossquareBond(double bend, double phi0);
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-  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,26 +72,30 @@ struct AngleCossquareBond {
  *  @return Forces on the second, first and third particles, in that order.
  */
 inline std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-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 {
 
   auto forceFactor = [this](double const cos_phi) {
     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.
  *  @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 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);
 }

src/core/bonded_interactions/angle_harmonic.hpp

@@ -51,9 +51,10 @@ struct AngleHarmonicBond {
   }
 
   std::tuple<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-  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<Utils::Vector3d, Utils::Vector3d, Utils::Vector3d>
-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,18 +84,21 @@ 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.
  *  @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 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);