ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
p3m.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2024 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22/** @file
23 *
24 * The corresponding header file is @ref p3m.hpp.
25 */
26
27#include "config/config.hpp"
28
29#ifdef P3M
30
33
36#ifdef CUDA
39#endif // CUDA
40
41#include "fft/fft.hpp"
43#include "p3m/TuningLogger.hpp"
44#include "p3m/for_each_3d.hpp"
46#include "p3m/math.hpp"
47
48#include "BoxGeometry.hpp"
49#include "LocalBox.hpp"
50#include "Particle.hpp"
52#include "ParticleRange.hpp"
53#include "PropagationMode.hpp"
54#include "actor/visitors.hpp"
57#include "communication.hpp"
58#include "errorhandling.hpp"
60#include "npt.hpp"
62#include "system/System.hpp"
63#include "tuning.hpp"
64
65#include <utils/Vector.hpp>
68#include <utils/math/sqr.hpp>
69
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>
76
77#include <algorithm>
78#include <array>
79#include <cassert>
80#include <complex>
81#include <cstddef>
82#include <functional>
83#include <initializer_list>
84#include <numbers>
85#include <optional>
86#include <span>
87#include <sstream>
88#include <stdexcept>
89#include <string>
90#include <tuple>
91#include <utility>
92
93#ifdef FFTW3_H
94#error "The FFTW3 library shouldn't be visible in this translation unit"
95#endif
96
97template <typename FloatType, Arch Architecture>
99 auto local_n = 0;
100 auto local_q2 = 0.0;
101 auto local_q = 0.0;
102
103 for (auto const &p : get_system().cell_structure->local_particles()) {
104 if (p.q() != 0.0) {
105 local_n++;
106 local_q2 += Utils::sqr(p.q());
107 local_q += p.q();
108 }
109 }
110
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);
115}
116
117/** Calculate the optimal influence function of @cite hockney88a.
118 * (optimised for force calculations)
119 *
120 * Each node calculates only the values for its domain in k-space.
121 *
122 * See also: @cite hockney88a eq. 8-22 (p. 275). Note the somewhat
123 * different convention for the prefactors, which is described in
124 * @cite deserno98a @cite deserno98b.
125 */
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());
132}
133
134/** Calculate the influence function optimized for the energy and the
135 * self energy correction.
136 */
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());
143}
144
145/** Aliasing sum used by @ref p3m_k_space_error. */
147 Utils::Vector3i const &mesh,
148 Utils::Vector3d const &mesh_i, int cao,
149 double alpha_L_i) {
150
151 auto constexpr mesh_start = Utils::Vector3i::broadcast(-P3M_BRILLOUIN);
152 auto constexpr mesh_stop = Utils::Vector3i::broadcast(P3M_BRILLOUIN + 1);
153 auto const factor1 = Utils::sqr(std::numbers::pi * alpha_L_i);
154 auto alias1 = 0.;
155 auto alias2 = 0.;
156
157 Utils::Vector3i indices{};
158 Utils::Vector3i nm{};
159 Utils::Vector3d fnm{};
161 mesh_start, mesh_stop, indices,
162 [&]() {
163 auto const norm_sq = nm.norm2();
164 auto const ex = exp(-factor1 * norm_sq);
165 auto const energy = std::pow(Utils::product(fnm), 2 * cao);
166 alias1 += Utils::sqr(ex) / norm_sq;
167 alias2 += energy * ex * (shift * nm) / norm_sq;
168 },
169 [&](unsigned dim, int n) {
170 nm[dim] = shift[dim] + n * mesh[dim];
171 fnm[dim] = math::sinc(nm[dim] * mesh_i[dim]);
172 });
173
174 return std::make_pair(alias1, alias2);
175}
176
177/** Calculate the real space contribution to the rms error in the force (as
178 * described by Kolafa and Perram).
179 * \param pref Prefactor of Coulomb interaction.
180 * \param r_cut_iL rescaled real space cutoff for p3m method.
181 * \param n_c_part number of charged particles in the system.
182 * \param sum_q2 sum of square of charges in the system
183 * \param alpha_L rescaled Ewald splitting parameter.
184 * \param box_l box dimensions.
185 * \return real space error
186 */
187static double p3m_real_space_error(double pref, double r_cut_iL, int n_c_part,
188 double sum_q2, double alpha_L,
189 Utils::Vector3d const &box_l) {
190 auto const volume = Utils::product(box_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);
193}
194
195/** Calculate the analytic expression of the error estimate for the
196 * P3M method in @cite hockney88a (eq. 8-23 p. 275) in
197 * order to obtain the rms error in the force for a system of N
198 * randomly distributed particles in a cubic box (k-space part).
199 * \param pref Prefactor of Coulomb interaction.
200 * \param mesh number of mesh points in one direction.
201 * \param cao charge assignment order.
202 * \param n_c_part number of charged particles in the system.
203 * \param sum_q2 sum of square of charges in the system
204 * \param alpha_L rescaled Ewald splitting parameter.
205 * \param box_l box dimensions.
206 * \return reciprocal (k) space error
207 */
208static double p3m_k_space_error(double pref, Utils::Vector3i const &mesh,
209 int cao, int n_c_part, double sum_q2,
210 double alpha_L, Utils::Vector3d const &box_l) {
211
212 auto const cotangent_sum = math::get_analytic_cotangent_sum_kernel(cao);
213 auto const mesh_i = 1. / Utils::Vector3d(mesh);
214 auto const alpha_L_i = 1. / alpha_L;
215 auto const mesh_stop = mesh / 2;
216 auto const mesh_start = -mesh_stop;
217 auto indices = Utils::Vector3i{};
218 auto values = Utils::Vector3d{};
219 auto he_q = 0.;
220
222 mesh_start, mesh_stop, indices,
223 [&]() {
224 if ((indices[0] != 0) or (indices[1] != 0) or (indices[2] != 0)) {
225 auto const n2 = indices.norm2();
226 auto const cs = Utils::product(values);
227 auto const [alias1, alias2] =
228 p3m_tune_aliasing_sums(indices, mesh, mesh_i, cao, alpha_L_i);
229 auto const d = alias1 - Utils::sqr(alias2 / cs) / n2;
230 /* at high precision, d can become negative due to extinction;
231 also, don't take values that have no significant digits left*/
232 if (d > 0. and std::fabs(d / alias1) > ROUND_ERROR_PREC) {
233 he_q += d;
234 }
235 }
236 },
237 [&values, &mesh_i, cotangent_sum](unsigned dim, int n) {
238 values[dim] = cotangent_sum(n, mesh_i[dim]);
239 });
240
241 return 2. * pref * sum_q2 * sqrt(he_q / static_cast<double>(n_c_part)) /
242 (box_l[1] * box_l[2]);
243}
244
245template <typename FloatType, Arch Architecture>
247 assert(p3m.params.mesh >= Utils::Vector3i::broadcast(1));
248 assert(p3m.params.cao >= 1 and p3m.params.cao <= 7);
249 assert(p3m.params.alpha > 0.);
250
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();
255
256 p3m.params.cao3 = Utils::int_pow<3>(p3m.params.cao);
257 p3m.params.recalc_a_ai_cao_cut(box_geo.length());
258
259 sanity_checks();
260
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;
265 }
266
267 assert(p3m.fft);
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();
275
276 /* fix box length dependent constants */
277 scaleby_box_l();
278
279 count_charged_particles();
280}
281
282namespace {
283template <int cao> struct AssignCharge {
284 void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos,
285 InterpolationWeights<cao> const &w) {
286 using value_type =
287 typename std::remove_reference_t<decltype(p3m)>::value_type;
288 p3m_interpolate(p3m.local_mesh, w, [q, &p3m](int ind, double w) {
289 p3m.mesh.rs_scalar[ind] += value_type(w * q);
290 });
291 }
292
293 void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos,
294 p3m_interpolation_cache &inter_weights) {
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);
299 }
300
301 void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos) {
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);
305 }
306
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);
312 if (p_q != 0.0) {
313 this->operator()(p3m, p_q, p_pos, p3m.inter_weights);
314 }
315 }
316 }
317};
318} // namespace
319
320template <typename FloatType, Arch Architecture>
322 ParticleRange const &particles) {
323 prepare_fft_mesh(true);
324
325 auto p_q_range = ParticlePropertyRange::charge_range(particles);
326 auto p_pos_range = ParticlePropertyRange::pos_range(particles);
327
328 Utils::integral_parameter<int, AssignCharge, 1, 7>(
329 p3m.params.cao, p3m, boost::combine(p_q_range, p_pos_range));
330}
331
332template <typename FloatType, Arch Architecture>
334 double q, Utils::Vector3d const &real_pos, bool skip_cache) {
335 if (skip_cache) {
336 Utils::integral_parameter<int, AssignCharge, 1, 7>(p3m.params.cao, p3m, q,
337 real_pos);
338 } else {
339 Utils::integral_parameter<int, AssignCharge, 1, 7>(
340 p3m.params.cao, p3m, q, real_pos, p3m.inter_weights);
341 }
342}
343
344template <int cao> struct AssignForces {
345 template <typename combined_ranges>
346 void operator()(auto &p3m, double force_prefac,
347 combined_ranges const &p_q_force_range) const {
348
349 assert(cao == p3m.inter_weights.cao());
350
351 /* charged particle counter */
352 auto p_index = std::size_t{0ul};
353
354 for (auto zipped : p_q_force_range) {
355 auto p_q = boost::get<0>(zipped);
356 auto &p_force = boost::get<1>(zipped);
357 if (p_q != 0.0) {
358 auto const pref = p_q * force_prefac;
359 auto const w = p3m.inter_weights.template load<cao>(p_index);
360
361 Utils::Vector3d force{};
362 p3m_interpolate(p3m.local_mesh, w, [&force, &p3m](int ind, double w) {
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]);
366 });
367
368 p_force -= pref * force;
369 ++p_index;
370 }
371 }
372 }
373};
374
375template <typename combined_ranges>
376static auto calc_dipole_moment(boost::mpi::communicator const &comm,
377 combined_ranges const &p_q_unfolded_pos_range) {
378 auto const local_dip =
379 boost::accumulate(p_q_unfolded_pos_range, Utils::Vector3d{},
380 [](Utils::Vector3d const &dip, auto const &q_pos) {
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;
384 });
385 return boost::mpi::all_reduce(comm, local_dip, std::plus<>());
386}
387
388/** @details Calculate the long range electrostatics part of the pressure
389 * tensor. This is part \f$\Pi_{\textrm{dir}, \alpha, \beta}\f$ eq. (2.6)
390 * in @cite essmann95a. The part \f$\Pi_{\textrm{corr}, \alpha, \beta}\f$
391 * eq. (2.8) is not present here since M is the empty set in our simulations.
392 */
393template <typename FloatType, Arch Architecture>
395 ParticleRange const &particles) {
396 auto const &box_geo = *get_system().box_geo;
397 Utils::Vector9d node_k_space_pressure_tensor{};
398
399 if (p3m.sum_q2 > 0.) {
400 charge_assign(particles);
401 p3m.fft_buffers->perform_scalar_halo_gather();
402 p3m.fft->forward_fft(p3m.fft_buffers->get_scalar_mesh());
403 p3m.update_mesh_views();
404
405 auto constexpr mesh_start = Utils::Vector3i::broadcast(0);
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();
411 auto indices = Utils::Vector3i{};
412 auto index = std::size_t(0u);
413 auto diagonal = 0.;
414
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];
420 auto const norm_sq = Utils::sqr(kx) + Utils::sqr(ky) + Utils::sqr(kz);
421
422 if (norm_sq != 0.) {
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; /* sigma_xx */
430 node_k_space_pressure_tensor[1u] += pref * kx * ky; /* sigma_xy */
431 node_k_space_pressure_tensor[2u] += pref * kx * kz; /* sigma_xz */
432 node_k_space_pressure_tensor[4u] += pref * ky * ky; /* sigma_yy */
433 node_k_space_pressure_tensor[5u] += pref * ky * kz; /* sigma_yz */
434 node_k_space_pressure_tensor[8u] += pref * kz * kz; /* sigma_zz */
435 diagonal += node_k_space_energy;
436 }
437 ++index;
438 });
439
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];
446 }
447
448 return node_k_space_pressure_tensor * prefactor / (2. * box_geo.volume());
449}
450
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;
456#ifdef NPT
457 auto const npt_flag =
458 force_flag and (system.propagation->integ_switch == INTEG_METHOD_NPT_ISO);
459#else
460 auto constexpr npt_flag = false;
461#endif
462
463 if (p3m.sum_q2 > 0.) {
464 if (not has_actor_of_type<ElectrostaticLayerCorrection>(
465 system.coulomb.impl->solver)) {
466 charge_assign(particles);
467 }
468 p3m.fft_buffers->perform_scalar_halo_gather();
469 p3m.fft->forward_fft(p3m.fft_buffers->get_scalar_mesh());
470 p3m.update_mesh_views();
471 }
472
473 auto p_q_range = ParticlePropertyRange::charge_range(particles);
474 auto p_force_range = ParticlePropertyRange::force_range(particles);
475 auto p_unfolded_pos_range =
477
478 // Note: after these calls, the grids are in the order yzx and not xyz
479 // anymore!!!
480 /* The dipole moment is only needed if we don't have metallic boundaries. */
481 auto const box_dipole =
482 (p3m.params.epsilon != P3M_EPSILON_METALLIC)
483 ? std::make_optional(calc_dipole_moment(
484 comm_cart, boost::combine(p_q_range, p_unfolded_pos_range)))
485 : std::nullopt;
486 auto const volume = box_geo.volume();
487 auto const pref =
488 4. * std::numbers::pi / volume / (2. * p3m.params.epsilon + 1.);
489
490 /* === k-space force calculation === */
491 if (force_flag) {
492 /* i*k differentiation */
493 auto constexpr mesh_start = Utils::Vector3i::broadcast(0);
494 auto const &mesh_stop = p3m.mesh.size;
495 auto const &offset = p3m.mesh.start;
496 auto const wavevector = Utils::Vector3<FloatType>((2. * std::numbers::pi) *
497 box_geo.length_inv());
498 auto indices = Utils::Vector3i{};
499 auto index = std::size_t(0u);
500
501 /* compute electric field */
502 // Eq. (3.49) @cite deserno00b
503 for_each_3d(mesh_start, mesh_stop, indices, [&]() {
504 auto const rho_hat =
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;
508
509 for (int d = 0; d < 3; d++) {
510 /* direction in r-space: */
511 int d_rs = (d + p3m.mesh.ks_pnum) % 3;
512 /* directions */
513 auto const k = FloatType(p3m.d_op[d_rs][indices[d] + offset[d]]) *
514 wavevector[d_rs];
515
516 /* i*k*(Re+i*Im) = - Im*k + i*Re*k (i=sqrt(-1)) */
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();
519 }
520
521 ++index;
522 });
523
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);
529 }
530 p3m.fft_buffers->perform_vector_halo_spread();
531 p3m.fft->check_complex_residuals = false;
532
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));
537
538 // add dipole forces
539 // Eq. (3.19) @cite deserno00b
540 if (box_dipole) {
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);
545 p_force -= p_q * dm;
546 }
547 }
548 }
549
550 /* === k-space energy calculation === */
551 if (energy_flag or npt_flag) {
552 auto node_energy = 0.;
553 auto const mesh_length = Utils::product(p3m.mesh.size);
554 for (int i = 0; i < mesh_length; i++) {
555 // Use the energy optimized influence function for energy!
556 // Eq. (3.40) @cite deserno00b
557 node_energy +=
558 double(p3m.g_energy[i] * (Utils::sqr(p3m.mesh.rs_scalar[2 * i + 0]) +
559 Utils::sqr(p3m.mesh.rs_scalar[2 * i + 1])));
560 }
561 node_energy /= 2. * volume;
562
563 auto energy = 0.;
564 boost::mpi::reduce(comm_cart, node_energy, energy, std::plus<>(), 0);
565 if (this_node == 0) {
566 /* self energy correction */
567 // Eq. (3.8) @cite deserno00b
568 energy -= p3m.sum_q2 * p3m.params.alpha * std::numbers::inv_sqrtpi;
569 /* net charge correction */
570 // Eq. (3.11) @cite deserno00b
571 energy -= p3m.square_sum_q * std::numbers::pi /
572 (2. * volume * Utils::sqr(p3m.params.alpha));
573 /* dipole correction */
574 // Eq. (3.9) @cite deserno00b
575 if (box_dipole) {
576 energy += pref * box_dipole.value().norm2();
577 }
578 }
579 energy *= prefactor;
580#ifdef NPT
581 if (npt_flag) {
583 }
584#endif
585 if (not energy_flag) {
586 energy = 0.;
587 }
588 return energy;
589 }
590
591 return 0.;
592}
593
594template <typename FloatType, Arch Architecture>
597 double m_mesh_density_min = -1., m_mesh_density_max = -1.;
598 // indicates if mesh should be tuned
599 bool m_tune_mesh = false;
600
601protected:
602 P3MParameters &get_params() override { return p3m.params; }
603
604public:
605 CoulombTuningAlgorithm(System::System &system, auto &input_p3m,
606 double prefactor, int timings)
607 : TuningAlgorithm(system, prefactor, timings), p3m{input_p3m} {}
608
609 void on_solver_change() const override { m_system.on_coulomb_change(); }
610
611 void setup_logger(bool verbose) override {
612 auto const &box_geo = *m_system.box_geo;
613#ifdef CUDA
614 auto const on_gpu = Architecture == Arch::GPU;
615#else
616 auto const on_gpu = false;
617#endif
618 m_logger = std::make_unique<TuningLogger>(
619 verbose and this_node == 0, (on_gpu) ? "CoulombP3MGPU" : "CoulombP3M",
621 m_logger->tuning_goals(p3m.params.accuracy, m_prefactor,
622 box_geo.length()[0], p3m.sum_qpart, p3m.sum_q2);
623 m_logger->log_tuning_start();
624 }
625
626 std::optional<std::string>
627 layer_correction_veto_r_cut(double r_cut) const override {
628 auto const &solver = m_system.coulomb.impl->solver;
629 if (auto actor = get_actor_by_type<ElectrostaticLayerCorrection>(solver)) {
630 return actor->veto_r_cut(r_cut);
631 }
632 return {};
633 }
634
635 std::tuple<double, double, double, double>
637 double r_cut_iL) const override {
638
639 auto const &box_geo = *m_system.box_geo;
640 double alpha_L, rs_err, ks_err;
641
642 /* calc maximal real space error for setting */
643 rs_err = p3m_real_space_error(m_prefactor, r_cut_iL, p3m.sum_qpart,
644 p3m.sum_q2, 0., box_geo.length());
645
646 if (std::numbers::sqrt2 * rs_err > p3m.params.accuracy) {
647 /* assume rs_err = ks_err -> rs_err = accuracy/sqrt(2.0) -> alpha_L */
648 alpha_L = sqrt(log(std::numbers::sqrt2 * rs_err / p3m.params.accuracy)) /
649 r_cut_iL;
650 } else {
651 /* even alpha=0 is ok, however, we cannot choose it since it kills the
652 k-space error formula.
653 Anyways, this very likely NOT the optimal solution */
654 alpha_L = 0.1;
655 }
656
657 /* calculate real-space and k-space error for this alpha_L */
658 rs_err = p3m_real_space_error(m_prefactor, r_cut_iL, p3m.sum_qpart,
659 p3m.sum_q2, alpha_L, box_geo.length());
660#ifdef CUDA
661 if constexpr (Architecture == Arch::GPU) {
662 if (this_node == 0) {
663 ks_err =
665 p3m.sum_q2, alpha_L, box_geo.length().data());
666 }
667 boost::mpi::broadcast(comm_cart, ks_err, 0);
668 } else
669#endif
670 ks_err = p3m_k_space_error(m_prefactor, mesh, cao, p3m.sum_qpart,
671 p3m.sum_q2, alpha_L, box_geo.length());
672
673 return {Utils::Vector2d{rs_err, ks_err}.norm(), rs_err, ks_err, alpha_L};
674 }
675
676 void determine_mesh_limits() override {
677 auto const &box_geo = *m_system.box_geo;
678 auto const mesh_density =
679 static_cast<double>(p3m.params.mesh[0]) * box_geo.length_inv()[0];
680
681 if (p3m.params.mesh == Utils::Vector3i::broadcast(-1)) {
682 /* avoid using more than 1 GB of FFT arrays */
683 auto const normalized_box_dim = std::cbrt(box_geo.volume());
684 auto const max_npart_per_dim = 512.;
685 /* simple heuristic to limit the tried meshes if the accuracy cannot
686 be obtained with smaller meshes, but normally not all these
687 meshes have to be tested */
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;
692 m_tune_mesh = true;
693 } else {
694 m_mesh_density_min = m_mesh_density_max = mesh_density;
695 assert(p3m.params.mesh[0] >= 1);
696 if (p3m.params.mesh[1] == -1 and p3m.params.mesh[2] == -1) {
697 // determine the two missing values by rescaling by the box length
698 for (auto i : {1, 2}) {
699 p3m.params.mesh[i] =
700 static_cast<int>(std::round(mesh_density * box_geo.length()[i]));
701 // make the mesh even in all directions
702 p3m.params.mesh[i] += p3m.params.mesh[i] % 2;
703 }
704 }
705 m_logger->report_fixed_mesh(p3m.params.mesh);
706 }
707 }
708
710 auto const &box_geo = *m_system.box_geo;
711 auto const &solver = m_system.coulomb.impl->solver;
712 auto tuned_params = TuningAlgorithm::Parameters{};
713 auto time_best = time_sentinel;
714 auto mesh_density = m_mesh_density_min;
715 while (mesh_density <= m_mesh_density_max) {
716 auto trial_params = TuningAlgorithm::Parameters{};
717 if (m_tune_mesh) {
718 for (auto i : {0, 1, 2}) {
719 trial_params.mesh[i] =
720 static_cast<int>(std::round(box_geo.length()[i] * mesh_density));
721 // make the mesh even in all directions
722 trial_params.mesh[i] += trial_params.mesh[i] % 2;
723 }
724 } else {
725 trial_params.mesh = p3m.params.mesh;
726 }
727 trial_params.cao = cao_best;
728
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);
732
733 if (trial_time >= 0.) {
734 /* the optimum r_cut for this mesh is the upper limit for higher meshes,
735 everything else is slower */
736 if (has_actor_of_type<CoulombP3M>(solver)) {
737 m_r_cut_iL_max = trial_params.r_cut_iL;
738 }
739
740 if (trial_time < time_best) {
741 /* new optimum */
743 tuned_params = trial_params;
744 time_best = tuned_params.time = trial_time;
745 } else if (trial_time > time_best + time_granularity or
747 /* no hope of further optimisation */
748 break;
749 }
750 }
751 mesh_density += 0.1;
752 }
753 return tuned_params;
754 }
755};
756
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];
763 }
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];
766 }
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");
772 }
773 try {
775 system, p3m, prefactor, tune_timings);
776 parameters.setup_logger(tune_verbose);
777 // parameter ranges
778 parameters.determine_mesh_limits();
779 parameters.determine_r_cut_limits();
780 parameters.determine_cao_limits(7);
781 // run tuning algorithm
782 parameters.tune();
783 m_is_tuned = true;
784 system.on_coulomb_change();
785 } catch (...) {
786 p3m.params.tuning = false;
787 throw;
788 }
789 }
790 init();
791}
792
794 auto const &system = get_system();
795 auto const &box_geo = *system.box_geo;
796 auto const &local_geo = *system.local_geo;
797 for (auto i = 0u; i < 3u; i++) {
798 /* check k-space cutoff */
799 if (p3m_params.cao_cut[i] >= box_geo.length_half()[i]) {
800 std::stringstream msg;
801 msg << "P3M_init: k-space cutoff " << p3m_params.cao_cut[i]
802 << " is larger than half of box dimension " << box_geo.length()[i];
803 throw std::runtime_error(msg.str());
804 }
805 if (p3m_params.cao_cut[i] >= local_geo.length()[i]) {
806 std::stringstream msg;
807 msg << "P3M_init: k-space cutoff " << p3m_params.cao_cut[i]
808 << " is larger than local box dimension " << local_geo.length()[i];
809 throw std::runtime_error(msg.str());
810 }
811 }
812
814 if ((box_geo.length()[0] != box_geo.length()[1]) or
815 (box_geo.length()[1] != box_geo.length()[2]) or
816 (p3m_params.mesh[0] != p3m_params.mesh[1]) or
817 (p3m_params.mesh[1] != p3m_params.mesh[2])) {
818 throw std::runtime_error(
819 "CoulombP3M: non-metallic epsilon requires cubic box");
820 }
821 }
822}
823
825 auto const &box_geo = *get_system().box_geo;
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)");
829 }
830}
831
833 auto const &local_geo = *get_system().local_geo;
834 if (local_geo.cell_structure_type() != CellStructureType::REGULAR and
835 local_geo.cell_structure_type() != CellStructureType::HYBRID) {
836 throw std::runtime_error(
837 "CoulombP3M: requires the regular or hybrid decomposition cell system");
838 }
839 if (::communicator.size > 1 and
840 local_geo.cell_structure_type() == CellStructureType::HYBRID) {
841 throw std::runtime_error(
842 "CoulombP3M: does not work with the hybrid decomposition cell system, "
843 "if using more than one MPI node");
844 }
845}
846
848 auto const &node_grid = ::communicator.node_grid;
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");
852 }
853}
854
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();
865}
866
867#ifdef CUDA
868template <typename FloatType, Arch Architecture>
870 ParticleRange const &particles) {
871 if constexpr (Architecture == Arch::GPU) {
872#ifdef NPT
873 if (get_system().propagation->integ_switch == INTEG_METHOD_NPT_ISO) {
874 auto const energy = long_range_energy(particles);
876 }
877#else
878 static_cast<void>(particles);
879#endif
880 if (this_node == 0) {
881 auto &gpu = get_system().gpu;
882 p3m_gpu_add_farfield_force(*m_gpu_data, gpu, prefactor,
883 gpu.n_particles());
884 }
885 }
886}
887
888/* Initialize the CPU kernels.
889 * This operation is time-consuming and sets up data members
890 * that are only relevant for ELC force corrections, since the
891 * GPU implementation uses CPU kernels to compute energies.
892 */
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)) {
899 init_cpu_kernels();
900 }
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());
903 }
904}
905
906template <typename FloatType, Arch Architecture>
908 if constexpr (Architecture == Arch::GPU) {
909 auto &gpu_particle_data = get_system().gpu;
911 gpu_particle_data.enable_property(GpuParticleData::prop::q);
912 gpu_particle_data.enable_property(GpuParticleData::prop::pos);
913 }
914}
915#endif // CUDA
916
921
922#endif // P3M
@ HYBRID
Hybrid decomposition.
@ REGULAR
Regular decomposition.
@ INTEG_METHOD_NPT_ISO
Vector implementation and trait types for boost qvm interoperability.
TuningAlgorithm::Parameters get_time() override
Definition p3m.cpp:709
P3MParameters & get_params() override
Definition p3m.cpp:602
std::tuple< double, double, double, double > calculate_accuracy(Utils::Vector3i const &mesh, int cao, double r_cut_iL) const override
Definition p3m.cpp:636
void setup_logger(bool verbose) override
Definition p3m.cpp:611
void on_solver_change() const override
Definition p3m.cpp:609
std::optional< std::string > layer_correction_veto_r_cut(double r_cut) const override
Definition p3m.cpp:627
CoulombTuningAlgorithm(System::System &system, auto &input_p3m, double prefactor, int timings)
Definition p3m.cpp:605
void determine_mesh_limits() override
Definition p3m.cpp:676
void enable_property(std::size_t property)
A range of particles.
Main system class.
GpuParticleData gpu
std::shared_ptr< Propagation > propagation
std::shared_ptr< CellStructure > cell_structure
Coulomb::Solver coulomb
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).
T norm() const
Definition Vector.hpp:138
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.
Definition Vector.hpp:110
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.
Definition config.hpp:66
#define P3M_BRILLOUIN
P3M: Number of Brillouin zones taken into account in the calculation of the optimal influence functio...
Definition config.hpp:53
void charge_assign(elc_data const &elc, CoulombP3M &solver, combined_ranges const &p_q_pos_range)
Definition elc.cpp:1111
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)
Definition Vector.hpp:374
VectorXd< 3 > Vector3d
Definition Vector.hpp:164
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
DEVICE_QUALIFIER auto sinc(T x)
Calculate the function .
Definition math.hpp:53
auto get_analytic_cotangent_sum_kernel(int cao)
Definition math.hpp:128
void npt_add_virial_contribution(double energy)
Definition npt.cpp:141
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.
Definition p3m.cpp:146
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.
Definition p3m.cpp:208
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...
Definition p3m.cpp:187
static auto calc_dipole_moment(boost::mpi::communicator const &comm, combined_ranges const &p_q_unfolded_pos_range)
Definition p3m.cpp:376
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
Definition p3m.cpp:346
Utils::Vector3i node_grid
void charge_assign(ParticleRange const &particles) override
Definition p3m.cpp:321
void tune() override
Definition p3m.cpp:758
void assign_charge(double q, Utils::Vector3d const &real_pos, bool skip_cache) override
Definition p3m.cpp:333
void calc_influence_function_energy() override
Calculate the influence function optimized for the energy and the self energy correction.
Definition p3m.cpp:138
void scaleby_box_l() override
Definition p3m.cpp:856
void calc_influence_function_force() override
Calculate the optimal influence function of .
Definition p3m.cpp:127
Utils::Vector9d long_range_pressure(ParticleRange const &particles) override
Definition p3m.cpp:394
void init_cpu_kernels()
Definition p3m.cpp:246
void count_charged_particles() override
Definition p3m.cpp:98
double long_range_kernel(bool force_flag, bool energy_flag, ParticleRange const &particles)
Compute the k-space part of forces and energies.
Definition p3m.cpp:452
void request_gpu() const
Definition p3m.cpp:907
void init_gpu_kernels()
Definition p3m.cpp:894
void add_long_range_forces_gpu(ParticleRange const &particles)
Definition p3m.cpp:869
void sanity_checks_periodicity() const
Definition p3m.cpp:824
void sanity_checks_boxl() const
Checks for correctness of the k-space cutoff.
Definition p3m.cpp:793
void sanity_checks_cell_structure() const
Definition p3m.cpp:832
P3MParameters const & p3m_params
Definition p3m.hpp:56
void sanity_checks_node_grid() const
Definition p3m.cpp:847
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
Definition Array.hpp:120
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos, p3m_interpolation_cache &inter_weights)
Definition p3m.cpp:293
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos, InterpolationWeights< cao > const &w)
Definition p3m.cpp:284
void operator()(auto &p3m, combined_ranges const &p_q_pos_range)
Definition p3m.cpp:308
void operator()(auto &p3m, double q, Utils::Vector3d const &real_pos)
Definition p3m.cpp:301
double sum_q2
Sum of square of charges (only on head node).
Definition p3m.impl.hpp:50
int sum_qpart
number of charged particles (only on head node).
Definition p3m.impl.hpp:48
P3MParameters params
P3M base parameters.