diff --git a/src/core/actor/mmm-common_cuda.hpp b/src/core/actor/mmm-common_cuda.hpp
index 9559bfa888f..1692b33c1c8 100644
--- a/src/core/actor/mmm-common_cuda.hpp
+++ b/src/core/actor/mmm-common_cuda.hpp
@@ -39,14 +39,12 @@ __constant__ int device_linModPsi_offsets[2 * modpsi_order],
 __constant__ mmm1dgpu_real device_linModPsi[modpsi_constant_size];
 
 int modpsi_init() {
-  if (n_modPsi < modpsi_order) {
-    create_mod_psi_up_to(modpsi_order);
-  }
+  create_mod_psi_up_to(modpsi_order);
 
   // linearize the coefficients array
-  linModPsi_offsets.resize(2 * n_modPsi);
-  linModPsi_lengths.resize(2 * n_modPsi);
-  for (int i = 0; i < 2 * n_modPsi; i++) {
+  linModPsi_offsets.resize(modPsi.size());
+  linModPsi_lengths.resize(modPsi.size());
+  for (int i = 0; i < modPsi.size(); i++) {
     if (i == 0)
       linModPsi_offsets[i] = 0;
     else
@@ -54,9 +52,9 @@ int modpsi_init() {
           linModPsi_offsets[i - 1] + linModPsi_lengths[i - 1];
     linModPsi_lengths[i] = modPsi[i].size();
   }
-  linModPsi.resize(linModPsi_offsets[2 * n_modPsi - 1] +
-                   linModPsi_lengths[2 * n_modPsi - 1]);
-  for (int i = 0; i < 2 * n_modPsi; i++) {
+  linModPsi.resize(linModPsi_offsets[modPsi.size() - 1] +
+                   linModPsi_lengths[modPsi.size() - 1]);
+  for (int i = 0; i < modPsi.size(); i++) {
     for (int j = 0; j < modPsi[i].size(); j++) {
       linModPsi[linModPsi_offsets[i] + j] =
           (mmm1dgpu_real)modPsi[i][j]; // cast to single-precision if necessary
@@ -67,21 +65,22 @@ int modpsi_init() {
     cudaSetDevice(d);
 
     // copy to GPU
-    int linModPsiSize = linModPsi_offsets[2 * n_modPsi - 1] +
-                        linModPsi_lengths[2 * n_modPsi - 1];
+    int linModPsiSize = linModPsi_offsets[modPsi.size() - 1] +
+                        linModPsi_lengths[modPsi.size() - 1];
     if (linModPsiSize > modpsi_constant_size) {
       printf("ERROR: __constant__ device_linModPsi[] is not large enough\n");
       std::abort();
     }
     cuda_safe_mem(cudaMemcpyToSymbol(HIP_SYMBOL(device_linModPsi_offsets),
                                      linModPsi_offsets.data(),
-                                     2 * n_modPsi * sizeof(int)));
+                                     modPsi.size() * sizeof(int)));
     cuda_safe_mem(cudaMemcpyToSymbol(HIP_SYMBOL(device_linModPsi_lengths),
                                      linModPsi_lengths.data(),
-                                     2 * n_modPsi * sizeof(int)));
+                                     modPsi.size() * sizeof(int)));
     cuda_safe_mem(cudaMemcpyToSymbol(HIP_SYMBOL(device_linModPsi),
                                      linModPsi.data(),
                                      linModPsiSize * sizeof(mmm1dgpu_real)));
+    auto const n_modPsi = static_cast<int>(modPsi.size() >> 1);
     cuda_safe_mem(cudaMemcpyToSymbol(HIP_SYMBOL(device_n_modPsi), &n_modPsi,
                                      sizeof(int)));
   }
diff --git a/src/core/electrostatics_magnetostatics/mmm-modpsi.cpp b/src/core/electrostatics_magnetostatics/mmm-modpsi.cpp
index 6ab43db4967..fe3a83d089b 100644
--- a/src/core/electrostatics_magnetostatics/mmm-modpsi.cpp
+++ b/src/core/electrostatics_magnetostatics/mmm-modpsi.cpp
@@ -27,7 +27,6 @@
 #include <cmath>
 
 std::vector<std::vector<double>> modPsi;
-int n_modPsi = 0;
 
 static void preparePolygammaEven(int n, double binom,
                                  std::vector<double> &series) {
@@ -102,17 +101,15 @@ static void preparePolygammaOdd(int n, double binom,
 void create_mod_psi_up_to(int new_n) {
   int n;
   double binom;
-
-  if (new_n > n_modPsi) {
-    int old = n_modPsi;
-    n_modPsi = new_n;
-    modPsi.resize(2 * n_modPsi);
+  auto const old_n = static_cast<int>(modPsi.size() >> 1);
+  if (new_n > old_n) {
+    modPsi.resize(2 * new_n);
 
     binom = 1.0;
-    for (n = 0; n < old; n++)
+    for (n = 0; n < old_n; n++)
       binom *= (-0.5 - n) / (double)(n + 1);
 
-    for (; n < n_modPsi; n++) {
+    for (; n < new_n; n++) {
       preparePolygammaEven(n, binom, modPsi[2 * n]);
       preparePolygammaOdd(n, binom, modPsi[2 * n + 1]);
       binom *= (-0.5 - n) / (double)(n + 1);
diff --git a/src/core/electrostatics_magnetostatics/mmm-modpsi.hpp b/src/core/electrostatics_magnetostatics/mmm-modpsi.hpp
index 447b0edee1d..ee4780622ef 100644
--- a/src/core/electrostatics_magnetostatics/mmm-modpsi.hpp
+++ b/src/core/electrostatics_magnetostatics/mmm-modpsi.hpp
@@ -37,7 +37,6 @@
 
 /** Table of the Taylor expansions of the modified polygamma functions */
 extern std::vector<std::vector<double>> modPsi;
-extern int n_modPsi;
 
 /** Create both the even and odd polygamma functions up to order 2*n */
 void create_mod_psi_up_to(int n);
diff --git a/src/core/electrostatics_magnetostatics/mmm1d.cpp b/src/core/electrostatics_magnetostatics/mmm1d.cpp
index 73acef7b72f..0c316d50f80 100644
--- a/src/core/electrostatics_magnetostatics/mmm1d.cpp
+++ b/src/core/electrostatics_magnetostatics/mmm1d.cpp
@@ -179,6 +179,7 @@ void MMM1D_init() {
 
 void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d,
                                   double r, Utils::Vector3d &force) {
+  auto const n_modPsi = static_cast<int>(modPsi.size() >> 1);
   int dim;
   Utils::Vector3d F;
   double rxy2, rxy2_d, z_d;
@@ -276,6 +277,7 @@ void add_mmm1d_coulomb_pair_force(double chpref, Utils::Vector3d const &d,
 
 double mmm1d_coulomb_pair_energy(double const chpref, Utils::Vector3d const &d,
                                  double r2, double r) {
+  auto const n_modPsi = static_cast<int>(modPsi.size() >> 1);
   double rxy2, rxy2_d, z_d;
   double E;