57#include "communication.hpp"
62#include "system/System.hpp"
70#include <boost/mpi/collectives/all_reduce.hpp>
71#include <boost/mpi/collectives/broadcast.hpp>
72#include <boost/mpi/collectives/reduce.hpp>
73#include <boost/mpi/communicator.hpp>
74#include <boost/range/combine.hpp>
75#include <boost/range/numeric.hpp>
83#include <initializer_list>
94#error "The FFTW3 library shouldn't be visible in this translation unit"
97template <
typename FloatType, Arch Architecture>
103 for (
auto const &p : get_system().
cell_structure->local_particles()) {
111 boost::mpi::all_reduce(
comm_cart, local_n, p3m.sum_qpart, std::plus<>());
112 boost::mpi::all_reduce(
comm_cart, local_q2, p3m.sum_q2, std::plus<>());
113 boost::mpi::all_reduce(
comm_cart, local_q, p3m.square_sum_q, std::plus<>());
114 p3m.square_sum_q =
Utils::sqr(p3m.square_sum_q);
126template <
typename FloatType, Arch Architecture>
128 auto const [KX, KY, KZ] = p3m.fft->get_permutations();
129 p3m.g_force = grid_influence_function<FloatType, 1>(
130 p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ,
131 get_system().box_geo->length_inv());
137template <
typename FloatType, Arch Architecture>
139 auto const [KX, KY, KZ] = p3m.fft->get_permutations();
140 p3m.g_energy = grid_influence_function<FloatType, 0>(
141 p3m.params, p3m.mesh.start, p3m.mesh.stop, KX, KY, KZ,
142 get_system().box_geo->length_inv());
153 auto const factor1 =
Utils::sqr(std::numbers::pi * alpha_L_i);
161 mesh_start, mesh_stop, indices,
163 auto const norm_sq = nm.norm2();
164 auto const ex = exp(-factor1 * norm_sq);
167 alias2 += energy * ex * (shift * nm) / norm_sq;
169 [&](
unsigned dim,
int n) {
170 nm[dim] = shift[dim] + n * mesh[dim];
174 return std::make_pair(alias1, alias2);
188 double sum_q2,
double alpha_L,
191 return (2. * pref * sum_q2 * exp(-
Utils::sqr(r_cut_iL * alpha_L))) /
192 sqrt(
static_cast<double>(n_c_part) * r_cut_iL * box_l[0] * volume);
209 int cao,
int n_c_part,
double sum_q2,
214 auto const alpha_L_i = 1. / alpha_L;
215 auto const mesh_stop = mesh / 2;
216 auto const mesh_start = -mesh_stop;
222 mesh_start, mesh_stop, indices,
224 if ((indices[0] != 0) or (indices[1] != 0) or (indices[2] != 0)) {
225 auto const n2 = indices.norm2();
227 auto const [alias1, alias2] =
229 auto const d = alias1 -
Utils::sqr(alias2 / cs) / n2;
237 [&values, &mesh_i, cotangent_sum](
unsigned dim,
int n) {
238 values[dim] = cotangent_sum(n, mesh_i[dim]);
241 return 2. * pref * sum_q2 * sqrt(he_q /
static_cast<double>(n_c_part)) /
242 (box_l[1] * box_l[2]);
245template <
typename FloatType, Arch Architecture>
248 assert(p3m.params.cao >= 1 and p3m.params.cao <= 7);
249 assert(p3m.params.alpha > 0.);
251 auto const &system = get_system();
252 auto const &box_geo = *system.box_geo;
253 auto const &local_geo = *system.local_geo;
254 auto const skin = system.cell_structure->get_verlet_skin();
256 p3m.params.cao3 = Utils::int_pow<3>(p3m.params.cao);
257 p3m.params.recalc_a_ai_cao_cut(box_geo.length());
261 auto const &solver = system.coulomb.impl->solver;
262 double elc_layer = 0.;
263 if (
auto actor = get_actor_by_type<ElectrostaticLayerCorrection>(solver)) {
264 elc_layer = actor->elc.space_layer;
268 p3m.local_mesh.calc_local_ca_mesh(p3m.params, local_geo, skin, elc_layer);
269 p3m.fft_buffers->init_halo();
270 p3m.fft->init(p3m.params);
271 p3m.mesh.ks_pnum = p3m.fft->get_ks_pnum();
272 p3m.fft_buffers->init_meshes(p3m.fft->get_ca_mesh_size());
273 p3m.update_mesh_views();
274 p3m.calc_differential_operator();
279 count_charged_particles();
287 typename std::remove_reference_t<
decltype(p3m)>::value_type;
289 p3m.mesh.rs_scalar[ind] += value_type(w * q);
295 auto const w = p3m_calculate_interpolation_weights<cao>(
296 real_pos, p3m.params.ai, p3m.local_mesh);
297 inter_weights.
store(w);
298 this->operator()(p3m, q, real_pos, w);
302 auto const w = p3m_calculate_interpolation_weights<cao>(
303 real_pos, p3m.params.ai, p3m.local_mesh);
304 this->operator()(p3m, q, real_pos, w);
307 template <
typename combined_ranges>
308 void operator()(
auto &p3m, combined_ranges
const &p_q_pos_range) {
309 for (
auto zipped : p_q_pos_range) {
310 auto const p_q = boost::get<0>(zipped);
311 auto const &p_pos = boost::get<1>(zipped);
313 this->operator()(p3m, p_q, p_pos, p3m.inter_weights);
320template <
typename FloatType, Arch Architecture>
323 prepare_fft_mesh(
true);
328 Utils::integral_parameter<int, AssignCharge, 1, 7>(
329 p3m.params.cao, p3m, boost::combine(p_q_range, p_pos_range));
332template <
typename FloatType, Arch Architecture>
336 Utils::integral_parameter<int, AssignCharge, 1, 7>(p3m.params.cao, p3m, q,
339 Utils::integral_parameter<int, AssignCharge, 1, 7>(
340 p3m.params.cao, p3m, q, real_pos, p3m.inter_weights);
345 template <
typename combined_ranges>
347 combined_ranges
const &p_q_force_range)
const {
349 assert(cao == p3m.inter_weights.cao());
352 auto p_index = std::size_t{0ul};
354 for (
auto zipped : p_q_force_range) {
355 auto p_q = boost::get<0>(zipped);
356 auto &p_force = boost::get<1>(zipped);
358 auto const pref = p_q * force_prefac;
359 auto const w = p3m.inter_weights.template load<cao>(p_index);
363 force[0u] += w * double(p3m.mesh.rs_fields[0u][ind]);
364 force[1u] += w * double(p3m.mesh.rs_fields[1u][ind]);
365 force[2u] += w * double(p3m.mesh.rs_fields[2u][ind]);
368 p_force -= pref * force;
375template <
typename combined_ranges>
377 combined_ranges
const &p_q_unfolded_pos_range) {
378 auto const local_dip =
381 auto const p_q = boost::get<0>(q_pos);
382 auto const &p_unfolded_pos = boost::get<1>(q_pos);
383 return dip + p_q * p_unfolded_pos;
385 return boost::mpi::all_reduce(comm, local_dip, std::plus<>());
393template <
typename FloatType, Arch Architecture>
396 auto const &box_geo = *get_system().
box_geo;
399 if (p3m.sum_q2 > 0.) {
401 p3m.fft_buffers->perform_scalar_halo_gather();
402 p3m.fft->forward_fft(p3m.fft_buffers->get_scalar_mesh());
403 p3m.update_mesh_views();
406 auto const &mesh_stop = p3m.mesh.size;
407 auto const &offset = p3m.mesh.start;
408 auto const half_alpha_inv_sq =
Utils::sqr(1. / 2. / p3m.params.alpha);
409 auto const wavevector = (2. * std::numbers::pi) * box_geo.length_inv();
410 auto const [KX, KY, KZ] = p3m.fft->get_permutations();
412 auto index = std::size_t(0u);
415 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
416 auto const shift = indices + offset;
417 auto const kx = p3m.d_op[0u][shift[KX]] * wavevector[0u];
418 auto const ky = p3m.d_op[1u][shift[KY]] * wavevector[1u];
419 auto const kz = p3m.d_op[2u][shift[KZ]] * wavevector[2u];
423 auto const node_k_space_energy =
424 double(p3m.g_energy[index] *
425 (
Utils::sqr(p3m.mesh.rs_scalar[2u * index + 0u]) +
426 Utils::sqr(p3m.mesh.rs_scalar[2u * index + 1u])));
427 auto const vterm = -2. * (1. / norm_sq + half_alpha_inv_sq);
428 auto const pref = node_k_space_energy * vterm;
429 node_k_space_pressure_tensor[0u] += pref * kx * kx;
430 node_k_space_pressure_tensor[1u] += pref * kx * ky;
431 node_k_space_pressure_tensor[2u] += pref * kx * kz;
432 node_k_space_pressure_tensor[4u] += pref * ky * ky;
433 node_k_space_pressure_tensor[5u] += pref * ky * kz;
434 node_k_space_pressure_tensor[8u] += pref * kz * kz;
435 diagonal += node_k_space_energy;
440 node_k_space_pressure_tensor[0u] += diagonal;
441 node_k_space_pressure_tensor[4u] += diagonal;
442 node_k_space_pressure_tensor[8u] += diagonal;
443 node_k_space_pressure_tensor[3u] = node_k_space_pressure_tensor[1u];
444 node_k_space_pressure_tensor[6u] = node_k_space_pressure_tensor[2u];
445 node_k_space_pressure_tensor[7u] = node_k_space_pressure_tensor[5u];
448 return node_k_space_pressure_tensor * prefactor / (2. * box_geo.volume());
451template <
typename FloatType, Arch Architecture>
453 bool force_flag,
bool energy_flag,
ParticleRange const &particles) {
454 auto const &system = get_system();
455 auto const &box_geo = *system.box_geo;
457 auto const npt_flag =
460 auto constexpr npt_flag =
false;
463 if (p3m.sum_q2 > 0.) {
464 if (not has_actor_of_type<ElectrostaticLayerCorrection>(
465 system.coulomb.impl->solver)) {
468 p3m.fft_buffers->perform_scalar_halo_gather();
469 p3m.fft->forward_fft(p3m.fft_buffers->get_scalar_mesh());
470 p3m.update_mesh_views();
475 auto p_unfolded_pos_range =
481 auto const box_dipole =
484 comm_cart, boost::combine(p_q_range, p_unfolded_pos_range)))
486 auto const volume = box_geo.volume();
488 4. * std::numbers::pi / volume / (2. * p3m.params.epsilon + 1.);
494 auto const &mesh_stop = p3m.mesh.size;
495 auto const &offset = p3m.mesh.start;
497 box_geo.length_inv());
499 auto index = std::size_t(0u);
503 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
505 std::complex<FloatType>(p3m.mesh.rs_scalar[2u * index + 0u],
506 p3m.mesh.rs_scalar[2u * index + 1u]);
507 auto const phi_hat = p3m.g_force[index] * rho_hat;
509 for (
int d = 0; d < 3; d++) {
511 int d_rs = (d + p3m.mesh.ks_pnum) % 3;
513 auto const k = FloatType(p3m.d_op[d_rs][indices[d] + offset[d]]) *
517 p3m.mesh.rs_fields[d_rs][2u * index + 0u] = -k * phi_hat.imag();
518 p3m.mesh.rs_fields[d_rs][2u * index + 1u] = +k * phi_hat.real();
524 auto const check_residuals =
525 not p3m.params.tuning and check_complex_residuals;
526 p3m.fft->check_complex_residuals = check_residuals;
527 for (
auto &rs_mesh : p3m.fft_buffers->get_vector_mesh()) {
528 p3m.fft->backward_fft(rs_mesh);
530 p3m.fft_buffers->perform_vector_halo_spread();
531 p3m.fft->check_complex_residuals =
false;
533 auto const force_prefac = prefactor / volume;
534 Utils::integral_parameter<int, AssignForces, 1, 7>(
535 p3m.params.cao, p3m, force_prefac,
536 boost::combine(p_q_range, p_force_range));
541 auto const dm = prefactor * pref * box_dipole.value();
542 for (
auto zipped : boost::combine(p_q_range, p_force_range)) {
543 auto p_q = boost::get<0>(zipped);
544 auto &p_force = boost::get<1>(zipped);
551 if (energy_flag or npt_flag) {
552 auto node_energy = 0.;
554 for (
int i = 0; i < mesh_length; i++) {
558 double(p3m.g_energy[i] * (
Utils::sqr(p3m.mesh.rs_scalar[2 * i + 0]) +
561 node_energy /= 2. * volume;
564 boost::mpi::reduce(
comm_cart, node_energy, energy, std::plus<>(), 0);
568 energy -= p3m.sum_q2 * p3m.params.alpha * std::numbers::inv_sqrtpi;
571 energy -= p3m.square_sum_q * std::numbers::pi /
576 energy += pref * box_dipole.value().norm2();
585 if (not energy_flag) {
594template <
typename FloatType, Arch Architecture>
597 double m_mesh_density_min = -1., m_mesh_density_max = -1.;
599 bool m_tune_mesh =
false;
606 double prefactor,
int timings)
614 auto const on_gpu = Architecture ==
Arch::GPU;
616 auto const on_gpu =
false;
618 m_logger = std::make_unique<TuningLogger>(
619 verbose and
this_node == 0, (on_gpu) ?
"CoulombP3MGPU" :
"CoulombP3M",
626 std::optional<std::string>
629 if (
auto actor = get_actor_by_type<ElectrostaticLayerCorrection>(solver)) {
630 return actor->veto_r_cut(r_cut);
635 std::tuple<double, double, double, double>
637 double r_cut_iL)
const override {
640 double alpha_L, rs_err, ks_err;
644 p3m.
sum_q2, 0., box_geo.length());
648 alpha_L = sqrt(log(std::numbers::sqrt2 * rs_err / p3m.
params.
accuracy)) /
659 p3m.
sum_q2, alpha_L, box_geo.length());
661 if constexpr (Architecture ==
Arch::GPU) {
665 p3m.
sum_q2, alpha_L, box_geo.length().data());
667 boost::mpi::broadcast(
comm_cart, ks_err, 0);
671 p3m.
sum_q2, alpha_L, box_geo.length());
678 auto const mesh_density =
679 static_cast<double>(p3m.
params.
mesh[0]) * box_geo.length_inv()[0];
683 auto const normalized_box_dim = std::cbrt(box_geo.volume());
684 auto const max_npart_per_dim = 512.;
688 auto const min_npart_per_dim = std::min(
689 max_npart_per_dim, std::cbrt(
static_cast<double>(p3m.
sum_qpart)));
690 m_mesh_density_min = min_npart_per_dim / normalized_box_dim;
691 m_mesh_density_max = max_npart_per_dim / normalized_box_dim;
694 m_mesh_density_min = m_mesh_density_max = mesh_density;
698 for (
auto i : {1, 2}) {
700 static_cast<int>(std::round(mesh_density * box_geo.length()[i]));
714 auto mesh_density = m_mesh_density_min;
715 while (mesh_density <= m_mesh_density_max) {
718 for (
auto i : {0, 1, 2}) {
719 trial_params.
mesh[i] =
720 static_cast<int>(std::round(box_geo.length()[i] * mesh_density));
722 trial_params.mesh[i] += trial_params.mesh[i] % 2;
729 auto const trial_time =
730 get_m_time(trial_params.mesh, trial_params.cao, trial_params.r_cut_iL,
731 trial_params.alpha_L, trial_params.accuracy);
733 if (trial_time >= 0.) {
736 if (has_actor_of_type<CoulombP3M>(solver)) {
740 if (trial_time < time_best) {
743 tuned_params = trial_params;
744 time_best = tuned_params.time = trial_time;
757template <
typename FloatType, Arch Architecture>
759 auto &system = get_system();
760 auto const &box_geo = *system.box_geo;
761 if (p3m.params.alpha_L == 0. and p3m.params.alpha != 0.) {
762 p3m.params.alpha_L = p3m.params.alpha * box_geo.length()[0];
764 if (p3m.params.r_cut_iL == 0. and p3m.params.r_cut != 0.) {
765 p3m.params.r_cut_iL = p3m.params.r_cut * box_geo.length_inv()[0];
767 if (not is_tuned()) {
768 count_charged_particles();
769 if (p3m.sum_qpart == 0) {
770 throw std::runtime_error(
771 "CoulombP3M: no charged particles in the system");
775 system, p3m, prefactor, tune_timings);
784 system.on_coulomb_change();
786 p3m.params.tuning =
false;
795 auto const &box_geo = *system.box_geo;
796 auto const &local_geo = *system.local_geo;
797 for (
auto i = 0u; i < 3u; i++) {
800 std::stringstream msg;
802 <<
" is larger than half of box dimension " << box_geo.length()[i];
803 throw std::runtime_error(msg.str());
806 std::stringstream msg;
808 <<
" is larger than local box dimension " << local_geo.length()[i];
809 throw std::runtime_error(msg.str());
814 if ((box_geo.length()[0] != box_geo.length()[1]) or
815 (box_geo.length()[1] != box_geo.length()[2]) or
818 throw std::runtime_error(
819 "CoulombP3M: non-metallic epsilon requires cubic box");
826 if (!box_geo.periodic(0) or !box_geo.periodic(1) or !box_geo.periodic(2)) {
827 throw std::runtime_error(
828 "CoulombP3M: requires periodicity (True, True, True)");
833 auto const &local_geo = *
get_system().local_geo;
836 throw std::runtime_error(
837 "CoulombP3M: requires the regular or hybrid decomposition cell system");
841 throw std::runtime_error(
842 "CoulombP3M: does not work with the hybrid decomposition cell system, "
843 "if using more than one MPI node");
849 if (node_grid[0] < node_grid[1] || node_grid[1] < node_grid[2]) {
850 throw std::runtime_error(
851 "CoulombP3M: node grid must be sorted, largest first");
855template <
typename FloatType, Arch Architecture>
857 auto const &box_geo = *get_system().
box_geo;
858 p3m.params.r_cut = p3m.params.r_cut_iL * box_geo.length()[0];
859 p3m.params.alpha = p3m.params.alpha_L * box_geo.length_inv()[0];
860 p3m.params.recalc_a_ai_cao_cut(box_geo.length());
861 p3m.local_mesh.recalc_ld_pos(p3m.params);
862 sanity_checks_boxl();
863 calc_influence_function_force();
864 calc_influence_function_energy();
868template <
typename FloatType, Arch Architecture>
871 if constexpr (Architecture ==
Arch::GPU) {
874 auto const energy = long_range_energy(particles);
878 static_cast<void>(particles);
881 auto &gpu = get_system().
gpu;
893template <
typename FloatType, Arch Architecture>
895 if constexpr (Architecture ==
Arch::GPU) {
896 auto &system = get_system();
897 if (has_actor_of_type<ElectrostaticLayerCorrection>(
898 system.coulomb.impl->solver)) {
901 p3m_gpu_init(m_gpu_data, p3m.params.cao, p3m.params.mesh, p3m.params.alpha,
902 system.box_geo->length(), system.gpu.n_particles());
906template <
typename FloatType, Arch Architecture>
908 if constexpr (Architecture ==
Arch::GPU) {
909 auto &gpu_particle_data = get_system().
gpu;
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
Vector implementation and trait types for boost qvm interoperability.
TuningAlgorithm::Parameters get_time() override
P3MParameters & get_params() override
std::tuple< double, double, double, double > calculate_accuracy(Utils::Vector3i const &mesh, int cao, double r_cut_iL) const override
void setup_logger(bool verbose) override
void on_solver_change() const override
std::optional< std::string > layer_correction_veto_r_cut(double r_cut) const override
CoulombTuningAlgorithm(System::System &system, auto &input_p3m, double prefactor, int timings)
void determine_mesh_limits() override
void enable_property(std::size_t property)
std::shared_ptr< Propagation > propagation
std::shared_ptr< CellStructure > cell_structure
std::shared_ptr< BoxGeometry > box_geo
Tuning algorithm for P3M.
double get_m_time(Utils::Vector3i const &mesh, int &tuned_cao, double &tuned_r_cut_iL, double &tuned_alpha_L, double &tuned_accuracy)
Get the optimal alpha and the corresponding computation time for a fixed mesh.
static auto constexpr time_sentinel
Value for invalid time measurements.
static auto constexpr max_n_consecutive_trials
Maximal number of consecutive trials that don't improve runtime.
System::System & m_system
void determine_cao_limits(int initial_cao)
Determine a sensible range for the charge assignment order.
void determine_r_cut_limits()
Determine a sensible range for the real-space cutoff.
std::unique_ptr< TuningLogger > m_logger
static auto constexpr time_granularity
Granularity of the time measurement (milliseconds).
static DEVICE_QUALIFIER constexpr Vector< T, N > broadcast(typename Base::value_type const &value)
Create a vector that has all entries set to the same value.
Cache for interpolation weights.
void store(const InterpolationWeights< cao > &w)
Push back weights for one point.
Communicator communicator
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.
#define P3M_BRILLOUIN
P3M: Number of Brillouin zones taken into account in the calculation of the optimal influence functio...
void charge_assign(elc_data const &elc, CoulombP3M &solver, combined_ranges const &p_q_pos_range)
ELC algorithm for long-range Coulomb interactions.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
Routines, row decomposition, data structures and communication for the 3D-FFT.
and std::invocable< Projector, unsigned, int > void for_each_3d(detail::IndexVectorConcept auto &&start, detail::IndexVectorConcept auto &&stop, detail::IndexVectorConcept auto &&counters, Kernel &&kernel, Projector &&projector=detail::noop_projector)
Repeat an operation on every element of a 3D grid.
void p3m_interpolate(P3MLocalMesh const &local_mesh, InterpolationWeights< cao > const &weights, Kernel kernel)
P3M grid interpolation.
auto charge_range(ParticleRange const &particles)
auto pos_range(ParticleRange const &particles)
auto force_range(ParticleRange const &particles)
auto unfolded_pos_range(ParticleRange const &particles, BoxGeometry const &box)
T product(Vector< T, N > const &v)
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
auto get_analytic_cotangent_sum_kernel(int cao)
void npt_add_virial_contribution(double energy)
Exports for the NpT code.
auto constexpr P3M_EPSILON_METALLIC
This value indicates metallic boundary conditions.
static auto p3m_tune_aliasing_sums(Utils::Vector3i const &shift, Utils::Vector3i const &mesh, Utils::Vector3d const &mesh_i, int cao, double alpha_L_i)
Aliasing sum used by p3m_k_space_error.
static double p3m_k_space_error(double pref, Utils::Vector3i const &mesh, int cao, int n_c_part, double sum_q2, double alpha_L, Utils::Vector3d const &box_l)
Calculate the analytic expression of the error estimate for the P3M method in (eq.
static double p3m_real_space_error(double pref, double r_cut_iL, int n_c_part, double sum_q2, double alpha_L, Utils::Vector3d const &box_l)
Calculate the real space contribution to the rms error in the force (as described by Kolafa and Perra...
static auto calc_dipole_moment(boost::mpi::communicator const &comm, combined_ranges const &p_q_unfolded_pos_range)
P3M algorithm for long-range Coulomb interaction.
void p3m_gpu_add_farfield_force(P3MGpuParams &data, GpuParticleData &gpu, double prefactor, std::size_t n_part)
The long-range part of the P3M algorithm.
void p3m_gpu_init(std::shared_ptr< P3MGpuParams > &data, int cao, Utils::Vector3i const &mesh, double alpha, Utils::Vector3d const &box_l, std::size_t n_part)
Initialize the internal data structure of the P3M GPU.
P3M electrostatics on GPU.
double p3m_k_space_error_gpu(double prefactor, const int *mesh, int cao, int npart, double sum_q2, double alpha_L, const double *box)
void operator()(auto &p3m, double force_prefac, combined_ranges const &p_q_force_range) const
Utils::Vector3i node_grid
void charge_assign(ParticleRange const &particles) override
void assign_charge(double q, Utils::Vector3d const &real_pos, bool skip_cache) override
void calc_influence_function_energy() override
Calculate the influence function optimized for the energy and the self energy correction.
void scaleby_box_l() override
void calc_influence_function_force() override
Calculate the optimal influence function of .
Utils::Vector9d long_range_pressure(ParticleRange const &particles) override
void count_charged_particles() override
double long_range_kernel(bool force_flag, bool energy_flag, ParticleRange const &particles)
Compute the k-space part of forces and energies.
void add_long_range_forces_gpu(ParticleRange const &particles)
void sanity_checks_periodicity() const
void sanity_checks_boxl() const
Checks for correctness of the k-space cutoff.
void sanity_checks_cell_structure() const
P3MParameters const & p3m_params
void sanity_checks_node_grid() const
std::unique_ptr< Implementation > impl
Pointer-to-implementation.
static constexpr std::size_t force
static constexpr std::size_t pos
static constexpr std::size_t q
Interpolation weights for one point.
Structure to hold P3M parameters and some dependent variables.
Utils::Vector3d cao_cut
cutoff for charge assignment.
double accuracy
accuracy of the actual parameter set.
Utils::Vector3i mesh
number of mesh points per coordinate direction (>0).
double epsilon
epsilon of the "surrounding dielectric".
DEVICE_QUALIFIER constexpr pointer data() noexcept
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos, p3m_interpolation_cache &inter_weights)
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos, InterpolationWeights< cao > const &w)
void operator()(auto &p3m, combined_ranges const &p_q_pos_range)
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos)
double sum_q2
Sum of square of charges (only on head node).
int sum_qpart
number of charged particles (only on head node).
P3MParameters params
P3M base parameters.