P3M interpolation cleanup (#3434)

Description of changes:
 - Removed interpolation function duplication from P3M/DP3M
 - Reduced redundancy from cached interpolation weights
 - Fixed syntax of p3m benchmark

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"]))
 

src/core/electrostatics_magnetostatics/elc.cpp

@@ -1301,35 +1301,32 @@ 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)) {
     auto const q_eff = elc_params.delta_mid_top * p.p.q;
     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) {

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) {

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);

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;