53#ifndef MMM1D_MACHINE_PREC
59 auto const wavenumber = 2. *
Utils::pi() * box_l_inv[2];
61 auto const rhores = wavenumber * minrad;
62 auto const pref = 4. * box_l_inv[2] * std::max(1., wavenumber);
64 return pref *
K1(rhores *
P) * exp(rhores) / rhores * (
P - 1. + 1. / rhores);
71 auto constexpr min_rad = 0.01;
72 auto const rgranularity = min_rad * box_l[2];
73 auto rmin = rgranularity;
74 auto rmax = std::min(box_l[0], box_l[1]);
75 auto const errmin =
far_error(
P, rmin, box_l_inv);
76 auto const errmax =
far_error(
P, rmax, box_l_inv);
83 return 2. * std::max(box_l[0], box_l[1]);
86 while (rmax - rmin > rgranularity) {
87 auto const c = 0.5 * (rmin + rmax);
95 return 0.5 * (rmin + rmax);
98void CoulombMMM1D::determine_bessel_radii() {
100 auto const &box_l = box_geo.length();
101 auto const &box_l_inv = box_geo.length_inv();
102 for (
int i = 0; i < MAXIMAL_B_CUT; ++i) {
107void CoulombMMM1D::prepare_polygamma_series() {
110 auto const rhomax2 = uz2 * far_switch_radius_sq;
113 auto rhomax2nm2 = 1.0;
118 err = 2. *
static_cast<double>(n) * fabs(
mod_psi_even(n, 0.5)) * rhomax2nm2;
119 rhomax2nm2 *= rhomax2;
125 double switch_rad,
int tune_timings,
128 tune_timings{tune_timings}, tune_verbose{tune_verbose}, m_is_tuned{false},
132 throw std::domain_error(
"Parameter 'maxPWerror' must be > 0");
135 throw std::domain_error(
"Parameter 'far_switch_radius' must be > 0");
141 throw std::domain_error(
"Parameter 'timings' must be > 0");
145void CoulombMMM1D::sanity_checks_periodicity()
const {
147 if (box_geo.periodic(0) || box_geo.periodic(1) || !box_geo.periodic(2)) {
148 throw std::runtime_error(
"MMM1D requires periodicity (False, False, True)");
152void CoulombMMM1D::sanity_checks_cell_structure()
const {
153 auto const &local_geo = *
get_system().local_geo;
155 throw std::runtime_error(
"MMM1D requires the N-square cellsystem");
159void CoulombMMM1D::recalc_boxl_parameters() {
162 if (far_switch_radius_sq >=
Utils::sqr(box_geo.length()[2]))
167 prefL3_i = prefuz2 * box_geo.length_inv()[2];
169 determine_bessel_radii();
170 prepare_polygamma_series();
177 auto const n_modPsi =
static_cast<int>(
modPsi.size()) >> 1;
178 auto const rxy2 = d[0] * d[0] + d[1] * d[1];
179 auto const rxy2_d = rxy2 * uz2;
180 auto const z_d = d[2] * box_geo.length_inv()[2];
188 for (
int n = 1; n < n_modPsi; n++) {
189 auto const deriv =
static_cast<double>(2 * n);
192 auto const r2n = r2nm1 * rxy2_d;
195 sr += deriv * r2nm1 * mpe;
203 double Fx = prefL3_i * sr * d[0];
204 double Fy = prefL3_i * sr * d[1];
205 double Fz = prefuz2 * sz;
209 double pref, rt, rt2, shift_z;
211 pref = 1. / Utils::int_pow<3>(dist);
216 shift_z = d[2] + box_geo.length()[2];
217 rt2 = rxy2 + shift_z * shift_z;
219 pref = 1. / (rt2 * rt);
222 Fz += pref * shift_z;
224 shift_z = d[2] - box_geo.length()[2];
225 rt2 = rxy2 + shift_z * shift_z;
227 pref = 1. / (rt2 * rt);
230 Fz += pref * shift_z;
232 force = {Fx, Fy, Fz};
235 auto const rxy = sqrt(rxy2);
236 auto const rxy_d = rxy * box_geo.length_inv()[2];
237 auto sr = 0., sz = 0.;
239 for (
int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
240 if (bessel_radii[bp - 1] < rxy)
243 auto const fq = c_2pi * bp;
244#ifdef MMM1D_MACHINE_PREC
245 auto const k0 =
K0(fq * rxy_d);
246 auto const k1 =
K1(fq * rxy_d);
248 auto const [k0, k1] =
LPK01(fq * rxy_d);
250 sr += bp * k1 * cos(fq * z_d);
251 sz += bp * k0 * sin(fq * z_d);
253 sr *= uz2 * 4. * c_2pi;
254 sz *= uz2 * 4. * c_2pi;
256 auto const pref = sr / rxy + 2. * box_geo.length_inv()[2] / rxy2;
258 force = {pref * d[0], pref * d[1], sz};
265 double const dist)
const {
271 auto const n_modPsi =
static_cast<int>(
modPsi.size()) >> 1;
272 auto const rxy2 = d[0] * d[0] + d[1] * d[1];
273 auto const rxy2_d = rxy2 * uz2;
274 auto const z_d = d[2] * box_geo.length_inv()[2];
283 for (
int n = 0; n < n_modPsi; n++) {
292 energy *= box_geo.length_inv()[2];
300 shift_z = d[2] + box_geo.length()[2];
301 rt = sqrt(rxy2 + shift_z * shift_z);
304 shift_z = d[2] - box_geo.length()[2];
305 rt = sqrt(rxy2 + shift_z * shift_z);
309 auto const rxy = sqrt(rxy2);
310 auto const rxy_d = rxy * box_geo.length_inv()[2];
314 for (
int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
315 if (bessel_radii[bp - 1] < rxy)
318 auto const fq = c_2pi * bp;
319 energy +=
K0(fq * rxy_d) * cos(fq * z_d);
321 energy *= 4. * box_geo.length_inv()[2];
331 recalc_boxl_parameters();
335 auto const &box_geo = *system.box_geo;
336 auto const maxrad = box_geo.length()[2];
337 auto min_time = std::numeric_limits<double>::infinity();
339 auto switch_radius = 0.2 * maxrad;
341 while (switch_radius < 0.4 * maxrad) {
342 if (switch_radius > bessel_radii.back()) {
345 system.on_coulomb_change();
351 std::printf(
"r= %f t= %f ms\n", switch_radius, int_time);
354 if (int_time < min_time) {
356 min_rad = switch_radius;
357 }
else if (int_time > 2. * min_time) {
362 switch_radius += 0.025 * maxrad;
364 switch_radius = min_rad;
368 throw std::runtime_error(
"MMM1D could not find a reasonable Bessel cutoff");
372 system.on_coulomb_change();
@ NSQUARE
Atom decomposition (N-square).
Vector implementation and trait types for boost qvm interoperability.
void set_prefactor(double new_prefactor)
double prefactor
Electrostatics prefactor.
This file contains the defaults for ESPResSo.
This file contains the errorhandling code for severe errors, like a broken bond or illegal parameter ...
Common parts of the MMM family of methods for the electrostatic interaction: MMM1D and ELC.
double mod_psi_even(int n, double x)
Modified polygamma for even order 2*n, n >= 0
double mod_psi_odd(int n, double x)
Modified polygamma for odd order 2*n+1, n>= 0
std::vector< std::vector< double > > modPsi
Table of the Taylor expansions of the modified polygamma functions.
void create_mod_psi_up_to(int new_n)
Create both the even and odd polygamma functions up to order 2*new_n
Common parts of the MMM family of methods for the electrostatic interaction: MMM1D and ELC.
static auto determine_minrad(double maxPWerror, int P, Utils::Vector3d const &box_l, Utils::Vector3d const &box_l_inv)
static auto far_error(int P, double minrad, Utils::Vector3d const &box_l_inv)
MMM1D algorithm for long-range Coulomb interactions on the CPU.
__constant__ float far_switch_radius_sq[1]
__constant__ float maxPWerror[1]
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.
DEVICE_QUALIFIER constexpr T ln_2()
Natural logarithm of 2.
DEVICE_QUALIFIER constexpr T gamma()
Euler-Mascheroni constant.
std::pair< double, double > LPK01(double x)
Modified Bessel functions of second kind, order 0 and 1, low precision.
double K1(double x)
Modified Bessel function of second kind, order 1.
double K0(double x)
Modified Bessel function of second kind, order 0.
This file contains implementations for some special functions which are needed by the MMM family of a...
CoulombMMM1D(double prefactor, double maxPWerror, double switch_rad, int tune_timings, bool tune_verbose)
double far_switch_radius
Far switch radius.
Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair force.
double pair_energy(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair energy.
double maxPWerror
Maximal allowed pairwise error for the potential and force.
double benchmark_integration_step(System::System &system, int int_steps)
Benchmark the integration loop.