51#ifndef MMM1D_MACHINE_PREC
57 auto const wavenumber = 2. * std::numbers::pi * box_l_inv[2];
59 auto const rhores = wavenumber * minrad;
60 auto const pref = 4. * box_l_inv[2] * std::max(1., wavenumber);
62 return pref *
K1(rhores *
P) * exp(rhores) / rhores * (
P - 1. + 1. / rhores);
69 auto constexpr min_rad = 0.01;
70 auto const rgranularity = min_rad * box_l[2];
71 auto rmin = rgranularity;
72 auto rmax = std::min(box_l[0], box_l[1]);
73 auto const errmin =
far_error(
P, rmin, box_l_inv);
74 auto const errmax =
far_error(
P, rmax, box_l_inv);
75 if (errmin < maxPWerror) {
79 if (errmax > maxPWerror) {
81 return 2. * std::max(box_l[0], box_l[1]);
84 while (rmax - rmin > rgranularity) {
85 auto const c = 0.5 * (rmin + rmax);
87 if (errc > maxPWerror) {
93 return 0.5 * (rmin + rmax);
106void CoulombMMM1D::determine_bessel_radii() {
108 auto const &box_l = box_geo.length();
109 auto const &box_l_inv = box_geo.length_inv();
110 for (
int i = 0; i < MAXIMAL_B_CUT; ++i) {
115void CoulombMMM1D::prepare_polygamma_series() {
118 auto const rhomax2 = uz2 * far_switch_radius_sq;
121 auto rhomax2nm2 = 1.0;
123 create_mod_psi_up_to(n + 1);
126 err = 2. *
static_cast<double>(n) * fabs(
mod_psi_even(modPsi, n, 0.5)) *
128 rhomax2nm2 *= rhomax2;
134 double switch_rad,
int tune_timings,
136 : maxPWerror{maxPWerror}, far_switch_radius{switch_rad},
137 tune_timings{tune_timings}, tune_verbose{tune_verbose}, m_is_tuned{false},
138 far_switch_radius_sq{-1.}, uz2{0.}, prefuz2{0.}, prefL3_i{0.} {
141 throw std::domain_error(
"Parameter 'maxPWerror' must be > 0");
144 throw std::domain_error(
"Parameter 'far_switch_radius' must be > 0");
150 throw std::domain_error(
"Parameter 'timings' must be > 0");
154void CoulombMMM1D::sanity_checks_periodicity()
const {
156 if (box_geo.periodic(0) || box_geo.periodic(1) || !box_geo.periodic(2)) {
157 throw std::runtime_error(
"MMM1D requires periodicity (False, False, True)");
161void CoulombMMM1D::sanity_checks_cell_structure()
const {
162 auto const &local_geo = *
get_system().local_geo;
164 throw std::runtime_error(
"MMM1D requires the N-square cellsystem");
168void CoulombMMM1D::recalc_boxl_parameters() {
171 if (far_switch_radius_sq >=
Utils::sqr(box_geo.length()[2]))
172 far_switch_radius_sq = 0.8 *
Utils::sqr(box_geo.length()[2]);
176 prefL3_i = prefuz2 * box_geo.length_inv()[2];
178 determine_bessel_radii();
179 prepare_polygamma_series();
184 auto constexpr c_2pi = 2. * std::numbers::pi;
186 auto const n_modPsi =
static_cast<int>(modPsi.size()) >> 1;
187 auto const rxy2 = d[0] * d[0] + d[1] * d[1];
188 auto const rxy2_d = rxy2 * uz2;
189 auto const z_d = d[2] * box_geo.length_inv()[2];
192 if (rxy2 <= far_switch_radius_sq) {
197 for (
int n = 1; n < n_modPsi; n++) {
198 auto const deriv =
static_cast<double>(2 * n);
201 auto const r2n = r2nm1 * rxy2_d;
204 sr += deriv * r2nm1 * mpe;
212 double Fx = prefL3_i * sr * d[0];
213 double Fy = prefL3_i * sr * d[1];
214 double Fz = prefuz2 * sz;
218 double pref, rt, rt2, shift_z;
220 pref = 1. / Utils::int_pow<3>(dist);
225 shift_z = d[2] + box_geo.length()[2];
226 rt2 = rxy2 + shift_z * shift_z;
228 pref = 1. / (rt2 * rt);
231 Fz += pref * shift_z;
233 shift_z = d[2] - box_geo.length()[2];
234 rt2 = rxy2 + shift_z * shift_z;
236 pref = 1. / (rt2 * rt);
239 Fz += pref * shift_z;
241 force = {Fx, Fy, Fz};
244 auto const rxy = sqrt(rxy2);
245 auto const rxy_d = rxy * box_geo.length_inv()[2];
246 auto sr = 0., sz = 0.;
248 for (
int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
249 if (bessel_radii[bp - 1] < rxy)
252 auto const fq = c_2pi * bp;
253#ifdef MMM1D_MACHINE_PREC
254 auto const k0 =
K0(fq * rxy_d);
255 auto const k1 =
K1(fq * rxy_d);
257 auto const [k0, k1] =
LPK01(fq * rxy_d);
259 sr += bp * k1 * cos(fq * z_d);
260 sz += bp * k0 * sin(fq * z_d);
262 sr *= uz2 * 4. * c_2pi;
263 sz *= uz2 * 4. * c_2pi;
265 auto const pref = sr / rxy + 2. * box_geo.length_inv()[2] / rxy2;
267 force = {pref * d[0], pref * d[1], sz};
274 double const dist)
const {
278 auto constexpr c_2pi = 2. * std::numbers::pi;
280 auto const n_modPsi =
static_cast<int>(modPsi.size()) >> 1;
281 auto const rxy2 = d[0] * d[0] + d[1] * d[1];
282 auto const rxy2_d = rxy2 * uz2;
283 auto const z_d = d[2] * box_geo.length_inv()[2];
286 if (rxy2 <= far_switch_radius_sq) {
288 energy = -2. * std::numbers::egamma;
292 for (
int n = 0; n < n_modPsi; n++) {
301 energy *= box_geo.length_inv()[2];
309 shift_z = d[2] + box_geo.length()[2];
310 rt = sqrt(rxy2 + shift_z * shift_z);
313 shift_z = d[2] - box_geo.length()[2];
314 rt = sqrt(rxy2 + shift_z * shift_z);
318 auto const rxy = sqrt(rxy2);
319 auto const rxy_d = rxy * box_geo.length_inv()[2];
323 -0.25 * log(rxy2_d) + 0.5 * (std::numbers::ln2 - std::numbers::egamma);
324 for (
int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
325 if (bessel_radii[bp - 1] < rxy)
328 auto const fq = c_2pi * bp;
329 energy +=
K0(fq * rxy_d) * cos(fq * z_d);
331 energy *= 4. * box_geo.length_inv()[2];
341 recalc_boxl_parameters();
344 if (far_switch_radius_sq < 0.) {
345 auto const &box_geo = *system.box_geo;
346 auto const maxrad = box_geo.length()[2];
347 auto min_time = std::numeric_limits<double>::infinity();
349 auto switch_radius = 0.2 * maxrad;
351 while (switch_radius < 0.4 * maxrad) {
352 if (switch_radius > bessel_radii.back()) {
354 far_switch_radius_sq =
Utils::sqr(switch_radius);
355 system.on_coulomb_change();
361 std::printf(
"r= %f t= %f ms\n", switch_radius, int_time);
364 if (int_time < min_time) {
366 min_rad = switch_radius;
367 }
else if (int_time > 2. * min_time) {
372 switch_radius += 0.025 * maxrad;
374 switch_radius = min_rad;
375 far_switch_radius_sq =
Utils::sqr(switch_radius);
376 }
else if (far_switch_radius_sq <=
Utils::sqr(bessel_radii.back())) {
378 throw std::runtime_error(
"MMM1D could not find a reasonable Bessel cutoff");
382 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 ...
static double mod_psi_odd(auto const &modPsi, int n, double x)
Modified polygamma for odd order 2*n+1, n>= 0
static auto determine_minrad(double maxPWerror, int P, Utils::Vector3d const &box_l, Utils::Vector3d const &box_l_inv)
static double mod_psi_even(auto const &modPsi, int n, double x)
Modified polygamma for even order 2*n, n >= 0
static auto far_error(int P, double minrad, Utils::Vector3d const &box_l_inv)
MMM1D algorithm for long-range Coulomb interactions on the CPU.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
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...
double evaluateAsTaylorSeriesAt(std::span< const double > series, double x)
Evaluate the polynomial interpreted as a Taylor series via the Horner scheme.
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.