diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp
index 98358dc5f7a..36a44813e53 100644
--- a/src/core/electrostatics_magnetostatics/elc.cpp
+++ b/src/core/electrostatics_magnetostatics/elc.cpp
@@ -50,11 +50,6 @@
 #include <cstddef>
 #include <vector>
 
-/** \name Inverse box dimensions */
-/**@{*/
-static double ux, uy, uz;
-/**@}*/
-
 ELC_struct elc_params = {1e100, 10,    1, 0, true, true, false, 1,
                          1,     false, 0, 0, 0,    0,    0.0};
 
@@ -107,12 +102,6 @@ static double dipole_energy(const ParticleRange &particles);
 static double z_energy(const ParticleRange &particles);
 static void add_z_force(const ParticleRange &particles);
 
-void ELC_setup_constants() {
-  ux = 1 / box_geo.length()[0];
-  uy = 1 / box_geo.length()[1];
-  uz = 1 / box_geo.length()[2];
-}
-
 /**
  * @brief Calculate cached sin/cos values for one direction.
  *
@@ -215,14 +204,16 @@ inline void check_gap_elc(const Particle &p) {
  *  See @cite yeh99a.
  */
 static void add_dipole_force(const ParticleRange &particles) {
-  double const pref = coulomb.prefactor * 4 * Utils::pi() * ux * uy * uz;
+  double const pref = coulomb.prefactor * 4 * Utils::pi() *
+                      box_geo.length_inv()[0] * box_geo.length_inv()[1] *
+                      box_geo.length_inv()[2];
   constexpr std::size_t size = 3;
 
   auto local_particles = particles;
 
   /* for nonneutral systems, this shift gives the background contribution
      (rsp. for this shift, the DM of the background is zero) */
-  double const shift = 0.5 * box_geo.length()[2];
+  double const shift = box_geo.length_half()[2];
 
   // collect moments
 
@@ -251,7 +242,7 @@ static void add_dipole_force(const ParticleRange &particles) {
   }
 
   gblcblk[0] *= pref;
-  gblcblk[1] *= pref / elc_params.h / uz;
+  gblcblk[1] *= pref / elc_params.h * box_geo.length()[2];
   gblcblk[2] *= pref;
 
   distribute(size);
@@ -280,11 +271,13 @@ static void add_dipole_force(const ParticleRange &particles) {
  *  See @cite yeh99a.
  */
 static double dipole_energy(const ParticleRange &particles) {
-  double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy * uz;
+  double const pref = coulomb.prefactor * 2 * Utils::pi() *
+                      box_geo.length_inv()[0] * box_geo.length_inv()[1] *
+                      box_geo.length_inv()[2];
   constexpr std::size_t size = 7;
   /* for nonneutral systems, this shift gives the background contribution
      (rsp. for this shift, the DM of the background is zero) */
-  double const shift = 0.5 * box_geo.length()[2];
+  double const shift = box_geo.length_half()[2];
 
   // collect moments
 
@@ -336,7 +329,8 @@ static double dipole_energy(const ParticleRange &particles) {
   if (elc_params.dielectric_contrast_on) {
     if (elc_params.const_pot) {
       // zero potential difference contribution
-      energy += pref / elc_params.h / uz * Utils::sqr(gblcblk[6]);
+      energy +=
+          pref / elc_params.h * box_geo.length()[2] * Utils::sqr(gblcblk[6]);
       // external potential shift contribution
       energy -= 2 * elc_params.pot_diff / elc_params.h * gblcblk[6];
     }
@@ -355,7 +349,7 @@ static double dipole_energy(const ParticleRange &particles) {
 /*****************************************************************/
 
 inline double image_sum_b(double q, double z) {
-  double const shift = 0.5 * box_geo.length()[2];
+  double const shift = box_geo.length_half()[2];
   double const fac = elc_params.delta_mid_top * elc_params.delta_mid_bot;
   double const image_sum =
       (q / (1.0 - fac) * (z - 2.0 * fac * box_geo.length()[2] / (1.0 - fac))) -
@@ -364,7 +358,7 @@ inline double image_sum_b(double q, double z) {
 }
 
 inline double image_sum_t(double q, double z) {
-  double const shift = 0.5 * box_geo.length()[2];
+  double const shift = box_geo.length_half()[2];
   double const fac = elc_params.delta_mid_top * elc_params.delta_mid_bot;
   double const image_sum =
       (q / (1.0 - fac) * (z + 2.0 * fac * box_geo.length()[2] / (1.0 - fac))) -
@@ -374,12 +368,13 @@ inline double image_sum_t(double q, double z) {
 
 /*****************************************************************/
 static double z_energy(const ParticleRange &particles) {
-  double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
+  double const pref = coulomb.prefactor * 2 * Utils::pi() *
+                      box_geo.length_inv()[0] * box_geo.length_inv()[1];
   constexpr std::size_t size = 4;
 
   /* for nonneutral systems, this shift gives the background contribution
      (rsp. for this shift, the DM of the background is zero) */
-  double const shift = 0.5 * box_geo.length()[2];
+  double const shift = box_geo.length_half()[2];
 
   if (elc_params.dielectric_contrast_on) {
     if (elc_params.const_pot) {
@@ -452,7 +447,8 @@ static double z_energy(const ParticleRange &particles) {
 
 /*****************************************************************/
 static void add_z_force(const ParticleRange &particles) {
-  double const pref = coulomb.prefactor * 2 * Utils::pi() * ux * uy;
+  double const pref = coulomb.prefactor * 2 * Utils::pi() *
+                      box_geo.length_inv()[0] * box_geo.length_inv()[1];
   constexpr std::size_t size = 1;
 
   if (elc_params.dielectric_contrast_on) {
@@ -512,7 +508,8 @@ template <PoQ axis>
 void setup_PoQ(std::size_t index, double omega,
                const ParticleRange &particles) {
   assert(index >= 1);
-  double const pref_di = coulomb.prefactor * 4 * Utils::pi() * ux * uy;
+  double const pref_di = coulomb.prefactor * 4 * Utils::pi() *
+                         box_geo.length_inv()[0] * box_geo.length_inv()[1];
   double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
   constexpr std::size_t size = 4;
   double lclimgebot[4], lclimgetop[4], lclimge[4];
@@ -659,7 +656,8 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega,
                      const ParticleRange &particles) {
   assert(index_p >= 1);
   assert(index_q >= 1);
-  double const pref_di = coulomb.prefactor * 8 * Utils::pi() * ux * uy;
+  double const pref_di = coulomb.prefactor * 8 * Utils::pi() *
+                         box_geo.length_inv()[0] * box_geo.length_inv()[1];
   double const pref = -pref_di / expm1(omega * box_geo.length()[2]);
   constexpr std::size_t size = 8;
   double lclimgebot[8], lclimgetop[8], lclimge[8];
@@ -789,8 +787,10 @@ static void setup_PQ(std::size_t index_p, std::size_t index_q, double omega,
 static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega,
                          const ParticleRange &particles) {
   constexpr double c_2pi = 2 * Utils::pi();
-  double const pref_x = c_2pi * ux * static_cast<double>(index_p) / omega;
-  double const pref_y = c_2pi * uy * static_cast<double>(index_q) / omega;
+  double const pref_x =
+      c_2pi * box_geo.length_inv()[0] * static_cast<double>(index_p) / omega;
+  double const pref_y =
+      c_2pi * box_geo.length_inv()[1] * static_cast<double>(index_q) / omega;
   constexpr std::size_t size = 8;
 
   std::size_t ic = 0;
@@ -847,45 +847,54 @@ static double PQ_energy(double omega, std::size_t n_part) {
 
 void ELC_add_force(const ParticleRange &particles) {
   constexpr double c_2pi = 2 * Utils::pi();
-  auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1);
-  auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1);
+  auto const n_scxcache =
+      std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1);
+  auto const n_scycache =
+      std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1);
 
-  prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy);
+  prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache,
+                   box_geo.length_inv()[1]);
   partblk.resize(particles.size() * 8);
 
   add_dipole_force(particles);
   add_z_force(particles);
 
   /* the second condition is just for the case of numerical accident */
-  for (std::size_t p = 1;
-       ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
+  for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
+                              elc_params.far_cut &&
+                          p <= n_scxcache;
        p++) {
-    auto const omega = c_2pi * ux * static_cast<double>(p);
+    auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
     setup_PoQ<PoQ::P>(p, omega, particles);
     distribute(4);
     add_PoQ_force<PoQ::P>(particles);
   }
 
-  for (std::size_t q = 1;
-       uy * static_cast<double>(q - 1) < elc_params.far_cut && q <= n_scycache;
+  for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast<double>(q - 1) <
+                              elc_params.far_cut &&
+                          q <= n_scycache;
        q++) {
-    auto const omega = c_2pi * uy * static_cast<double>(q);
+    auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
     setup_PoQ<PoQ::Q>(q, omega, particles);
     distribute(4);
     add_PoQ_force<PoQ::Q>(particles);
   }
 
-  for (std::size_t p = 1;
-       ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
+  for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
+                              elc_params.far_cut &&
+                          p <= n_scxcache;
        p++) {
     for (std::size_t q = 1;
-         Utils::sqr(ux * static_cast<double>(p - 1)) +
-                 Utils::sqr(uy * static_cast<double>(q - 1)) <
+         Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
+                 Utils::sqr(box_geo.length_inv()[1] *
+                            static_cast<double>(q - 1)) <
              elc_params.far_cut2 &&
          q <= n_scycache;
          q++) {
-      auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast<double>(p)) +
-                                      Utils::sqr(uy * static_cast<double>(q)));
+      auto const omega =
+          c_2pi *
+          sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
+               Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
       setup_PQ(p, q, omega, particles);
       distribute(8);
       add_PQ_force(p, q, omega, particles);
@@ -898,43 +907,52 @@ double ELC_energy(const ParticleRange &particles) {
   auto energy = dipole_energy(particles);
   energy += z_energy(particles);
 
-  auto const n_scxcache = std::size_t(ceil(elc_params.far_cut / ux) + 1);
-  auto const n_scycache = std::size_t(ceil(elc_params.far_cut / uy) + 1);
-  prepare_sc_cache(particles, n_scxcache, ux, n_scycache, uy);
+  auto const n_scxcache =
+      std::size_t(ceil(elc_params.far_cut * box_geo.length()[0]) + 1);
+  auto const n_scycache =
+      std::size_t(ceil(elc_params.far_cut * box_geo.length()[1]) + 1);
+  prepare_sc_cache(particles, n_scxcache, box_geo.length_inv()[0], n_scycache,
+                   box_geo.length_inv()[1]);
 
   auto const n_localpart = particles.size();
   partblk.resize(n_localpart * 8);
 
   /* the second condition is just for the case of numerical accident */
-  for (std::size_t p = 1;
-       ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
+  for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
+                              elc_params.far_cut &&
+                          p <= n_scxcache;
        p++) {
-    auto const omega = c_2pi * ux * static_cast<double>(p);
+    auto const omega = c_2pi * box_geo.length_inv()[0] * static_cast<double>(p);
     setup_PoQ<PoQ::P>(p, omega, particles);
     distribute(4);
     energy += PoQ_energy(omega, n_localpart);
   }
 
-  for (std::size_t q = 1;
-       uy * static_cast<double>(q - 1) < elc_params.far_cut && q <= n_scycache;
+  for (std::size_t q = 1; box_geo.length_inv()[1] * static_cast<double>(q - 1) <
+                              elc_params.far_cut &&
+                          q <= n_scycache;
        q++) {
-    auto const omega = c_2pi * uy * static_cast<double>(q);
+    auto const omega = c_2pi * box_geo.length_inv()[1] * static_cast<double>(q);
     setup_PoQ<PoQ::Q>(q, omega, particles);
     distribute(4);
     energy += PoQ_energy(omega, n_localpart);
   }
 
-  for (std::size_t p = 1;
-       ux * static_cast<double>(p - 1) < elc_params.far_cut && p <= n_scxcache;
+  for (std::size_t p = 1; box_geo.length_inv()[0] * static_cast<double>(p - 1) <
+                              elc_params.far_cut &&
+                          p <= n_scxcache;
        p++) {
     for (std::size_t q = 1;
-         Utils::sqr(ux * static_cast<double>(p - 1)) +
-                 Utils::sqr(uy * static_cast<double>(q - 1)) <
+         Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p - 1)) +
+                 Utils::sqr(box_geo.length_inv()[1] *
+                            static_cast<double>(q - 1)) <
              elc_params.far_cut2 &&
          q <= n_scycache;
          q++) {
-      auto const omega = c_2pi * sqrt(Utils::sqr(ux * static_cast<double>(p)) +
-                                      Utils::sqr(uy * static_cast<double>(q)));
+      auto const omega =
+          c_2pi *
+          sqrt(Utils::sqr(box_geo.length_inv()[0] * static_cast<double>(p)) +
+               Utils::sqr(box_geo.length_inv()[1] * static_cast<double>(q)));
       setup_PQ(p, q, omega, particles);
       distribute(8);
       energy += PQ_energy(omega, n_localpart);
@@ -949,7 +967,8 @@ double ELC_tune_far_cut(ELC_struct const &params) {
   constexpr auto maximal_far_cut = 50.;
   double const h = params.h;
   double lz = box_geo.length()[2];
-  double const min_inv_boxl = std::min(ux, uy);
+  double const min_inv_boxl =
+      std::min(box_geo.length_inv()[0], box_geo.length_inv()[1]);
 
   if (params.dielectric_contrast_on) {
     // adjust lz according to dielectric layer method
@@ -965,7 +984,8 @@ double ELC_tune_far_cut(ELC_struct const &params) {
   do {
     const auto prefactor = 2 * Utils::pi() * far_cut;
 
-    const auto sum = prefactor + 2 * (ux + uy);
+    const auto sum =
+        prefactor + 2 * (box_geo.length_inv()[0] + box_geo.length_inv()[1]);
     const auto den = -expm1(-prefactor * lz);
     const auto num1 = exp(prefactor * (h - lz));
     const auto num2 = exp(-prefactor * (h + lz));
@@ -1018,8 +1038,6 @@ void ELC_sanity_checks(ELC_struct const &params) {
 }
 
 void ELC_init() {
-
-  ELC_setup_constants();
   elc_params.h = box_geo.length()[2] - elc_params.gap_size;
 
   if (elc_params.dielectric_contrast_on) {
@@ -1107,7 +1125,6 @@ void ELC_set_params(double maxPWerror, double gap_size, double far_cut,
 
   ELC_sanity_checks(new_elc_params);
 
-  ELC_setup_constants();
   if (new_elc_params.far_calculated) {
     new_elc_params.far_cut = ELC_tune_far_cut(new_elc_params);
   }
diff --git a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
index 6e7f4bde1b3..50b7ba5bad4 100644
--- a/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
+++ b/src/core/electrostatics_magnetostatics/mdlc_correction.cpp
@@ -126,8 +126,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &fs,
   double s1z, s2z, s3z, s4z;
   double ss;
 
-  auto const facux = 2.0 * Utils::pi() / box_geo.length()[0];
-  auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1];
+  auto const facux = 2.0 * Utils::pi() * box_geo.length_inv()[0];
+  auto const facuy = 2.0 * Utils::pi() * box_geo.length_inv()[1];
   double energy = 0.0;
 
   for (int ix = -kcut; ix <= +kcut; ix++) {
@@ -243,7 +243,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &fs,
 
   // Multiply by the factors we have left during the loops
 
-  auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]);
+  auto const piarea =
+      Utils::pi() * box_geo.length_inv()[0] * box_geo.length_inv()[1];
 
   for (int j = 0; j < n_local_particles; j++) {
     fs[j] *= piarea;
@@ -259,8 +260,8 @@ double get_DLC_dipolar(int kcut, std::vector<Utils::Vector3d> &fs,
  *  %Algorithm implemented accordingly to @cite brodka04a.
  */
 double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) {
-  auto const facux = 2.0 * Utils::pi() / box_geo.length()[0];
-  auto const facuy = 2.0 * Utils::pi() / box_geo.length()[1];
+  auto const facux = 2.0 * Utils::pi() * box_geo.length_inv()[0];
+  auto const facuy = 2.0 * Utils::pi() * box_geo.length_inv()[1];
 
   double energy = 0.0;
   for (int ix = -kcut; ix <= +kcut; ix++) {
@@ -314,7 +315,8 @@ double get_DLC_energy_dipolar(int kcut, const ParticleRange &particles) {
 
   // Multiply by the factors we have left during the loops
 
-  auto const piarea = Utils::pi() / (box_geo.length()[0] * box_geo.length()[1]);
+  auto const piarea =
+      Utils::pi() * box_geo.length_inv()[0] * box_geo.length_inv()[1];
   energy *= (-piarea);
   return (this_node == 0) ? energy : 0.0;
 }
diff --git a/src/core/electrostatics_magnetostatics/mmm1d.cpp b/src/core/electrostatics_magnetostatics/mmm1d.cpp
index 3204db5018d..1453fcec6d5 100644
--- a/src/core/electrostatics_magnetostatics/mmm1d.cpp
+++ b/src/core/electrostatics_magnetostatics/mmm1d.cpp
@@ -75,7 +75,7 @@
 
 /** @name inverse box dimensions and other constants */
 /**@{*/
-static double uz, L2, uz2, prefuz2, prefL3_i;
+static double uz2, prefuz2, prefL3_i;
 /**@}*/
 
 MMM1D_struct mmm1d_params = {0.05, 1e-5, 0};
@@ -85,8 +85,9 @@ static std::vector<double> bessel_radii;
 
 static double far_error(int P, double minrad) {
   // this uses an upper bound to all force components and the potential
-  auto const rhores = 2 * Utils::pi() * uz * minrad;
-  auto const pref = 4 * uz * std::max(1.0, 2 * Utils::pi() * uz);
+  auto const rhores = 2 * Utils::pi() * box_geo.length_inv()[2] * minrad;
+  auto const pref = 4 * box_geo.length_inv()[2] *
+                    std::max(1.0, 2 * Utils::pi() * box_geo.length_inv()[2]);
 
   return pref * K1(rhores * P) * exp(rhores) / rhores * (P - 1 + 1 / rhores);
 }
@@ -171,11 +172,9 @@ int MMM1D_init() {
   if (mmm1d_params.far_switch_radius_2 >= Utils::sqr(box_geo.length()[2]))
     mmm1d_params.far_switch_radius_2 = 0.8 * Utils::sqr(box_geo.length()[2]);
 
-  uz = 1 / box_geo.length()[2];
-  L2 = box_geo.length()[2] * box_geo.length()[2];
-  uz2 = uz * uz;
+  uz2 = Utils::sqr(box_geo.length_inv()[2]);
   prefuz2 = coulomb.prefactor * uz2;
-  prefL3_i = prefuz2 * uz;
+  prefL3_i = prefuz2 * box_geo.length_inv()[2];
 
   determine_bessel_radii(mmm1d_params.maxPWerror, MAXIMAL_B_CUT);
   prepare_polygamma_series(mmm1d_params.maxPWerror,
@@ -189,7 +188,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d,
   auto const n_modPsi = static_cast<int>(modPsi.size() >> 1);
   auto const rxy2 = d[0] * d[0] + d[1] * d[1];
   auto const rxy2_d = rxy2 * uz2;
-  auto const z_d = d[2] * uz;
+  auto const z_d = d[2] * box_geo.length_inv()[2];
   Utils::Vector3d F;
 
   if (rxy2 <= mmm1d_params.far_switch_radius_2) {
@@ -245,7 +244,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d,
   } else {
     /* far range formula */
     auto const rxy = sqrt(rxy2);
-    auto const rxy_d = rxy * uz;
+    auto const rxy_d = rxy * box_geo.length_inv()[2];
     double sr = 0, sz = 0;
 
     for (int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
@@ -266,7 +265,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d,
     sr *= uz2 * 4 * c_2pi;
     sz *= uz2 * 4 * c_2pi;
 
-    auto const pref = sr / rxy + 2 * uz / rxy2;
+    auto const pref = sr / rxy + 2 * box_geo.length_inv()[2] / rxy2;
 
     F = {pref * d[0], pref * d[1], sz};
   }
@@ -283,7 +282,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
   auto const n_modPsi = static_cast<int>(modPsi.size() >> 1);
   auto const rxy2 = d[0] * d[0] + d[1] * d[1];
   auto const rxy2_d = rxy2 * uz2;
-  auto const z_d = d[2] * uz;
+  auto const z_d = d[2] * box_geo.length_inv()[2];
   double E;
 
   if (rxy2 <= mmm1d_params.far_switch_radius_2) {
@@ -301,7 +300,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
 
       r2n *= rxy2_d;
     }
-    E *= uz;
+    E *= box_geo.length_inv()[2];
 
     /* real space parts */
 
@@ -319,7 +318,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
   } else {
     /* far range formula */
     auto const rxy = sqrt(rxy2);
-    auto const rxy_d = rxy * uz;
+    auto const rxy_d = rxy * box_geo.length_inv()[2];
     /* The first Bessel term will compensate a little bit the
        log term, so add them close together */
     E = -0.25 * log(rxy2_d) + 0.5 * (Utils::ln_2() - Utils::gamma());
@@ -330,7 +329,7 @@ double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
       auto const fq = c_2pi * bp;
       E += K0(fq * rxy_d) * cos(fq * z_d);
     }
-    E *= 4 * uz;
+    E *= 4 * box_geo.length_inv()[2];
   }
 
   return chpref * E;
diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
index d2fe0cef479..79b9cdee07d 100644
--- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
@@ -229,7 +229,7 @@ void dp3m_init() {
 void dp3m_set_tune_params(double r_cut, int mesh, int cao, double accuracy) {
   if (r_cut >= 0) {
     dp3m.params.r_cut = r_cut;
-    dp3m.params.r_cut_iL = r_cut / box_geo.length()[0];
+    dp3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0];
   }
 
   if (mesh >= 0)
@@ -268,7 +268,7 @@ void dp3m_set_params(double r_cut, int mesh, int cao, double alpha,
     Dipole::set_method_local(DIPOLAR_P3M);
 
   dp3m.params.r_cut = r_cut;
-  dp3m.params.r_cut_iL = r_cut / box_geo.length()[0];
+  dp3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0];
   dp3m.params.mesh[2] = dp3m.params.mesh[1] = dp3m.params.mesh[0] = mesh;
   dp3m.params.cao = cao;
   dp3m.params.alpha = alpha;
@@ -443,7 +443,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         }
       }
       node_k_space_energy_dip *=
-          dipole_prefac * Utils::pi() / box_geo.length()[0];
+          dipole_prefac * Utils::pi() * box_geo.length_inv()[0];
       boost::mpi::reduce(comm_cart, node_k_space_energy_dip, k_space_energy_dip,
                          std::plus<>(), 0);
 
@@ -454,7 +454,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         k_space_energy_dip -=
             dipole.prefactor *
             (dp3m.sum_mu2 * 2 *
-             pow(dp3m.params.alpha_L / box_geo.length()[0], 3) *
+             pow(dp3m.params.alpha_L * box_geo.length_inv()[0], 3) *
              Utils::sqrt_pi_i() / 3.0);
 
         auto const volume = box_geo.volume();
@@ -534,7 +534,7 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         /* Assign force component from mesh to particle */
         Utils::integral_parameter<AssignTorques, 1, 7>(
             dp3m.params.cao,
-            dipole_prefac * (2 * Utils::pi() / box_geo.length()[0]), d_rs,
+            dipole_prefac * (2 * Utils::pi() * box_geo.length_inv()[0]), d_rs,
             particles);
       }
 
@@ -623,8 +623,8 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         /* Assign force component from mesh to particle */
         Utils::integral_parameter<AssignForces, 1, 7>(
             dp3m.params.cao,
-            dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), d_rs,
-            particles);
+            dipole_prefac * pow(2 * Utils::pi() * box_geo.length_inv()[0], 2),
+            d_rs, particles);
       }
     } /* if (dp3m.sum_mu2 > 0) */
   }   /* if (force_flag) */
@@ -1143,8 +1143,8 @@ int dp3m_adaptive_tune(bool verbose) {
 
     r_cut_iL_min = 0;
     r_cut_iL_max = std::min(min_local_box_l, min_box_l / 2) - skin;
-    r_cut_iL_min *= 1. / box_geo.length()[0];
-    r_cut_iL_max *= 1. / box_geo.length()[0];
+    r_cut_iL_min *= box_geo.length_inv()[0];
+    r_cut_iL_max *= box_geo.length_inv()[0];
   } else {
     r_cut_iL_min = r_cut_iL_max = dp3m.params.r_cut_iL;
 
@@ -1382,7 +1382,7 @@ double dp3m_rtbisection(double box_size, double prefac, double r_cut_iL,
 
 void dp3m_init_a_ai_cao_cut() {
   for (int i = 0; i < 3; i++) {
-    dp3m.params.ai[i] = (double)dp3m.params.mesh[i] / box_geo.length()[i];
+    dp3m.params.ai[i] = dp3m.params.mesh[i] * box_geo.length_inv()[i];
     dp3m.params.a[i] = 1.0 / dp3m.params.ai[i];
     dp3m.params.cao_cut[i] = 0.5 * dp3m.params.a[i] * dp3m.params.cao;
   }
@@ -1394,7 +1394,7 @@ bool dp3m_sanity_checks_boxl() {
   bool ret = false;
   for (int i = 0; i < 3; i++) {
     /* check k-space cutoff */
-    if (dp3m.params.cao_cut[i] >= 0.5 * box_geo.length()[i]) {
+    if (dp3m.params.cao_cut[i] >= box_geo.length_half()[i]) {
       runtimeErrorMsg() << "dipolar P3M_init: k-space cutoff "
                         << dp3m.params.cao_cut[i]
                         << " is larger than half of box dimension "
@@ -1469,7 +1469,7 @@ void dp3m_scaleby_box_l() {
   }
 
   dp3m.params.r_cut = dp3m.params.r_cut_iL * box_geo.length()[0];
-  dp3m.params.alpha = dp3m.params.alpha_L / box_geo.length()[0];
+  dp3m.params.alpha = dp3m.params.alpha_L * box_geo.length_inv()[0];
   dp3m_init_a_ai_cao_cut();
   p3m_calc_lm_ld_pos(dp3m.local_mesh, dp3m.params);
   dp3m_sanity_checks_boxl();
diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp
index 58eb290e686..25758ec84c9 100644
--- a/src/core/electrostatics_magnetostatics/p3m.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m.cpp
@@ -216,7 +216,7 @@ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao,
                          double accuracy) {
   if (r_cut >= 0) {
     p3m.params.r_cut = r_cut;
-    p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]);
+    p3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0];
   } else if (r_cut == -1.0) {
     p3m.params.r_cut = 0;
     p3m.params.r_cut_iL = 0;
@@ -260,7 +260,7 @@ void p3m_set_params(double r_cut, const int *mesh, int cao, double alpha,
     coulomb.method = COULOMB_P3M;
 
   p3m.params.r_cut = r_cut;
-  p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]);
+  p3m.params.r_cut_iL = r_cut * box_geo.length_inv()[0];
   p3m.params.mesh[2] = mesh[2];
   p3m.params.mesh[1] = mesh[1];
   p3m.params.mesh[0] = mesh[0];
@@ -434,14 +434,14 @@ Utils::Vector9d p3m_calc_kspace_pressure_tensor() {
       for (j[1] = 0; j[1] < p3m.fft.plan[3].new_mesh[RY]; j[1]++) {
         for (j[2] = 0; j[2] < p3m.fft.plan[3].new_mesh[RZ]; j[2]++) {
           auto const kx = 2.0 * Utils::pi() *
-                          p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] /
-                          box_geo.length()[RX];
+                          p3m.d_op[RX][j[KX] + p3m.fft.plan[3].start[KX]] *
+                          box_geo.length_inv()[RX];
           auto const ky = 2.0 * Utils::pi() *
-                          p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] /
-                          box_geo.length()[RY];
+                          p3m.d_op[RY][j[KY] + p3m.fft.plan[3].start[KY]] *
+                          box_geo.length_inv()[RY];
           auto const kz = 2.0 * Utils::pi() *
-                          p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] /
-                          box_geo.length()[RZ];
+                          p3m.d_op[RZ][j[KZ] + p3m.fft.plan[3].start[KZ]] *
+                          box_geo.length_inv()[RZ];
           auto const sqk = Utils::sqr(kx) + Utils::sqr(ky) + Utils::sqr(kz);
 
           auto const node_k_space_energy =
@@ -509,8 +509,8 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
             int d_rs = (d + p3m.ks_pnum) % 3;
             /* directions */
             auto const k = 2.0 * Utils::pi() *
-                           p3m.d_op[d_rs][j[d] + p3m.fft.plan[3].start[d]] /
-                           box_geo.length()[d_rs];
+                           p3m.d_op[d_rs][j[d] + p3m.fft.plan[3].start[d]] *
+                           box_geo.length_inv()[d_rs];
 
             /* i*k*(Re+i*Im) = - Im*k + i*Re*k     (i=sqrt(-1)) */
             p3m.E_mesh[d_rs][2 * ind + 0] = -k * phi_hat.imag();
@@ -675,7 +675,7 @@ static double p3m_mcr_time(const int mesh[3], int cao, double r_cut_iL,
   p3m.params.mesh[2] = mesh[2];
   p3m.params.cao = cao;
   p3m.params.alpha_L = alpha_L;
-  p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]);
+  p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0];
 
   /* initialize p3m structures */
   mpi_bcast_coulomb_params();
@@ -1004,7 +1004,7 @@ int p3m_adaptive_tune(bool verbose) {
      */
   } else if (p3m.params.mesh[1] == -1 && p3m.params.mesh[2] == -1) {
     double mesh_density = mesh_density_min = mesh_density_max =
-        p3m.params.mesh[0] / box_geo.length()[0];
+        p3m.params.mesh[0] * box_geo.length_inv()[0];
     p3m.params.mesh[1] =
         static_cast<int>(std::round(mesh_density * box_geo.length()[1]));
     p3m.params.mesh[2] =
@@ -1020,7 +1020,7 @@ int p3m_adaptive_tune(bool verbose) {
     }
   } else {
     mesh_density_min = mesh_density_max =
-        p3m.params.mesh[0] / box_geo.length()[0];
+        p3m.params.mesh[0] * box_geo.length_inv()[0];
 
     if (verbose) {
       std::printf("fixed mesh %d %d %d\n", p3m.params.mesh[0],
@@ -1034,8 +1034,8 @@ int p3m_adaptive_tune(bool verbose) {
 
     r_cut_iL_min = 0;
     r_cut_iL_max = std::min(min_local_box_l, min_box_l / 2.0) - skin;
-    r_cut_iL_min *= (1. / box_geo.length()[0]);
-    r_cut_iL_max *= (1. / box_geo.length()[0]);
+    r_cut_iL_min *= box_geo.length_inv()[0];
+    r_cut_iL_max *= box_geo.length_inv()[0];
   } else {
     r_cut_iL_min = r_cut_iL_max = p3m.params.r_cut_iL;
 
@@ -1136,7 +1136,7 @@ int p3m_adaptive_tune(bool verbose) {
   p3m.params.mesh[2] = mesh[2];
   p3m.params.cao = cao;
   p3m.params.alpha_L = alpha_L;
-  p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]);
+  p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0];
   p3m.params.accuracy = accuracy;
   /* broadcast tuned p3m parameters */
   mpi_bcast_coulomb_params();
@@ -1244,7 +1244,7 @@ void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3],
 
 void p3m_init_a_ai_cao_cut() {
   for (int i = 0; i < 3; i++) {
-    p3m.params.ai[i] = (double)p3m.params.mesh[i] / box_geo.length()[i];
+    p3m.params.ai[i] = p3m.params.mesh[i] * box_geo.length_inv()[i];
     p3m.params.a[i] = 1.0 / p3m.params.ai[i];
     p3m.params.cao_cut[i] = 0.5 * p3m.params.a[i] * p3m.params.cao;
   }
@@ -1254,7 +1254,7 @@ bool p3m_sanity_checks_boxl() {
   bool ret = false;
   for (int i = 0; i < 3; i++) {
     /* check k-space cutoff */
-    if (p3m.params.cao_cut[i] >= 0.5 * box_geo.length()[i]) {
+    if (p3m.params.cao_cut[i] >= box_geo.length_half()[i]) {
       runtimeErrorMsg() << "P3M_init: k-space cutoff " << p3m.params.cao_cut[i]
                         << " is larger than half of box dimension "
                         << box_geo.length()[i];
@@ -1346,7 +1346,7 @@ void p3m_scaleby_box_l() {
   }
 
   p3m.params.r_cut = p3m.params.r_cut_iL * box_geo.length()[0];
-  p3m.params.alpha = p3m.params.alpha_L * (1. / box_geo.length()[0]);
+  p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0];
   p3m_init_a_ai_cao_cut();
   p3m_calc_lm_ld_pos(p3m.local_mesh, p3m.params);
   p3m_sanity_checks_boxl();
diff --git a/src/core/grid.cpp b/src/core/grid.cpp
index 843e909fc95..593fba1f562 100644
--- a/src/core/grid.cpp
+++ b/src/core/grid.cpp
@@ -111,8 +111,8 @@ void grid_changed_n_nodes() {
 }
 
 void rescale_boxl(int dir, double d_new) {
-  double scale =
-      (dir - 3) ? d_new / box_geo.length()[dir] : d_new / box_geo.length()[0];
+  double scale = (dir - 3) ? d_new * box_geo.length_inv()[dir]
+                           : d_new * box_geo.length_inv()[0];
 
   /* If shrinking, rescale the particles first. */
   if (scale <= 1.) {
diff --git a/src/core/integrators/velocity_verlet_npt.cpp b/src/core/integrators/velocity_verlet_npt.cpp
index 34790cb375c..f73b0ea123c 100644
--- a/src/core/integrators/velocity_verlet_npt.cpp
+++ b/src/core/integrators/velocity_verlet_npt.cpp
@@ -115,7 +115,7 @@ void velocity_verlet_npt_propagate_pos(const ParticleRange &particles) {
 
     L_new = pow(nptiso.volume, 1.0 / nptiso.dimension);
 
-    scal[1] = L_new / box_geo.length()[nptiso.non_const_dim];
+    scal[1] = L_new * box_geo.length_inv()[nptiso.non_const_dim];
     scal[0] = 1 / scal[1];
   }
   MPI_Bcast(scal, 3, MPI_DOUBLE, 0, comm_cart);
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index d4f8571e818..e34ac40b269 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -274,7 +274,7 @@ void calc_structurefactor(PartCfg &partCfg, std::vector<int> const &p_types,
 
   auto const order_sq = order * order;
   std::vector<double> ff(2 * order_sq + 1);
-  auto const twoPI_L = 2 * Utils::pi() / box_geo.length()[0];
+  auto const twoPI_L = 2 * Utils::pi() * box_geo.length_inv()[0];
 
   for (int i = 0; i <= order; i++) {
     for (int j = -order; j <= order; j++) {