diff --git a/maintainer/benchmarks/p3m.py b/maintainer/benchmarks/p3m.py
index 43fe781596d..7c75363f921 100644
--- a/maintainer/benchmarks/p3m.py
+++ b/maintainer/benchmarks/p3m.py
@@ -56,7 +56,7 @@
 
 
 import espressomd
-from espressomd import electrostatics
+from espressomd import electrostatics, minimize_energy
 if args.visualizer:
     from espressomd import visualization
 
@@ -140,10 +140,10 @@
 
 energy = system.analysis.energy()
 print("Before Minimization: E_total = {}".format(energy["total"]))
-system.minimize_energy.init(f_max=1000, gamma=30.0,
-                            max_steps=1000, max_displacement=0.05)
-system.minimize_energy.minimize()
-system.minimize_energy.minimize()
+minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0,
+                                 max_steps=1000, max_displacement=0.05)
+minimize_energy.steepest_descent(system, f_max=1000, gamma=30.0,
+                                 max_steps=1000, max_displacement=0.05)
 energy = system.analysis.energy()
 print("After Minimization: E_total = {}".format(energy["total"]))
 
diff --git a/src/core/electrostatics_magnetostatics/elc.cpp b/src/core/electrostatics_magnetostatics/elc.cpp
index a98126bbad6..4beb7f02c73 100644
--- a/src/core/electrostatics_magnetostatics/elc.cpp
+++ b/src/core/electrostatics_magnetostatics/elc.cpp
@@ -1301,7 +1301,7 @@ void assign_image_charge(const Particle &p) {
     auto const q_eff = elc_params.delta_mid_bot * p.p.q;
     auto const pos = Utils::Vector3d{p.r.p[0], p.r.p[1], -p.r.p[2]};
 
-    p3m_assign_charge(q_eff, pos, -1);
+    p3m_assign_charge(q_eff, pos);
   }
 
   if (p.r.p[2] > (elc_params.h - elc_params.space_layer)) {
@@ -1309,27 +1309,24 @@ void assign_image_charge(const Particle &p) {
     auto const pos =
         Utils::Vector3d{p.r.p[0], p.r.p[1], 2 * elc_params.h - p.r.p[2]};
 
-    p3m_assign_charge(q_eff, pos, -1);
+    p3m_assign_charge(q_eff, pos);
   }
 }
 } // namespace
 
 void ELC_p3m_charge_assign_both(const ParticleRange &particles) {
-  /* charged particle counter, charge fraction counter */
-  int cp_cnt = 0;
+  p3m.inter_weights.reset(p3m.params.cao);
+
   /* prepare local FFT mesh */
   for (int i = 0; i < p3m.local_mesh.size; i++)
     p3m.rs_mesh[i] = 0.0;
 
   for (auto &p : particles) {
     if (p.p.q != 0.0) {
-      p3m_assign_charge(p.p.q, p.r.p, cp_cnt);
+      p3m_assign_charge(p.p.q, p.r.p, p3m.inter_weights);
       assign_image_charge(p);
-
-      cp_cnt++;
     }
   }
-  p3m_shrink_wrap_charge_grid(cp_cnt);
 }
 
 void ELC_p3m_charge_assign_image(const ParticleRange &particles) {
diff --git a/src/core/electrostatics_magnetostatics/fft.cpp b/src/core/electrostatics_magnetostatics/fft.cpp
index 45a5620f958..7499281a655 100644
--- a/src/core/electrostatics_magnetostatics/fft.cpp
+++ b/src/core/electrostatics_magnetostatics/fft.cpp
@@ -489,7 +489,7 @@ void calc_2d_grid(int n, int grid[3]) {
 }
 } // namespace
 
