This is inspired by the ESPResSo tuning code. In their tuning code, they take the running time into account (which we have not yet). Their code combines grid-search (for the mesh density and the interpolation order) and handcraft (for the smearing and the real space cutoff) together. Basically, what they have done is:

1. Search the mesh_density iteratively from small to large in a certain range (m_mesh_density_min, m_mesh_density_max) .
   https://github.com/espressomd/espresso/blob/a8a00a06e4eb8b0726a318739d069e379619e400/src/core/electrostatics/p3m.cpp#L715-L754
2. In the iteration of 1, for each mesh_density, iterate cao (charge assignment order, our interpolation_nodes) from low to high. https://github.com/espressomd/espresso/blob/a8a00a06e4eb8b0726a318739d069e379619e400/src/core/p3m/TuningAlgorithm.cpp#L228-L326
3. For each pair of (mesh_density, cao), bisectionly search r_cut_iL
   https://github.com/espressomd/espresso/blob/a8a00a06e4eb8b0726a318739d069e379619e400/src/core/p3m/TuningAlgorithm.cpp#L155-L168
   and calculate an alpha_L (our smearing) according to the real space error and the desired accuracy
   https://github.com/espressomd/espresso/blob/a8a00a06e4eb8b0726a318739d069e379619e400/src/core/electrostatics/p3m.cpp#L646-L655
4. For each pair of (mesh_density, cao, r_cut_iL, alpha_L), initialize P3M, run the calculation and record the calculation time.
5. Return a pair of parameters with the shortest calculation time and satisfying the accuracy requirement.

There are some complicated exit conditions for the loops in 1 and 2, in order to shrink the size of the search space. The methods in 3 look like heuristic, and perhaps we can replace them with gradient-based optimization. For our case, we can grid-search each possible (mesh_spacing, interpolation_nodes), and use Adam optimizer to get (smearing, rc_cutoff). And for each possible parameter, we can do the calculation, record the time, and finally return a parameter pair with the shortest time.

I have already tried this for PME and here is the result, but I forgot to record the calculation time. One other problem is that the tuning time is very long. It took the computer 76 min to get all the results below (7 * 8 = 56 cases, avg. 1.4 min/case), so we may introduce some rules to help reduce the scale of searching.