CoulombGalore
coulombgalore.h
1 /*
2  * CoulombGalore - A Library for Electrostatic Interactions
3  *
4  * MIT License
5  * Copyright (c) 2019 Björn Stenqvist and Mikael Lund
6  *
7  * Permission is hereby granted, free of charge, to any person obtaining a copy
8  * of this software and associated documentation files (the "Software"), to deal
9  * in the Software without restriction, including without limitation the rights
10  * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
11  * copies of the Software, and to permit persons to whom the Software is
12  * furnished to do so, subject to the following conditions:
13  *
14  * The above copyright notice and this permission notice shall be included in all
15  * copies or substantial portions of the Software.
16  *
17  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
18  * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
19  * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
20  * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
21  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
22  * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
23  * SOFTWARE.
24  */
25 #pragma once
26 
27 #include <string>
28 #include <limits>
29 #include <cmath>
30 #include <iostream>
31 #include <vector>
32 #include <array>
33 #include <Eigen/Core>
34 
36 namespace CoulombGalore {
37 
39 typedef Eigen::Vector3d vec3;
40 
41 constexpr double infinity = std::numeric_limits<double>::infinity();
42 
44 enum class Scheme {
45  plain,
46  ewald,
47  reactionfield,
48  wolf,
49  poisson,
50  qpotential,
51  fanourgakis,
52  zerodipole,
53  zahn,
54  fennell,
55  qpotential5,
56  spline
57 };
58 
64 inline double powi(double x, int n) {
65 #if defined(__GNUG__)
66  return __builtin_powi(x, n);
67 #else
68  return std::pow(x, n);
69 #endif
70 }
71 
76 inline constexpr unsigned int factorial(unsigned int n) { return n <= 1 ? 1 : n * factorial(n - 1); }
77 
78 constexpr unsigned int binomial(signed int n, signed int k) {
79  return factorial(n) / (factorial(k) * factorial(n - k));
80 }
81 
101 inline double qPochhammerSymbol(double q, int l = 0, int P = 300) {
102  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
103  for (int n = 1; n < P + 1; n++) {
104  double val = 0.0;
105  for (int k = 1; k < n + l + 1; k++)
106  val += powi(q, k - 1);
107  Ct *= val;
108  }
109  double Dt = powi(1.0 - q, P); // (1-q)^P
110  return (Ct * Dt);
111 }
112 
116 inline double qPochhammerSymbolDerivative(double q, int l = 0, int P = 300) {
117  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
118  for (int n = 1; n < P + 1; n++) {
119  double val = 0.0;
120  for (int k = 1; k < n + l + 1; k++)
121  val += powi(q, k - 1);
122  Ct *= val;
123  }
124  double dCt = 0.0; // evaluates to derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
125  for (int n = 1; n < P + 1; n++) {
126  double nom = 0.0;
127  double denom = 1.0;
128  for (int k = 2; k < n + l + 1; k++) {
129  nom += (k - 1) * powi(q, k - 2);
130  denom += powi(q, k - 1);
131  }
132  dCt += nom / denom;
133  }
134  dCt *= Ct;
135  double Dt = powi(1.0 - q, P); // (1-q)^P
136  double dDt = 0.0;
137  if (P > 0)
138  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
139  return (Ct * dDt + dCt * Dt);
140 }
141 
145 inline double qPochhammerSymbolSecondDerivative(double q, int l = 0, int P = 300) {
146  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
147  double DS = 0.0;
148  double dDS = 0.0;
149  for (int n = 1; n < P + 1; n++) {
150  double tmp = 0.0;
151  for (int k = 1; k < n + l + 1; k++)
152  tmp += powi(q, k - 1);
153  Ct *= tmp;
154  double nom = 0.0;
155  double denom = 1.0;
156  for (int k = 2; k < n + l + 1; k++) {
157  nom += (k - 1) * powi(q, k - 2);
158  denom += powi(q, k - 1);
159  }
160  DS += nom / denom;
161  double diffNom = 0.0;
162  double diffDenom = 1.0;
163  for (int k = 3; k < n + l + 1; k++) {
164  diffNom += (k - 1) * (k - 2) * powi(q, k - 3);
165  diffDenom += (k - 1) * powi(q, k - 2);
166  }
167  dDS += (diffNom * denom - nom * diffDenom) / denom / denom;
168  }
169  double dCt = Ct * DS; // derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
170  double ddCt = dCt * DS + Ct * dDS; // second derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
171  double Dt = powi(1.0 - q, P); // (1-q)^P
172  double dDt = 0.0;
173  if (P > 0)
174  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
175  double ddDt = 0.0;
176  if (P > 1)
177  ddDt = P * (P - 1) * powi(1 - q, P - 2); // second derivative of (1-q)^P
178  return (Ct * ddDt + 2 * dCt * dDt + ddCt * Dt);
179 }
180 
184 inline double qPochhammerSymbolThirdDerivative(double q, int l = 0, int P = 300) {
185  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
186  double DS = 0.0;
187  double dDS = 0.0;
188  double ddDS = 0.0;
189  for (int n = 1; n < P + 1; n++) {
190  double tmp = 0.0;
191  for (int k = 1; k < n + l + 1; k++)
192  tmp += powi(q, k - 1);
193  Ct *= tmp;
194  double f = 0.0;
195  double g = 1.0;
196  for (int k = 2; k < n + l + 1; k++) {
197  f += (k - 1) * powi(q, k - 2);
198  g += powi(q, k - 1);
199  }
200  DS += f / g;
201  double df = 0.0;
202  double dg = 0.0;
203  if (n + l > 1)
204  dg = 1.0;
205  for (int k = 3; k < n + l + 1; k++) {
206  df += (k - 1) * (k - 2) * powi(q, k - 3);
207  dg += (k - 1) * powi(q, k - 2);
208  }
209  dDS += (df * g - f * dg) / g / g;
210  double ddf = 0.0;
211  double ddg = 0.0;
212  if (n + l > 2)
213  ddg = 2.0;
214  for (int k = 4; k < n + l + 1; k++) {
215  ddf += (k - 1) * (k - 2) * (k - 3) * powi(q, k - 4);
216  ddg += (k - 1) * (k - 2) * powi(q, k - 3);
217  }
218  ddDS += (ddf * g * g - 2.0 * df * dg * g + 2.0 * f * dg * dg - f * ddg * g) / g / g / g;
219  }
220  double dCt = Ct * DS; // derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
221  double ddCt = dCt * DS + Ct * dDS; // second derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
222  double dddCt = ddCt * DS + 2.0 * dCt * dDS + Ct * ddDS; // third derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
223  double Dt = powi(1.0 - q, P); // (1-q)^P
224  double dDt = 0.0;
225  if (P > 0)
226  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
227  double ddDt = 0.0;
228  if (P > 1)
229  ddDt = P * (P - 1) * powi(1 - q, P - 2); // second derivative of (1-q)^P
230  double dddDt = 0.0;
231  if (P > 2)
232  dddDt = -P * (P - 1) * (P - 2) * powi(1 - q, P - 3); // third derivative of (1-q)^P
233  return (dddCt * Dt + 3.0 * ddCt * dDt + 3 * dCt * ddDt + Ct * dddDt);
234 }
235 
236 namespace Tabulate {
237 
238 /* base class for all tabulators - no dependencies */
239 template <typename T = double> class TabulatorBase {
240  protected:
241  T utol = 1e-5, ftol = -1, umaxtol = -1, fmaxtol = -1;
242  T numdr = 0.0001; // dr for derivative evaluation
243 
244  // First derivative with respect to x
245  T f1(std::function<T(T)> f, T x) const { return (f(x + numdr * 0.5) - f(x - numdr * 0.5)) / (numdr); }
246 
247  // Second derivative with respect to x
248  T f2(std::function<T(T)> f, T x) const { return (f1(f, x + numdr * 0.5) - f1(f, x - numdr * 0.5)) / (numdr); }
249 
250  void check() const {
251  if (ftol != -1 && ftol <= 0.0) {
252  std::cerr << "ftol=" << ftol << " too small\n" << std::endl;
253  abort();
254  }
255  if (umaxtol != -1 && umaxtol <= 0.0) {
256  std::cerr << "umaxtol=" << umaxtol << " too small\n" << std::endl;
257  abort();
258  }
259  if (fmaxtol != -1 && fmaxtol <= 0.0) {
260  std::cerr << "fmaxtol=" << fmaxtol << " too small\n" << std::endl;
261  abort();
262  }
263  }
264 
265  public:
266  struct data {
267  std::vector<T> r2; // r2 for intervals
268  std::vector<T> c; // c for coefficents
269  T rmin2 = 0, rmax2 = 0; // useful to save these with table
270  bool empty() const { return r2.empty() && c.empty(); }
271  };
272 
273  void setTolerance(T _utol, T _ftol = -1, T _umaxtol = -1, T _fmaxtol = -1) {
274  utol = _utol;
275  ftol = _ftol;
276  umaxtol = _umaxtol;
277  fmaxtol = _fmaxtol;
278  }
279 
280  void setNumdr(T _numdr) { numdr = _numdr; }
281 };
282 
283 /*
284  * @brief Andrea table with logarithmic search
285  *
286  * Tabulator with logarithmic search.
287  * Code mainly from MolSim (Per Linse) with some upgrades
288  * Reference: doi:10/frzp4d
289  *
290  * @note Slow on Intel compiler
291  */
292 template <typename T = double> class Andrea : public TabulatorBase<T> {
293  private:
294  typedef TabulatorBase<T> base; // for convenience
295  int mngrid = 1200; // Max number of controlpoints
296  int ndr = 100; // Max number of trials to decr dr
297  T drfrac = 0.9; // Multiplicative factor to decr dr
298 
299  std::vector<T> SetUBuffer(T, T zlow, T, T zupp, T u0low, T u1low, T u2low, T u0upp, T u1upp, T u2upp) {
300 
301  // Zero potential and force return no coefficients
302  if (std::fabs(u0low) < 1e-9)
303  if (std::fabs(u1low) < 1e-9)
304  return {0, 0, 0, 0, 0, 0, 0};
305 
306  T dz1 = zupp - zlow;
307  T dz2 = dz1 * dz1;
308  T dz3 = dz2 * dz1;
309  T w0low = u0low;
310  T w1low = u1low;
311  T w2low = u2low;
312  T w0upp = u0upp;
313  T w1upp = u1upp;
314  T w2upp = u2upp;
315  T c0 = w0low;
316  T c1 = w1low;
317  T c2 = w2low * 0.5;
318  T a = 6 * (w0upp - c0 - c1 * dz1 - c2 * dz2) / dz3;
319  T b = 2 * (w1upp - c1 - 2 * c2 * dz1) / dz2;
320  T c = (w2upp - 2 * c2) / dz1;
321  T c3 = (10 * a - 12 * b + 3 * c) / 6;
322  T c4 = (-15 * a + 21 * b - 6 * c) / (6 * dz1);
323  T c5 = (2 * a - 3 * b + c) / (2 * dz2);
324 
325  return {zlow, c0, c1, c2, c3, c4, c5};
326  }
327 
328  /*
329  * @returns boolean vector.
330  * - `[0]==true`: tolerance is approved,
331  * - `[1]==true` Repulsive part is found.
332  */
333  std::vector<bool> CheckUBuffer(std::vector<T> &ubuft, T rlow, T rupp, std::function<T(T)> f) const {
334 
335  // Number of points to control
336  int ncheck = 11;
337  T dr = (rupp - rlow) / (ncheck - 1);
338  std::vector<bool> vb(2, false);
339 
340  for (int i = 0; i < ncheck; i++) {
341  T r1 = rlow + dr * ((T)i);
342  T r2 = r1 * r1;
343  T u0 = f(r2);
344  T u1 = base::f1(f, r2);
345  T dz = r2 - rlow * rlow;
346  T usum =
347  ubuft.at(1) +
348  dz * (ubuft.at(2) + dz * (ubuft.at(3) + dz * (ubuft.at(4) + dz * (ubuft.at(5) + dz * ubuft.at(6)))));
349 
350  T fsum = ubuft.at(2) +
351  dz * (2 * ubuft.at(3) + dz * (3 * ubuft.at(4) + dz * (4 * ubuft.at(5) + dz * (5 * ubuft.at(6)))));
352 
353  if (std::fabs(usum - u0) > base::utol)
354  return vb;
355  if (base::ftol != -1 && std::fabs(fsum - u1) > base::ftol)
356  return vb;
357  if (base::umaxtol != -1 && std::fabs(usum) > base::umaxtol)
358  vb[1] = true;
359  if (base::fmaxtol != -1 && std::fabs(usum) > base::fmaxtol)
360  vb[1] = true;
361  }
362  vb[0] = true;
363  return vb;
364  }
365 
366  public:
367  /*
368  * @brief Get tabulated value at f(x)
369  * @param d Table data
370  * @param r2 value
371  */
372  inline T eval(const typename base::data &d, T r2) const {
373  size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1;
374  size_t pos6 = 6 * pos;
375  assert((pos6 + 5) < d.c.size() && "out of bounds");
376  T dz = r2 - d.r2[pos];
377  return d.c[pos6] +
378  dz * (d.c[pos6 + 1] +
379  dz * (d.c[pos6 + 2] + dz * (d.c[pos6 + 3] + dz * (d.c[pos6 + 4] + dz * (d.c[pos6 + 5])))));
380  }
381 
382  /*
383  * @brief Get tabulated value at df(x)/dx
384  * @param d Table data
385  * @param r2 value
386  */
387  T evalDer(const typename base::data &d, T r2) const {
388  size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1;
389  size_t pos6 = 6 * pos;
390  T dz = r2 - d.r2[pos];
391  return (d.c[pos6 + 1] +
392  dz * (2.0 * d.c[pos6 + 2] +
393  dz * (3.0 * d.c[pos6 + 3] + dz * (4.0 * d.c[pos6 + 4] + dz * (5.0 * d.c[pos6 + 5])))));
394  }
395 
399  typename base::data generate(std::function<T(T)> f, double rmin, double rmax) {
400  rmin = std::sqrt(rmin);
401  rmax = std::sqrt(rmax);
402  base::check();
403  typename base::data td;
404  td.rmin2 = rmin * rmin;
405  td.rmax2 = rmax * rmax;
406 
407  T rumin = rmin;
408  T rmax2 = rmax * rmax;
409  T dr = rmax - rmin;
410  T rupp = rmax;
411  T zupp = rmax2;
412  bool repul = false; // Stop tabulation if repul is true
413 
414  td.r2.push_back(zupp);
415 
416  int i;
417  for (i = 0; i < mngrid; i++) {
418  T rlow = rupp;
419  T zlow;
420  std::vector<T> ubuft;
421  int j;
422 
423  dr = (rupp - rmin);
424 
425  for (j = 0; j < ndr; j++) {
426  zupp = rupp * rupp;
427  rlow = rupp - dr;
428  if (rumin > rlow)
429  rlow = rumin;
430 
431  zlow = rlow * rlow;
432 
433  T u0low = f(zlow);
434  T u1low = base::f1(f, zlow);
435  T u2low = base::f2(f, zlow);
436  T u0upp = f(zupp);
437  T u1upp = base::f1(f, zupp);
438  T u2upp = base::f2(f, zupp);
439 
440  ubuft = SetUBuffer(rlow, zlow, rupp, zupp, u0low, u1low, u2low, u0upp, u1upp, u2upp);
441  std::vector<bool> vb = CheckUBuffer(ubuft, rlow, rupp, f);
442  repul = vb[1];
443  if (vb[0]) {
444  rupp = rlow;
445  break;
446  }
447  dr *= drfrac;
448  }
449 
450  if (j >= ndr)
451  throw std::runtime_error("Andrea spline: try to increase utol/ftol");
452  if (ubuft.size() != 7)
453  throw std::runtime_error("Andrea spline: wrong size of ubuft, min value + 6 coefficients");
454 
455  td.r2.push_back(zlow);
456  for (size_t k = 1; k < ubuft.size(); k++)
457  td.c.push_back(ubuft.at(k));
458 
459  // Entered a highly repulsive part, stop tabulation
460  if (repul) {
461  rumin = rlow;
462  td.rmin2 = rlow * rlow;
463  }
464  if (rlow <= rumin || repul)
465  break;
466  }
467 
468  if (i >= mngrid)
469  throw std::runtime_error("Andrea spline: try to increase utol/ftol");
470 
471  // create final reversed c and r2
472 #if __cplusplus >= 201703L
473  // C++17 only code
474  assert(td.c.size() % 6 == 0);
475  assert(td.c.size() / (td.r2.size() - 1) == 6);
476  assert(std::is_sorted(td.r2.rbegin(), td.r2.rend()));
477  std::reverse(td.r2.begin(), td.r2.end()); // reverse all elements
478  for (size_t i = 0; i < td.c.size() / 2; i += 6) // reverse knot order in packets of six
479  std::swap_ranges(td.c.begin() + i, td.c.begin() + i + 6, td.c.end() - i - 6); // c++17 only
480  return td;
481 #else
482  typename base::data tdsort;
483  tdsort.rmax2 = td.rmax2;
484  tdsort.rmin2 = td.rmin2;
485 
486  // reverse copy all elements in r2
487  tdsort.r2.resize(td.r2.size());
488  std::reverse_copy(td.r2.begin(), td.r2.end(), tdsort.r2.begin());
489 
490  // sanity check before reverse knot copy
491  assert(std::is_sorted(td.r2.rbegin(), td.r2.rend()));
492  assert(td.c.size() % 6 == 0);
493  assert(td.c.size() / (td.r2.size() - 1) == 6);
494 
495  // reverse copy knots
496  tdsort.c.resize(td.c.size());
497  auto dst = tdsort.c.end();
498  for (auto src = td.c.begin(); src != td.c.end(); src += 6)
499  std::copy(src, src + 6, dst -= 6);
500  return tdsort;
501 #endif
502  }
503 };
504 
505 } // namespace Tabulate
506 
518 class SchemeBase {
519  protected:
520  double invcutoff; // inverse cutoff distance, UNIT: [ ( input length )^-1 ]
521  double cutoff2; // square cutoff distance, UNIT: [ ( input length )^2 ]
522  double kappa; // inverse Debye-length, UNIT: [ ( input length )^-1 ]
523  double T0; // Spatial Fourier transformed modified interaction tensor, used to calculate the dielectric constant, UNIT: [ 1 ]
524  double chi; // Negative integrated volume potential, used to neutralize charged system, UNIT: [ ( input length )^2 ]
525  bool dipolar_selfenergy; // is there a valid dipolar self-energy?
526  std::array<double, 2> self_energy_prefactor; // Prefactor for self-energies, UNIT: [ 1 ]
527 
528  public:
529  std::string doi;
530  std::string name;
532  double cutoff;
533  double debye_length;
534 
535  inline SchemeBase(Scheme scheme, double cutoff, double debye_length = infinity)
536  : scheme(scheme), cutoff(cutoff), debye_length(debye_length) {}
537 
538  virtual ~SchemeBase() = default;
539 
540  virtual double neutralization_energy(const std::vector<double> &, double) const { return 0.0; }
541 
556  double calc_dielectric(double M2V) { return (M2V * T0 + 2.0 * M2V + 1.0) / (M2V * T0 - M2V + 1.0); }
557 
558  virtual double short_range_function(double q) const = 0;
559  virtual double short_range_function_derivative(double q) const = 0;
560  virtual double short_range_function_second_derivative(double q) const = 0;
561  virtual double short_range_function_third_derivative(double q) const = 0;
562 
563  virtual double ion_ion_energy(double, double, double) const = 0;
564  virtual double ion_dipole_energy(double, const vec3 &, const vec3 &) const = 0;
565  virtual double dipole_dipole_energy(const vec3 &, const vec3 &, const vec3 &) const = 0;
566  virtual double multipole_multipole_energy(double, double, const vec3 &, const vec3 &, const vec3 &) const = 0;
567 
568  virtual vec3 ion_field(double, const vec3 &) const = 0;
569  virtual vec3 dipole_field(const vec3 &, const vec3 &) const = 0;
570  virtual vec3 multipole_field(double, const vec3 &, const vec3 &) const = 0;
571 
572  virtual vec3 ion_ion_force(double, double, const vec3 &) const = 0;
573  virtual vec3 ion_dipole_force(double, const vec3 &, const vec3 &) const = 0;
574  virtual vec3 dipole_dipole_force(const vec3 &, const vec3 &, const vec3 &) const = 0;
575  virtual vec3 multipole_multipole_force(double, double, const vec3 &, const vec3 &, const vec3 &) const = 0;
576 
577  virtual double ion_potential(double, double) const = 0;
578  virtual double dipole_potential(const vec3 &, const vec3 &) const = 0;
579 
580  // add remaining funtions here...
581 
582  // virtual double surface_energy(const std::vector<vec3> &, const std::vector<double> &, const std::vector<vec3> &,
583  // double) const { return 0.0; } virtual vec3 surface_force(const std::vector<vec3> &, const std::vector<double> &,
584  // const std::vector<vec3> &, int, double) const { return {0.0,0.0,0.0}; }
585 
586 #ifdef NLOHMANN_JSON_HPP
587  private:
588  virtual void _to_json(nlohmann::json &) const = 0;
589 
590  public:
591  inline void to_json(nlohmann::json &j) const {
592  _to_json(j);
593  if (std::isfinite(cutoff))
594  j["cutoff"] = cutoff;
595  if (not doi.empty())
596  j["doi"] = doi;
597  if (not name.empty())
598  j["type"] = name;
599  if (std::isfinite(debye_length))
600  j["debyelength"] = debye_length;
601  }
602 #endif
603 };
604 
610 template <class T, bool debyehuckel = true> class EnergyImplementation : public SchemeBase {
611  public:
612  EnergyImplementation(Scheme type, double cutoff, double debyelength = infinity)
613  : SchemeBase(type, cutoff, debyelength) {
614  invcutoff = 1.0 / cutoff;
615  cutoff2 = cutoff * cutoff;
616  kappa = 1.0 / debyelength;
617  }
618 
630  inline double ion_potential(double z, double r) const override {
631  if (r < cutoff) {
632  double q = r * invcutoff;
633  if (debyehuckel) // determined at compile time
634  return z / r * static_cast<const T *>(this)->short_range_function(q) * std::exp(-kappa * r);
635  else
636  return z / r * static_cast<const T *>(this)->short_range_function(q);
637  } else {
638  return 0.0;
639  }
640  }
641 
654  inline double dipole_potential(const vec3 &mu, const vec3 &r) const override {
655  double r2 = r.squaredNorm();
656  if (r2 < cutoff2) {
657  double r1 = std::sqrt(r2);
658  double q = r1 * invcutoff;
659  return mu.dot(r) / (r2 * r1) *
660  (static_cast<const T *>(this)->short_range_function(q) * (1.0 + kappa * r1) -
661  q * static_cast<const T *>(this)->short_range_function_derivative(q)) *
662  std::exp(-kappa * r1);
663  } else {
664  return 0.0;
665  }
666  }
667 
679  inline vec3 ion_field(double z, const vec3 &r) const {
680  double r2 = r.squaredNorm();
681  if (r2 < cutoff2) {
682  double r1 = std::sqrt(r2);
683  double q = r1 * invcutoff;
684  return z * r / (r2 * r1) *
685  (static_cast<const T *>(this)->short_range_function(q) * (1.0 + kappa * r1) -
686  q * static_cast<const T *>(this)->short_range_function_derivative(q)) *
687  std::exp(-kappa * r1);
688  } else {
689  return {0, 0, 0};
690  }
691  }
692 
706  inline vec3 dipole_field(const vec3 &mu, const vec3 &r) const {
707  double r2 = r.squaredNorm();
708  if (r2 < cutoff2) {
709  double r1 = std::sqrt(r2);
710  double r3 = r1 * r2;
711  double q = r1 * invcutoff;
712  double q2 = q * q;
713  double kappa2 = kappa * kappa;
714  double kappa_x_r1 = kappa * r1;
715  double srf = static_cast<const T *>(this)->short_range_function(q);
716  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
717  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
718  vec3 fieldD = (3.0 * mu.dot(r) * r / r2 - mu) / r3;
719  fieldD *= (srf * (1.0 + kappa_x_r1 + kappa2 * r2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kappa_x_r1) +
720  q2 / 3.0 * ddsrf);
721  vec3 fieldI = mu / r3 * (srf * kappa2 * r2 - 2.0 * kappa_x_r1 * q * dsrf + ddsrf * q2) / 3.0;
722  return (fieldD + fieldI) * std::exp(-kappa_x_r1);
723  } else {
724  return {0, 0, 0};
725  }
726  }
727 
743  inline vec3 multipole_field(double z, const vec3 &mu, const vec3 &r) const {
744  double r2 = r.squaredNorm();
745  if (r2 < cutoff2) {
746  double r1 = std::sqrt(r2);
747  double q = r1 * invcutoff;
748  double r3 = r1 * r2;
749  double q2 = q * q;
750  double kappa2 = kappa * kappa;
751  double kappa_x_r1 = kappa * r1;
752  double srf = static_cast<const T *>(this)->short_range_function(q);
753  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
754  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
755  vec3 fieldIon = z * r / r3 * ( srf * (1.0 + kappa_x_r1) - q * dsrf ); // field from ion
756  vec3 fieldD = (3.0 * mu.dot(r) * r / r2 - mu) / r3;
757  fieldD *= (srf * (1.0 + kappa_x_r1 + kappa2 * r2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kappa_x_r1) + q2 / 3.0 * ddsrf);
758  vec3 fieldI = mu / r3 * (srf * kappa2 * r2 - 2.0 * kappa_x_r1 * q * dsrf + ddsrf * q2) / 3.0;
759  return (fieldD + fieldI + fieldIon) * std::exp(-kappa_x_r1);
760  } else {
761  return {0, 0, 0};
762  }
763  }
764 
778  inline double ion_ion_energy(double zA, double zB, double r) const override { return zB * ion_potential(zA, r); }
779 
798  inline double ion_dipole_energy(double z, const vec3 &mu, const vec3 &r) const override {
799  // Both expressions below gives same answer. Keep for possible optimization in future.
800  // return -mu.dot(ion_field(z,r)); // field from charge interacting with dipole
801  return z * dipole_potential(mu, -r); // potential of dipole interacting with charge
802  }
803 
818  inline double dipole_dipole_energy(const vec3 &muA, const vec3 &muB, const vec3 &r) const override {
819  return -muA.dot(dipole_field(muB, r));
820  }
821 
833  inline double multipole_multipole_energy(double zA, double zB, const vec3 &muA, const vec3 &muB,
834  const vec3 &r) const override {
835  double r2 = r.squaredNorm();
836  if (r2 < cutoff2) {
837  double r1 = std::sqrt(r2);
838  double r3 = r1 * r2;
839  double q = r1 / cutoff;
840  double kappa_r1 = kappa * r1;
841 
842  double srf = static_cast<const T *>(this)->short_range_function(q);
843  double dsrfq = static_cast<const T *>(this)->short_range_function_derivative(q) * q;
844  double ddsrfq2 = static_cast<const T *>(this)->short_range_function_second_derivative(q) * q * q / 3.0;
845 
846  double tmp1 = (srf * (1.0 + kappa_r1) - dsrfq) / r3;
847  double tmp2 = (srf * kappa_r1 * kappa_r1 / 3.0 - dsrfq * (2.0 / 3.0 * kappa_r1) + ddsrfq2) / r3;
848 
849  vec3 field_dipoleB = (3.0 * muB.dot(r) * r / r2 - muB) * (tmp1 + tmp2) + muB * tmp2;
850 
851  double ion_ion = zA * zB / r1 * srf;
852  double ion_dipole = (zB * muA.dot(r) + zA * muB.dot(-r)) * tmp1;
853  double dipole_dipole = -muA.dot(field_dipoleB);
854  return (ion_ion + ion_dipole + dipole_dipole) * std::exp(-kappa_r1);
855  } else {
856  return 0.0;
857  }
858  }
859 
873  inline vec3 ion_ion_force(double zA, double zB, const vec3 &r) const override { return zB * ion_field(zA, r); }
874 
888  inline vec3 ion_dipole_force(double z, const vec3 &mu, const vec3 &r) const override {
889  return z * dipole_field(mu, r);
890  }
891 
918  inline vec3 dipole_dipole_force(const vec3 &muA, const vec3 &muB, const vec3 &r) const override {
919  double r2 = r.squaredNorm();
920  if (r2 < cutoff2) {
921  double r1 = std::sqrt(r2);
922  vec3 rh = r / r1;
923  double q = r1 * invcutoff;
924  double q2 = q * q;
925  double r4 = r2 * r2;
926  double muAdotRh = muA.dot(rh);
927  double muBdotRh = muB.dot(rh);
928  vec3 forceD =
929  3.0 * ((5.0 * muAdotRh * muBdotRh - muA.dot(muB)) * rh - muBdotRh * muA - muAdotRh * muB) / r4;
930  double srf = static_cast<const T *>(this)->short_range_function(q);
931  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
932  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
933  double dddsrf = static_cast<const T *>(this)->short_range_function_third_derivative(q);
934  forceD *= (srf * (1.0 + kappa * r1 + kappa * kappa * r2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kappa * r1) +
935  q2 / 3.0 * ddsrf);
936  vec3 forceI = muAdotRh * muBdotRh * rh / r4;
937  forceI *=
938  (srf * (1.0 + kappa * r1) * kappa * kappa * r2 - q * dsrf * (3.0 * kappa * r1 + 2.0) * kappa * r1 +
939  ddsrf * (1.0 + 3.0 * kappa * r1) * q2 - q2 * q * dddsrf);
940  return (forceD + forceI) * std::exp(-kappa * r1);
941  } else {
942  return {0, 0, 0};
943  }
944  }
945 
957  inline vec3 multipole_multipole_force(double zA, double zB, const vec3 &muA, const vec3 &muB,
958  const vec3 &r) const override {
959  double r2 = r.squaredNorm();
960  if (r2 < cutoff2) {
961  double r1 = std::sqrt(r2);
962  double q = r1 * invcutoff;
963  double q2 = q * q;
964  double kr = kappa * r1;
965  vec3 rh = r / r1;
966  double muAdotRh = muA.dot(rh);
967  double muBdotRh = muB.dot(rh);
968 
969  double srf = static_cast<const T *>(this)->short_range_function(q);
970  double dsrfq = static_cast<const T *>(this)->short_range_function_derivative(q) * q;
971  double ddsrfq2 = static_cast<const T *>(this)->short_range_function_second_derivative(q) * q2;
972  double dddsrfq3 = static_cast<const T *>(this)->short_range_function_third_derivative(q) * q2 * q;
973 
974  double tmp0 = (srf * (1.0 + kr) - dsrfq);
975  double tmp1 = (srf * kr * kr - 2.0 * kr * dsrfq + ddsrfq2) / 3.0;
976  double tmp2 = tmp1 + tmp0;
977  double tmp3 = (srf * (1.0 + kr) * kr * kr - dsrfq * (2.0 * (1.0 + kr) + kr) * kr +
978  ddsrfq2 * (1.0 + 3.0 * kr) - dddsrfq3) /
979  r1;
980 
981  vec3 ion_ion = zB * zA * r * (srf * (1.0 + kr) - dsrfq);
982  // vec3 ion_ion = zB * zA * r * tmp0;
983 
984  vec3 ion_dipoleA = zA * ((3.0 * muBdotRh * rh - muB) * tmp2 + muB * tmp1);
985  vec3 ion_dipoleB = zB * ((3.0 * muAdotRh * rh - muA) * tmp2 + muA * tmp1);
986 
987  vec3 forceD =
988  3.0 * ((5.0 * muAdotRh * muBdotRh - muA.dot(muB)) * rh - muBdotRh * muA - muAdotRh * muB) * tmp2 / r1;
989  vec3 dipole_dipole = (forceD + muAdotRh * muBdotRh * rh * tmp3);
990  return (ion_ion + ion_dipoleA + ion_dipoleB + dipole_dipole) * std::exp(-kr) / r2 / r1;
991  } else {
992  return {0, 0, 0};
993  }
994  }
995 
1009  inline vec3 dipole_torque(const vec3 &mu, const vec3 &E) const { return mu.cross(E); }
1010 
1023  inline double self_energy(const std::array<double, 2> &m2) const {
1024  static_assert(decltype(m2){0}.size() == decltype(self_energy_prefactor){0}.size(), "static assert");
1025  double e_self = 0.0;
1026  for (int i = 0; i < (int)m2.size(); i++)
1027  e_self += self_energy_prefactor[i] * m2[i] * powi(invcutoff, 2 * i + 1);
1028  return e_self;
1029  }
1030 
1039  inline double neutralization_energy(const std::vector<double> &charges, double volume) const override {
1040  double squaredSumQ = 0.0;
1041  for (unsigned int i = 0; i < charges.size(); i++)
1042  squaredSumQ += charges.at(i);
1043  return ((this)->chi / 2.0 / volume * squaredSumQ);
1044  }
1045 };
1046 
1047 // -------------- Plain ---------------
1048 
1053 class Plain : public EnergyImplementation<Plain> {
1054  public:
1055  inline Plain(double debye_length = infinity)
1056  : EnergyImplementation(Scheme::plain, std::numeric_limits<double>::max(), debye_length) {
1057  name = "plain";
1058  dipolar_selfenergy = true;
1059  doi = "Premier mémoire sur l’électricité et le magnétisme by Charles-Augustin de Coulomb"; // :P
1060  self_energy_prefactor = {0.0, 0.0};
1061  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1062  chi = -2.0 * std::acos(-1.0) * cutoff2; // should not be used!
1063  };
1064  inline double short_range_function(double) const override { return 1.0; };
1065  inline double short_range_function_derivative(double) const override { return 0.0; }
1066  inline double short_range_function_second_derivative(double) const override { return 0.0; }
1067  inline double short_range_function_third_derivative(double) const override { return 0.0; }
1068 #ifdef NLOHMANN_JSON_HPP
1069  inline Plain(const nlohmann::json &j) : Plain(j.value("debyelength", infinity)) {}
1070 
1071  private:
1072  inline void _to_json(nlohmann::json &) const override {}
1073 #endif
1074 };
1075 
1076 // -------------- Ewald real-space ---------------
1077 
1085 class Ewald : public EnergyImplementation<Ewald> {
1086  double alpha, alpha2;
1087  double alphaRed, alphaRed2, alphaRed3;
1088  double eps_sur;
1089  double debye_length;
1090  double kappa, kappa2;
1091  double beta, beta2, beta3;
1092  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1093  const double pi = 4.0 * std::atan(1.0);
1094 
1095  public:
1100  inline Ewald(double cutoff, double alpha, double eps_sur = infinity, double debye_length = infinity)
1101  : EnergyImplementation(Scheme::ewald, cutoff), alpha(alpha), eps_sur(eps_sur) {
1102  name = "Ewald real-space";
1103  dipolar_selfenergy = true;
1104  doi = "10.1002/andp.19213690304";
1105  alpha2 = alpha * alpha;
1106  alphaRed = alpha * cutoff;
1107  alphaRed2 = alphaRed * alphaRed;
1108  alphaRed3 = alphaRed2 * alphaRed;
1109  if (eps_sur < 1.0)
1110  eps_sur = infinity;
1111  T0 = (std::isinf(eps_sur)) ? 1.0 : 2.0 * (eps_sur - 1.0) / (2.0 * eps_sur + 1.0);
1112  chi =
1113  (2.0 * alphaRed * exp(-alphaRed2) / pi_sqrt + (-2.0 * alphaRed2 + 1) * erfc(alphaRed) - 1.0) * pi / alpha2;
1114  kappa = 1.0 / debye_length;
1115  kappa2 = kappa * kappa;
1116  beta = kappa / (2.0 * alpha);
1117  beta2 = beta * beta;
1118  beta3 = beta2 * beta;
1119  self_energy_prefactor = {
1120  -alphaRed / pi_sqrt * (std::exp(-beta2) - pi_sqrt * beta * std::erfc(beta)),
1121  -alphaRed3 * 2.0 / 3.0 / pi_sqrt *
1122  (2.0 * pi_sqrt * beta3 * std::erfc(beta) + (1.0 - 2.0 * beta2) * std::exp(-beta2))};
1123  }
1124 
1125  inline double short_range_function(double q) const override {
1126  return 0.5 *
1127  (std::erfc(alphaRed * q + beta) * std::exp(4.0 * alphaRed * beta * q) + std::erfc(alphaRed * q - beta));
1128  }
1129  inline double short_range_function_derivative(double q) const override {
1130  double expC = std::exp(-powi(alphaRed * q - beta, 2));
1131  double erfcC = std::erfc(alphaRed * q + beta);
1132  return (-2.0 * alphaRed / pi_sqrt * expC + 2.0 * alphaRed * beta * erfcC * std::exp(4.0 * alphaRed * beta * q));
1133  }
1134  inline double short_range_function_second_derivative(double q) const override {
1135  double expC = std::exp(-powi(alphaRed * q - beta, 2));
1136  double erfcC = std::erfc(alphaRed * q + beta);
1137  return (4.0 * alphaRed2 / pi_sqrt * (alphaRed * q - 2.0 * beta) * expC +
1138  8.0 * alphaRed2 * beta2 * erfcC * std::exp(4.0 * alphaRed * beta * q));
1139  }
1140  inline double short_range_function_third_derivative(double q) const override {
1141  double expC = std::exp(-powi(alphaRed * q - beta, 2));
1142  double erfcC = std::erfc(alphaRed * q + beta);
1143  return (4.0 * alphaRed3 / pi_sqrt *
1144  (1.0 - 2.0 * (alphaRed * q - 2.0 * beta) * (alphaRed * q - beta) - 4.0 * beta2) * expC +
1145  32.0 * alphaRed3 * beta3 * erfcC * std::exp(4.0 * alphaRed * beta * q));
1146  }
1147 
1154  inline double neutralization_energy(const std::vector<double> &charges, double volume) const override {
1155  double squaredSumQ = 0.0;
1156  for (size_t i = 0; i < charges.size(); i++)
1157  squaredSumQ += charges[i] * charges[i]; // this should be squared, right?!
1158  return (-pi / (2.0 * alpha2 * volume) * squaredSumQ);
1159  }
1160 
1170  inline double reciprocal_energy(const std::vector<vec3> &positions, const std::vector<double> &charges,
1171  const std::vector<vec3> &dipoles, const vec3 &L, int nmax) const {
1172 
1173  assert(positions.size() == charges.size());
1174  assert(positions.size() == dipoles.size());
1175 
1176  double volume = L[0] * L[1] * L[2];
1177  std::vector<vec3> kvec;
1178  std::vector<double> Ak;
1179  // kvec.reserve( expected_size_of_kvec ); // speeds up push_back below
1180  // Ak.reserve( expected_size_of_Ak ); // speeds up push_back below
1181  for (int nx = -nmax; nx < nmax + 1; nx++) {
1182  for (int ny = -nmax; ny < nmax + 1; ny++) {
1183  for (int nz = -nmax; nz < nmax + 1; nz++) {
1184  vec3 kv = {2.0 * pi * nx / L[0], 2.0 * pi * ny / L[1], 2.0 * pi * nz / L[2]};
1185  double k2 = kv.squaredNorm() + kappa2;
1186  vec3 nv = {double(nx), double(ny), double(nz)};
1187  double nv1 = nv.squaredNorm();
1188  if (nv1 > 0) {
1189  if (nv1 <= nmax * nmax) {
1190  kvec.push_back(kv);
1191  Ak.push_back(std::exp(-k2 / (4.0 * alpha2) - beta2) / k2);
1192  }
1193  }
1194  }
1195  }
1196  }
1197 
1198  assert(kvec.size() == Ak.size());
1199 
1200  double E = 0.0;
1201  for (size_t k = 0; k < kvec.size(); k++) {
1202  std::complex<double> Qq(0.0, 0.0);
1203  std::complex<double> Qmu(0.0, 0.0);
1204  for (size_t i = 0; i < positions.size(); i++) {
1205  double kDotR = kvec[k].dot(positions[i]);
1206  double coskDotR = std::cos(kDotR);
1207  double sinkDotR = std::sin(kDotR);
1208  Qq += charges[i] * std::complex<double>(coskDotR, sinkDotR);
1209  Qmu += dipoles[i].dot(kvec[k]) * std::complex<double>(-sinkDotR, coskDotR);
1210  }
1211  std::complex<double> Q = Qq + Qmu;
1212  E += (powi(std::abs(Q), 2) * Ak[k]);
1213  }
1214  return (E * 2.0 * pi / volume);
1215  }
1216 
1224  inline double surface_energy(const std::vector<vec3> &positions, const std::vector<double> &charges,
1225  const std::vector<vec3> &dipoles, double volume) const {
1226  assert(positions.size() == charges.size());
1227  assert(positions.size() == dipoles.size());
1228  vec3 sum_r_charges = {0.0, 0.0, 0.0};
1229  vec3 sum_dipoles = {0.0, 0.0, 0.0};
1230  for (size_t i = 0; i < positions.size(); i++) {
1231  sum_r_charges += positions[i] * charges[i];
1232  sum_dipoles += dipoles[i];
1233  }
1234  double sqDipoles =
1235  sum_r_charges.dot(sum_r_charges) + 2.0 * sum_r_charges.dot(sum_dipoles) + sum_dipoles.dot(sum_dipoles);
1236 
1237  return (2.0 * pi / (2.0 * eps_sur + 1.0) / volume * sqDipoles);
1238  }
1239 
1240 #ifdef NLOHMANN_JSON_HPP
1241  inline Ewald(const nlohmann::json &j)
1242  : Ewald(j.at("cutoff").get<double>(), j.at("alpha").get<double>(), j.value("epss", infinity),
1243  j.value("debyelength", infinity)) {}
1244 
1245  private:
1246  inline void _to_json(nlohmann::json &j) const override {
1247  j = {{ "alpha", alpha }};
1248  if (std::isinf(eps_sur))
1249  j["epss"] = "inf";
1250  else
1251  j["epss"] = eps_sur;
1252  }
1253 #endif
1254 };
1255 
1256 // -------------- Reaction-field ---------------
1257 
1261 class ReactionField : public EnergyImplementation<ReactionField> {
1262  private:
1263  double epsRF;
1264  double epsr;
1265  bool shifted;
1266  const double pi = 4.0 * std::atan(1.0);
1267 
1268  public:
1275  inline ReactionField(double cutoff, double epsRF, double epsr, bool shifted)
1276  : EnergyImplementation(Scheme::reactionfield, cutoff), epsRF(epsRF), epsr(epsr), shifted(shifted) {
1277  name = "Reaction-field";
1278  dipolar_selfenergy = true;
1279  doi = "10.1080/00268977300102101";
1280  epsRF = epsRF;
1281  epsr = epsr;
1282  shifted = shifted;
1283  self_energy_prefactor = {-3.0 * epsRF * double(shifted) / (4.0 * epsRF + 2.0 * epsr),
1284  -(2.0 * epsRF - 2.0 * epsr) / (2.0 * (2.0 * epsRF + epsr))};
1285  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1286  chi = -6.0 * cutoff * cutoff * pi * ((-10.0 * double(shifted) / 3.0 + 4.0) * epsRF + epsr) /
1287  ((5.0 * (2.0 * epsRF + epsr)));
1288  }
1289 
1290  inline double short_range_function(double q) const override {
1291  return (1.0 + (epsRF - epsr) * q * q * q / (2.0 * epsRF + epsr) -
1292  3.0 * epsRF * q / (2.0 * epsRF + epsr) * double(shifted));
1293  }
1294  inline double short_range_function_derivative(double q) const override {
1295  return (3.0 * (epsRF - epsr) * q * q / (2.0 * epsRF + epsr) -
1296  3.0 * epsRF * double(shifted) / (2.0 * epsRF + epsr));
1297  }
1298  inline double short_range_function_second_derivative(double q) const override {
1299  return 6.0 * (epsRF - epsr) * q / (2.0 * epsRF + epsr);
1300  }
1301  inline double short_range_function_third_derivative(double) const override {
1302  return 6.0 * (epsRF - epsr) / (2.0 * epsRF + epsr);
1303  }
1304 
1305 #ifdef NLOHMANN_JSON_HPP
1306 
1307  inline ReactionField(const nlohmann::json &j)
1308  : ReactionField(j.at("cutoff").get<double>(), j.at("epsRF").get<double>(), j.at("epsr").get<double>(),
1309  j.at("shifted").get<bool>()) {}
1310 
1311  private:
1312  inline void _to_json(nlohmann::json &j) const override {
1313  j = {{"epsr", epsr}, {"epsRF", epsRF}, { "shifted", shifted }};
1314  }
1315 #endif
1316 };
1317 
1318 // -------------- Zahn ---------------
1319 
1323 class Zahn : public EnergyImplementation<Zahn> {
1324  private:
1325  double alpha;
1326  double alphaRed, alphaRed2;
1327  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1328  const double pi = 4.0 * std::atan(1.0);
1329 
1330  public:
1336  inline Zahn(double cutoff, double alpha) : EnergyImplementation(Scheme::zahn, cutoff), alpha(alpha) {
1337  name = "Zahn";
1338  dipolar_selfenergy = false;
1339  doi = "10.1021/jp025949h";
1340  alphaRed = alpha * cutoff;
1341  alphaRed2 = alphaRed * alphaRed;
1342  self_energy_prefactor = {-alphaRed * (1.0 - std::exp(-alphaRed2)) / pi_sqrt + 0.5 * std::erfc(alphaRed),
1343  0.0}; // Dipole self-energy undefined!
1344  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1345  chi = -(2.0 * (alphaRed * (alphaRed2 - 3.0) * std::exp(-alphaRed2) -
1346  0.5 * std::sqrt(pi) * ((7.0 - 3.0 / alphaRed2) * std::erf(alphaRed) - 7.0) * alphaRed2)) *
1347  cutoff * cutoff * std::sqrt(pi) / (3.0 * alphaRed2);
1348  }
1349 
1350  inline double short_range_function(double q) const override {
1351  return (std::erfc(alphaRed * q) -
1352  (q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1353  }
1354  inline double short_range_function_derivative(double q) const override {
1355  return (-(4.0 * (0.5 * std::exp(-alphaRed2 * q * q) * alphaRed +
1356  (alphaRed * std::exp(-alphaRed2) + 0.5 * pi_sqrt * std::erfc(alphaRed)) * (q - 0.5))) /
1357  pi_sqrt);
1358  }
1359  inline double short_range_function_second_derivative(double q) const override {
1360  return (4.0 * (alphaRed2 * alphaRed * q * std::exp(-alphaRed2 * q * q) - alphaRed * std::exp(-alphaRed2) -
1361  0.5 * pi_sqrt * std::erfc(alphaRed))) /
1362  pi_sqrt;
1363  }
1364  inline double short_range_function_third_derivative(double q) const override {
1365  return (-8.0 * std::exp(-alphaRed2 * q * q) * (alphaRed2 * q * q - 0.5) * alphaRed2 * alphaRed / pi_sqrt);
1366  }
1367 
1368 #ifdef NLOHMANN_JSON_HPP
1369 
1370  inline Zahn(const nlohmann::json &j) : Zahn(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1371 
1372  private:
1373  inline void _to_json(nlohmann::json &j) const override {
1374  j = {{ "alpha", alpha }};
1375  }
1376 #endif
1377 };
1378 
1379 // -------------- Fennell ---------------
1380 
1384 class Fennell : public EnergyImplementation<Fennell> {
1385  private:
1386  double alpha;
1387  double alphaRed, alphaRed2;
1388  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1389  const double pi = 4.0 * std::atan(1.0);
1390 
1391  public:
1396  inline Fennell(double cutoff, double alpha) : EnergyImplementation(Scheme::fennell, cutoff), alpha(alpha) {
1397  name = "Fennell";
1398  dipolar_selfenergy = false;
1399  doi = "10.1063/1.2206581";
1400  alphaRed = alpha * cutoff;
1401  alphaRed2 = alphaRed * alphaRed;
1402  self_energy_prefactor = {-alphaRed * (1.0 + std::exp(-alphaRed2)) / pi_sqrt - std::erfc(alphaRed),
1403  0.0}; // Dipole self-energy undefined!
1404  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1405  chi = (2.0 * ((alphaRed2 + 3.0) * alphaRed * std::exp(-alphaRed2) +
1406  0.5 * (std::erf(alphaRed) * alphaRed2 - alphaRed2 - 3.0 * std::erf(alphaRed)) * std::sqrt(pi))) *
1407  std::sqrt(pi) * cutoff * cutoff / (3.0 * alphaRed2);
1408  }
1409 
1410  inline double short_range_function(double q) const override {
1411  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed) +
1412  (q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1413  }
1414  inline double short_range_function_derivative(double q) const override {
1415  return (2.0 * alphaRed * (2.0 * (q - 0.5) * std::exp(-alphaRed2) - std::exp(-alphaRed2 * q * q)) / pi_sqrt +
1416  2.0 * std::erfc(alphaRed) * (q - 1.0));
1417  }
1418  inline double short_range_function_second_derivative(double q) const override {
1419  return (4.0 * alphaRed * (alphaRed2 * q * std::exp(-alphaRed2 * q * q) + std::exp(-alphaRed2)) / pi_sqrt +
1420  2.0 * std::erfc(alphaRed));
1421  }
1422  inline double short_range_function_third_derivative(double q) const override {
1423  return 4.0 * alphaRed2 * alphaRed * (1.0 - 2.0 * alphaRed2 * q * q) * std::exp(-alphaRed2 * q * q) / pi_sqrt;
1424  }
1425 
1426 #ifdef NLOHMANN_JSON_HPP
1427 
1428  inline Fennell(const nlohmann::json &j) : Fennell(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1429 
1430  private:
1431  inline void _to_json(nlohmann::json &j) const override {
1432  j = {{ "alpha", alpha }};
1433  }
1434 #endif
1435 };
1436 
1437 // -------------- Zero-dipole ---------------
1438 
1442 class ZeroDipole : public EnergyImplementation<ZeroDipole> {
1443  private:
1444  double alpha;
1445  double alphaRed, alphaRed2;
1446  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1447  const double pi = 4.0 * std::atan(1.0);
1448 
1449  public:
1454  inline ZeroDipole(double cutoff, double alpha) : EnergyImplementation(Scheme::zerodipole, cutoff), alpha(alpha) {
1455  name = "ZeroDipole";
1456  dipolar_selfenergy = true;
1457  doi = "10.1063/1.3582791";
1458  alphaRed = alpha * cutoff;
1459  alphaRed2 = alphaRed * alphaRed;
1460  self_energy_prefactor = {-alphaRed * (1.0 + 0.5 * std::exp(-alphaRed2)) / pi_sqrt - 0.75 * std::erfc(alphaRed),
1461  -alphaRed * (2.0 * alphaRed2 * (1.0 / 3.0) + std::exp(-alphaRed2)) / pi_sqrt -
1462  0.5 * std::erfc(alphaRed)};
1463  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1464  chi = cutoff * cutoff *
1465  ((6.0 * alphaRed2 - 15.0) * std::erf(alphaRed) * pi - 6.0 * pi * alphaRed2 +
1466  (8.0 * alphaRed2 + 30.0) * alphaRed * std::exp(-alphaRed2) * std::sqrt(pi)) /
1467  (15.0 * alphaRed2);
1468  }
1469 
1470  inline double short_range_function(double q) const override {
1471  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed) +
1472  0.5 * (q * q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1473  }
1474  inline double short_range_function_derivative(double q) const override {
1475  return (alphaRed * ((3.0 * q * q - 1.0) * std::exp(-alphaRed2) - 2.0 * std::exp(-alphaRed2 * q * q)) / pi_sqrt +
1476  1.5 * std::erfc(alphaRed) * (q * q - 1.0));
1477  }
1478  inline double short_range_function_second_derivative(double q) const override {
1479  return (2.0 * alphaRed * q * (2.0 * alphaRed2 * std::exp(-alphaRed2 * q * q) + 3.0 * std::exp(-alphaRed2)) /
1480  pi_sqrt +
1481  3.0 * q * std::erfc(alphaRed));
1482  }
1483  inline double short_range_function_third_derivative(double q) const override {
1484  return (2.0 * alphaRed *
1485  (2.0 * alphaRed2 * (1.0 - 2.0 * alphaRed2 * q * q) * std::exp(-alphaRed2 * q * q) +
1486  3.0 * std::exp(-alphaRed2)) /
1487  pi_sqrt +
1488  3.0 * std::erfc(alphaRed));
1489  }
1490 
1491 #ifdef NLOHMANN_JSON_HPP
1492 
1493  inline ZeroDipole(const nlohmann::json &j)
1494  : ZeroDipole(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1495 
1496  private:
1497  inline void _to_json(nlohmann::json &j) const override {
1498  j = {{ "alpha", alpha }};
1499  }
1500 #endif
1501 };
1502 
1503 // -------------- Wolf ---------------
1504 
1508 class Wolf : public EnergyImplementation<Wolf> {
1509  private:
1510  double alpha;
1511  double alphaRed, alphaRed2;
1512  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1513  const double pi = 4.0 * std::atan(1.0);
1514 
1515  public:
1520  inline Wolf(double cutoff, double alpha) : EnergyImplementation(Scheme::wolf, cutoff), alpha(alpha) {
1521  name = "Wolf";
1522  dipolar_selfenergy = true;
1523  doi = "10.1063/1.478738";
1524  alphaRed = alpha * cutoff;
1525  alphaRed2 = alphaRed * alphaRed;
1526  self_energy_prefactor = {-alphaRed / pi_sqrt - std::erfc(alphaRed) / 2.0,
1527  -powi(alphaRed, 3) * 2.0 / 3.0 / pi_sqrt};
1528  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1529  chi = 2.0 * std::sqrt(pi) * cutoff * cutoff *
1530  (3.0 * std::exp(-alphaRed2) * alphaRed -
1531  std::sqrt(pi) * (std::erfc(alphaRed) * alphaRed2 + 3.0 * std::erf(alphaRed) * 0.5)) /
1532  (3.0 * alphaRed2);
1533  }
1534 
1535  inline double short_range_function(double q) const override {
1536  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed));
1537  }
1538  inline double short_range_function_derivative(double q) const override {
1539  return (-2.0 * std::exp(-alphaRed2 * q * q) * alphaRed / pi_sqrt - std::erfc(alphaRed));
1540  }
1541  inline double short_range_function_second_derivative(double q) const override {
1542  return 4.0 * std::exp(-alphaRed2 * q * q) * alphaRed2 * alphaRed * q / pi_sqrt;
1543  }
1544  inline double short_range_function_third_derivative(double q) const override {
1545  return -8.0 * std::exp(-alphaRed2 * q * q) * alphaRed2 * alphaRed * (alphaRed2 * q * q - 0.5) / pi_sqrt;
1546  }
1547 
1548 #ifdef NLOHMANN_JSON_HPP
1549 
1550  inline Wolf(const nlohmann::json &j) : Wolf(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1551 
1552  private:
1553  inline void _to_json(nlohmann::json &j) const override {
1554  j = {{ "alpha", alpha }};
1555  }
1556 #endif
1557 };
1558 
1559 // -------------- qPotential ---------------
1560 template <int order> class qPotentialFixedOrder : public EnergyImplementation<qPotentialFixedOrder<order>> {
1561  public:
1563  using base::chi;
1564  using base::name;
1565  using base::self_energy_prefactor;
1566  using base::T0;
1571  inline qPotentialFixedOrder(double cutoff) : base(Scheme::qpotential, cutoff) {
1572  name = "qpotential";
1573  this->dipolar_selfenergy = true;
1574  this->doi = "10.1039/c9cp03875b";
1575  self_energy_prefactor = {-0.5, -0.5};
1576  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1577  chi = -86459.0 * std::acos(-1.0) * cutoff * cutoff / 235620.0;
1578  }
1579 
1580  inline double short_range_function(double q) const override { return qPochhammerSymbol(q, 0, order); }
1581  inline double short_range_function_derivative(double q) const override {
1582  return qPochhammerSymbolDerivative(q, 0, order);
1583  }
1584  inline double short_range_function_second_derivative(double q) const override {
1585  return qPochhammerSymbolSecondDerivative(q, 0, order);
1586  }
1587  inline double short_range_function_third_derivative(double q) const override {
1588  return qPochhammerSymbolThirdDerivative(q, 0, order);
1589  }
1590 
1591 #ifdef NLOHMANN_JSON_HPP
1592  inline qPotentialFixedOrder(const nlohmann::json &j) : qPotentialFixedOrder(j.at("cutoff").get<double>()) {}
1593 
1594  private:
1595  inline void _to_json(nlohmann::json &j) const override {
1596  j = {{ "order", order }};
1597  }
1598 #endif
1599 };
1600 
1611 class qPotential : public EnergyImplementation<qPotential> {
1612  private:
1613  int order;
1614 
1615  public:
1620  inline qPotential(double cutoff, int order) : EnergyImplementation(Scheme::qpotential, cutoff), order(order) {
1621  name = "qpotential";
1622  dipolar_selfenergy = true;
1623  doi = "10.1039/c9cp03875b";
1624  self_energy_prefactor = {-0.5, -0.5};
1625  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1626  chi = 0.0; // -Pi*Rc^2 * [ 2/3 7/15 17/42 146/385 86459/235620 ]
1627  }
1628 
1629  inline double short_range_function(double q) const override { return qPochhammerSymbol(q, 0, order); }
1630  inline double short_range_function_derivative(double q) const override {
1631  return qPochhammerSymbolDerivative(q, 0, order);
1632  }
1633  inline double short_range_function_second_derivative(double q) const override {
1634  return qPochhammerSymbolSecondDerivative(q, 0, order);
1635  }
1636  inline double short_range_function_third_derivative(double q) const override {
1637  return qPochhammerSymbolThirdDerivative(q, 0, order);
1638  }
1639 
1640 #ifdef NLOHMANN_JSON_HPP
1641 
1642  inline qPotential(const nlohmann::json &j)
1643  : qPotential(j.at("cutoff").get<double>(), j.at("order").get<double>()) {}
1644 
1645  private:
1646  inline void _to_json(nlohmann::json &j) const override {
1647  j = {{ "order", order }};
1648  }
1649 #endif
1650 };
1651 
1684 class Poisson : public EnergyImplementation<Poisson> {
1685  private:
1686  signed int C, D;
1687  double kappaRed, kappaRed2;
1688  double yukawa_denom, binomCDC;
1689  bool yukawa;
1690 
1691  public:
1698  inline Poisson(double cutoff, signed int C, signed int D, double debye_length = infinity)
1699  : EnergyImplementation(Scheme::poisson, cutoff, debye_length), C(C), D(D) {
1700  if ( C < 1 )
1701  throw std::runtime_error("`C` must be larger than zero");
1702  if ( ( D < -1 ) && ( D != -C ) )
1703  throw std::runtime_error("If `D` is less than negative one, then it has to equal negative `C`");
1704  if ( ( D == 0 ) && ( C != 1 ) )
1705  throw std::runtime_error("If `D` is zero, then `C` has to equal one ");
1706  name = "poisson";
1707  dipolar_selfenergy = true;
1708  if(C < 2)
1709  dipolar_selfenergy = false;
1710  doi = "10/c5fr";
1711  double a1 = -double(C + D) / double(C);
1712  kappaRed = 0.0;
1713  yukawa = false;
1714  if( !std::isinf(debye_length) ) {
1715  kappaRed = cutoff / debye_length;
1716  if (std::fabs(kappaRed) > 1e-6) {
1717  yukawa = true;
1718  kappaRed2 = kappaRed * kappaRed;
1719  yukawa_denom = 1.0 / (1.0 - std::exp(2.0 * kappaRed));
1720  a1 *= -2.0 * kappaRed * yukawa_denom;
1721  }
1722  }
1723  binomCDC = 0.0;
1724  if( D != -C )
1725  binomCDC = double(binomial(C + D, C) * D);
1726  self_energy_prefactor = {0.5 * a1, 0.0}; // Dipole self-energy seems to be 0 for C >= 2
1727  if(C == 2)
1728  self_energy_prefactor.at(2) = -double(D) * ( double(D*D) + 3.0 * double(D) + 2.0 ) / 12.0;
1729  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) +
1730  short_range_function(0.0); // Is this OK for Yukawa-interactions?
1731  chi = -2.0 * std::acos(-1.0) * cutoff * cutoff * (1.0 + double(C)) * (2.0 + double(C)) /
1732  (3.0 * double(D + 1 + C) *
1733  double(D + 2 + C)); // not confirmed, but have worked for all tested values of 'C' and 'D'
1734  }
1735 
1736  inline double short_range_function(double q) const override {
1737  if( D == -C )
1738  return 1.0;
1739  double tmp = 0;
1740  double qp = q;
1741  if (yukawa)
1742  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
1743  if( ( D == 0 ) && ( C == 1 ) )
1744  return ( 1.0 - qp );
1745  for (signed int c = 0; c < C; c++)
1746  tmp += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
1747  return powi(1.0 - qp, D + 1) * tmp;
1748  }
1749 
1750  inline double short_range_function_derivative(double q) const override {
1751  if( D == -C )
1752  return 0.0;
1753  if( ( D == 0 ) && ( C == 1 ) )
1754  return 0.0;
1755  double qp = q;
1756  double dqpdq = 1.0;
1757  if (yukawa) {
1758  double exp2kq = std::exp(2.0 * kappaRed * q);
1759  qp = (1.0 - exp2kq) * yukawa_denom;
1760  dqpdq = -2.0 * kappaRed * exp2kq * yukawa_denom;
1761  }
1762  double tmp1 = 1.0;
1763  double tmp2 = 0.0;
1764  for (int c = 1; c < C; c++) {
1765  double _fact = double(binomial(D - 1 + c, c)) * double(C - c) / double(C);
1766  tmp1 += _fact * powi(qp, c);
1767  tmp2 += _fact * double(c) * powi(qp, c - 1);
1768  }
1769  double dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
1770  return dSdqp * dqpdq;
1771  }
1772 
1773  inline double short_range_function_second_derivative(double q) const override {
1774  if( D == -C )
1775  return 0.0;
1776  if( ( D == 0 ) && ( C == 1 ) )
1777  return 0.0;
1778  double qp = q;
1779  double dqpdq = 1.0;
1780  double d2qpdq2 = 0.0;
1781  double dSdqp = 0.0;
1782  if (yukawa) {
1783  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
1784  dqpdq = -2.0 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1785  d2qpdq2 = -4.0 * kappaRed2 * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1786  double tmp1 = 1.0;
1787  double tmp2 = 0.0;
1788  for (signed int c = 1; c < C; c++) {
1789  tmp1 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
1790  tmp2 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * double(c) * powi(qp, c - 1);
1791  }
1792  dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
1793  }
1794  double d2Sdqp2 = binomCDC * powi(1.0 - qp, D - 1) * powi(qp, C - 1);
1795  return (d2Sdqp2 * dqpdq * dqpdq + dSdqp * d2qpdq2);
1796  };
1797 
1798  inline double short_range_function_third_derivative(double q) const override {
1799  if( D == -C )
1800  return 0.0;
1801  if( ( D == 0 ) && ( C == 1 ) )
1802  return 0.0;
1803  double qp = q;
1804  double dqpdq = 1.0;
1805  double d2qpdq2 = 0.0;
1806  double d3qpdq3 = 0.0;
1807  double d2Sdqp2 = 0.0;
1808  double dSdqp = 0.0;
1809  if (yukawa) {
1810  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
1811  dqpdq = -2.0 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1812  d2qpdq2 = -4.0 * kappaRed2 * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1813  d3qpdq3 = -8.0 * kappaRed2 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1814  d2Sdqp2 = binomCDC * powi(1.0 - qp, D - 1) * powi(qp, C - 1);
1815  double tmp1 = 1.0;
1816  double tmp2 = 0.0;
1817  for (signed int c = 1; c < C; c++) {
1818  tmp1 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
1819  tmp2 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * double(c) * powi(qp, c - 1);
1820  }
1821  dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
1822  }
1823  double d3Sdqp3 =
1824  binomCDC * powi(1.0 - qp, D - 2) * powi(qp, C - 2) * ((2.0 - double(C + D)) * qp + double(C) - 1.0);
1825  return (d3Sdqp3 * dqpdq * dqpdq * dqpdq + 3.0 * d2Sdqp2 * dqpdq * d2qpdq2 + dSdqp * d3qpdq3);
1826  };
1827 
1828 #ifdef NLOHMANN_JSON_HPP
1829 
1831  inline Poisson(const nlohmann::json &j)
1832  : Poisson(j.at("cutoff").get<double>(), j.at("C").get<int>(), j.at("D").get<int>(),
1833  j.value("debyelength", infinity)) {}
1834 
1835  private:
1836  inline void _to_json(nlohmann::json &j) const override {
1837  j = {{"C", C}, { "D", D }};
1838  }
1839 #endif
1840 };
1841 
1842 // -------------- Fanourgakis ---------------
1843 
1848 class Fanourgakis : public EnergyImplementation<Fanourgakis> {
1849  public:
1853  inline Fanourgakis(double cutoff) : EnergyImplementation(Scheme::fanourgakis, cutoff) {
1854  name = "fanourgakis";
1855  dipolar_selfenergy = true;
1856  doi = "10.1063/1.3216520";
1857  self_energy_prefactor = {-0.875, 0.0};
1858  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1859  chi = -5.0 * std::acos(-1.0) * cutoff * cutoff / 18.0;
1860  }
1861 
1862  inline double short_range_function(double q) const override {
1863  double q2 = q * q;
1864  return powi(1.0 - q, 4) * (1.0 + 2.25 * q + 3.0 * q2 + 2.5 * q2 * q);
1865  }
1866  inline double short_range_function_derivative(double q) const override {
1867  return (-1.75 + 26.25 * powi(q, 4) - 42.0 * powi(q, 5) + 17.5 * powi(q, 6));
1868  }
1869  inline double short_range_function_second_derivative(double q) const override {
1870  return 105.0 * powi(q, 3) * powi(q - 1.0, 2);
1871  };
1872  inline double short_range_function_third_derivative(double q) const override {
1873  return 525.0 * powi(q, 2) * (q - 0.6) * (q - 1.0);
1874  };
1875 #ifdef NLOHMANN_JSON_HPP
1876 
1877  inline Fanourgakis(const nlohmann::json &j) : Fanourgakis(j.at("cutoff").get<double>()) {}
1878 
1879  private:
1880  inline void _to_json(nlohmann::json &) const override {}
1881 #endif
1882 };
1883 
1884 #ifdef NLOHMANN_JSON_HPP
1885 inline std::shared_ptr<SchemeBase> createScheme(const nlohmann::json &j) {
1886  const std::map<std::string, Scheme> m = {{"plain", Scheme::plain},
1887  {"qpotential", Scheme::qpotential},
1888  {"wolf", Scheme::wolf},
1889  {"poisson", Scheme::poisson},
1890  {"reactionfield", Scheme::reactionfield},
1891  {"spline", Scheme::spline},
1892  {"fanourgakis", Scheme::fanourgakis},
1893  {"fennell", Scheme::fennell},
1894  {"zahn", Scheme::zahn},
1895  {"zerodipole", Scheme::zerodipole},
1896  {"ewald", Scheme::ewald}}; // map string keyword to scheme type
1897 
1898  std::string name = j.at("type").get<std::string>();
1899  auto it = m.find(name);
1900  if (it == m.end())
1901  throw std::runtime_error("unknown coulomb scheme " + name);
1902 
1903  std::shared_ptr<SchemeBase> scheme;
1904  switch (it->second) {
1905  case Scheme::plain:
1906  scheme = std::make_shared<Plain>(j);
1907  break;
1908  case Scheme::wolf:
1909  scheme = std::make_shared<Wolf>(j);
1910  break;
1911  case Scheme::zahn:
1912  scheme = std::make_shared<Zahn>(j);
1913  break;
1914  case Scheme::fennell:
1915  scheme = std::make_shared<Fennell>(j);
1916  break;
1917  case Scheme::zerodipole:
1918  scheme = std::make_shared<ZeroDipole>(j);
1919  break;
1920  case Scheme::fanourgakis:
1921  scheme = std::make_shared<Fanourgakis>(j);
1922  break;
1923  case Scheme::qpotential5:
1924  scheme = std::make_shared<qPotentialFixedOrder<5>>(j);
1925  break;
1926  case Scheme::qpotential:
1927  scheme = std::make_shared<qPotential>(j);
1928  break;
1929  case Scheme::ewald:
1930  scheme = std::make_shared<Ewald>(j);
1931  break;
1932  case Scheme::poisson:
1933  scheme = std::make_shared<Poisson>(j);
1934  break;
1935  case Scheme::reactionfield:
1936  scheme = std::make_shared<ReactionField>(j);
1937  break;
1938  default:
1939  break;
1940  }
1941  return scheme;
1942 }
1943 #endif
1944 
1945 // -------------- Splined ---------------
1946 
1963 class Splined : public EnergyImplementation<Splined> {
1964  private:
1965  std::shared_ptr<SchemeBase> pot;
1966  Tabulate::Andrea<double> splined_srf; // spline class
1967  std::array<Tabulate::TabulatorBase<double>::data, 4> splinedata; // 0=original, 1=first derivative, ...
1968 
1969  inline void generate_spline_data() {
1970  assert(pot);
1971  SchemeBase::operator=(*pot); // copy base data from pot -> Splined
1972  splined_srf.setTolerance(1e-3, 1e-1);
1973  splinedata[0] = splined_srf.generate([pot = pot](double q) { return pot->short_range_function(q); }, 0, 1);
1974  splinedata[1] =
1975  splined_srf.generate([pot = pot](double q) { return pot->short_range_function_derivative(q); }, 0, 1);
1976  splinedata[2] = splined_srf.generate(
1977  [pot = pot](double q) { return pot->short_range_function_second_derivative(q); }, 0, 1);
1978  splinedata[3] =
1979  splined_srf.generate([pot = pot](double q) { return pot->short_range_function_third_derivative(q); }, 0, 1);
1980  }
1981 
1982  public:
1983  inline Splined() : EnergyImplementation<Splined>(Scheme::spline, infinity) {}
1984 
1991  template <class T, class... Args> void spline(Args &&... args) {
1992  pot = std::make_shared<T>(args...);
1993  generate_spline_data();
1994  }
1995  inline double short_range_function(double q) const override { return splined_srf.eval(splinedata[0], q); };
1996 
1997  inline double short_range_function_derivative(double q) const override {
1998  return splined_srf.eval(splinedata[1], q);
1999  }
2000  inline double short_range_function_second_derivative(double q) const override {
2001  return splined_srf.eval(splinedata[2], q);
2002  }
2003  inline double short_range_function_third_derivative(double q) const override {
2004  return splined_srf.eval(splinedata[3], q);
2005  }
2006 #ifdef NLOHMANN_JSON_HPP
2007  public:
2008  inline void to_json(nlohmann::json &j) const { pot->to_json(j); }
2009 
2010  private:
2011  inline void _to_json(nlohmann::json &) const override {}
2012 #endif
2013 };
2014 
2015 } // namespace CoulombGalore
constexpr double infinity
Numerical infinity.
Definition: coulombgalore.h:41
qPotential scheme
double ion_dipole_energy(double z, const vec3 &mu, const vec3 &r) const override
interaction energy between a point charges and a point dipole
Wolf(const nlohmann::json &j)
double self_energy(const std::array< double, 2 > &m2) const
Self-energy for all type of interactions.
vec3 ion_ion_force(double zA, double zB, const vec3 &r) const override
interaction force between two point charges
constexpr unsigned int factorial(unsigned int n)
Returns the factorial of &#39;n&#39;. Note that &#39;n&#39; must be positive semidefinite.
Definition: coulombgalore.h:76
Fanourgakis(const nlohmann::json &j)
ReactionField(const nlohmann::json &j)
vec3 multipole_multipole_force(double zA, double zB, const vec3 &muA, const vec3 &muB, const vec3 &r) const override
interaction force between two point multipoles
ReactionField(double cutoff, double epsRF, double epsr, bool shifted)
vec3 ion_dipole_force(double z, const vec3 &mu, const vec3 &r) const override
interaction force between a point charges and a point dipole
Base class for truncation schemes.
ZeroDipole(const nlohmann::json &j)
vec3 multipole_field(double z, const vec3 &mu, const vec3 &r) const
electrostatic field from point multipole
std::string name
Descriptive name.
Ewald real-space scheme.
Wolf(double cutoff, double alpha)
double cutoff
Cut-off distance, UNIT: [ input length ].
double debye_length
Debye-length, UNIT: [ input length ].
double qPochhammerSymbolThirdDerivative(double q, int l=0, int P=300)
Gives the third derivative of the q-Pochhammer Symbol.
double dipole_dipole_energy(const vec3 &muA, const vec3 &muB, const vec3 &r) const override
interaction energy between two point dipoles
Poisson(double cutoff, signed int C, signed int D, double debye_length=infinity)
vec3 dipole_torque(const vec3 &mu, const vec3 &E) const
torque exerted on point dipole due to field
void spline(Args &&... args)
Spline given potential type.
Fennell(const nlohmann::json &j)
double qPochhammerSymbol(double q, int l=0, int P=300)
Help-function for the q-potential scheme.
Fanourgakis scheme.
Scheme scheme
Truncation scheme.
Zahn(double cutoff, double alpha)
Contructor.
Intermediate base class that implements the interaction energies.
Poisson scheme with and without specified Debye-length.
double neutralization_energy(const std::vector< double > &charges, double volume) const override
Compensating term for non-neutral systems.
qPotential(double cutoff, int order)
double calc_dielectric(double M2V)
Calculate dielectric constant.
double ion_potential(double z, double r) const override
electrostatic potential from point charge
ZeroDipole(double cutoff, double alpha)
double surface_energy(const std::vector< vec3 > &positions, const std::vector< double > &charges, const std::vector< vec3 > &dipoles, double volume) const
Surface-term.
vec3 dipole_field(const vec3 &mu, const vec3 &r) const
electrostatic field from point dipole
Fennell scheme.
Eigen::Vector3d vec3
Definition: coulombgalore.h:39
double qPochhammerSymbolDerivative(double q, int l=0, int P=300)
Gives the derivative of the q-Pochhammer Symbol.
double qPochhammerSymbolSecondDerivative(double q, int l=0, int P=300)
Gives the second derivative of the q-Pochhammer Symbol.
double dipole_potential(const vec3 &mu, const vec3 &r) const override
electrostatic potential from point dipole
Zahn(const nlohmann::json &j)
double powi(double x, int n)
n&#39;th integer power of float
Definition: coulombgalore.h:64
vec3 dipole_dipole_force(const vec3 &muA, const vec3 &muB, const vec3 &r) const override
interaction force between two point dipoles
qPotential(const nlohmann::json &j)
Ewald(double cutoff, double alpha, double eps_sur=infinity, double debye_length=infinity)
double multipole_multipole_energy(double zA, double zB, const vec3 &muA, const vec3 &muB, const vec3 &r) const override
interaction energy between two multipoles with charges and dipole moments
std::string doi
DOI for original citation.
Dynamic scheme where all short ranged functions are splined.
double ion_ion_energy(double zA, double zB, double r) const override
interaction energy between two point charges
Fennell(double cutoff, double alpha)
Poisson(const nlohmann::json &j)
No truncation scheme, cutoff = infinity.
vec3 ion_field(double z, const vec3 &r) const
electrostatic field from point charge
Reaction-field scheme.
Zero-dipole scheme.
double reciprocal_energy(const std::vector< vec3 > &positions, const std::vector< double > &charges, const std::vector< vec3 > &dipoles, const vec3 &L, int nmax) const
Reciprocal-space energy.
double neutralization_energy(const std::vector< double > &charges, double volume) const override
Compensating term for non-neutral systems.