-int fft_init(int const *ca_mesh_dim, int const *ca_mesh_margin,
+int fft_init(const Utils::Vector3i &ca_mesh_dim, int const *ca_mesh_margin,
              int *global_mesh_dim, double *global_mesh_off, int *ks_pnum,
              fft_data_struct &fft, const Utils::Vector3i &grid,
              const boost::mpi::communicator &comm) {
diff --git a/src/core/electrostatics_magnetostatics/fft.hpp b/src/core/electrostatics_magnetostatics/fft.hpp
index 8694974487d..4bf414768c8 100644
--- a/src/core/electrostatics_magnetostatics/fft.hpp
+++ b/src/core/electrostatics_magnetostatics/fft.hpp
@@ -186,7 +186,7 @@ struct fft_data_struct {
  *  \param comm            MPI communicator.
  *  \return Maximal size of local fft mesh (needed for allocation of ca_mesh).
  */
-int fft_init(int const *ca_mesh_dim, int const *ca_mesh_margin,
+int fft_init(const Utils::Vector3i &ca_mesh_dim, int const *ca_mesh_margin,
              int *global_mesh_dim, double *global_mesh_off, int *ks_pnum,
              fft_data_struct &fft, const Utils::Vector3i &grid,
              const boost::mpi::communicator &comm);
diff --git a/src/core/electrostatics_magnetostatics/p3m-common.hpp b/src/core/electrostatics_magnetostatics/p3m-common.hpp
index 3eddddf3400..87da9c0665f 100644
--- a/src/core/electrostatics_magnetostatics/p3m-common.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m-common.hpp
@@ -38,6 +38,8 @@
 
 #include "LocalBox.hpp"
 
+#include <utils/Vector.hpp>
+
 /** Error Codes for p3m tuning (version 2) */
 enum P3M_TUNE_ERROR {
   /** force evaluation failed */
@@ -71,7 +73,7 @@ enum P3M_TUNE_ERROR {
 typedef struct {
   /* local mesh characterization. */
   /** dimension (size) of local mesh. */
-  int dim[3];
+  Utils::Vector3i dim;
   /** number of local mesh points. */
   int size;
   /** index of lower left corner of the
@@ -112,8 +114,6 @@ typedef struct {
   double mesh_off[3] = {P3M_MESHOFF, P3M_MESHOFF, P3M_MESHOFF};
   /** charge assignment order ([0,7]). */
   int cao = 0;
-  /** number of interpolation points for charge assignment function */
-  int inter = P3M_N_INTERPOL;
   /** accuracy of the actual parameter set. */
   double accuracy = 0.0;
 
@@ -124,15 +124,13 @@ typedef struct {
   /** mesh constant. */
   double a[3] = {};
   /** inverse mesh constant. */
-  double ai[3] = {};
+  Utils::Vector3d ai = {};
   /** unscaled @ref P3MParameters::alpha_L "alpha_L" for use with fast
    *  inline functions only */
   double alpha = 0.0;
   /** unscaled @ref P3MParameters::r_cut_iL "r_cut_iL" for use with fast
    *  inline functions only */
   double r_cut = -1.;
-  /** full size of the interpolated assignment function */
-  int inter2 = 0;
   /** number of points unto which a single charge is interpolated, i.e.
    *  p3m.cao^3 */
   int cao3 = 0;
@@ -142,8 +140,8 @@ typedef struct {
 
   template <typename Archive> void serialize(Archive &ar, long int) {
     ar &tuning &alpha_L &r_cut_iL &mesh;
-    ar &mesh_off &cao &inter &accuracy &epsilon &cao_cut;
-    ar &a &ai &alpha &r_cut &inter2 &cao3 &additional_mesh;
+    ar &mesh_off &cao &accuracy &epsilon &cao_cut;
+    ar &a &ai &alpha &r_cut &cao3 &additional_mesh;
   }
 
 } P3MParameters;
diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
index 4d80c4d3dbb..2ac93b6afca 100644
--- a/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.cpp
@@ -51,7 +51,7 @@ using Utils::strcat_alloc;
 #include <utils/math/sinc.hpp>
 using Utils::sinc;
 #include <utils/constants.hpp>
-#include <utils/math/bspline.hpp>
+#include <utils/integral_parameter.hpp>
 #include <utils/math/sqr.hpp>
 
 #include <boost/range/algorithm/min_element.hpp>
@@ -82,20 +82,11 @@ dp3m_data_struct dp3m;
  */
 static void dp3m_init_a_ai_cao_cut();
 
-/** realloc charge assignment fields. */
-static void dp3m_realloc_ca_fields(int newsize);
-
 /** Checks for correctness for magnetic dipoles in P3M of the cao_cut,
  *  necessary when the box length changes
  */
 static bool dp3m_sanity_checks_boxl();
 
-/** Interpolate the P-th order charge assignment function from
- *  @cite hockney88a 5-189 (or 8-61). The following charge fractions
- *  are also tabulated in @cite deserno98a @cite deserno98b.
- */
-static void dp3m_interpolate_dipole_assignment_function();
-
 /** Shift the mesh points by mesh/2 */
 static void dp3m_calc_meshift();
 
@@ -172,72 +163,6 @@ static double dp3m_average_dipolar_self_energy(double box_l, int mesh);
 static double dp3m_perform_aliasing_sums_dipolar_self_energy(const int n[3]);
 
 /************************************************************/
-/* functions related to the correction of the dipolar p3m-energy */
-
-/*
-
-// Do the sum over k<>0 where k={kx,ky,kz} with kx,ky,kz INTEGERS, of
-// exp(-Utils::pi()**2*k**2/alpha**2/L**2)
-static double dp3m_sumi1(double alpha_L){
-       int k2,kx,ky,kz,kx2,ky2,limit=60;
-       double suma,alpha_L2;
-
-       alpha_L2= alpha_L* alpha_L;
-
-       //fprintf(stderr,"alpha_L=%le\n",alpha_L);
-       //fprintf(stderr,"Utils::pi()=%le\n",Utils::pi());
-
-
-       suma=0.0;
-       for(kx=-limit;kx<=limit;kx++){
-         kx2=kx*kx;
-       for(ky=-limit;ky<=limit;ky++){
-         ky2=ky*ky;
-       for(kz=-limit;kz<=limit;kz++){
-           k2=kx2+ky2+kz*kz;
-           suma+=exp(-Utils::pi()*Utils::pi()*k2/(alpha_L*alpha_L));
-       }}}
-       suma-=1; //It's easier to subtract the term k=0 later than put an if
-inside the loops
-
-
-         //fprintf(stderr,"suma=%le\n",suma);
-
-
-   return suma;
-}
-
-*/
-
-/************************************************************/
-
-/*
-// Do the sum over n<>0 where n={nx*L,ny*L,nz*L} with nx,ny,nz INTEGERS, of
-// exp(-alpha_iL**2*n**2)
-static double dp3m_sumi2(double alpha_L){
-       int n2,nx,ny,nz,nx2,ny2,limit=60;
-       double suma;
-
-
-
-       suma=0.0;
-       for(nx=-limit;nx<=limit;nx++){
-         nx2=nx*nx;
-       for(ny=-limit;ny<=limit;ny++){
-         ny2=ny*ny;
-       for(nz=-limit;nz<=limit;nz++){
-           n2=nx2+ny2+nz*nz;
-           suma+=exp(-alpha_L*alpha_L*n2);
-       }}}
-       suma-=1; //It's easier to subtract the term n=0 later than put an if
-inside the loops
-
-
-
-   return suma;
-}
-
-*/
 
 dp3m_data_struct::dp3m_data_struct() {
   params.epsilon = P3M_EPSILON_MAGNETIC;
@@ -248,8 +173,6 @@ dp3m_data_struct::dp3m_data_struct() {
   sum_mu2 = 0.0;
 
   pos_shift = 0.0;
-
-  ca_num = 0;
   ks_pnum = 0;
 
   energy_correction = 0.0;
@@ -267,8 +190,6 @@ void dp3m_deactivate() {
 }
 
 void dp3m_init() {
-  int n;
-
   if (dipole.prefactor <= 0.0) {
     // dipolar prefactor is zero: magnetostatics switched off
     dp3m.params.r_cut = 0.0;
@@ -283,14 +204,6 @@ void dp3m_init() {
      * and the cutoff for charge assignment dp3m.params.cao_cut */
     dp3m_init_a_ai_cao_cut();
 
-    /* initialize ca fields to size CA_INCREMENT: dp3m.ca_frac and dp3m.ca_fmp
-     */
-    dp3m.ca_num = 0;
-    if (dp3m.ca_num < CA_INCREMENT) {
-      dp3m.ca_num = 0;
-      dp3m_realloc_ca_fields(CA_INCREMENT);
-    }
-
     p3m_calc_local_ca_mesh(dp3m.local_mesh, dp3m.params, local_geo, skin);
 
     dp3m.sm.resize(comm_cart, dp3m.local_mesh);
@@ -298,12 +211,6 @@ void dp3m_init() {
     /* fix box length dependent constants */
     dp3m_scaleby_box_l();
 
-    if (dp3m.params.inter > 0)
-      dp3m_interpolate_dipole_assignment_function();
-
-    dp3m.pos_shift =
-        std::floor((dp3m.params.cao - 1) / 2.0) - (dp3m.params.cao % 2) / 2.0;
-
     int ca_mesh_size = fft_init(dp3m.local_mesh.dim, dp3m.local_mesh.margin,
                                 dp3m.params.mesh, dp3m.params.mesh_off,
                                 &dp3m.ks_pnum, dp3m.fft, node_grid, comm_cart);
@@ -419,9 +326,6 @@ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha,
 
   if (accuracy >= 0)
     dp3m.params.accuracy = accuracy;
-
-  if (n_interpol != -1)
-    dp3m.params.inter = n_interpol;
 }
 
 /*****************************************************************************/
@@ -486,43 +390,25 @@ int dp3m_set_eps(double eps) {
   return ES_OK;
 }
 
-int dp3m_set_ninterpol(int n) {
-  if (n < 0)
-    return ES_ERROR;
-
-  dp3m.params.inter = n;
-
-  mpi_bcast_coulomb_params();
-
-  return ES_OK;
-}
-
-/*****************************************************************************/
-
-void dp3m_interpolate_dipole_assignment_function() {
-  double dInterpol = 0.5 / (double)dp3m.params.inter;
-  int i;
-  long j;
-
-  if (dp3m.params.inter == 0)
-    return;
-
-  dp3m.params.inter2 = 2 * dp3m.params.inter + 1;
-
-  for (i = 0; i < dp3m.params.cao; i++) {
-    /* allocate memory for interpolation array */
-    dp3m.int_caf[i].resize(2 * dp3m.params.inter + 1);
-
-    /* loop over all interpolation points */
-    for (j = -dp3m.params.inter; j <= dp3m.params.inter; j++)
-      dp3m.int_caf[i][j + dp3m.params.inter] =
-          Utils::bspline(i, j * dInterpol, dp3m.params.cao);
+namespace {
+template <size_t cao> struct AssignDipole {
+  void operator()(Utils::Vector3d const &real_pos,
+                  Utils::Vector3d const &dip) const {
+    auto const weights = p3m_calculate_interpolation_weights<cao>(
+        real_pos, dp3m.params.ai, dp3m.local_mesh);
+    p3m_interpolate<cao>(dp3m.local_mesh, weights, [&dip](int ind, double w) {
+      dp3m.rs_mesh_dip[0][ind] += w * dip[0];
+      dp3m.rs_mesh_dip[1][ind] += w * dip[1];
+      dp3m.rs_mesh_dip[2][ind] += w * dip[2];
+    });
+
+    dp3m.inter_weights.store<cao>(weights);
   }
-}
+};
+} // namespace
 
 void dp3m_dipole_assign(const ParticleRange &particles) {
-  /* magnetic particle counter, dipole fraction counter */
-  int cp_cnt = 0;
+  dp3m.inter_weights.reset(dp3m.params.cao);
 
   /* prepare local FFT mesh */
   for (auto &i : dp3m.rs_mesh_dip)
@@ -531,210 +417,52 @@ void dp3m_dipole_assign(const ParticleRange &particles) {
 
   for (auto const &p : particles) {
     if (p.p.dipm != 0.0) {
-      dp3m_assign_dipole(p.r.p.data(), p.p.dipm, p.calc_dip().data(), cp_cnt);
-      cp_cnt++;
+      Utils::integral_parameter<AssignDipole, 1, 7>(dp3m.params.cao, p.r.p,
+                                                    p.calc_dip());
     }
   }
-
-  dp3m_shrink_wrap_dipole_grid(cp_cnt);
 }
 
-void dp3m_assign_dipole(double const real_pos[3], double mu,
-                        double const dip[3], int cp_cnt) {
-  /* we do not really want to export these, but this function should be inlined
-   */
-  int d, i0, i1, i2;
-  double tmp0, tmp1;
-  /* position of a particle in local mesh units */
-  double pos;
-  /* 1d-index of nearest mesh point */
-  int nmp;
-  /* distance to nearest mesh point */
-  double dist[3];
-  /* index for caf interpolation grid */
-  int arg[3];
-  /* index, index jumps for dp3m.rs_mesh array */
-  int q_ind = 0;
-  double cur_ca_frac_val, *cur_ca_frac;
-
-  // make sure we have enough space
-  if (cp_cnt >= dp3m.ca_num)
-    dp3m_realloc_ca_fields(cp_cnt + 1);
-  // do it here, since p3m_realloc_ca_fields may change the address of
-  // dp3m.ca_frac
-  cur_ca_frac = dp3m.ca_frac.data() + dp3m.params.cao3 * cp_cnt;
-
-  if (dp3m.params.inter == 0) {
-    for (d = 0; d < 3; d++) {
-      /* particle position in mesh coordinates */
-      pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) -
-            dp3m.pos_shift;
-      /* nearest mesh point */
-      nmp = (int)pos;
-      /* distance to nearest mesh point */
-      dist[d] = (pos - nmp) - 0.5;
-      /* 3d-array index of nearest mesh point */
-      q_ind = (d == 0) ? nmp : nmp + dp3m.local_mesh.dim[d] * q_ind;
-    }
+namespace {
+template <size_t cao> struct AssignTorques {
+  void operator()(double prefac, int d_rs,
+                  const ParticleRange &particles) const {
+    /* particle counter */
+    int cp_cnt = 0;
+    for (auto &p : particles) {
+      auto const w = dp3m.inter_weights.load<cao>(cp_cnt++);
 
-    if (cp_cnt >= 0)
-      dp3m.ca_fmp[cp_cnt] = q_ind;
-
-    for (i0 = 0; i0 < dp3m.params.cao; i0++) {
-      tmp0 = Utils::bspline(i0, dist[0], dp3m.params.cao);
-      for (i1 = 0; i1 < dp3m.params.cao; i1++) {
-        tmp1 = tmp0 * Utils::bspline(i1, dist[1], dp3m.params.cao);
-        for (i2 = 0; i2 < dp3m.params.cao; i2++) {
-          cur_ca_frac_val = tmp1 * Utils::bspline(i2, dist[2], dp3m.params.cao);
-          if (cp_cnt >= 0)
-            *(cur_ca_frac++) = cur_ca_frac_val;
-          if (mu != 0.0) {
-            dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val;
-            dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val;
-            dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val;
-          }
-          q_ind++;
-        }
-        q_ind += dp3m.local_mesh.q_2_off;
-      }
-      q_ind += dp3m.local_mesh.q_21_off;
-    }
-  } else {
-    /* particle position in mesh coordinates */
-    for (d = 0; d < 3; d++) {
-      pos = ((real_pos[d] - dp3m.local_mesh.ld_pos[d]) * dp3m.params.ai[d]) -
-            dp3m.pos_shift;
-      nmp = (int)pos;
-      arg[d] = (int)((pos - nmp) * dp3m.params.inter2);
-      /* for the first dimension, q_ind is always zero, so this shifts correctly
-       */
-      q_ind = nmp + dp3m.local_mesh.dim[d] * q_ind;
-    }
-    if (cp_cnt >= 0)
-      dp3m.ca_fmp[cp_cnt] = q_ind;
-
-    for (i0 = 0; i0 < dp3m.params.cao; i0++) {
-      tmp0 = dp3m.int_caf[i0][arg[0]];
-      for (i1 = 0; i1 < dp3m.params.cao; i1++) {
-        tmp1 = tmp0 * dp3m.int_caf[i1][arg[1]];
-        for (i2 = 0; i2 < dp3m.params.cao; i2++) {
-          cur_ca_frac_val = tmp1 * dp3m.int_caf[i2][arg[2]];
-          if (cp_cnt >= 0)
-            *(cur_ca_frac++) = cur_ca_frac_val;
-          if (mu != 0.0) {
-            dp3m.rs_mesh_dip[0][q_ind] += dip[0] * cur_ca_frac_val;
-            dp3m.rs_mesh_dip[1][q_ind] += dip[1] * cur_ca_frac_val;
-            dp3m.rs_mesh_dip[2][q_ind] += dip[2] * cur_ca_frac_val;
-          }
-          q_ind++;
-        }
-        q_ind += dp3m.local_mesh.q_2_off;
-      }
-      q_ind += dp3m.local_mesh.q_21_off;
+      Utils::Vector3d E{};
+      p3m_interpolate(dp3m.local_mesh, w, [&E, d_rs](int ind, double w) {
+        E[d_rs] += w * dp3m.rs_mesh[ind];
+      });
+
+      p.f.torque -= vector_product(p.calc_dip(), prefac * E);
     }
   }
-}
+};
 
-void dp3m_shrink_wrap_dipole_grid(int n_dipoles) {
-  if (n_dipoles < dp3m.ca_num)
-    dp3m_realloc_ca_fields(n_dipoles);
-}
+template <size_t cao> struct AssignForces {
+  void operator()(double prefac, int d_rs,
+                  const ParticleRange &particles) const {
+    /* particle counter */
+    int cp_cnt = 0;
+    for (auto &p : particles) {
+      auto const w = dp3m.inter_weights.load<cao>(cp_cnt++);
 
-#ifdef ROTATION
-/** Assign the torques obtained from k-space */
-static void P3M_assign_torques(double prefac, int d_rs,
-                               const ParticleRange &particles) {
-  /* particle counter, charge fraction counter */
-  int cp_cnt = 0, cf_cnt = 0;
-  /* index, index jumps for dp3m.rs_mesh array */
-  int q_ind;
-  auto const q_m_off = (dp3m.local_mesh.dim[2] - dp3m.params.cao);
-  auto const q_s_off =
-      dp3m.local_mesh.dim[2] * (dp3m.local_mesh.dim[1] - dp3m.params.cao);
-
-  for (auto &p : particles) {
-    if ((p.p.dipm) != 0.0) {
-      const Utils::Vector3d dip = p.calc_dip();
-      q_ind = dp3m.ca_fmp[cp_cnt];
-      for (int i0 = 0; i0 < dp3m.params.cao; i0++) {
-        for (int i1 = 0; i1 < dp3m.params.cao; i1++) {
-          for (int i2 = 0; i2 < dp3m.params.cao; i2++) {
-            /*
-            The following line would fill the torque with the k-space electric
-            field (without the self-field term) [notice the minus sign!]:
-                                p.f.torque[d_rs] -=
-                prefac*dp3m.ca_frac[cf_cnt]*dp3m.rs_mesh[q_ind];
-            Since the torque is the dipole moment cross-product with E, we
-            have:
-            */
-            switch (d_rs) {
-            case 0: // E_x
-              p.f.torque[1] -=
-                  dip[2] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-              p.f.torque[2] +=
-                  dip[1] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-              break;
-            case 1: // E_y
-              p.f.torque[0] +=
-                  dip[2] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-              p.f.torque[2] -=
-                  dip[0] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-              break;
-            case 2: // E_z
-              p.f.torque[0] -=
-                  dip[1] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-              p.f.torque[1] +=
-                  dip[0] * prefac * dp3m.ca_frac[cf_cnt] * dp3m.rs_mesh[q_ind];
-            }
-            q_ind++;
-            cf_cnt++;
-          }
-          q_ind += q_m_off;
-        }
-        q_ind += q_s_off;
-      }
-      cp_cnt++;
-    }
-  }
-}
-#endif
+      Utils::Vector3d E{};
 
-/** Assign the dipolar forces obtained from k-space */
-static void dp3m_assign_forces_dip(double prefac, int d_rs,
-                                   const ParticleRange &particles) {
-  /* particle counter, charge fraction counter */
-  int cp_cnt = 0, cf_cnt = 0;
-  /* index, index jumps for dp3m.rs_mesh array */
-  int q_ind;
-  int q_m_off = (dp3m.local_mesh.dim[2] - dp3m.params.cao);
-  int q_s_off =
-      dp3m.local_mesh.dim[2] * (dp3m.local_mesh.dim[1] - dp3m.params.cao);
-
-  cp_cnt = 0;
-  cf_cnt = 0;
-
-  for (auto &p : particles) {
-    if ((p.p.dipm) != 0.0) {
-      const Utils::Vector3d dip = p.calc_dip();
-      q_ind = dp3m.ca_fmp[cp_cnt];
-      for (int i0 = 0; i0 < dp3m.params.cao; i0++) {
-        for (int i1 = 0; i1 < dp3m.params.cao; i1++) {
-          for (int i2 = 0; i2 < dp3m.params.cao; i2++) {
-            p.f.f[d_rs] += prefac * dp3m.ca_frac[cf_cnt] *
-                           (dp3m.rs_mesh_dip[0][q_ind] * dip[0] +
-                            dp3m.rs_mesh_dip[1][q_ind] * dip[1] +
-                            dp3m.rs_mesh_dip[2][q_ind] * dip[2]);
-            q_ind++;
-            cf_cnt++;
-          }
-          q_ind += q_m_off;
-        }
-        q_ind += q_s_off;
-      }
-      cp_cnt++;
+      p3m_interpolate(dp3m.local_mesh, w, [&E](int ind, double w) {
+        E[0] += w * dp3m.rs_mesh_dip[0][ind];
+        E[1] += w * dp3m.rs_mesh_dip[1][ind];
+        E[2] += w * dp3m.rs_mesh_dip[2][ind];
+      });
+
+      p.f.f[d_rs] += p.calc_dip() * prefac * E;
     }
   }
-}
+};
+} // namespace
 
 /*****************************************************************************/
 
@@ -890,9 +618,10 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         dp3m.sm.spread_grid(dp3m.rs_mesh.data(), comm_cart,
                             dp3m.local_mesh.dim);
         /* Assign force component from mesh to particle */
-        P3M_assign_torques(dipole_prefac *
-                               (2 * Utils::pi() / box_geo.length()[0]),
-                           d_rs, particles);
+        Utils::integral_parameter<AssignTorques, 1, 7>(
+            dp3m.params.cao,
+            dipole_prefac * (2 * Utils::pi() / box_geo.length()[0]), d_rs,
+            particles);
       }
 #endif /*ifdef ROTATION */
 
@@ -979,7 +708,8 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
         dp3m.sm.spread_grid(Utils::make_span(meshes), comm_cart,
                             dp3m.local_mesh.dim);
         /* Assign force component from mesh to particle */
-        dp3m_assign_forces_dip(
+        Utils::integral_parameter<AssignForces, 1, 7>(
+            dp3m.params.cao,
             dipole_prefac * pow(2 * Utils::pi() / box_geo.length()[0], 2), d_rs,
             particles);
       }
@@ -1074,20 +804,6 @@ double calc_surface_term(bool force_flag, bool energy_flag,
 
 /*****************************************************************************/
 
-void dp3m_realloc_ca_fields(int newsize) {
-  newsize = ((newsize + CA_INCREMENT - 1) / CA_INCREMENT) * CA_INCREMENT;
-  if (newsize == dp3m.ca_num)
-    return;
-  if (newsize < CA_INCREMENT)
-    newsize = CA_INCREMENT;
-
-  dp3m.ca_num = newsize;
-  dp3m.ca_frac.resize(dp3m.params.cao3 * dp3m.ca_num);
-  dp3m.ca_fmp.resize(dp3m.ca_num);
-}
-
-/*****************************************************************************/
-
 void dp3m_calc_meshift() {
   int i;
   double dmesh;
@@ -1987,11 +1703,7 @@ bool dp3m_sanity_checks(const Utils::Vector3i &grid) {
     runtimeErrorMsg() << "dipolar P3M requires periodicity 1 1 1";
     ret = true;
   }
-  /*
-  if (n_nodes != 1) {
-      runtimeErrorMsg() <<"dipolar P3M does not run in parallel";
-    ret = true;
-  } */
+
   if (cell_structure.type != CELL_STRUCTURE_DOMDEC) {
     runtimeErrorMsg() << "dipolar P3M at present requires the domain "
                          "decomposition cell system";
diff --git a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
index 4cffb17a23a..6e7b565eb74 100644
--- a/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m-dipolar.hpp
@@ -40,6 +40,7 @@
 #include "electrostatics_magnetostatics/dipole.hpp"
 #include "fft.hpp"
 #include "p3m-common.hpp"
+#include "p3m_interpolation.hpp"
 #include "p3m_send_mesh.hpp"
 
 #include <ParticleRange.hpp>
@@ -65,9 +66,6 @@ struct dp3m_data_struct {
   /** Sum of square of magnetic dipoles (only on master node). */
   double sum_mu2;
 
-  /** interpolation of the charge assignment function. */
-  std::vector<std::vector<double>> int_caf{7};
-
   /** position shift for calc. of first assignment mesh point. */
   double pos_shift;
   /** help variable for calculation of aliasing sums */
@@ -81,13 +79,8 @@ struct dp3m_data_struct {
   /** Energy optimised influence function (k-space) */
   std::vector<double> g_energy;
 
-  /** number of charged particles on the node. */
-  int ca_num;
+  p3m_interpolation_cache inter_weights;
 
-  /** Charge fractions for mesh assignment. */
-  std::vector<double> ca_frac;
-  /** index of first mesh point for charge assignment. */
-  std::vector<int> ca_fmp;
   /** number of permutations in k_space */
   int ks_pnum;
 
@@ -115,9 +108,6 @@ void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha,
 int dp3m_set_params(double r_cut, int mesh, int cao, double alpha,
                     double accuracy);
 
-/** @copydoc p3m_set_ninterpol */
-int dp3m_set_ninterpol(int n);
-
 /** @copydoc p3m_set_mesh_offset */
 int dp3m_set_mesh_offset(double x, double y, double z);
 
@@ -195,22 +185,6 @@ double dp3m_calc_kspace_forces(bool force_flag, bool energy_flag,
  */
 void dp3m_count_magnetic_particles();
 
-/** Assign a single dipole into the current dipole grid.
- *
- *  @param[in] real_pos   %Particle position in real space
- *  @param[in] mu         %Particle magnetic dipole magnitude
- *  @param[in] dip        %Particle magnetic dipole vector
- *  @param[in] cp_cnt     The running index, which may be smaller than 0, in
- *                        which case the dipole is assumed to be virtual and
- *                        is not stored in the @ref dp3m_data_struct::ca_frac
- *                        "ca_frac" arrays
- */
-void dp3m_assign_dipole(double const real_pos[3], double mu,
-                        double const dip[3], int cp_cnt);
-
-/** Shrink wrap the dipoles grid */
-void dp3m_shrink_wrap_dipole_grid(int n_dipoles);
-
 /** Calculate real space contribution of p3m dipolar pair forces and torques.
  *  If NPT is compiled in, it returns the energy, which is needed for NPT.
  */
diff --git a/src/core/electrostatics_magnetostatics/p3m.cpp b/src/core/electrostatics_magnetostatics/p3m.cpp
index 4f028eae22d..45fec45185d 100644
--- a/src/core/electrostatics_magnetostatics/p3m.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m.cpp
@@ -50,7 +50,7 @@ using Utils::sinc;
 #include <utils/strcat_alloc.hpp>
 using Utils::strcat_alloc;
 #include <utils/constants.hpp>
-#include <utils/math/bspline.hpp>
+#include <utils/integral_parameter.hpp>
 #include <utils/math/sqr.hpp>
 
 #include <boost/optional.hpp>
@@ -95,9 +95,6 @@ p3m_data_struct p3m;
  */
 static void p3m_init_a_ai_cao_cut();
 
-/** Realloc charge assignment fields. */
-static void p3m_realloc_ca_fields(int newsize);
-
 static bool p3m_sanity_checks_system(const Utils::Vector3i &grid);
 
 /** Checks for correctness for charges in P3M of the cao_cut,
@@ -105,12 +102,6 @@ static bool p3m_sanity_checks_system(const Utils::Vector3i &grid);
  */
 static bool p3m_sanity_checks_boxl();
 
-/** Interpolate the P-th order charge assignment function from
- *  @cite hockney88a 5-189 (or 8-61). The following charge fractions
- *  are also tabulated in @cite deserno98a @cite deserno98b.
- */
-static void p3m_interpolate_charge_assignment_function();
-
 /** Shift the mesh points by mesh/2 */
 static void p3m_calc_meshift();
 
@@ -137,19 +128,6 @@ static void p3m_calc_influence_function_force();
  */
 static void p3m_calc_influence_function_energy();
 
-/** Calculate the aliasing sums for the optimal influence function.
- *
- *  Calculate the aliasing sums in the nominator and denominator of
- *  the expression for the optimal influence function (see
- *  @cite hockney88a : 8-22, p. 275).
- *
- *  \param  n           n-vector for which the aliasing sum is to be performed.
- *  \param  nominator   aliasing sums in the nominator.
- *  \retval denominator aliasing sum in the denominator
- */
-double p3m_perform_aliasing_sums_force(int n[3], double nominator[3]);
-double p3m_perform_aliasing_sums_energy(int n[3]);
-
 /*@}*/
 
 /** @name P3M tuning helper functions */
@@ -188,17 +166,6 @@ static void p3m_tune_aliasing_sums(int nx, int ny, int nz, const int mesh[3],
                                    double alpha_L_i, double *alias1,
                                    double *alias2);
 
-/** Template parameterized calculation of the charge assignment to be called by
- *  wrapper.
- *  \tparam cao      charge assignment order.
- */
-template <int cao>
-static void p3m_do_charge_assign(const ParticleRange &particles);
-
-template <int cao>
-void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos,
-                          int cp_cnt);
-
 p3m_data_struct::p3m_data_struct() {
   /* local_mesh is uninitialized */
   /* sm is uninitialized */
@@ -207,9 +174,6 @@ p3m_data_struct::p3m_data_struct() {
   sum_q2 = 0.0;
   square_sum_q = 0.0;
 
-  pos_shift = 0.0;
-
-  ca_num = 0;
   ks_pnum = 0;
 }
 
@@ -228,10 +192,6 @@ void p3m_init() {
      * the cutoff for charge assignment p3m.params.cao_cut */
     p3m_init_a_ai_cao_cut();
 
-    /* initialize ca fields to size CA_INCREMENT: p3m.ca_frac and p3m.ca_fmp */
-    p3m.ca_num = 0;
-    p3m_realloc_ca_fields(CA_INCREMENT);
-
     p3m_calc_local_ca_mesh(p3m.local_mesh, p3m.params, local_geo, skin);
 
     p3m.sm.resize(comm_cart, p3m.local_mesh);
@@ -250,19 +210,12 @@ void p3m_init() {
     /* fix box length dependent constants */
     p3m_scaleby_box_l();
 
-    if (p3m.params.inter > 0)
-      p3m_interpolate_charge_assignment_function();
-
-    /* position offset for calc. of first meshpoint */
-    p3m.pos_shift =
-        std::floor((p3m.params.cao - 1) / 2.0) - (p3m.params.cao % 2) / 2.0;
-
     p3m_count_charged_particles();
   }
 }
 
 void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha,
-                         double accuracy, int n_interpol) {
+                         double accuracy) {
   if (r_cut >= 0) {
     p3m.params.r_cut = r_cut;
     p3m.params.r_cut_iL = r_cut * (1. / box_geo.length()[0]);
@@ -284,9 +237,6 @@ void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha,
 
   if (accuracy >= 0)
     p3m.params.accuracy = accuracy;
-
-  if (n_interpol != -1)
-    p3m.params.inter = n_interpol;
 }
 
 /*@}*/
@@ -350,239 +300,91 @@ int p3m_set_eps(double eps) {
   return ES_OK;
 }
 
-int p3m_set_ninterpol(int n) {
-  if (n < 0)
-    return ES_ERROR;
-
-  p3m.params.inter = n;
-
-  mpi_bcast_coulomb_params();
-
-  return ES_OK;
-}
-
-void p3m_interpolate_charge_assignment_function() {
-  double dInterpol = 0.5 / (double)p3m.params.inter;
-  int i;
-  long j;
+namespace {
+template <size_t cao> struct AssignCharge {
+  void operator()(double q, const Utils::Vector3d &real_pos,
+                  const Utils::Vector3d &ai, p3m_local_mesh const &local_mesh,
+                  p3m_interpolation_cache &inter_weights) {
+    auto const w =
+        p3m_calculate_interpolation_weights<cao>(real_pos, ai, local_mesh);
 
-  if (p3m.params.inter == 0)
-    return;
+    inter_weights.store(w);
 
-  p3m.params.inter2 = 2 * p3m.params.inter + 1;
+    p3m_interpolate(local_mesh, w,
+                    [q](int ind, double w) { p3m.rs_mesh[ind] += w * q; });
+  }
 
-  for (i = 0; i < p3m.params.cao; i++) {
-    /* allocate memory for interpolation array */
-    p3m.int_caf[i].resize(2 * p3m.params.inter + 1);
+  void operator()(double q, const Utils::Vector3d &real_pos,
+                  const Utils::Vector3d &ai, p3m_local_mesh const &local_mesh) {
+    p3m_interpolate(
+        local_mesh,
+        p3m_calculate_interpolation_weights<cao>(real_pos, ai, local_mesh),
+        [q](int ind, double w) { p3m.rs_mesh[ind] += w * q; });
+  }
 
-    /* loop over all interpolation points */
-    for (j = -p3m.params.inter; j <= p3m.params.inter; j++)
-      p3m.int_caf[i][j + p3m.params.inter] =
-          Utils::bspline(i, j * dInterpol, p3m.params.cao);
+  void operator()(const ParticleRange &particles) {
+    for (auto &p : particles) {
+      if (p.p.q != 0.0) {
+        this->operator()(p.p.q, p.r.p, p3m.params.ai, p3m.local_mesh,
+                         p3m.inter_weights);
+      }
+    }
   }
-}
+};
+} // namespace
 
-/* Template wrapper for p3m_do_charge_assign() */
 void p3m_charge_assign(const ParticleRange &particles) {
-  switch (p3m.params.cao) {
-  case 1:
-    p3m_do_charge_assign<1>(particles);
-    break;
-  case 2:
-    p3m_do_charge_assign<2>(particles);
-    break;
-  case 3:
-    p3m_do_charge_assign<3>(particles);
-    break;
-  case 4:
-    p3m_do_charge_assign<4>(particles);
-    break;
-  case 5:
-    p3m_do_charge_assign<5>(particles);
-    break;
-  case 6:
-    p3m_do_charge_assign<6>(particles);
-    break;
-  case 7:
-    p3m_do_charge_assign<7>(particles);
-    break;
-  }
-}
+  p3m.inter_weights.reset(p3m.params.cao);
 
-/** Assign the charges */
-template <int cao> void p3m_do_charge_assign(const ParticleRange &particles) {
-  /* charged particle counter, charge fraction counter */
-  int cp_cnt = 0;
   /* prepare local FFT mesh */
   for (int i = 0; i < p3m.local_mesh.size; i++)
     p3m.rs_mesh[i] = 0.0;
 
-  for (auto &p : particles) {
-    if (p.p.q != 0.0) {
-      p3m_do_assign_charge<cao>(p.p.q, p.r.p, cp_cnt);
-      cp_cnt++;
-    }
-  }
-
-  p3m_shrink_wrap_charge_grid(cp_cnt);
-}
-
-/* Template wrapper for p3m_do_assign_charge() */
-void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt) {
-  switch (p3m.params.cao) {
-  case 1:
-    p3m_do_assign_charge<1>(q, real_pos, cp_cnt);
-    break;
-  case 2:
-    p3m_do_assign_charge<2>(q, real_pos, cp_cnt);
-    break;
-  case 3:
-    p3m_do_assign_charge<3>(q, real_pos, cp_cnt);
-    break;
-  case 4:
-    p3m_do_assign_charge<4>(q, real_pos, cp_cnt);
-    break;
-  case 5:
-    p3m_do_assign_charge<5>(q, real_pos, cp_cnt);
-    break;
-  case 6:
-    p3m_do_assign_charge<6>(q, real_pos, cp_cnt);
-    break;
-  case 7:
-    p3m_do_assign_charge<7>(q, real_pos, cp_cnt);
-    break;
-  }
+  Utils::integral_parameter<AssignCharge, 1, 7>(p3m.params.cao, particles);
 }
 
-template <int cao>
-void p3m_do_assign_charge(double q, const Utils::Vector3d &real_pos,
-                          int cp_cnt) {
-  auto const inter = not(p3m.params.inter == 0);
-  /* distance to nearest mesh point */
-  double dist[3];
-  /* index for caf interpolation grid */
-  int arg[3];
-  /* index, index jumps for rs_mesh array */
-  int q_ind = 0;
-  // make sure we have enough space
-  if (cp_cnt >= p3m.ca_num)
-    p3m_realloc_ca_fields(cp_cnt + 1);
-  // do it here, since p3m_realloc_ca_fields may change the address of
-  // p3m.ca_frac
-  double *cur_ca_frac = p3m.ca_frac.data() + cao * cao * cao * cp_cnt;
-
-  for (int d = 0; d < 3; d++) {
-    /* particle position in mesh coordinates */
-    auto const pos =
-        ((real_pos[d] - p3m.local_mesh.ld_pos[d]) * p3m.params.ai[d]) -
-        p3m.pos_shift;
-    /* nearest mesh point */
-    auto const nmp = (int)pos;
-    /* 3d-array index of nearest mesh point */
-    q_ind = (d == 0) ? nmp : nmp + p3m.local_mesh.dim[d] * q_ind;
-
-    if (!inter)
-      /* distance to nearest mesh point */
-      dist[d] = (pos - nmp) - 0.5;
-    else
-      /* distance to nearest mesh point for interpolation */
-      arg[d] = (int)((pos - nmp) * p3m.params.inter2);
-
-#ifdef ADDITIONAL_CHECKS
-    if (pos < -skin * p3m.params.ai[d]) {
-      fprintf(stderr, "%d: rs_mesh underflow! (pos %f)\n", this_node,
-              real_pos[d]);
-      fprintf(stderr, "%d: allowed coordinates: %f - %f\n", this_node,
-              local_geo.my_left()[d] - skin, local_geo.my_right()[d] + skin);
-    }
-    if ((nmp + cao) > p3m.local_mesh.dim[d]) {
-      fprintf(stderr, "%d: rs_mesh overflow! (pos %f, nmp=%d)\n", this_node,
-              real_pos[d], nmp);
-      fprintf(stderr, "%d: allowed coordinates: %f - %f\n", this_node,
-              local_geo.my_left()[d] - skin, local_geo.my_right()[d] + skin);
-    }
-#endif
-  }
-
-  if (cp_cnt >= 0)
-    p3m.ca_fmp[cp_cnt] = q_ind;
-
-  Utils::Array<double, cao> w_x, w_y, w_z;
-
-  if (!inter) {
-    for (int i = 0; i < cao; i++) {
-      using Utils::bspline;
-
-      w_x[i] = bspline<cao>(i, dist[0]);
-      w_y[i] = bspline<cao>(i, dist[1]);
-      w_z[i] = bspline<cao>(i, dist[2]);
-    }
-  } else {
-    for (int i = 0; i < cao; i++) {
-      w_x[i] = p3m.int_caf[i][arg[0]];
-      w_y[i] = p3m.int_caf[i][arg[1]];
-      w_z[i] = p3m.int_caf[i][arg[2]];
-    }
-  }
-
-  for (int i0 = 0; i0 < cao; i0++) {
-    auto const tmp0 = w_x[i0];
-    for (int i1 = 0; i1 < cao; i1++) {
-      auto const tmp1 = tmp0 * w_y[i1];
-      for (int i2 = 0; i2 < cao; i2++) {
-        auto const cur_ca_frac_val = q * tmp1 * w_z[i2];
-        p3m.rs_mesh[q_ind] += cur_ca_frac_val;
-        /* store current ca frac */
-        if (cp_cnt >= 0)
-          *(cur_ca_frac++) = cur_ca_frac_val;
-        q_ind++;
-      }
-      q_ind += p3m.local_mesh.q_2_off;
-    }
-    q_ind += p3m.local_mesh.q_21_off;
-  }
+void p3m_assign_charge(double q, const Utils::Vector3d &real_pos,
+                       p3m_interpolation_cache &inter_weights) {
+  Utils::integral_parameter<AssignCharge, 1, 7>(p3m.params.cao, q, real_pos,
+                                                p3m.params.ai, p3m.local_mesh,
+                                                inter_weights);
 }
 
-void p3m_shrink_wrap_charge_grid(int n_charges) {
-  /* we do not really want to export these */
-  if (n_charges < p3m.ca_num)
-    p3m_realloc_ca_fields(n_charges);
+void p3m_assign_charge(double q, const Utils::Vector3d &real_pos) {
+  Utils::integral_parameter<AssignCharge, 1, 7>(p3m.params.cao, q, real_pos,
+                                                p3m.params.ai, p3m.local_mesh);
 }
 
-/* Assign the forces obtained from k-space */
-template <int cao>
-static void P3M_assign_forces(double force_prefac,
-                              const ParticleRange &particles) {
-  /* charged particle counter, charge fraction counter */
-  int cp_cnt = 0;
-  int cf_cnt = 0;
-  /* index, index jumps for rs_mesh array */
-  int q_ind = 0;
-
-  for (auto &p : particles) {
-    auto const q = p.p.q;
-    if (q != 0.0) {
-      q_ind = p3m.ca_fmp[cp_cnt];
-      for (int i0 = 0; i0 < cao; i0++) {
-        for (int i1 = 0; i1 < cao; i1++) {
-          for (int i2 = 0; i2 < cao; i2++) {
-            p.f.f -= force_prefac * p3m.ca_frac[cf_cnt] *
-                     Utils::Vector3d{p3m.E_mesh[0][q_ind], p3m.E_mesh[1][q_ind],
-                                     p3m.E_mesh[2][q_ind]};
-
-            q_ind++;
-            cf_cnt++;
-          }
-          q_ind += p3m.local_mesh.q_2_off;
-        }
-        q_ind += p3m.local_mesh.q_21_off;
+namespace {
+template <size_t cao> struct AssignForces {
+  void operator()(double force_prefac, const ParticleRange &particles) const {
+    using Utils::make_const_span;
+    using Utils::Span;
+    using Utils::Vector;
+
+    assert(cao == p3m.inter_weights.cao());
+
+    /* charged particle counter, charge fraction counter */
+    int cp_cnt = 0;
+
+    for (auto &p : particles) {
+      auto const q = p.p.q;
+      if (q != 0.0) {
+        auto const pref = q * force_prefac;
+        auto const w = p3m.inter_weights.load<cao>(cp_cnt++);
+
+        Utils::Vector3d E{};
+        p3m_interpolate(p3m.local_mesh, w, [&E](int ind, double w) {
+          E += w * Utils::Vector3d{p3m.E_mesh[0][ind], p3m.E_mesh[1][ind],
+                                   p3m.E_mesh[2][ind]};
+        });
+
+        p.f.f -= pref * E;
       }
-      cp_cnt++;
     }
   }
-}
+};
 
-namespace {
 auto dipole_moment(Particle const &p, BoxGeometry const &box) {
   return p.p.q * unfolded_position(p.r.p, p.l.i, box.length());
 }
@@ -699,7 +501,7 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
 
   // Note: after these calls, the grids are in the order yzx and not xyz
   // anymore!!!
-  /* The dipol moment is only needed if we don't have metallic boundaries. */
+  /* The dipole moment is only needed if we don't have metallic boundaries. */
   auto const box_dipole = (p3m.params.epsilon != P3M_EPSILON_METALLIC)
                               ? boost::make_optional(calc_dipole_moment(
                                     comm_cart, particles, box_geo))
@@ -754,30 +556,8 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
                          p3m.local_mesh.dim);
     }
 
-    /* Assign force component from mesh to particle */
-    switch (p3m.params.cao) {
-    case 1:
-      P3M_assign_forces<1>(force_prefac, particles);
-      break;
-    case 2:
-      P3M_assign_forces<2>(force_prefac, particles);
-      break;
-    case 3:
-      P3M_assign_forces<3>(force_prefac, particles);
-      break;
-    case 4:
-      P3M_assign_forces<4>(force_prefac, particles);
-      break;
-    case 5:
-      P3M_assign_forces<5>(force_prefac, particles);
-      break;
-    case 6:
-      P3M_assign_forces<6>(force_prefac, particles);
-      break;
-    case 7:
-      P3M_assign_forces<7>(force_prefac, particles);
-      break;
-    }
+    Utils::integral_parameter<AssignForces, 1, 7>(p3m.params.cao, force_prefac,
+                                                  particles);
 
     if (p3m.params.epsilon != P3M_EPSILON_METALLIC) {
       add_dipole_correction(box_dipole.value(), particles);
@@ -817,20 +597,6 @@ double p3m_calc_kspace_forces(bool force_flag, bool energy_flag,
   return 0.0;
 }
 
-/************************************************************/
-
-void p3m_realloc_ca_fields(int newsize) {
-  newsize = ((newsize + CA_INCREMENT - 1) / CA_INCREMENT) * CA_INCREMENT;
-  if (newsize == p3m.ca_num)
-    return;
-  if (newsize < CA_INCREMENT)
-    newsize = CA_INCREMENT;
-
-  p3m.ca_num = newsize;
-  p3m.ca_frac.resize(p3m.params.cao3 * p3m.ca_num);
-  p3m.ca_fmp.resize(p3m.ca_num);
-}
-
 void p3m_calc_meshift() {
   p3m.meshift_x.resize(p3m.params.mesh[0]);
   p3m.meshift_y.resize(p3m.params.mesh[1]);
diff --git a/src/core/electrostatics_magnetostatics/p3m.hpp b/src/core/electrostatics_magnetostatics/p3m.hpp
index fd733a52071..5c72dac27ed 100644
--- a/src/core/electrostatics_magnetostatics/p3m.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m.hpp
@@ -39,6 +39,7 @@
 
 #include "fft.hpp"
 #include "p3m-common.hpp"
+#include "p3m_interpolation.hpp"
 #include "p3m_send_mesh.hpp"
 
 #include <ParticleRange.hpp>
@@ -68,11 +69,6 @@ struct p3m_data_struct {
   /** square of sum of charges (only on master node). */
   double square_sum_q;
 
-  /** interpolation of the charge assignment function. */
-  std::array<std::vector<double>, 7> int_caf;
-
-  /** position shift for calc. of first assignment mesh point. */
-  double pos_shift;
   /** help variable for calculation of aliasing sums */
   std::vector<double> meshift_x;
   std::vector<double> meshift_y;
@@ -86,12 +82,7 @@ struct p3m_data_struct {
   /** Energy optimised influence function (k-space) */
   std::vector<double> g_energy;
 
-  /** number of charged particles on the node. */
-  int ca_num;
-  /** Charge fractions for mesh assignment. */
-  std::vector<double> ca_frac;
-  /** index of first mesh point for charge assignment. */
-  std::vector<int> ca_fmp;
+  p3m_interpolation_cache inter_weights;
 
   /** number of permutations in k_space */
   int ks_pnum;
@@ -173,9 +164,6 @@ bool p3m_sanity_checks();
 void p3m_count_charged_particles();
 
 /** Assign the physical charges using the tabulated charge assignment function.
- *  The charge fractions are buffered
- *  in @ref p3m_data_struct::ca_fmp "ca_fmp" and @ref p3m_data_struct::ca_frac
- *  "ca_frac".
  */
 void p3m_charge_assign(const ParticleRange &particles);
 
@@ -183,15 +171,12 @@ void p3m_charge_assign(const ParticleRange &particles);
  *
  *  @param[in] q          %Particle charge
  *  @param[in] real_pos   %Particle position in real space
- *  @param[in] cp_cnt     The running index, which may be smaller than 0, in
- *                        which case the charge is assumed to be virtual and
- *                        is not stored in the @ref p3m_data_struct::ca_frac
- *                        "ca_frac" arrays
+ *  @param[in] inter_weights Cached interpolation weights to be used.
  */
-void p3m_assign_charge(double q, const Utils::Vector3d &real_pos, int cp_cnt);
-
-/** Shrink wrap the charge grid */
-void p3m_shrink_wrap_charge_grid(int n_charges);
+void p3m_assign_charge(double q, const Utils::Vector3d &real_pos,
+                       p3m_interpolation_cache &inter_weights);
+/** @overload */
+void p3m_assign_charge(double q, const Utils::Vector3d &real_pos);
 
 /** Calculate real space contribution of Coulomb pair forces. */
 inline void p3m_add_pair_force(double q1q2, Utils::Vector3d const &d,
@@ -226,10 +211,9 @@ inline void p3m_add_pair_force(double q1q2, Utils::Vector3d const &d,
  *  @param[in]  cao          @copybrief P3MParameters::cao
  *  @param[in]  alpha        @copybrief P3MParameters::alpha
  *  @param[in]  accuracy     @copybrief P3MParameters::accuracy
- *  @param[in]  n_interpol   @copybrief P3MParameters::inter
  */
 void p3m_set_tune_params(double r_cut, const int mesh[3], int cao, double alpha,
-                         double accuracy, int n_interpol);
+                         double accuracy);
 
 /** Set custom parameters
  *
@@ -256,12 +240,6 @@ int p3m_set_mesh_offset(double x, double y, double z);
  */
 int p3m_set_eps(double eps);
 
-/** Set @ref P3MParameters::inter "inter" parameter
- *
- *  @param[in]  n            @copybrief P3MParameters::inter
- */
-int p3m_set_ninterpol(int n);
-
 /** Calculate real space contribution of Coulomb pair energy. */
 inline double p3m_pair_energy(double chgfac, double dist) {
   if (dist < p3m.params.r_cut && dist != 0) {
diff --git a/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp
new file mode 100644
index 00000000000..56e638fa155
--- /dev/null
+++ b/src/core/electrostatics_magnetostatics/p3m_interpolation.hpp
@@ -0,0 +1,207 @@
+/*
+ * Copyright (C) 2010-2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef ESPRESSO_P3M_INTERPOLATION_HPP
+#define ESPRESSO_P3M_INTERPOLATION_HPP
+
+#include <utils/Span.hpp>
+#include <utils/index.hpp>
+#include <utils/math/bspline.hpp>
+
+#include <boost/range/algorithm/copy.hpp>
+
+#include <cassert>
+#include <tuple>
+#include <vector>
+
+/**
+ * @brief Interpolation weights for one point.
+ *
+ * Interpolation weights and  grid offset for one point.
+ *
+ * @tparam cao Interpolation order.
+ */
+template <int cao> struct InterpolationWeights {
+  /** Linear index of the corner of the interpolation cube. */
+  int ind;
+  /** Weights for the directions */
+  Utils::Array<double, cao> w_x, w_y, w_z;
+};
+
+/**
+ * @brief Cache for interpolation weights.
+ *
+ * This is a storage container for interpolation weights of
+ * type InterpolationWeights.
+ */
+class p3m_interpolation_cache {
+  size_t m_cao = 0;
+  /** Charge fractions for mesh assignment. */
+  std::vector<double> ca_frac;
+  /** index of first mesh point for charge assignment. */
+  std::vector<int> ca_fmp;
+
+public:
+  /**
+   * @brief Number of points in the cache.
+   * @return Number of points currently in the cache.
+   */
+  auto size() const { return ca_fmp.size(); }
+
+  /**
+   * @brief Charge assignment order the weights are for.
+   * @return The charge assigment order.
+   */
+  auto cao() const { return m_cao; }
+
+  /**
+   * @brief Push back weights for one point.
+   *
+   * @tparam cao Interpolation order has to match the order
+   *         set at last call to @ref p3m_interpolation_cache::reset.
+   * @param w Interpolation weights to store.
+   */
+  template <int cao> void store(const InterpolationWeights<cao> &w) {
+    assert(cao == m_cao);
+
+    ca_fmp.push_back(w.ind);
+    auto it = std::back_inserter(ca_frac);
+    boost::copy(w.w_x, it);
+    boost::copy(w.w_y, it);
+    boost::copy(w.w_z, it);
+  }
+
+  /**
+   * @brief Load entry from the cache.
+   *
+   * This loads an entry at an index from the cache,
+   * the entries are indexed by the order they were stored.
+   *
+   * @tparam cao Interpolation order has to match the order
+   *         set at last call to @ref p3m_interpolation_cache::reset.
+   * @param i Index of the entry to load.
+   * @return i-it interpolation weights.
+   */
+  template <int cao> InterpolationWeights<cao> load(size_t i) const {
+    assert(cao == m_cao);
+
+    using Utils::make_const_span;
+    assert(i < size());
+
+    InterpolationWeights<cao> ret;
+    ret.ind = ca_fmp[i];
+
+    auto const offset = ca_frac.data() + 3 * i * m_cao;
+    boost::copy(make_const_span(offset + 0 * m_cao, m_cao), ret.w_x.begin());
+    boost::copy(make_const_span(offset + 1 * m_cao, m_cao), ret.w_y.begin());
+    boost::copy(make_const_span(offset + 2 * m_cao, m_cao), ret.w_z.begin());
+
+    return ret;
+  }
+
+  /**
+   * @brief Reset the cache.
+   *
+   * @param cao Interpolation order.
+   */
+  void reset(int cao) {
+    m_cao = cao;
+    ca_frac.clear();
+    ca_fmp.clear();
+  }
+};
+
+/**
+ * @brief Calculate the P-th order interpolation weights.
+ *
+ * As described in from @cite hockney88a 5-189 (or 8-61).
+ * The weights are also tabulated in @cite deserno98a @cite deserno98b.
+ */
+template <int cao>
+InterpolationWeights<cao>
+p3m_calculate_interpolation_weights(const Utils::Vector3d &position,
+                                    const Utils::Vector3d &ai,
+                                    p3m_local_mesh const &local_mesh) {
+  /** position shift for calc. of first assignment mesh point. */
+  static auto const pos_shift = std::floor((cao - 1) / 2.0) - (cao % 2) / 2.0;
+
+  /* distance to nearest mesh point */
+  double dist[3];
+
+  /* nearest mesh point */
+  Utils::Vector3i nmp;
+
+  for (int d = 0; d < 3; d++) {
+    /* particle position in mesh coordinates */
+    auto const pos = ((position[d] - local_mesh.ld_pos[d]) * ai[d]) - pos_shift;
+
+    nmp[d] = (int)pos;
+
+    /* distance to nearest mesh point */
+    dist[d] = (pos - nmp[d]) - 0.5;
+  }
+
+  InterpolationWeights<cao> ret;
+
+  /* 3d-array index of nearest mesh point */
+  ret.ind = Utils::get_linear_index(nmp, local_mesh.dim,
+                                    Utils::MemoryOrder::ROW_MAJOR);
+  for (int i = 0; i < cao; i++) {
+    using Utils::bspline;
+
+    ret.w_x[i] = bspline<cao>(i, dist[0]);
+    ret.w_y[i] = bspline<cao>(i, dist[1]);
+    ret.w_z[i] = bspline<cao>(i, dist[2]);
+  }
+
+  return ret;
+}
+
+/**
+ * @brief P3M grid interpolation.
+ *
+ * This runs an kernel for every interpolation point
+ * in a set of interpolation weights with the linear
+ * grid index and the weight of the point as arguments.
+ *
+ * @param local_mesh Mesh info.
+ * @param weights Set of weights
+ * @param kernel The kernel to run.
+ */
+template <int cao, class Kernel>
+void p3m_interpolate(p3m_local_mesh const &local_mesh,
+                     InterpolationWeights<cao> const &weights, Kernel kernel) {
+  auto q_ind = weights.ind;
+  for (int i0 = 0; i0 < cao; i0++) {
+    auto const tmp0 = weights.w_x[i0];
+    for (int i1 = 0; i1 < cao; i1++) {
+      auto const tmp1 = tmp0 * weights.w_y[i1];
+      for (int i2 = 0; i2 < cao; i2++) {
+        kernel(q_ind, tmp1 * weights.w_z[i2]);
+
+        q_ind++;
+      }
+      q_ind += local_mesh.q_2_off;
+    }
+    q_ind += local_mesh.q_21_off;
+  }
+}
+
+#endif // ESPRESSO_P3M_INTERPOLATION_HPP
diff --git a/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp b/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp
index 674a83c5a2a..8984a53989e 100644
--- a/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp
+++ b/src/core/electrostatics_magnetostatics/p3m_send_mesh.cpp
@@ -103,7 +103,7 @@ void p3m_send_mesh::resize(const boost::mpi::communicator &comm,
 
 void p3m_send_mesh::gather_grid(Utils::Span<double *> meshes,
                                 const boost::mpi::communicator &comm,
-                                const int *dim) {
+                                const Utils::Vector3i &dim) {
   auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm);
   send_grid.resize(max * meshes.size());
   recv_grid.resize(max * meshes.size());
@@ -116,7 +116,7 @@ void p3m_send_mesh::gather_grid(Utils::Span<double *> meshes,
     if (s_size[s_dir] > 0)
       for (size_t i = 0; i < meshes.size(); i++) {
         fft_pack_block(meshes[i], send_grid.data() + i * s_size[s_dir],
-                       s_ld[s_dir], s_dim[s_dir], dim, 1);
+                       s_ld[s_dir], s_dim[s_dir], dim.data(), 1);
       }
 
     /* communication */
@@ -133,7 +133,7 @@ void p3m_send_mesh::gather_grid(Utils::Span<double *> meshes,
     if (r_size[r_dir] > 0) {
       for (size_t i = 0; i < meshes.size(); i++) {
         p3m_add_block(recv_grid.data() + i * r_size[r_dir], meshes[i],
-                      r_ld[r_dir], r_dim[r_dir], dim);
+                      r_ld[r_dir], r_dim[r_dir], dim.data());
       }
     }
   }
@@ -141,7 +141,7 @@ void p3m_send_mesh::gather_grid(Utils::Span<double *> meshes,
 
 void p3m_send_mesh::spread_grid(Utils::Span<double *> meshes,
                                 const boost::mpi::communicator &comm,
-                                const int *dim) {
+                                const Utils::Vector3i &dim) {
   auto const node_neighbors = Utils::Mpi::cart_neighbors<3>(comm);
   send_grid.resize(max * meshes.size());
   recv_grid.resize(max * meshes.size());
@@ -154,7 +154,7 @@ void p3m_send_mesh::spread_grid(Utils::Span<double *> meshes,
     if (r_size[r_dir] > 0)
       for (size_t i = 0; i < meshes.size(); i++) {
         fft_pack_block(meshes[i], send_grid.data() + i * r_size[r_dir],
-                       r_ld[r_dir], r_dim[r_dir], dim, 1);
+                       r_ld[r_dir], r_dim[r_dir], dim.data(), 1);
       }
     /* communication */
     if (node_neighbors[r_dir] != comm.rank()) {
@@ -170,7 +170,7 @@ void p3m_send_mesh::spread_grid(Utils::Span<double *> meshes,
     if (s_size[s_dir] > 0) {
       for (size_t i = 0; i < meshes.size(); i++) {
         fft_unpack_block(recv_grid.data() + i * s_size[s_dir], meshes[i],
-                         s_ld[s_dir], s_dim[s_dir], dim, 1);
+                         s_ld[s_dir], s_dim[s_dir], dim.data(), 1);
       }
     }
   }
diff --git a/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp b/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp
index 6197d90cd8c..d4e1106d1a0 100644
--- a/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp
+++ b/src/core/electrostatics_magnetostatics/p3m_send_mesh.hpp
@@ -66,15 +66,17 @@ class p3m_send_mesh {
   void resize(const boost::mpi::communicator &comm,
               const p3m_local_mesh &local_mesh);
   void gather_grid(Utils::Span<double *> meshes,
-                   const boost::mpi::communicator &comm, const int *dim);
+                   const boost::mpi::communicator &comm,
+                   const Utils::Vector3i &dim);
   void gather_grid(double *mesh, const boost::mpi::communicator &comm,
-                   const int *dim) {
+                   const Utils::Vector3i &dim) {
     gather_grid(Utils::make_span(&mesh, 1), comm, dim);
   }
   void spread_grid(Utils::Span<double *> meshes,
-                   const boost::mpi::communicator &comm, const int *dim);
+                   const boost::mpi::communicator &comm,
+                   const Utils::Vector3i &dim);
   void spread_grid(double *mesh, const boost::mpi::communicator &comm,
-                   const int *dim) {
+                   const Utils::Vector3i &dim) {
     spread_grid(Utils::make_span(&mesh, 1), comm, dim);
   }
 };
diff --git a/src/python/espressomd/electrostatics.pxd b/src/python/espressomd/electrostatics.pxd
index ba0c0ba02d8..0dab50df9c1 100644
--- a/src/python/espressomd/electrostatics.pxd
+++ b/src/python/espressomd/electrostatics.pxd
@@ -65,31 +65,13 @@ IF ELECTROSTATICS:
         void deactivate_method()
 
     IF P3M:
-        cdef extern from "electrostatics_magnetostatics/p3m-common.hpp":
-            ctypedef struct P3MParameters:
-                double alpha_L
-                double r_cut_iL
-                int    mesh[3]
-                double mesh_off[3]
-                int    cao
-                int    inter
-                double accuracy
-                double epsilon
-                double cao_cut[3]
-                double a[3]
-                double ai[3]
-                double alpha
-                double r_cut
-                int    inter2
-                int    cao3
-                double additional_mesh[3]
+        from p3m_common cimport P3MParameters
 
         cdef extern from "electrostatics_magnetostatics/p3m.hpp":
             int p3m_set_params(double r_cut, int * mesh, int cao, double alpha, double accuracy)
-            void p3m_set_tune_params(double r_cut, int mesh[3], int cao, double alpha, double accuracy, int n_interpol)
+            void p3m_set_tune_params(double r_cut, int mesh[3], int cao, double alpha, double accuracy)
             int p3m_set_mesh_offset(double x, double y, double z)
             int p3m_set_eps(double eps)
-            int p3m_set_ninterpol(int n)
             int p3m_adaptive_tune(char ** log)
 
             ctypedef struct p3m_data_struct:
@@ -153,19 +135,17 @@ IF ELECTROSTATICS:
 
             return p3m_set_params(r_cut, mesh, cao, alpha, accuracy)
 
-        cdef inline python_p3m_set_tune_params(p_r_cut, p_mesh, p_cao, p_alpha, p_accuracy, p_n_interpol):
-            # cdef inline python_p3m_set_tune_params():
+        cdef inline python_p3m_set_tune_params(p_r_cut, p_mesh, p_cao, p_alpha, p_accuracy):
             cdef int mesh[3]
             cdef double r_cut
             cdef int cao
             cdef double alpha
             cdef double accuracy
-            cdef int n_interpol
             r_cut = p_r_cut
             cao = p_cao
             alpha = p_alpha
             accuracy = p_accuracy
-            n_interpol = p_n_interpol
+
             if is_valid_type(p_mesh, int):
                 mesh[0] = p_mesh
                 mesh[1] = p_mesh
@@ -173,7 +153,7 @@ IF ELECTROSTATICS:
             else:
                 mesh = p_mesh
 
-            p3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy, n_interpol)
+            p3m_set_tune_params(r_cut, mesh, cao, alpha, accuracy)
 
     cdef extern from "electrostatics_magnetostatics/debye_hueckel.hpp":
         ctypedef struct Debye_hueckel_params:
diff --git a/src/python/espressomd/electrostatics.pyx b/src/python/espressomd/electrostatics.pyx
index ceb0bf59717..884d931740f 100644
--- a/src/python/espressomd/electrostatics.pyx
+++ b/src/python/espressomd/electrostatics.pyx
@@ -267,10 +267,6 @@ IF P3M == 1:
                     or self._params["epsilon"] == "metallic"):
                 raise ValueError("epsilon should be a double or 'metallic'")
 
-            if not (self._params["inter"] == default_params["inter"]
-                    or self._params["inter"] >= 0):
-                raise ValueError("inter should be a positive integer")
-
             if not (self._params["mesh_off"] == default_params["mesh_off"]
                     or len(self._params) != 3):
                 raise ValueError(
@@ -283,14 +279,13 @@ IF P3M == 1:
 
         def valid_keys(self):
             return ["mesh", "cao", "accuracy", "epsilon", "alpha", "r_cut",
-                    "prefactor", "tune", "check_neutrality", "inter"]
+                    "prefactor", "tune", "check_neutrality"]
 
         def required_keys(self):
             return ["prefactor", "accuracy"]
 
         def default_params(self):
             return {"cao": 0,
-                    "inter": -1,
                     "r_cut": -1,
                     "alpha": 0,
                     "accuracy": 0,
@@ -319,8 +314,6 @@ IF P3M == 1:
             #         which resets r_cut if lb is zero. OK.
             # Sets eps, bcast
             p3m_set_eps(self._params["epsilon"])
-            # Sets ninterpol, bcast
-            p3m_set_ninterpol(self._params["inter"])
             python_p3m_set_mesh_offset(self._params["mesh_off"])
 
         def _tune(self):
@@ -329,8 +322,7 @@ IF P3M == 1:
                                        self._params["mesh"],
                                        self._params["cao"],
                                        -1.0,
-                                       self._params["accuracy"],
-                                       self._params["inter"])
+                                       self._params["accuracy"])
             resp = python_p3m_adaptive_tune()
             if resp:
                 raise Exception(
@@ -419,10 +411,6 @@ IF P3M == 1:
                     raise ValueError(
                         "epsilon should be a double or 'metallic'")
 
-                if not (self._params["inter"] == default_params["inter"]
-                        or self._params["inter"] > 0):
-                    raise ValueError("inter should be a positive integer")
-
                 if not (self._params["mesh_off"] == default_params["mesh_off"]
                         or len(self._params) != 3):
                     raise ValueError(
@@ -437,7 +425,6 @@ IF P3M == 1:
 
             def default_params(self):
                 return {"cao": 0,
-                        "inter": -1,
                         "r_cut": -1,
                         "alpha": 0,
                         "accuracy": 0,
@@ -460,8 +447,7 @@ IF P3M == 1:
                                            self._params["mesh"],
                                            self._params["cao"],
                                            -1.0,
-                                           self._params["accuracy"],
-                                           self._params["inter"])
+                                           self._params["accuracy"])
                 resp = python_p3m_adaptive_tune()
                 if resp:
                     raise Exception(
@@ -485,7 +471,6 @@ IF P3M == 1:
                                       self._params["alpha"],
                                       self._params["accuracy"])
                 p3m_set_eps(self._params["epsilon"])
-                p3m_set_ninterpol(self._params["inter"])
                 python_p3m_set_mesh_offset(self._params["mesh_off"])
                 handle_errors("p3m gpu init")
 
diff --git a/src/python/espressomd/magnetostatics.pxd b/src/python/espressomd/magnetostatics.pxd
index 4373c0a509e..14ef94f12a0 100644
--- a/src/python/espressomd/magnetostatics.pxd
+++ b/src/python/espressomd/magnetostatics.pxd
@@ -63,7 +63,6 @@ IF DP3M == 1:
         void dp3m_set_tune_params(double r_cut, int mesh, int cao, double alpha, double accuracy, int n_interpol)
         int dp3m_set_mesh_offset(double x, double y, double z)
         int dp3m_set_eps(double eps)
-        int dp3m_set_ninterpol(int n)
         int dp3m_adaptive_tune(char ** log)
         int dp3m_deactivate()
 
diff --git a/src/python/espressomd/magnetostatics.pyx b/src/python/espressomd/magnetostatics.pyx
index a5c18babd14..d662e16eac0 100644
--- a/src/python/espressomd/magnetostatics.pyx
+++ b/src/python/espressomd/magnetostatics.pyx
@@ -165,7 +165,6 @@ IF DP3M == 1:
         def _set_params_in_es_core(self):
             self.set_magnetostatics_prefactor()
             dp3m_set_eps(self._params["epsilon"])
-            dp3m_set_ninterpol(self._params["inter"])
             self.python_dp3m_set_mesh_offset(self._params["mesh_off"])
             self.python_dp3m_set_params(
                 self._params["r_cut"], self._params["mesh"],
@@ -230,7 +229,6 @@ IF DP3M == 1:
 
         def python_dp3m_set_tune_params(self, p_r_cut, p_mesh, p_cao, p_alpha,
                                         p_accuracy, p_n_interpol):
-            # cdef  python_p3m_set_tune_params():
             cdef int mesh
             cdef double r_cut
             cdef int cao
diff --git a/src/python/espressomd/p3m_common.pxd b/src/python/espressomd/p3m_common.pxd
index c6c99b4be71..3f7b685038c 100644
--- a/src/python/espressomd/p3m_common.pxd
+++ b/src/python/espressomd/p3m_common.pxd
@@ -24,14 +24,10 @@ IF P3M == 1 or DP3M == 1:
             int    mesh[3]
             double mesh_off[3]
             int    cao
-            int    inter
             double accuracy
             double epsilon
             double cao_cut[3]
             double a[3]
-            double ai[3]
             double alpha
             double r_cut
-            int    inter2
-            int    cao3
             double additional_mesh[3]
diff --git a/src/utils/include/utils/Vector.hpp b/src/utils/include/utils/Vector.hpp
index bdebb05c061..fa4b0158ee7 100644
--- a/src/utils/include/utils/Vector.hpp
+++ b/src/utils/include/utils/Vector.hpp
@@ -83,7 +83,7 @@ template <typename T, std::size_t N> class Vector : public Array<T, N> {
   }
 
   template <typename InputIterator>
-  Vector(InputIterator first, InputIterator last) {
+  Vector(InputIterator first, InputIterator last) : Base() {
     if (std::distance(first, last) == N) {
       std::copy_n(first, N, begin());
     } else {
diff --git a/src/utils/include/utils/integral_parameter.hpp b/src/utils/include/utils/integral_parameter.hpp
new file mode 100644
index 00000000000..1df9595b0c2
--- /dev/null
+++ b/src/utils/include/utils/integral_parameter.hpp
@@ -0,0 +1,63 @@
+/*
+ * Copyright (C) 2020 The ESPResSo project
+ * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
+ *   Max-Planck-Institute for Polymer Research, Theory Group
+ *
+ * This file is part of ESPResSo.
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+#ifndef ESPRESSO_INTEGRAL_PARAMETER_HPP
+#define ESPRESSO_INTEGRAL_PARAMETER_HPP
+
+#include <stdexcept>
+
+namespace Utils {
+namespace detail {
+template <template <size_t> class F, size_t I, size_t N>
+struct integral_parameter_impl {
+  template <class... Args>
+  static decltype(auto) eval(size_t i, Args &&... args) {
+    if (i == I)
+      return F<I>{}(std::forward<Args>(args)...);
+
+    return integral_parameter_impl<F, I + 1, N>::eval(
+        i, std::forward<Args>(args)...);
+  }
+};
+
+template <template <size_t> class F, size_t N>
+struct integral_parameter_impl<F, N, N> {
+  template <class... Args>
+  static decltype(auto) eval(size_t i, Args &&... args) {
+    if (i == N)
+      return F<N>{}(std::forward<Args>(args)...);
+
+    throw std::runtime_error("Invalid parameter value");
+  }
+};
+} // namespace detail
+
+/**
+ * @brief Generate a call table for a integral non-type template parameter.
+ */
+template <template <size_t> class F, size_t min, size_t max, class... Args>
+decltype(auto) integral_parameter(size_t i, Args &&... args) {
+  return detail::integral_parameter_impl<F, min, max>::eval(
+      i, std::forward<Args>(args)...);
+}
+
+} // namespace Utils
+
+#endif // ESPRESSO_INTEGRAL_PARAMETER_HPP
diff --git a/src/utils/tests/CMakeLists.txt b/src/utils/tests/CMakeLists.txt
index 4f6226a2775..4066270c62c 100644
--- a/src/utils/tests/CMakeLists.txt
+++ b/src/utils/tests/CMakeLists.txt
@@ -46,6 +46,8 @@ unit_test(NAME memcpy_archive_test SRC memcpy_archive_test.cpp DEPENDS utils)
 unit_test(NAME triangle_functions_test SRC triangle_functions_test.cpp DEPENDS
           utils)
 unit_test(NAME Bag_test SRC Bag_test.cpp DEPENDS utils Boost::serialization)
+unit_test(NAME integral_parameter_test SRC integral_parameter_test.cpp DEPENDS
+          utils)
 
 unit_test(NAME gather_buffer_test SRC gather_buffer_test.cpp DEPENDS utils
           Boost::mpi MPI::MPI_CXX NUM_PROC 4)
diff --git a/src/utils/tests/integral_parameter_test.cpp b/src/utils/tests/integral_parameter_test.cpp
new file mode 100644
index 00000000000..01c70b79ea3
--- /dev/null
+++ b/src/utils/tests/integral_parameter_test.cpp
@@ -0,0 +1,46 @@
+/*
+ * Copyright (C) 2017-2019 The ESPResSo project
+ *
+ * ESPResSo is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * ESPResSo is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#define BOOST_TEST_MODULE integral_parameter test
+#define BOOST_TEST_DYN_LINK
+#include <boost/test/unit_test.hpp>
+
+#include <utils/integral_parameter.hpp>
+
+#include <utility>
+
+template <size_t I> struct F {
+  template <class T> auto operator()(T arg) const {
+    return std::make_pair(I, arg);
+  }
+};
+
+BOOST_AUTO_TEST_CASE(integral_parameter_) {
+  static_assert(
+      std::is_same<decltype(Utils::integral_parameter<F, 1, 5>(5, 13)),
+                   std::pair<size_t, int>>::value,
+      "");
+
+  BOOST_CHECK(std::make_pair(1ul, 13) ==
+              (Utils::integral_parameter<F, 1, 5>(1, 13)));
+  BOOST_CHECK(std::make_pair(3ul, 13) ==
+              (Utils::integral_parameter<F, 1, 5>(3, 13)));
+  BOOST_CHECK(std::make_pair(5ul, 13) ==
+              (Utils::integral_parameter<F, 1, 5>(5, 13)));
+  BOOST_CHECK_THROW((Utils::integral_parameter<F, 1, 5>(6, 13)),
+                    std::exception);
+}
\ No newline at end of file
diff --git a/testsuite/python/analyze_energy.py b/testsuite/python/analyze_energy.py
index 03000fb49ed..a9768b078fe 100644
--- a/testsuite/python/analyze_energy.py
+++ b/testsuite/python/analyze_energy.py
@@ -180,7 +180,7 @@ def test_electrostatics(self):
         self.assertAlmostEqual(energy["kinetic"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["bonded"], 0., delta=1e-7)
         self.assertAlmostEqual(energy["non_bonded"], 0, delta=1e-7)
-        self.assertAlmostEqual(energy["coulomb"], u_p3m, delta=1e-7)
+        self.assertAlmostEqual(energy["coulomb"], u_p3m, delta=1e-5)
         self.system.part[0].q = 0
         self.system.part[1].q = 0
         self.system.part[0].pos = [1, 2, 2]
diff --git a/testsuite/python/coulomb_cloud_wall.py b/testsuite/python/coulomb_cloud_wall.py
index e2fa89e9053..ab72e11b1b2 100644
--- a/testsuite/python/coulomb_cloud_wall.py
+++ b/testsuite/python/coulomb_cloud_wall.py
@@ -88,26 +88,9 @@ def compare(self, method_name, energy=True, prefactor=None):
     # Tests for individual methods
 
     @utx.skipIfMissingFeatures(["P3M"])
-    def test_p3m_direct_caf(self):
+    def test_p3m_direct(self):
         """
-        This checks P3M with using the charge assignment
-        function (window function) directly by setting the
-        `inter` parameter to zero.
-
-        """
-
-        self.S.actors.add(
-            espressomd.electrostatics.P3M(
-                prefactor=3, r_cut=1.001, accuracy=1e-3,
-                mesh=64, cao=7, alpha=2.70746, tune=False, inter=0))
-        self.S.integrator.run(0)
-        self.compare("p3m", energy=True, prefactor=3)
-
-    @utx.skipIfMissingFeatures(["P3M"])
-    def test_p3m_interpolated_caf(self):
-        """
-        This checks P3M with using an interpolated charge assignment
-        function (window function), which is the default.
+        This checks P3M.
 
         """
 
diff --git a/testsuite/python/electrostaticInteractions.py b/testsuite/python/electrostaticInteractions.py
index 13eb6c5047c..7dcbf6301dc 100644
--- a/testsuite/python/electrostaticInteractions.py
+++ b/testsuite/python/electrostaticInteractions.py
@@ -93,13 +93,13 @@ def test_p3m(self):
                                             tune=False)
         self.system.actors.add(p3m)
         self.assertAlmostEqual(self.system.analysis.energy()['coulomb'],
-                               p3m_energy)
+                               p3m_energy, places=5)
         # need to update forces
         self.system.integrator.run(0)
         np.testing.assert_allclose(np.copy(self.system.part[0].f),
-                                   [p3m_force, 0, 0], atol=1E-5)
+                                   [p3m_force, 0, 0], atol=1E-4)
         np.testing.assert_allclose(np.copy(self.system.part[1].f),
-                                   [-p3m_force, 0, 0], atol=1E-10)
+                                   [-p3m_force, 0, 0], atol=1E-5)
         self.system.actors.remove(p3m)
 
     def test_dh(self):