38#include "communication.hpp"
40#include "system/System.hpp"
46#include <boost/mpi/collectives/all_reduce.hpp>
47#include <boost/range/combine.hpp>
108template <std::
size_t dir>
110 std::size_t n_freq,
double u) {
112 auto const n_part = particles.
size();
113 std::vector<SCCache> ret(n_freq * n_part);
115 for (std::size_t freq = 1; freq <= n_freq; freq++) {
116 auto const pref = c_2pi *
u *
static_cast<double>(freq);
118 std::size_t o = (freq - 1) * n_part;
119 for (
auto const &p : particles) {
120 auto const arg = pref * p.pos()[dir];
121 ret[o++] = {sin(arg), cos(arg)};
128static std::pair<std::size_t, std::size_t>
131 assert(far_cut >= 0.);
132 auto const n_freq_x =
133 static_cast<std::size_t
>(std::ceil(far_cut * box_geo.
length()[0]) + 1.);
134 auto const n_freq_y =
135 static_cast<std::size_t
>(std::ceil(far_cut * box_geo.
length()[1]) + 1.);
138 scxcache = calc_sc_cache<0>(particles, n_freq_x, u_x);
139 scycache = calc_sc_cache<1>(particles, n_freq_y, u_y);
140 return {n_freq_x, n_freq_y};
148 for (std::size_t i = 0; i < size; i++)
152static void copy_vec(
double *pdc_d,
double const *pdc_s, std::size_t size) {
153 for (std::size_t i = 0; i < size; i++)
157static void add_vec(
double *pdc_d,
double const *pdc_s1,
double const *pdc_s2,
159 for (std::size_t i = 0; i < size; i++)
160 pdc_d[i] = pdc_s1[i] + pdc_s2[i];
163static void addscale_vec(
double *pdc_d,
double scale,
double const *pdc_s1,
164 double const *pdc_s2, std::size_t size) {
165 for (std::size_t i = 0; i < size; i++)
166 pdc_d[i] = scale * pdc_s1[i] + pdc_s2[i];
169static void scale_vec(
double scale,
double *pdc, std::size_t size) {
170 for (std::size_t i = 0; i < size; i++)
174static double *
block(
double *p, std::size_t index, std::size_t size) {
175 return &p[index * size];
182 boost::mpi::all_reduce(
comm_cart, send_buf,
static_cast<int>(size),
gblcblk,
186void ElectrostaticLayerCorrection::check_gap(
Particle const &p)
const {
188 auto const z = p.
pos()[2];
191 <<
"region by " << ((z < 0.) ? z : z -
elc.box_h);
203void ElectrostaticLayerCorrection::add_dipole_force(
205 constexpr std::size_t size = 3;
211 auto const shift = box_geo.length_half()[2];
221 auto const q = p.
q();
222 auto const z = p.
pos()[2];
251 auto const field_induced =
gblcblk[1];
253 field_tot -= field_applied + field_induced;
257 p.
force()[2] -= field_tot * p.
q();
269double ElectrostaticLayerCorrection::dipole_energy(
271 constexpr std::size_t size = 7;
274 auto const lz = box_geo.length()[2];
277 auto const shift = box_geo.length_half()[2];
291 auto const q = p.
q();
292 auto const z = p.
pos()[2];
321 energy += 2. * pref *
357 double b(
double q,
double z)
const {
362 double t(
double q,
double z)
const {
368ElectrostaticLayerCorrection::z_energy(
ParticleRange const &particles)
const {
369 constexpr std::size_t size = 4;
371 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
376 auto const fac_delta = delta / (1. - delta);
380 auto const shift = box_geo.length_half()[2];
381 auto const lz = box_geo.length()[2];
387 auto const z = p.
pos()[2];
388 auto const q = p.
q();
406 auto const z = p.
pos()[2];
407 auto const q = p.
q();
440 return (
this_node == 0) ? -pref * energy : 0.;
443void ElectrostaticLayerCorrection::add_z_force(
445 constexpr std::size_t size = 1;
447 auto const xy_area_inv = box_geo.length_inv()[0] * box_geo.length_inv()[1];
452 auto const fac_delta = delta / (1. - delta);
459 auto const z = p.
pos()[2];
460 auto const q = p.
q();
469 auto const z = p.
pos()[2];
470 auto const q = p.
q();
507 constexpr std::size_t size = 4;
509 auto const pref_di = prefactor * 4. *
Utils::pi() * xy_area_inv;
510 auto const pref = -pref_di / expm1(omega * box_geo.
length()[2]);
511 double lclimgebot[4], lclimgetop[4], lclimge[4];
512 double fac_delta_mid_bot = 1., fac_delta_mid_top = 1., fac_delta = 1.;
516 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.
box_h));
527 auto const o = (index - 1) * particles.
size();
528 for (
auto const &p : particles) {
529 auto const z = p.
pos()[2];
530 auto const q = p.
q();
531 auto e = exp(omega * z);
548 lclimgebot[
POQESM] = sc_cache[o + ic].s / e;
549 lclimgebot[
POQESP] = sc_cache[o + ic].s * e;
550 lclimgebot[
POQECM] = sc_cache[o + ic].c / e;
551 lclimgebot[
POQECP] = sc_cache[o + ic].c * e;
556 exp(omega * (+z - 2. * elc.
box_h))) *
559 e = (exp(-omega * z) +
564 lclimge[
POQESP] += q * sc_cache[o + ic].s * e;
565 lclimge[
POQECP] += q * sc_cache[o + ic].c * e;
568 e = exp(omega * (2. * elc.
box_h - z));
572 lclimgetop[
POQESM] = sc_cache[o + ic].s / e;
573 lclimgetop[
POQESP] = sc_cache[o + ic].s * e;
574 lclimgetop[
POQECM] = sc_cache[o + ic].c / e;
575 lclimgetop[
POQECP] = sc_cache[o + ic].c * e;
580 exp(omega * (-z - 2. * elc.
box_h))) *
583 e = (exp(omega * (+z - 2. * elc.
box_h)) +
588 lclimge[
POQESM] += q * sc_cache[o + ic].s * e;
589 lclimge[
POQECM] += q * sc_cache[o + ic].c * e;
604 constexpr auto i =
static_cast<int>(axis);
605 constexpr std::size_t size = 4;
608 for (
auto &p : particles) {
623 constexpr std::size_t size = 4;
626 for (std::size_t ic = 0; ic < n_part; ic++) {
633 return energy / omega;
644 std::size_t index_q,
double omega,
647 assert(index_p >= 1);
648 assert(index_q >= 1);
649 constexpr std::size_t size = 8;
651 auto const pref_di = prefactor * 8 *
Utils::pi() * xy_area_inv;
652 auto const pref = -pref_di / expm1(omega * box_geo.
length()[2]);
653 double lclimgebot[8], lclimgetop[8], lclimge[8];
654 double fac_delta_mid_bot = 1, fac_delta_mid_top = 1, fac_delta = 1;
657 auto const fac_elc = 1. / (1. - delta * exp(-omega * 2. * elc.
box_h));
667 auto const ox = (index_p - 1) * particles.
size();
668 auto const oy = (index_q - 1) * particles.
size();
669 for (
auto const &p : particles) {
670 auto const z = p.
pos()[2];
671 auto const q = p.
q();
672 auto e = exp(omega * z);
714 exp(omega * (+z - 2. * elc.
box_h))) *
719 e = (exp(-omega * z) +
721 fac_delta_mid_bot * q;
731 e = exp(omega * (2. * elc.
box_h - z));
747 exp(omega * (-z - 2. * elc.
box_h))) *
752 e = (exp(omega * (+z - 2. * elc.
box_h)) +
754 fac_delta_mid_top * q;
773static void add_PQ_force(std::size_t index_p, std::size_t index_q,
double omega,
778 c_2pi * box_geo.
length_inv()[0] *
static_cast<double>(index_p) / omega;
780 c_2pi * box_geo.
length_inv()[1] *
static_cast<double>(index_q) / omega;
781 constexpr std::size_t size = 8;
784 for (
auto &p : particles) {
814static double PQ_energy(
double omega, std::size_t n_part) {
815 constexpr std::size_t size = 8;
818 for (std::size_t ic = 0; ic < n_part; ic++) {
828 return energy / omega;
832void ElectrostaticLayerCorrection::add_force(
837 auto const n_scxcache = std::get<0>(n_freqs);
838 auto const n_scycache = std::get<1>(n_freqs);
841 add_dipole_force(particles);
842 add_z_force(particles);
845 for (std::size_t p = 1;
846 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
849 auto const omega = c_2pi * box_geo.length_inv()[0] *
static_cast<double>(p);
850 setup_PoQ<PoQ::P>(
elc,
prefactor, p, omega, particles, box_geo);
852 add_PoQ_force<PoQ::P>(particles);
855 for (std::size_t q = 1;
856 box_geo.length_inv()[1] *
static_cast<double>(q - 1) <
elc.
far_cut &&
859 auto const omega = c_2pi * box_geo.length_inv()[1] *
static_cast<double>(q);
860 setup_PoQ<PoQ::Q>(
elc,
prefactor, q, omega, particles, box_geo);
862 add_PoQ_force<PoQ::Q>(particles);
865 for (std::size_t p = 1;
866 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
869 for (std::size_t q = 1;
870 Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p - 1)) +
872 static_cast<double>(q - 1)) <
878 sqrt(
Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p)) +
879 Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q)));
887double ElectrostaticLayerCorrection::calc_energy(
891 auto energy = dipole_energy(particles) + z_energy(particles);
893 auto const n_scxcache = std::get<0>(n_freqs);
894 auto const n_scycache = std::get<1>(n_freqs);
897 partblk.resize(n_localpart * 8);
900 for (std::size_t p = 1;
901 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
904 auto const omega = c_2pi * box_geo.length_inv()[0] *
static_cast<double>(p);
905 setup_PoQ<PoQ::P>(
elc,
prefactor, p, omega, particles, box_geo);
910 for (std::size_t q = 1;
911 box_geo.length_inv()[1] *
static_cast<double>(q - 1) <
elc.
far_cut &&
914 auto const omega = c_2pi * box_geo.length_inv()[1] *
static_cast<double>(q);
915 setup_PoQ<PoQ::Q>(
elc,
prefactor, q, omega, particles, box_geo);
920 for (std::size_t p = 1;
921 box_geo.length_inv()[0] *
static_cast<double>(p - 1) <
elc.
far_cut &&
924 for (std::size_t q = 1;
925 Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p - 1)) +
927 static_cast<double>(q - 1)) <
933 sqrt(
Utils::sqr(box_geo.length_inv()[0] *
static_cast<double>(p)) +
934 Utils::sqr(box_geo.length_inv()[1] *
static_cast<double>(q)));
944double ElectrostaticLayerCorrection::tune_far_cut()
const {
946 auto constexpr maximal_far_cut = 50.;
948 auto const box_l_x_inv = box_geo.length_inv()[0];
949 auto const box_l_y_inv = box_geo.length_inv()[1];
950 auto const min_inv_boxl = std::min(box_l_x_inv, box_l_y_inv);
951 auto const box_l_z = box_geo.length()[2];
956 auto tuned_far_cut = min_inv_boxl;
959 auto const pref = 2. *
Utils::pi() * tuned_far_cut;
960 auto const sum = pref + 2. * (box_l_x_inv + box_l_y_inv);
961 auto const den = -expm1(-pref * lz);
962 auto const num1 = exp(pref * (
elc.
box_h - lz));
963 auto const num2 = exp(-pref * (
elc.
box_h + lz));
969 tuned_far_cut += min_inv_boxl;
970 }
while (err >
elc.
maxPWerror and tuned_far_cut < maximal_far_cut);
971 if (tuned_far_cut >= maximal_far_cut) {
972 throw std::runtime_error(
"ELC tuning failed: maxPWerror too small");
974 return tuned_far_cut - min_inv_boxl;
982 return boost::mpi::all_reduce(
comm_cart, local_q, std::plus<>());
985void ElectrostaticLayerCorrection::sanity_checks_periodicity()
const {
987 if (!box_geo.periodic(0) || !box_geo.periodic(1) || !box_geo.periodic(2)) {
988 throw std::runtime_error(
"ELC: requires periodicity (True, True, True)");
992void ElectrostaticLayerCorrection::sanity_checks_dielectric_contrasts()
const {
994 auto const &cell_structure = *
get_system().cell_structure;
997 if (total_charge >= precision_threshold) {
1001 throw std::runtime_error(
"ELC does not currently support non-neutral "
1002 "systems with a dielectric contrast.");
1006 throw std::runtime_error(
"ELC does not work for non-neutral systems and "
1007 "non-metallic dielectric contrast.");
1012void ElectrostaticLayerCorrection::adapt_solver() {
1014 [
this](
auto &solver) {
1021void ElectrostaticLayerCorrection::recalc_box_h() {
1025 if (new_box_h < 0.) {
1026 throw std::runtime_error(
"ELC gap size (" + std::to_string(
elc.
gap_size) +
1027 ") larger than box length in z-direction (" +
1028 std::to_string(box_z) +
")");
1033void ElectrostaticLayerCorrection::recalc_space_layer() {
1035 auto const p3m_r_cut = std::visit(
1036 [](
auto &solver) {
return solver->p3m.params.r_cut; },
base_solver);
1044 auto const half_box_h =
elc.
box_h / 2.;
1045 auto const max_space_layer = std::min(free_space, half_box_h);
1047 if (max_space_layer <= 0.) {
1048 throw std::runtime_error(
"P3M real-space cutoff too large for ELC w/ "
1049 "dielectric contrast");
1058 bool neutralize,
double delta_top,
double delta_bot,
1059 bool with_const_pot,
double potential_diff)
1061 far_cut2{-1.}, far_calculated{far_cut == -1.},
1062 dielectric_contrast_on{delta_top != 0. or delta_bot != 0.},
1063 const_pot{with_const_pot and dielectric_contrast_on},
1064 neutralize{neutralize and !dielectric_contrast_on},
1065 delta_mid_top{std::clamp(delta_top, -1., +1.)},
1066 delta_mid_bot{std::clamp(delta_bot, -1., +1.)},
1067 pot_diff{(with_const_pot) ? potential_diff : 0.},
1070 space_layer{(dielectric_contrast_on) ? gap_size / 3. : 0.},
1071 space_box{gap_size - ((dielectric_contrast_on) ? 2. * space_layer : 0.)} {
1075 throw std::domain_error(
"Parameter 'far_cut' must be > 0");
1078 throw std::domain_error(
"Parameter 'maxPWerror' must be > 0");
1081 throw std::domain_error(
"Parameter 'gap_size' must be > 0");
1083 if (potential_diff != 0. and not with_const_pot) {
1084 throw std::invalid_argument(
1085 "Parameter 'const_pot' must be True when 'pot_diff' is non-zero");
1087 if (delta_top < -delta_range or delta_top > delta_range) {
1088 throw std::domain_error(
1089 "Parameter 'delta_mid_top' must be >= -1 and <= +1");
1091 if (delta_bot < -delta_range or delta_bot > delta_range) {
1092 throw std::domain_error(
1093 "Parameter 'delta_mid_bot' must be >= -1 and <= +1");
1100 throw std::domain_error(
"ELC with two parallel metallic boundaries "
1101 "requires the const_pot option");
1107 : elc{parameters}, base_solver{solver} {
1123template <ChargeProtocol protocol,
typename combined_ranges>
1125 combined_ranges
const &p_q_pos_range) {
1133 for (
auto zipped : p_q_pos_range) {
1134 auto const p_q = boost::get<0>(zipped);
1135 auto const &p_pos = boost::get<1>(zipped);
1149template <ChargeProtocol protocol,
typename combined_range>
1151 combined_range
const &p_q_pos_range) {
1154 for (
auto zipped : p_q_pos_range) {
1155 auto const p_q = boost::get<0>(zipped);
1156 auto const &p_pos = boost::get<1>(zipped);
1158 auto const p_z = p_pos[2];
1164 node_sums[2] += p_q;
1184 auto const tot_sums =
1185 boost::mpi::all_reduce(
comm_cart, node_sums, std::plus<>());
1186 solver.
p3m.
sum_qpart =
static_cast<int>(tot_sums[0] + 0.1);
1193 auto const energy = std::visit(
1194 [
this, &particles](
auto const &solver_ptr) {
1195 auto &solver = *solver_ptr;
1200 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1203 solver.charge_assign(particles);
1206 return solver.long_range_energy(particles);
1210 energy += 0.5 * solver.long_range_energy(particles);
1215 charge_assign<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1216 modify_p3m_sums<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1217 energy += 0.5 * solver.long_range_energy(particles);
1220 charge_assign<ChargeProtocol::IMAGE>(
elc, solver, p_q_pos_range);
1221 modify_p3m_sums<ChargeProtocol::IMAGE>(
elc, solver, p_q_pos_range);
1222 energy -= 0.5 * solver.long_range_energy(particles);
1225 modify_p3m_sums<ChargeProtocol::REAL>(
elc, solver, p_q_pos_range);
1230 return energy + calc_energy(particles);
1236 [
this, &particles](
auto const &solver_ptr) {
1237 auto &solver = *solver_ptr;
1240 auto p_q_pos_range = boost::combine(p_q_range, p_pos_range);
1243 modify_p3m_sums<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1244 charge_assign<ChargeProtocol::BOTH>(
elc, solver, p_q_pos_range);
1247 solver.charge_assign(particles);
1249 solver.add_long_range_forces(particles);
1251 modify_p3m_sums<ChargeProtocol::REAL>(
elc, solver, p_q_pos_range);
1255 add_force(particles);
Vector implementation and trait types for boost qvm interoperability.
__shared__ int pos[MAXDEPTH *THREADS5/WARPSIZE]
Utils::Vector3d const & length() const
Box length.
Utils::Vector3d const & length_inv() const
Inverse box length.
void set_prefactor(double new_prefactor)
double prefactor
Electrostatics prefactor.
base_type::size_type size() const
void reset(int cao)
Reset the cache.
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
boost::mpi::communicator comm_cart
The communicator.
int this_node
The number of this node.
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
T image_sum(InputIterator begin, InputIterator end, InputIterator it, bool with_replicas, Utils::Vector3i const &ncut, BoxGeometry const &box_geo, T init, F f)
Sum over all pairs with periodic images.
static void addscale_vec(double *pdc_d, double scale, double const *pdc_s1, double const *pdc_s2, std::size_t size)
static std::pair< std::size_t, std::size_t > prepare_sc_cache(ParticleRange const &particles, BoxGeometry const &box_geo, double far_cut)
static void add_PQ_force(std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static auto calc_total_charge(CellStructure const &cell_structure)
static std::vector< double > partblk
temporary buffers for product decomposition
static void clear_vec(double *pdc, std::size_t size)
static double * block(double *p, std::size_t index, std::size_t size)
void setup_PoQ(elc_data const &elc, double prefactor, std::size_t index, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static void distribute(std::size_t size)
void modify_p3m_sums(elc_data const &elc, CoulombP3M &solver, combined_range const &p_q_pos_range)
static double PoQ_energy(double omega, std::size_t n_part)
static std::vector< SCCache > scxcache
Cached sin/cos values along the x-axis and y-axis.
static void p3m_assign_image_charge(elc_data const &elc, CoulombP3M &p3m, double q, Utils::Vector3d const &pos)
static double PQ_energy(double omega, std::size_t n_part)
void charge_assign(elc_data const &elc, CoulombP3M &solver, combined_ranges const &p_q_pos_range)
static std::vector< SCCache > scycache
static void setup_PQ(elc_data const &elc, double prefactor, std::size_t index_p, std::size_t index_q, double omega, ParticleRange const &particles, BoxGeometry const &box_geo)
static void copy_vec(double *pdc_d, double const *pdc_s, std::size_t size)
static void add_vec(double *pdc_d, double const *pdc_s1, double const *pdc_s2, std::size_t size)
static void scale_vec(double scale, double *pdc, std::size_t size)
void add_PoQ_force(ParticleRange const &particles)
ChargeProtocol
ELC charge sum/assign protocol: real charges, image charges, or both.
static std::vector< SCCache > calc_sc_cache(ParticleRange const &particles, std::size_t n_freq, double u)
Calculate cached sin/cos values for one direction.
static double gblcblk[8]
collected data from the other cells
PoQ
ELC axes (x and y directions)
ELC algorithm for long-range Coulomb interactions.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
#define runtimeErrorMsg()
Common parts of the MMM family of methods for the electrostatic interaction: MMM1D and ELC.
__constant__ float maxPWerror[1]
ParticleRange particles(Utils::Span< Cell *const > cells)
auto charge_range(ParticleRange const &particles)
auto pos_range(ParticleRange const &particles)
DEVICE_QUALIFIER constexpr T pi()
Ratio of diameter and circumference of a circle.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Vector< T, N > sqrt(Vector< T, N > const &a)
P3M algorithm for long-range Coulomb interaction.
P3M electrostatics on GPU.
Describes a cell structure / cell system.
ParticleRange local_particles() const
void assign_charge(double q, Utils::Vector3d const &real_pos, p3m_interpolation_cache &inter_weights)
Assign a single charge into the current charge grid.
p3m_data_struct p3m
P3M parameters.
BaseSolver base_solver
Electrostatics solver that is adapted.
ElectrostaticLayerCorrection(elc_data &¶meters, BaseSolver &&solver)
std::variant< std::shared_ptr< CoulombP3MGPU >, std::shared_ptr< CoulombP3M > > BaseSolver
void add_long_range_forces(ParticleRange const &particles) const
double long_range_energy(ParticleRange const &particles) const
double b(double q, double z) const
Image sum from the bottom layer.
ImageSum(double delta, double shift, double lz)
double t(double q, double z) const
Image sum from the top layer.
int size
number of local mesh points.
int cao
charge assignment order ([0,7]).
Struct holding all information for one particle.
auto const & force() const
structure for caching sin and cos values
Parameters for the ELC method.
double dielectric_layers_self_energy(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
self energies of top and bottom layers with their virtual images
double maxPWerror
Maximal allowed pairwise error for the potential and force.
double pot_diff
Constant potential difference.
double box_h
Up to where particles can be found.
bool dielectric_contrast_on
Flag whether there is any dielectric contrast in the system.
elc_data(double maxPWerror, double gap_size, double far_cut, bool neutralize, double delta_top, double delta_bot, bool const_pot, double pot_diff)
double space_box
The space that is finally left.
bool neutralize
Flag whether the box is neutralized by a homogeneous background.
double far_cut
Cutoff of the exponential sum.
double space_layer
Layer around the dielectric contrast in which we trick around.
bool far_calculated
Flag whether far_cut was set by the user, or calculated by ESPResSo.
double gap_size
Size of the empty gap.
double delta_mid_bot
dielectric contrast in the lower part of the simulation cell.
void dielectric_layers_self_forces(CoulombP3M const &p3m, BoxGeometry const &box_geo, ParticleRange const &particles) const
forces of particles in border layers with themselves
bool const_pot
Flag whether a constant potential difference is applied.
double far_cut2
Squared value of far_cut.
double delta_mid_top
dielectric contrast in the upper part of the simulation cell.
P3MLocalMesh local_mesh
local mesh.
double sum_q2
Sum of square of charges (only on head node).
fft_vector< double > rs_mesh
real space mesh (local) for CA/FFT.
double square_sum_q
square of sum of charges (only on head node).
int sum_qpart
number of charged particles (only on head node).
p3m_interpolation_cache inter_weights