ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
mmm1d.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 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#include "config/config.hpp"
23
24#ifdef ELECTROSTATICS
25
27
29
30#include "BoxGeometry.hpp"
31#include "LocalBox.hpp"
32#include "Particle.hpp"
34#include "errorhandling.hpp"
35#include "specfunc.hpp"
36#include "tuning.hpp"
37
38#include <utils/Vector.hpp>
39#include <utils/math/sqr.hpp>
40
41#include <algorithm>
42#include <cmath>
43#include <cstdio>
44#include <limits>
45#include <numbers>
46#include <vector>
47
48/* if you define this feature, the Bessel functions are calculated up
49 * to machine precision, otherwise 10^-14, which should be
50 * definitely enough for daily life. */
51#ifndef MMM1D_MACHINE_PREC
52#define K0 LPK0
53#define K1 LPK1
54#endif
55
56static auto far_error(int P, double minrad, Utils::Vector3d const &box_l_inv) {
57 auto const wavenumber = 2. * std::numbers::pi * box_l_inv[2];
58 // this uses an upper bound to all force components and the potential
59 auto const rhores = wavenumber * minrad;
60 auto const pref = 4. * box_l_inv[2] * std::max(1., wavenumber);
61
62 return pref * K1(rhores * P) * exp(rhores) / rhores * (P - 1. + 1. / rhores);
63}
64
65static auto determine_minrad(double maxPWerror, int P,
66 Utils::Vector3d const &box_l,
67 Utils::Vector3d const &box_l_inv) {
68 // bisection to search for where the error is maxPWerror
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) {
76 // we can do almost all radii with this P
77 return rmin;
78 }
79 if (errmax > maxPWerror) {
80 // make sure that this switching radius cannot be reached
81 return 2. * std::max(box_l[0], box_l[1]);
82 }
83
84 while (rmax - rmin > rgranularity) {
85 auto const c = 0.5 * (rmin + rmax);
86 auto const errc = far_error(P, c, box_l_inv);
87 if (errc > maxPWerror) {
88 rmin = c;
89 } else {
90 rmax = c;
91 }
92 }
93 return 0.5 * (rmin + rmax);
94}
95
96/** Modified polygamma for even order <tt>2*n, n >= 0</tt> */
97static double mod_psi_even(auto const &modPsi, int n, double x) {
98 return evaluateAsTaylorSeriesAt(modPsi[2 * n], x * x);
99}
100
101/** Modified polygamma for odd order <tt>2*n+1, n>= 0</tt> */
102static double mod_psi_odd(auto const &modPsi, int n, double x) {
103 return x * evaluateAsTaylorSeriesAt(modPsi[2 * n + 1], x * x);
104}
105
106void CoulombMMM1D::determine_bessel_radii() {
107 auto const &box_geo = *get_system().box_geo;
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) {
111 bessel_radii[i] = determine_minrad(maxPWerror, i + 1, box_l, box_l_inv);
112 }
113}
114
115void CoulombMMM1D::prepare_polygamma_series() {
116 /* polygamma, determine order */
117 double err;
118 auto const rhomax2 = uz2 * far_switch_radius_sq;
119 /* rhomax2 < 1, so rhomax2m2 falls monotonously */
120 int n = 1;
121 auto rhomax2nm2 = 1.0;
122 do {
123 create_mod_psi_up_to(n + 1);
124
125 /* |uz*z| <= 0.5 */
126 err = 2. * static_cast<double>(n) * fabs(mod_psi_even(modPsi, n, 0.5)) *
127 rhomax2nm2;
128 rhomax2nm2 *= rhomax2;
129 n++;
130 } while (err > 0.1 * maxPWerror);
131}
132
133CoulombMMM1D::CoulombMMM1D(double prefactor, double maxPWerror,
134 double switch_rad, int tune_timings,
135 bool tune_verbose)
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.} {
140 if (maxPWerror <= 0.) {
141 throw std::domain_error("Parameter 'maxPWerror' must be > 0");
142 }
143 if (far_switch_radius <= 0. and far_switch_radius != -1.) {
144 throw std::domain_error("Parameter 'far_switch_radius' must be > 0");
145 }
146 if (far_switch_radius > 0.) {
147 far_switch_radius_sq = Utils::sqr(far_switch_radius);
148 }
149 if (tune_timings <= 0) {
150 throw std::domain_error("Parameter 'timings' must be > 0");
151 }
152}
153
154void CoulombMMM1D::sanity_checks_periodicity() const {
155 auto const &box_geo = *get_system().box_geo;
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)");
158 }
159}
160
161void CoulombMMM1D::sanity_checks_cell_structure() const {
162 auto const &local_geo = *get_system().local_geo;
163 if (local_geo.cell_structure_type() != CellStructureType::NSQUARE) {
164 throw std::runtime_error("MMM1D requires the N-square cellsystem");
165 }
166}
167
168void CoulombMMM1D::recalc_boxl_parameters() {
169 auto const &box_geo = *get_system().box_geo;
170
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]);
173
174 uz2 = Utils::sqr(box_geo.length_inv()[2]);
175 prefuz2 = prefactor * uz2;
176 prefL3_i = prefuz2 * box_geo.length_inv()[2];
177
178 determine_bessel_radii();
179 prepare_polygamma_series();
180}
181
183 double dist) const {
184 auto constexpr c_2pi = 2. * std::numbers::pi;
185 auto const &box_geo = *get_system().box_geo;
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];
190 Utils::Vector3d force;
191
192 if (rxy2 <= far_switch_radius_sq) {
193 /* polygamma summation */
194 auto sr = 0.;
195 auto sz = mod_psi_odd(modPsi, 0, z_d);
196 auto r2nm1 = 1.;
197 for (int n = 1; n < n_modPsi; n++) {
198 auto const deriv = static_cast<double>(2 * n);
199 auto const mpe = mod_psi_even(modPsi, n, z_d);
200 auto const mpo = mod_psi_odd(modPsi, n, z_d);
201 auto const r2n = r2nm1 * rxy2_d;
202
203 sz += r2n * mpo;
204 sr += deriv * r2nm1 * mpe;
205
206 if (fabs(deriv * r2nm1 * mpe) < maxPWerror)
207 break;
208
209 r2nm1 = r2n;
210 }
211
212 double Fx = prefL3_i * sr * d[0];
213 double Fy = prefL3_i * sr * d[1];
214 double Fz = prefuz2 * sz;
215
216 /* real space parts */
217
218 double pref, rt, rt2, shift_z;
219
220 pref = 1. / Utils::int_pow<3>(dist);
221 Fx += pref * d[0];
222 Fy += pref * d[1];
223 Fz += pref * d[2];
224
225 shift_z = d[2] + box_geo.length()[2];
226 rt2 = rxy2 + shift_z * shift_z;
227 rt = sqrt(rt2);
228 pref = 1. / (rt2 * rt);
229 Fx += pref * d[0];
230 Fy += pref * d[1];
231 Fz += pref * shift_z;
232
233 shift_z = d[2] - box_geo.length()[2];
234 rt2 = rxy2 + shift_z * shift_z;
235 rt = sqrt(rt2);
236 pref = 1. / (rt2 * rt);
237 Fx += pref * d[0];
238 Fy += pref * d[1];
239 Fz += pref * shift_z;
240
241 force = {Fx, Fy, Fz};
242 } else {
243 /* far range formula */
244 auto const rxy = sqrt(rxy2);
245 auto const rxy_d = rxy * box_geo.length_inv()[2];
246 auto sr = 0., sz = 0.;
247
248 for (int bp = 1; bp < MAXIMAL_B_CUT; bp++) {
249 if (bessel_radii[bp - 1] < rxy)
250 break;
251
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);
256#else
257 auto const [k0, k1] = LPK01(fq * rxy_d);
258#endif
259 sr += bp * k1 * cos(fq * z_d);
260 sz += bp * k0 * sin(fq * z_d);
261 }
262 sr *= uz2 * 4. * c_2pi;
263 sz *= uz2 * 4. * c_2pi;
264
265 auto const pref = sr / rxy + 2. * box_geo.length_inv()[2] / rxy2;
266
267 force = {pref * d[0], pref * d[1], sz};
268 }
269
270 return (prefactor * q1q2) * force;
271}
272
273double CoulombMMM1D::pair_energy(double const q1q2, Utils::Vector3d const &d,
274 double const dist) const {
275 if (q1q2 == 0.)
276 return 0.;
277
278 auto constexpr c_2pi = 2. * std::numbers::pi;
279 auto const &box_geo = *get_system().box_geo;
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];
284 double energy;
285
286 if (rxy2 <= far_switch_radius_sq) {
287 /* near range formula */
288 energy = -2. * std::numbers::egamma;
289
290 /* polygamma summation */
291 double r2n = 1.;
292 for (int n = 0; n < n_modPsi; n++) {
293 auto const add = mod_psi_even(modPsi, n, z_d) * r2n;
294 energy -= add;
295
296 if (fabs(add) < maxPWerror)
297 break;
298
299 r2n *= rxy2_d;
300 }
301 energy *= box_geo.length_inv()[2];
302
303 /* real space parts */
304
305 double rt, shift_z;
306
307 energy += 1. / dist;
308
309 shift_z = d[2] + box_geo.length()[2];
310 rt = sqrt(rxy2 + shift_z * shift_z);
311 energy += 1. / rt;
312
313 shift_z = d[2] - box_geo.length()[2];
314 rt = sqrt(rxy2 + shift_z * shift_z);
315 energy += 1. / rt;
316 } else {
317 /* far range formula */
318 auto const rxy = sqrt(rxy2);
319 auto const rxy_d = rxy * box_geo.length_inv()[2];
320 /* The first Bessel term will compensate a little bit the
321 log term, so add them close together */
322 energy =
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)
326 break;
327
328 auto const fq = c_2pi * bp;
329 energy += K0(fq * rxy_d) * cos(fq * z_d);
330 }
331 energy *= 4. * box_geo.length_inv()[2];
332 }
333
334 return prefactor * q1q2 * energy;
335}
336
338 if (is_tuned()) {
339 return;
340 }
341 recalc_boxl_parameters();
342 auto &system = get_system();
343
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();
348 auto min_rad = -1.;
349 auto switch_radius = 0.2 * maxrad;
350 /* determine optimal switching radius. Should be around 0.33 */
351 while (switch_radius < 0.4 * maxrad) {
352 if (switch_radius > bessel_radii.back()) {
353 // this switching radius is large enough for our Bessel series
354 far_switch_radius_sq = Utils::sqr(switch_radius);
355 system.on_coulomb_change();
356
357 /* perform force calculation test */
358 auto const int_time = benchmark_integration_step(system, tune_timings);
359
360 if (tune_verbose) {
361 std::printf("r= %f t= %f ms\n", switch_radius, int_time);
362 }
363
364 if (int_time < min_time) {
365 min_time = int_time;
366 min_rad = switch_radius;
367 } else if (int_time > 2. * min_time) {
368 // simple heuristic to skip remaining radii when performance drops
369 break;
370 }
371 }
372 switch_radius += 0.025 * maxrad;
373 }
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())) {
377 // this switching radius is too small for our Bessel series
378 throw std::runtime_error("MMM1D could not find a reasonable Bessel cutoff");
379 }
380
381 m_is_tuned = true;
382 system.on_coulomb_change();
383}
384
385#endif // ELECTROSTATICS
@ 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
Definition mmm1d.cpp:102
static auto determine_minrad(double maxPWerror, int P, Utils::Vector3d const &box_l, Utils::Vector3d const &box_l_inv)
Definition mmm1d.cpp:65
static double mod_psi_even(auto const &modPsi, int n, double x)
Modified polygamma for even order 2*n, n >= 0
Definition mmm1d.cpp:97
static auto far_error(int P, double minrad, Utils::Vector3d const &box_l_inv)
Definition mmm1d.cpp:56
MMM1D algorithm for long-range Coulomb interactions on the CPU.
DEVICE_QUALIFIER constexpr T sqr(T x)
Calculates the SQuaRe of x.
Definition sqr.hpp:28
std::pair< double, double > LPK01(double x)
Modified Bessel functions of second kind, order 0 and 1, low precision.
Definition specfunc.cpp:408
double K1(double x)
Modified Bessel function of second kind, order 1.
Definition specfunc.cpp:262
double K0(double x)
Modified Bessel function of second kind, order 0.
Definition specfunc.cpp:250
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.
Definition specfunc.hpp:91
bool tune_verbose
Definition mmm1d.hpp:63
CoulombMMM1D(double prefactor, double maxPWerror, double switch_rad, int tune_timings, bool tune_verbose)
Definition mmm1d.cpp:133
double far_switch_radius
Far switch radius.
Definition mmm1d.hpp:61
Utils::Vector3d pair_force(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair force.
Definition mmm1d.cpp:182
double pair_energy(double q1q2, Utils::Vector3d const &d, double dist) const
Compute the pair energy.
Definition mmm1d.cpp:273
bool is_tuned() const
Definition mmm1d.hpp:84
double maxPWerror
Maximal allowed pairwise error for the potential and force.
Definition mmm1d.hpp:56
int tune_timings
Definition mmm1d.hpp:62
void tune()
Definition mmm1d.cpp:337
double benchmark_integration_step(System::System &system, int int_steps)
Benchmark the integration loop.
Definition tuning.cpp:73