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 <functional>
34 #include <memory>
35 #include <Eigen/Core>
36 
38 #ifdef NLOHMANN_JSON_HPP_
39 #define NLOHMANN_JSON_HPP
40 #endif
41 
43 namespace CoulombGalore {
44 
46 typedef Eigen::Vector3d vec3;
47 typedef Eigen::Matrix3d mat33;
48 
49 constexpr double infinity = std::numeric_limits<double>::infinity();
50 
52 enum class Scheme {
53  plain,
54  ewald,
55  ewaldt,
56  reactionfield,
57  wolf,
58  poisson,
59  qpotential,
60  fanourgakis,
61  zerodipole,
62  zahn,
63  fennell,
64  qpotential5,
65  spline
66 };
67 
73 inline double powi(double x, int n) {
74 #if defined(__GNUG__)
75  return __builtin_powi(x, n);
76 #else
77  return std::pow(x, n);
78 #endif
79 }
80 
85 inline constexpr unsigned int factorial(unsigned int n) { return n <= 1 ? 1 : n * factorial(n - 1); }
86 
87 constexpr unsigned int binomial(signed int n, signed int k) {
88  return factorial(n) / (factorial(k) * factorial(n - k));
89 }
90 
110 inline double qPochhammerSymbol(double q, int l = 0, int P = 300) {
111  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
112  for (int n = 1; n < P + 1; n++) {
113  double val = 0.0;
114  for (int k = 1; k < n + l + 1; k++)
115  val += powi(q, k - 1);
116  Ct *= val;
117  }
118  double Dt = powi(1.0 - q, P); // (1-q)^P
119  return (Ct * Dt);
120 }
121 
125 inline double qPochhammerSymbolDerivative(double q, int l = 0, int P = 300) {
126  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
127  for (int n = 1; n < P + 1; n++) {
128  double val = 0.0;
129  for (int k = 1; k < n + l + 1; k++)
130  val += powi(q, k - 1);
131  Ct *= val;
132  }
133  double dCt = 0.0; // evaluates to derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
134  for (int n = 1; n < P + 1; n++) {
135  double nom = 0.0;
136  double denom = 1.0;
137  for (int k = 2; k < n + l + 1; k++) {
138  nom += (k - 1) * powi(q, k - 2);
139  denom += powi(q, k - 1);
140  }
141  dCt += nom / denom;
142  }
143  dCt *= Ct;
144  double Dt = powi(1.0 - q, P); // (1-q)^P
145  double dDt = 0.0;
146  if (P > 0)
147  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
148  return (Ct * dDt + dCt * Dt);
149 }
150 
154 inline double qPochhammerSymbolSecondDerivative(double q, int l = 0, int P = 300) {
155  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
156  double DS = 0.0;
157  double dDS = 0.0;
158  for (int n = 1; n < P + 1; n++) {
159  double tmp = 0.0;
160  for (int k = 1; k < n + l + 1; k++)
161  tmp += powi(q, k - 1);
162  Ct *= tmp;
163  double nom = 0.0;
164  double denom = 1.0;
165  for (int k = 2; k < n + l + 1; k++) {
166  nom += (k - 1) * powi(q, k - 2);
167  denom += powi(q, k - 1);
168  }
169  DS += nom / denom;
170  double diffNom = 0.0;
171  double diffDenom = 1.0;
172  for (int k = 3; k < n + l + 1; k++) {
173  diffNom += (k - 1) * (k - 2) * powi(q, k - 3);
174  diffDenom += (k - 1) * powi(q, k - 2);
175  }
176  dDS += (diffNom * denom - nom * diffDenom) / denom / denom;
177  }
178  double dCt = Ct * DS; // derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
179  double ddCt = dCt * DS + Ct * dDS; // second derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
180  double Dt = powi(1.0 - q, P); // (1-q)^P
181  double dDt = 0.0;
182  if (P > 0)
183  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
184  double ddDt = 0.0;
185  if (P > 1)
186  ddDt = P * (P - 1) * powi(1 - q, P - 2); // second derivative of (1-q)^P
187  return (Ct * ddDt + 2 * dCt * dDt + ddCt * Dt);
188 }
189 
193 inline double qPochhammerSymbolThirdDerivative(double q, int l = 0, int P = 300) {
194  double Ct = 1.0; // evaluates to \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
195  double DS = 0.0;
196  double dDS = 0.0;
197  double ddDS = 0.0;
198  for (int n = 1; n < P + 1; n++) {
199  double tmp = 0.0;
200  for (int k = 1; k < n + l + 1; k++)
201  tmp += powi(q, k - 1);
202  Ct *= tmp;
203  double f = 0.0;
204  double g = 1.0;
205  for (int k = 2; k < n + l + 1; k++) {
206  f += (k - 1) * powi(q, k - 2);
207  g += powi(q, k - 1);
208  }
209  DS += f / g;
210  double df = 0.0;
211  double dg = 0.0;
212  if (n + l > 1)
213  dg = 1.0;
214  for (int k = 3; k < n + l + 1; k++) {
215  df += (k - 1) * (k - 2) * powi(q, k - 3);
216  dg += (k - 1) * powi(q, k - 2);
217  }
218  dDS += (df * g - f * dg) / g / g;
219  double ddf = 0.0;
220  double ddg = 0.0;
221  if (n + l > 2)
222  ddg = 2.0;
223  for (int k = 4; k < n + l + 1; k++) {
224  ddf += (k - 1) * (k - 2) * (k - 3) * powi(q, k - 4);
225  ddg += (k - 1) * (k - 2) * powi(q, k - 3);
226  }
227  ddDS += (ddf * g * g - 2.0 * df * dg * g + 2.0 * f * dg * dg - f * ddg * g) / g / g / g;
228  }
229  double dCt = Ct * DS; // derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
230  double ddCt = dCt * DS + Ct * dDS; // second derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
231  double dddCt = ddCt * DS + 2.0 * dCt * dDS + Ct * ddDS; // third derivative of \prod_{n=1}^P\sum_{k=0}^{n+l}q^k
232  double Dt = powi(1.0 - q, P); // (1-q)^P
233  double dDt = 0.0;
234  if (P > 0)
235  dDt = -P * powi(1 - q, P - 1); // derivative of (1-q)^P
236  double ddDt = 0.0;
237  if (P > 1)
238  ddDt = P * (P - 1) * powi(1 - q, P - 2); // second derivative of (1-q)^P
239  double dddDt = 0.0;
240  if (P > 2)
241  dddDt = -P * (P - 1) * (P - 2) * powi(1 - q, P - 3); // third derivative of (1-q)^P
242  return (dddCt * Dt + 3.0 * ddCt * dDt + 3 * dCt * ddDt + Ct * dddDt);
243 }
244 
245 namespace Tabulate {
246 
247 /* base class for all tabulators - no dependencies */
248 template <typename T = double> class TabulatorBase {
249  protected:
250  T utol = 1e-5, ftol = -1, umaxtol = -1, fmaxtol = -1;
251  T numdr = 0.0001; // dr for derivative evaluation
252 
253  // First derivative with respect to x
254  T f1(std::function<T(T)> f, T x) const { return (f(x + numdr * 0.5) - f(x - numdr * 0.5)) / (numdr); }
255 
256  // Second derivative with respect to x
257  T f2(std::function<T(T)> f, T x) const { return (f1(f, x + numdr * 0.5) - f1(f, x - numdr * 0.5)) / (numdr); }
258 
259  void check() const {
260  if (ftol != -1 && ftol <= 0.0) {
261  std::cerr << "ftol=" << ftol << " too small\n" << std::endl;
262  abort();
263  }
264  if (umaxtol != -1 && umaxtol <= 0.0) {
265  std::cerr << "umaxtol=" << umaxtol << " too small\n" << std::endl;
266  abort();
267  }
268  if (fmaxtol != -1 && fmaxtol <= 0.0) {
269  std::cerr << "fmaxtol=" << fmaxtol << " too small\n" << std::endl;
270  abort();
271  }
272  }
273 
274  public:
275  struct data {
276  std::vector<T> r2; // r2 for intervals
277  std::vector<T> c; // c for coefficents
278  T rmin2 = 0, rmax2 = 0; // useful to save these with table
279  bool empty() const { return r2.empty() && c.empty(); }
280  inline size_t numKnots() const { return r2.size(); }
281  };
282 
283  void setTolerance(T _utol, T _ftol = -1, T _umaxtol = -1, T _fmaxtol = -1) {
284  utol = _utol;
285  ftol = _ftol;
286  umaxtol = _umaxtol;
287  fmaxtol = _fmaxtol;
288  }
289 
290  void setNumdr(T _numdr) { numdr = _numdr; }
291 };
292 
293 /*
294  * @brief Andrea table with logarithmic search
295  *
296  * Tabulator with logarithmic search.
297  * Code mainly from MolSim (Per Linse) with some upgrades
298  * Reference: doi:10/frzp4d
299  *
300  * @note Slow on Intel compiler
301  */
302 template <typename T = double> class Andrea : public TabulatorBase<T> {
303  private:
304  typedef TabulatorBase<T> base; // for convenience
305  int mngrid = 1200; // Max number of controlpoints
306  int ndr = 100; // Max number of trials to decr dr
307  T drfrac = 0.9; // Multiplicative factor to decr dr
308 
309  std::vector<T> SetUBuffer(T, T zlow, T, T zupp, T u0low, T u1low, T u2low, T u0upp, T u1upp, T u2upp) {
310 
311  // Zero potential and force return no coefficients
312  if (std::fabs(u0low) < 1e-9)
313  if (std::fabs(u1low) < 1e-9)
314  return {0, 0, 0, 0, 0, 0, 0};
315 
316  T dz1 = zupp - zlow;
317  T dz2 = dz1 * dz1;
318  T dz3 = dz2 * dz1;
319  T w0low = u0low;
320  T w1low = u1low;
321  T w2low = u2low;
322  T w0upp = u0upp;
323  T w1upp = u1upp;
324  T w2upp = u2upp;
325  T c0 = w0low;
326  T c1 = w1low;
327  T c2 = w2low * 0.5;
328  T a = 6 * (w0upp - c0 - c1 * dz1 - c2 * dz2) / dz3;
329  T b = 2 * (w1upp - c1 - 2 * c2 * dz1) / dz2;
330  T c = (w2upp - 2 * c2) / dz1;
331  T c3 = (10 * a - 12 * b + 3 * c) / 6;
332  T c4 = (-15 * a + 21 * b - 6 * c) / (6 * dz1);
333  T c5 = (2 * a - 3 * b + c) / (2 * dz2);
334 
335  return {zlow, c0, c1, c2, c3, c4, c5};
336  }
337 
338  /*
339  * @returns boolean vector.
340  * - `[0]==true`: tolerance is approved,
341  * - `[1]==true` Repulsive part is found.
342  */
343  std::vector<bool> CheckUBuffer(std::vector<T> &ubuft, T rlow, T rupp, std::function<T(T)> f) const {
344 
345  // Number of points to control
346  int ncheck = 11;
347  T dr = (rupp - rlow) / (ncheck - 1);
348  std::vector<bool> vb(2, false);
349 
350  for (int i = 0; i < ncheck; i++) {
351  T r1 = rlow + dr * ((T)i);
352  T r2 = r1 * r1;
353  T u0 = f(r2);
354  T u1 = base::f1(f, r2);
355  T dz = r2 - rlow * rlow;
356  T usum =
357  ubuft.at(1) +
358  dz * (ubuft.at(2) + dz * (ubuft.at(3) + dz * (ubuft.at(4) + dz * (ubuft.at(5) + dz * ubuft.at(6)))));
359 
360  T fsum = ubuft.at(2) +
361  dz * (2 * ubuft.at(3) + dz * (3 * ubuft.at(4) + dz * (4 * ubuft.at(5) + dz * (5 * ubuft.at(6)))));
362 
363  if (std::fabs(usum - u0) > base::utol)
364  return vb;
365  if (base::ftol != -1 && std::fabs(fsum - u1) > base::ftol)
366  return vb;
367  if (base::umaxtol != -1 && std::fabs(usum) > base::umaxtol)
368  vb[1] = true;
369  if (base::fmaxtol != -1 && std::fabs(usum) > base::fmaxtol)
370  vb[1] = true;
371  }
372  vb[0] = true;
373  return vb;
374  }
375 
376  public:
377  /*
378  * @brief Get tabulated value at f(x)
379  * @param d Table data
380  * @param r2 value
381  */
382  inline T eval(const typename base::data &d, T r2) const {
383  assert(r2!=0); // r2 cannot be *exactly* zero
384  size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1;
385  size_t pos6 = 6 * pos;
386  assert((pos6 + 5) < d.c.size() && "out of bounds");
387  T dz = r2 - d.r2[pos];
388  return d.c[pos6] +
389  dz * (d.c[pos6 + 1] +
390  dz * (d.c[pos6 + 2] + dz * (d.c[pos6 + 3] + dz * (d.c[pos6 + 4] + dz * (d.c[pos6 + 5])))));
391  }
392 
393  /*
394  * @brief Get tabulated value at df(x)/dx
395  * @param d Table data
396  * @param r2 value
397  */
398  T evalDer(const typename base::data &d, T r2) const {
399  size_t pos = std::lower_bound(d.r2.begin(), d.r2.end(), r2) - d.r2.begin() - 1;
400  size_t pos6 = 6 * pos;
401  T dz = r2 - d.r2[pos];
402  return (d.c[pos6 + 1] +
403  dz * (2.0 * d.c[pos6 + 2] +
404  dz * (3.0 * d.c[pos6 + 3] + dz * (4.0 * d.c[pos6 + 4] + dz * (5.0 * d.c[pos6 + 5])))));
405  }
406 
410  typename base::data generate(std::function<T(T)> f, double rmin, double rmax) {
411  rmin = std::sqrt(rmin);
412  rmax = std::sqrt(rmax);
413  base::check();
414  typename base::data td;
415  td.rmin2 = rmin * rmin;
416  td.rmax2 = rmax * rmax;
417 
418  T rumin = rmin;
419  T rmax2 = rmax * rmax;
420  T dr = rmax - rmin;
421  T rupp = rmax;
422  T zupp = rmax2;
423  bool repul = false; // Stop tabulation if repul is true
424 
425  td.r2.push_back(zupp);
426 
427  int i;
428  for (i = 0; i < mngrid; i++) {
429  T rlow = rupp;
430  T zlow;
431  std::vector<T> ubuft;
432  int j;
433 
434  dr = (rupp - rmin);
435 
436  for (j = 0; j < ndr; j++) {
437  zupp = rupp * rupp;
438  rlow = rupp - dr;
439  if (rumin > rlow)
440  rlow = rumin;
441 
442  zlow = rlow * rlow;
443 
444  T u0low = f(zlow);
445  T u1low = base::f1(f, zlow);
446  T u2low = base::f2(f, zlow);
447  T u0upp = f(zupp);
448  T u1upp = base::f1(f, zupp);
449  T u2upp = base::f2(f, zupp);
450 
451  ubuft = SetUBuffer(rlow, zlow, rupp, zupp, u0low, u1low, u2low, u0upp, u1upp, u2upp);
452  std::vector<bool> vb = CheckUBuffer(ubuft, rlow, rupp, f);
453  repul = vb[1];
454  if (vb[0]) {
455  rupp = rlow;
456  break;
457  }
458  dr *= drfrac;
459  }
460 
461  if (j >= ndr)
462  throw std::runtime_error("Andrea spline: try to increase utol/ftol");
463  if (ubuft.size() != 7)
464  throw std::runtime_error("Andrea spline: wrong size of ubuft, min value + 6 coefficients");
465 
466  td.r2.push_back(zlow);
467  for (size_t k = 1; k < ubuft.size(); k++)
468  td.c.push_back(ubuft.at(k));
469 
470  // Entered a highly repulsive part, stop tabulation
471  if (repul) {
472  rumin = rlow;
473  td.rmin2 = rlow * rlow;
474  }
475  if (rlow <= rumin || repul)
476  break;
477  }
478 
479  if (i >= mngrid)
480  throw std::runtime_error("Andrea spline: try to increase utol/ftol");
481 
482  // create final reversed c and r2
483 #if __cplusplus >= 201703L
484  // C++17 only code
485  assert(td.c.size() % 6 == 0);
486  assert(td.c.size() / (td.r2.size() - 1) == 6);
487  assert(std::is_sorted(td.r2.rbegin(), td.r2.rend()));
488  std::reverse(td.r2.begin(), td.r2.end()); // reverse all elements
489  for (size_t i = 0; i < td.c.size() / 2; i += 6) // reverse knot order in packets of six
490  std::swap_ranges(td.c.begin() + i, td.c.begin() + i + 6, td.c.end() - i - 6); // c++17 only
491  return td;
492 #else
493  typename base::data tdsort;
494  tdsort.rmax2 = td.rmax2;
495  tdsort.rmin2 = td.rmin2;
496 
497  // reverse copy all elements in r2
498  tdsort.r2.resize(td.r2.size());
499  std::reverse_copy(td.r2.begin(), td.r2.end(), tdsort.r2.begin());
500 
501  // sanity check before reverse knot copy
502  assert(std::is_sorted(td.r2.rbegin(), td.r2.rend()));
503  assert(td.c.size() % 6 == 0);
504  assert(td.c.size() / (td.r2.size() - 1) == 6);
505 
506  // reverse copy knots
507  tdsort.c.resize(td.c.size());
508  auto dst = tdsort.c.end();
509  for (auto src = td.c.begin(); src != td.c.end(); src += 6)
510  std::copy(src, src + 6, dst -= 6);
511  return tdsort;
512 #endif
513  }
514 };
515 
516 } // namespace Tabulate
517 
529 class SchemeBase {
530  private:
531  std::array<double, 2> self_energy_prefactor; // Prefactor for self-energies, UNIT: [ 1 ]
532 
533  protected:
534  double invcutoff = 0; // inverse cutoff distance, UNIT: [ ( input length )^-1 ]
535  double cutoff2 = 0; // square cutoff distance, UNIT: [ ( input length )^2 ]
536  double kappa = 0; // inverse Debye-length, UNIT: [ ( input length )^-1 ]
537  double T0 = 0; // Spatial Fourier transformed modified interaction tensor, used to calculate the dielectric
538  // constant, UNIT: [ 1 ]
539  double chi = 0; // Negative integrated volume potential to neutralize charged system, UNIT: [ ( input length )^2 ]
540  bool dipolar_selfenergy = false; // is there a valid dipolar self-energy?
541 
542  void setSelfEnergyPrefactor(const std::array<double, 2> &factor) {
543  self_energy_prefactor = factor;
544  selfEnergyFunctor = [invcutoff=invcutoff, factor = factor](const std::array<double, 2> &squared_moments) {
545  double e_self = 0.0;
546  for (int i = 0; i < (int)squared_moments.size(); i++)
547  e_self += factor[i] * squared_moments[i] * powi(invcutoff, 2 * i + 1);
548  return e_self;
549  };
550  }
551 
552  public:
553  std::string doi;
554  std::string name;
556  double cutoff;
557  double debye_length;
558 
559  std::function<double(const std::array<double, 2> &)> selfEnergyFunctor = nullptr;
560 
561  inline SchemeBase(Scheme scheme, double cutoff, double debye_length = infinity)
562  : invcutoff(1.0/cutoff), cutoff2(cutoff*cutoff), kappa(1.0/debye_length), scheme(scheme), cutoff(cutoff), debye_length(debye_length) {}
563 
564  virtual ~SchemeBase() = default;
565 
566  virtual double neutralization_energy(const std::vector<double> &, double) const { return 0.0; }
567 
582  double calc_dielectric(double M2V) { return (M2V * T0 + 2.0 * M2V + 1.0) / (M2V * T0 - M2V + 1.0); }
583 
584  virtual double short_range_function(double q) const = 0;
585  virtual double short_range_function_derivative(double q) const = 0;
586  virtual double short_range_function_second_derivative(double q) const = 0;
587  virtual double short_range_function_third_derivative(double q) const = 0;
588 
589  virtual double ion_potential(double, double) const = 0;
590  virtual double dipole_potential(const vec3 &, const vec3 &) const = 0;
591  virtual double quadrupole_potential(const mat33 &, const vec3 &) const = 0;
592 
593  virtual double ion_ion_energy(double, double, double) const = 0;
594  virtual double ion_dipole_energy(double, const vec3 &, const vec3 &) const = 0;
595  virtual double dipole_dipole_energy(const vec3 &, const vec3 &, const vec3 &) const = 0;
596  virtual double ion_quadrupole_energy(double, const mat33 &, const vec3 &) const = 0;
597  virtual double multipole_multipole_energy(double, double, const vec3 &, const vec3 &, const mat33 &, const mat33 &, const vec3 &) const = 0;
598 
599  virtual vec3 ion_field(double, const vec3 &) const = 0;
600  virtual vec3 dipole_field(const vec3 &, const vec3 &) const = 0;
601  virtual vec3 quadrupole_field(const mat33 &, const vec3 &) const = 0;
602  virtual vec3 multipole_field(double, const vec3 &, const mat33 &, const vec3 &) const = 0;
603 
604  virtual vec3 ion_ion_force(double, double, const vec3 &) const = 0;
605  virtual vec3 ion_dipole_force(double, const vec3 &, const vec3 &) const = 0;
606  virtual vec3 dipole_dipole_force(const vec3 &, const vec3 &, const vec3 &) const = 0;
607  virtual vec3 ion_quadrupole_force(double, const mat33 &, const vec3 &) const = 0;
608  virtual vec3 multipole_multipole_force(double, double, const vec3 &, const vec3 &, const mat33 &, const mat33 &, const vec3 &) const = 0;
609 
610  // add remaining funtions here...
611 
612  // virtual double surface_energy(const std::vector<vec3> &, const std::vector<double> &, const std::vector<vec3> &,
613  // double) const { return 0.0; } virtual vec3 surface_force(const std::vector<vec3> &, const std::vector<double> &,
614  // const std::vector<vec3> &, int, double) const { return {0.0,0.0,0.0}; }
615 
616 #ifdef NLOHMANN_JSON_HPP
617  private:
618  virtual void _to_json(nlohmann::json &) const = 0;
619 
620  public:
621  inline void to_json(nlohmann::json &j) const {
622  _to_json(j);
623  if (std::isfinite(cutoff))
624  j["cutoff"] = cutoff;
625  if (not doi.empty())
626  j["doi"] = doi;
627  if (not name.empty())
628  j["type"] = name;
629  if (std::isfinite(debye_length))
630  j["debyelength"] = debye_length;
631  }
632 #endif
633 };
634 
640 template <class T, bool debyehuckel = true> class EnergyImplementation : public SchemeBase {
641 
642  public:
643 
644  EnergyImplementation(Scheme type, double cutoff, double debyelength = infinity)
645  : SchemeBase(type, cutoff, debyelength) {
646  }
647 
659  inline double ion_potential(double z, double r) const override {
660  if (r < cutoff) {
661  double q = r * invcutoff;
662  if (debyehuckel) // determined at compile time
663  return z / r * static_cast<const T *>(this)->short_range_function(q) * std::exp(-kappa * r);
664  else
665  return z / r * static_cast<const T *>(this)->short_range_function(q);
666  } else {
667  return 0.0;
668  }
669  }
670 
683  inline double dipole_potential(const vec3 &mu, const vec3 &r) const override {
684  double r2 = r.squaredNorm();
685  if (r2 < cutoff2) {
686  double r1 = std::sqrt(r2);
687  double q = r1 * invcutoff;
688  double kr = kappa * r1;
689  return mu.dot(r) / (r2 * r1) *
690  (static_cast<const T *>(this)->short_range_function(q) * (1.0 + kr) -
691  q * static_cast<const T *>(this)->short_range_function_derivative(q)) *
692  std::exp(-kr);
693  } else {
694  return 0.0;
695  }
696  }
697 
709  inline double quadrupole_potential(const mat33 &quad, const vec3 &r) const override {
710  double r2 = r.squaredNorm();
711  if (r2 < cutoff2) {
712  double r1 = std::sqrt(r2);
713  double q = r1 * invcutoff;
714  double q2 = q * q;
715  double kr = kappa * r1;
716  double kr2 = kr * kr;
717  double srf = static_cast<const T *>(this)->short_range_function(q);
718  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
719  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
720 
721  double a = (srf * (1.0 + kr + kr2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kr) + q2 / 3.0 * ddsrf);
722  double b = (srf * kr2 - 2.0 * kr * q * dsrf + ddsrf * q2) / 3.0;
723  return 0.5 * ( ( 3.0/r2*r.transpose()*quad*r - quad.trace() ) * a + quad.trace() * b ) / r2 / r1 * std::exp(-kappa * r1);
724  } else {
725  return 0.0;
726  }
727  }
728 
740  inline vec3 ion_field(double z, const vec3 &r) const override {
741  double r2 = r.squaredNorm();
742  if (r2 < cutoff2) {
743  double r1 = std::sqrt(r2);
744  double q = r1 * invcutoff;
745  double kr = kappa * r1;
746  return z * r / (r2 * r1) *
747  (static_cast<const T *>(this)->short_range_function(q) * (1.0 + kr) -
748  q * static_cast<const T *>(this)->short_range_function_derivative(q)) *
749  std::exp(-kr);
750  } else {
751  return {0, 0, 0};
752  }
753  }
754 
768  inline vec3 dipole_field(const vec3 &mu, const vec3 &r) const override {
769  double r2 = r.squaredNorm();
770  if (r2 < cutoff2) {
771  double r1 = std::sqrt(r2);
772  double r3 = r1 * r2;
773  double q = r1 * invcutoff;
774  double q2 = q * q;
775  double kr = kappa * r1;
776  double kr2 = kr * kr;
777  double srf = static_cast<const T *>(this)->short_range_function(q);
778  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
779  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
780  vec3 fieldD = (3.0 * mu.dot(r) * r / r2 - mu) / r3;
781  fieldD *= (srf * (1.0 + kr + kr2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kr) +
782  q2 / 3.0 * ddsrf);
783  vec3 fieldI = mu / r3 * (srf * kr2 - 2.0 * kr * q * dsrf + ddsrf * q2) / 3.0;
784  return (fieldD + fieldI) * std::exp(-kr);
785  } else {
786  return {0, 0, 0};
787  }
788  }
789 
801  inline vec3 quadrupole_field(const mat33 &quad, const vec3 &r) const {
802  double r2 = r.squaredNorm();
803  if (r2 < cutoff2) {
804  double r1 = std::sqrt(r2);
805  vec3 rh = r / r1;
806  double q = r1 * invcutoff;
807  double q2 = q * q;
808  double kr = kappa * r1;
809  double kr2 = kr * kr;
810  double r4 = r2 * r2;
811  vec3 quadrh = quad*rh;
812  vec3 quadTrh = quad.transpose()*rh;
813  double quadfactor = 1.0/r2*r.transpose()*quad*r;
814  vec3 fieldD =
815  3.0 * ((5.0 * quadfactor - quad.trace()) * rh - quadrh - quadTrh) / r4;
816  double srf = static_cast<const T *>(this)->short_range_function(q);
817  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
818  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
819  double dddsrf = static_cast<const T *>(this)->short_range_function_third_derivative(q);
820  fieldD *= (srf * (1.0 + kr + kr2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kr) + q2 / 3.0 * ddsrf);
821  vec3 fieldI = quadfactor * rh / r4;
822  fieldI *= (srf * (1.0 + kr) * kr2 - q * dsrf * (3.0 * kr + 2.0) * kr + ddsrf * (1.0 + 3.0 * kr) * q2 - q2 * q * dddsrf);
823  return 0.5 * (fieldD + fieldI) * std::exp(-kr);
824  } else {
825  return {0, 0, 0};
826  }
827  }
828 
845  inline vec3 multipole_field(double z, const vec3 &mu, const mat33 &quad, const vec3 &r) const {
846  double r2 = r.squaredNorm();
847  if (r2 < cutoff2) {
848  double r1 = std::sqrt(r2);
849  vec3 rh = r / r1;
850  double q = r1 * invcutoff;
851  double q2 = q * q;
852  double r3 = r1 * r2;
853  double kr = kappa * r1;
854  double kr2 = kr * kr;
855  double quadfactor = 1.0/r2*r.transpose()*quad*r;
856  double srf = static_cast<const T *>(this)->short_range_function(q);
857  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
858  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
859  double dddsrf = static_cast<const T *>(this)->short_range_function_third_derivative(q);
860  vec3 fieldIon = z * r / r3 * ( srf * (1.0 + kr) - q * dsrf ); // field from ion
861  double postfactor = (srf * (1.0 + kr + kr2 / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kr) + q2 / 3.0 * ddsrf);
862  vec3 fieldDd = (3.0 * mu.dot(r) * r / r2 - mu) / r3 * postfactor;
863  vec3 fieldId = mu / r3 * (srf * kr2 - 2.0 * kr * q * dsrf + ddsrf * q2) / 3.0;
864  vec3 fieldDq = 3.0 * ((5.0 * quadfactor - quad.trace()) * rh - quad * rh - quad.transpose() * rh) / r3 / r1 * postfactor;
865  vec3 fieldIq = quadfactor * rh / r3 / r1;
866  fieldIq *= (srf * (1.0 + kr) * kr2 - q * dsrf * (3.0 * kr + 2.0) * kr + ddsrf * (1.0 + 3.0 * kr) * q2 - q2 * q * dddsrf);
867  return ( fieldIon + fieldDd + fieldId + 0.5 * (fieldDq + fieldIq) ) * std::exp(-kr);
868  } else {
869  return {0, 0, 0};
870  }
871  }
872 
886  inline double ion_ion_energy(double zA, double zB, double r) const override { return zB * ion_potential(zA, r); }
887 
906  inline double ion_dipole_energy(double z, const vec3 &mu, const vec3 &r) const override {
907  // Both expressions below gives same answer. Keep for possible optimization in future.
908  // return -mu.dot(ion_field(z,r)); // field from charge interacting with dipole
909  return z * dipole_potential(mu, -r); // potential of dipole interacting with charge
910  }
911 
926  inline double dipole_dipole_energy(const vec3 &muA, const vec3 &muB, const vec3 &r) const override {
927  return -muA.dot(dipole_field(muB, r));
928  }
929 
943  inline double ion_quadrupole_energy(double z, const mat33 &quad, const vec3 &r) const override {
944  return z * quadrupole_potential(quad, -r); // potential of quadrupole interacting with charge
945  }
946 
960  inline double multipole_multipole_energy(double zA, double zB, const vec3 &muA, const vec3 &muB, const mat33 &quadA, const mat33 &quadB,
961  const vec3 &r) const override {
962  double r2 = r.squaredNorm();
963  if (r2 < cutoff2) {
964  double r1 = std::sqrt(r2);
965  double q = r1 / cutoff;
966  double kr = kappa * r1;
967  double quadAtrace = quadA.trace();
968  double quadBtrace = quadB.trace();
969 
970  double srf = static_cast<const T *>(this)->short_range_function(q);
971  double dsrfq = static_cast<const T *>(this)->short_range_function_derivative(q) * q;
972  double ddsrfq2 = static_cast<const T *>(this)->short_range_function_second_derivative(q) * q * q / 3.0;
973 
974  double angcor = (srf * (1.0 + kr) - dsrfq);
975  double unicor = (srf * kr * kr / 3.0 - 2.0 / 3.0 * dsrfq * kr + ddsrfq2);
976  double muBdotr = muB.dot(r);
977  vec3 field_dipoleB = (3.0 * muBdotr * r / r2 - muB) * (angcor + unicor) + muB * unicor;
978 
979  double ion_ion = zA * zB * srf * r2; // will later be divided by r3
980  double ion_dipole = (zB * muA.dot(r) - zA * muBdotr) * angcor; // will later be divided by r3
981  double dipole_dipole = -muA.dot(field_dipoleB); // will later be divided by r3
982  double ion_quadrupole = zA * 0.5 * ( ( 3.0/r2*r.transpose()*quadB*r - quadBtrace ) * (angcor + unicor) + quadBtrace * unicor ); // will later be divided by r3
983  ion_quadrupole += zB * 0.5 * ( ( 3.0/r2*r.transpose()*quadA*r - quadAtrace ) * (angcor + unicor) + quadAtrace * unicor );
984 
985  return (ion_ion + ion_dipole + dipole_dipole + ion_quadrupole) * std::exp(-kr) / r2 / r1;
986  } else {
987  return 0.0;
988  }
989  }
990 
1004  inline vec3 ion_ion_force(double zA, double zB, const vec3 &r) const override { return zB * ion_field(zA, r); }
1005 
1019  inline vec3 ion_dipole_force(double z, const vec3 &mu, const vec3 &r) const override {
1020  return z * dipole_field(mu, r);
1021  }
1022 
1049  inline vec3 dipole_dipole_force(const vec3 &muA, const vec3 &muB, const vec3 &r) const override {
1050  double r2 = r.squaredNorm();
1051  if (r2 < cutoff2) {
1052  double r1 = std::sqrt(r2);
1053  vec3 rh = r / r1;
1054  double q = r1 * invcutoff;
1055  double q2 = q * q;
1056  double kr = kappa * r1;
1057  double r4 = r2 * r2;
1058  double muAdotRh = muA.dot(rh);
1059  double muBdotRh = muB.dot(rh);
1060  vec3 forceD =
1061  3.0 * ((5.0 * muAdotRh * muBdotRh - muA.dot(muB)) * rh - muBdotRh * muA - muAdotRh * muB) / r4;
1062  double srf = static_cast<const T *>(this)->short_range_function(q);
1063  double dsrf = static_cast<const T *>(this)->short_range_function_derivative(q);
1064  double ddsrf = static_cast<const T *>(this)->short_range_function_second_derivative(q);
1065  double dddsrf = static_cast<const T *>(this)->short_range_function_third_derivative(q);
1066  forceD *= (srf * (1.0 + kr + kr * kr / 3.0) - q * dsrf * (1.0 + 2.0 / 3.0 * kr) + q2 / 3.0 * ddsrf);
1067  vec3 forceI = muAdotRh * muBdotRh * rh / r4;
1068  forceI *= (srf * (1.0 + kr) * kr * kr - q * dsrf * (3.0 * kr + 2.0) * kr + ddsrf * (1.0 + 3.0 * kr) * q2 - q2 * q * dddsrf);
1069  return (forceD + forceI) * std::exp(-kr);
1070  } else {
1071  return {0, 0, 0};
1072  }
1073  }
1074 
1088  inline vec3 ion_quadrupole_force(double z, const mat33 &quad, const vec3 &r) const override { return z * quadrupole_field(quad, r); }
1089 
1103  inline vec3 multipole_multipole_force(double zA, double zB, const vec3 &muA, const vec3 &muB, const mat33 &quadA, const mat33 &quadB,
1104  const vec3 &r) const override {
1105  double r2 = r.squaredNorm();
1106  if (r2 < cutoff2) {
1107  double r1 = std::sqrt(r2);
1108  double q = r1 * invcutoff;
1109  double q2 = q * q;
1110  double kr = kappa * r1;
1111  vec3 rh = r / r1;
1112  double muAdotRh = muA.dot(rh);
1113  double muBdotRh = muB.dot(rh);
1114 
1115  double srf = static_cast<const T *>(this)->short_range_function(q);
1116  double dsrfq = static_cast<const T *>(this)->short_range_function_derivative(q) * q;
1117  double ddsrfq2 = static_cast<const T *>(this)->short_range_function_second_derivative(q) * q2 / 3.0;
1118  double dddsrfq3 = static_cast<const T *>(this)->short_range_function_third_derivative(q) * q2 * q;
1119 
1120  double angcor = (srf * (1.0 + kr) - dsrfq);
1121  double unicor = (srf * kr - 2.0 * dsrfq) * kr / 3.0 + ddsrfq2;
1122  double totcor = unicor + angcor;
1123  double r3corr = ( angcor * kr * kr - dsrfq * 2.0 * (1.0 + kr) * kr + 3.0 * ddsrfq2 * (1.0 + 3.0 * kr) - dddsrfq3);
1124 
1125  vec3 ion_ion = zB * zA * r * angcor * r1;
1126  vec3 ion_dipole = zA * ((3.0 * muBdotRh * rh - muB) * totcor + muB * unicor);
1127  ion_dipole += zB * ((3.0 * muAdotRh * rh - muA) * totcor + muA * unicor);
1128  ion_dipole *= r1;
1129  vec3 forceD = 3.0 * ((5.0 * muAdotRh * muBdotRh - muA.dot(muB)) * rh - muBdotRh * muA - muAdotRh * muB) * totcor;
1130  vec3 dipole_dipole = (forceD + muAdotRh * muBdotRh * rh * r3corr);
1131  double quadfactor = 1.0/r2*r.transpose()*quadB*r;
1132  vec3 fieldD = 3.0 * (-(5.0 * quadfactor - quadB.trace()) * rh + quadB*rh + quadB.transpose()*rh) * totcor;
1133  vec3 ion_quadrupole = zA * 0.5 * (fieldD - quadfactor * rh * r3corr);
1134  quadfactor = 1.0/r2*r.transpose()*quadA*r;
1135  fieldD = 3.0 * ((5.0 * quadfactor - quadA.trace()) * rh - quadA*rh - quadA.transpose()*rh) * totcor;
1136  ion_quadrupole += zB * 0.5 * (fieldD + quadfactor * rh * r3corr);
1137 
1138  return (ion_ion + ion_dipole + dipole_dipole + ion_quadrupole) * std::exp(-kr) / r2 / r2;
1139  } else {
1140  return {0, 0, 0};
1141  }
1142  }
1143 
1157  inline vec3 dipole_torque(const vec3 &mu, const vec3 &E) const { return mu.cross(E); }
1158 
1171  inline double self_energy(const std::array<double, 2> &squared_moments) const {
1172  assert(selfEnergyFunctor != nullptr);
1173  return selfEnergyFunctor(squared_moments);
1174  }
1175 
1183  inline double neutralization_energy(const std::vector<double> &charges, double volume) const override {
1184  double charge_total = 0.0;
1185  for (unsigned int i = 0; i < charges.size(); i++)
1186  charge_total += charges.at(i);
1187  return ((this)->chi / 2.0 / volume * charge_total * charge_total);
1188  }
1189 };
1190 
1191 // -------------- Plain ---------------
1192 
1197 class Plain : public EnergyImplementation<Plain> {
1198  public:
1199  inline Plain(double debye_length = infinity)
1200  : EnergyImplementation(Scheme::plain, std::sqrt(std::numeric_limits<double>::max()), debye_length) {
1201  name = "plain";
1202  dipolar_selfenergy = true;
1203  doi = "Premier mémoire sur l’électricité et le magnétisme by Charles-Augustin de Coulomb"; // :P
1204  setSelfEnergyPrefactor({0.0, 0.0});
1205  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1206  chi = -2.0 * std::acos(-1.0) * cutoff2; // should not be used!
1207  };
1208  inline double short_range_function(double) const override { return 1.0; };
1209  inline double short_range_function_derivative(double) const override { return 0.0; }
1210  inline double short_range_function_second_derivative(double) const override { return 0.0; }
1211  inline double short_range_function_third_derivative(double) const override { return 0.0; }
1212 #ifdef NLOHMANN_JSON_HPP
1213  inline Plain(const nlohmann::json &j) : Plain(j.value("debyelength", infinity)) {}
1214 
1215  private:
1216  inline void _to_json(nlohmann::json &) const override {}
1217 #endif
1218 };
1219 
1220 // -------------- Ewald real-space (using Gaussian) ---------------
1221 
1230 class Ewald : public EnergyImplementation<Ewald> {
1231  double eta, eta2, eta3;
1232  double zeta, zeta2, zeta3;
1233  double eps_sur;
1234  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1235  const double pi = 4.0 * std::atan(1.0);
1236 
1237  public:
1242  inline Ewald(double cutoff, double alpha, double eps_sur = infinity, double debye_length = infinity)
1243  : EnergyImplementation(Scheme::ewald, cutoff), eps_sur(eps_sur) {
1244  name = "Ewald real-space";
1245  dipolar_selfenergy = true;
1246  doi = "10.1002/andp.19213690304";
1247  eta = alpha * cutoff;
1248  eta2 = eta * eta;
1249  eta3 = eta2 * eta;
1250  if (eps_sur < 1.0)
1251  eps_sur = infinity;
1252  double Q = 1.0 - std::erfc(eta) - 2.0 * eta / pi_sqrt * std::exp(-eta2); // Eq. 12 in DOI: 10.1016/0009-2614(83)80585-5 using 'K = cutoff region'
1253  T0 = (std::isinf(eps_sur)) ? Q : ( Q - 1.0 + 2.0 * (eps_sur - 1.0) / (2.0 * eps_sur + 1.0) ); // Eq. 17 in DOI: 10.1016/0009-2614(83)80585-5
1254  zeta = cutoff / debye_length;
1255  zeta2 = zeta * zeta;
1256  zeta3 = zeta2 * zeta;
1257  if(zeta < 1e-6) { // if close to zero the general expression numerically diverges, and this expresion is used instead
1258  chi = -pi * cutoff2 * ( 1.0 - std::erfc( eta ) * ( 1.0 - 2.0 * eta2 ) - 2.0 * eta * std::exp( -eta2 ) / pi_sqrt ) / eta2;
1259  } else {
1260  chi = 4.0 * ( 0.5 * ( 1.0 - zeta ) * std::erfc( eta + zeta / ( 2.0 * eta ) ) * std::exp( zeta ) + std::erf( eta ) * std::exp(-zeta2 / ( 4.0 * eta2 ) ) +
1261  0.5 * ( 1.0 + zeta ) * std::erfc( eta - zeta / ( 2.0 * eta ) ) * std::exp( -zeta ) - 1.0 ) * pi * cutoff2 / zeta2;
1262  }
1263  // chi = -pi * cutoff2 / eta2 according to DOI:10.1021/ct400626b, for uncscreened system
1264 
1265  setSelfEnergyPrefactor({
1266  -eta / pi_sqrt * (std::exp(-zeta2 / 4.0 / eta2) - pi_sqrt * zeta / (2.0 * eta) * std::erfc(zeta / (2.0 * eta) ) ),
1267  -eta3 / pi_sqrt * 2.0 / 3.0 *
1268  (pi_sqrt * zeta3 / 4.0 / eta3 * std::erfc(zeta / (2.0 * eta) ) + (1.0 - zeta2 / 2.0 / eta2) * std::exp(-zeta2 / 4.0 / eta2))}); // ion-quadrupole self-energy term: XYZ
1269  }
1270 
1271  inline double short_range_function(double q) const override {
1272  return 0.5 * (std::erfc(eta * q + zeta / (2.0 * eta)) * std::exp(2.0 * zeta * q) + std::erfc(eta * q - zeta / (2.0 * eta)));
1273  }
1274  inline double short_range_function_derivative(double q) const override {
1275  double expC = std::exp(-powi(eta * q - zeta / (2.0 * eta), 2));
1276  double erfcC = std::erfc(eta * q + zeta / (2.0 * eta));
1277  return (-2.0 * eta / pi_sqrt * expC + zeta * erfcC * std::exp(2.0 * zeta * q));
1278  }
1279  inline double short_range_function_second_derivative(double q) const override {
1280  double expC = std::exp(-powi(eta * q - zeta / (2.0 * eta), 2));
1281  double erfcC = std::erfc(eta * q + zeta / (2.0 * eta));
1282  return (4.0 * eta2 / pi_sqrt * (eta * q - zeta / eta) * expC + 2.0 * zeta2 * erfcC * std::exp(2.0 * zeta * q));
1283  }
1284  inline double short_range_function_third_derivative(double q) const override {
1285  double expC = std::exp(-powi(eta * q - zeta / (2.0 * eta), 2));
1286  double erfcC = std::erfc(eta * q + zeta / (2.0 * eta));
1287  return (4.0 * eta3 / pi_sqrt * (1.0 - 2.0 * (eta * q - zeta / eta) * (eta * q - zeta / (2.0 * eta) ) - zeta2 / eta2) * expC + 4.0 * zeta3 * erfcC * std::exp(2.0 * zeta * q));
1288  }
1289 
1299  inline double reciprocal_energy(const std::vector<vec3> &positions, const std::vector<double> &charges,
1300  const std::vector<vec3> &dipoles, const vec3 &L, int nmax) const {
1301 
1302  assert(positions.size() == charges.size());
1303  assert(positions.size() == dipoles.size());
1304 
1305  double volume = L[0] * L[1] * L[2];
1306  std::vector<vec3> kvec;
1307  std::vector<double> Ak;
1308  // kvec.reserve( expected_size_of_kvec ); // speeds up push_back below
1309  // Ak.reserve( expected_size_of_Ak ); // speeds up push_back below
1310  for (int nx = -nmax; nx < nmax + 1; nx++) {
1311  for (int ny = -nmax; ny < nmax + 1; ny++) {
1312  for (int nz = -nmax; nz < nmax + 1; nz++) {
1313  vec3 kv = {2.0 * pi * nx / L[0], 2.0 * pi * ny / L[1], 2.0 * pi * nz / L[2]};
1314  double k2 = kv.squaredNorm() + zeta2 / cutoff2;
1315  vec3 nv = {double(nx), double(ny), double(nz)};
1316  double nv1 = nv.squaredNorm();
1317  if (nv1 > 0) {
1318  if (nv1 <= nmax * nmax) {
1319  kvec.push_back(kv);
1320  Ak.push_back(std::exp(-( k2 * cutoff2 + zeta2 ) / 4.0 / eta2 ) / k2);
1321  }
1322  }
1323  }
1324  }
1325  }
1326 
1327  assert(kvec.size() == Ak.size());
1328 
1329  double E = 0.0;
1330  for (size_t k = 0; k < kvec.size(); k++) {
1331  std::complex<double> Qq(0.0, 0.0);
1332  std::complex<double> Qmu(0.0, 0.0);
1333  for (size_t i = 0; i < positions.size(); i++) {
1334  double kDotR = kvec[k].dot(positions[i]);
1335  double coskDotR = std::cos(kDotR);
1336  double sinkDotR = std::sin(kDotR);
1337  Qq += charges[i] * std::complex<double>(coskDotR, sinkDotR);
1338  Qmu += dipoles[i].dot(kvec[k]) * std::complex<double>(-sinkDotR, coskDotR);
1339  }
1340  std::complex<double> Q = Qq + Qmu;
1341  E += (powi(std::abs(Q), 2) * Ak[k]);
1342  }
1343  return (E * 2.0 * pi / volume);
1344  }
1345 
1353  inline double surface_energy(const std::vector<vec3> &positions, const std::vector<double> &charges,
1354  const std::vector<vec3> &dipoles, double volume) const {
1355  assert(positions.size() == charges.size());
1356  assert(positions.size() == dipoles.size());
1357  vec3 sum_r_charges = {0.0, 0.0, 0.0};
1358  vec3 sum_dipoles = {0.0, 0.0, 0.0};
1359  for (size_t i = 0; i < positions.size(); i++) {
1360  sum_r_charges += positions[i] * charges[i];
1361  sum_dipoles += dipoles[i];
1362  }
1363  double sqDipoles =
1364  sum_r_charges.dot(sum_r_charges) + 2.0 * sum_r_charges.dot(sum_dipoles) + sum_dipoles.dot(sum_dipoles);
1365 
1366  return (2.0 * pi / (2.0 * eps_sur + 1.0) / volume * sqDipoles);
1367  }
1368 
1369 #ifdef NLOHMANN_JSON_HPP
1370  inline Ewald(const nlohmann::json &j)
1371  : Ewald(j.at("cutoff").get<double>(), j.at("alpha").get<double>(), j.value("epss", infinity),
1372  j.value("debyelength", infinity)) {}
1373 
1374  private:
1375  inline void _to_json(nlohmann::json &j) const override {
1376  j = {{ "alpha", eta/cutoff }};
1377  if (std::isinf(eps_sur))
1378  j["epss"] = "inf";
1379  else
1380  j["epss"] = eps_sur;
1381  }
1382 #endif
1383 };
1384 
1385 // -------------- Ewald real-space (using truncated Gaussian) ---------------
1386 
1390 class EwaldT : public EnergyImplementation<EwaldT> {
1391  double eta, eta2, eta3;
1392  double zeta, zeta2, zeta3;
1393  double eps_sur;
1394  double F0;
1395  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1396  const double pi = 4.0 * std::atan(1.0);
1397 
1398  public:
1403  inline EwaldT(double cutoff, double alpha, double eps_sur = infinity, double debye_length = infinity)
1404  : EnergyImplementation(Scheme::ewaldt, cutoff), eps_sur(eps_sur) {
1405  name = "EwaldT real-space";
1406  dipolar_selfenergy = true;
1407  doi = "XYZ";
1408  eta = alpha * cutoff;
1409  eta2 = eta * eta;
1410  eta3 = eta2 * eta;
1411  if (eps_sur < 1.0)
1412  eps_sur = infinity;
1413  F0 = 1.0 - std::erfc(eta) - 2.0 * eta / pi_sqrt * std::exp(-eta2);
1414  T0 = (std::isinf(eps_sur)) ? 1.0 : 2.0 * (eps_sur - 1.0) / (2.0 * eps_sur + 1.0);
1415  chi = -( 1.0 - 4.0 * eta3 * std::exp( -eta2 ) / ( 3.0 * pi_sqrt * F0 ) ) * cutoff2 * pi / eta2;
1416  setSelfEnergyPrefactor({-eta / pi_sqrt * (1.0 - std::exp( -eta2 ) ) / F0, -eta3 * 2.0 / 3.0 / ( std::erf( eta ) * pi_sqrt - 2.0 * eta * std::exp( -eta2 ) ) });
1417  }
1418 
1419  inline double short_range_function(double q) const override {
1420  return ( std::erfc(eta * q) - std::erfc(eta) - (1.0 - q) * 2.0 * eta / pi_sqrt * std::exp(-eta2) ) / F0;
1421  }
1422  inline double short_range_function_derivative(double q) const override {
1423  return - 2.0 * eta * ( std::exp( -eta2 * q * q ) - std::exp( -eta2 ) ) / pi_sqrt / F0;
1424  }
1425  inline double short_range_function_second_derivative(double q) const override {
1426  return 4.0 * eta3 * q * std::exp( -eta2 * q * q ) / pi_sqrt / F0;
1427  }
1428  inline double short_range_function_third_derivative(double q) const override {
1429  return - 8.0 * ( eta2 * q * q - 0.5 ) * eta3 * std::exp( -eta2 * q * q ) / pi_sqrt / F0;
1430  }
1431 
1439  inline double surface_energy(const std::vector<vec3> &positions, const std::vector<double> &charges,
1440  const std::vector<vec3> &dipoles, double volume) const {
1441  assert(positions.size() == charges.size());
1442  assert(positions.size() == dipoles.size());
1443  vec3 sum_r_charges = {0.0, 0.0, 0.0};
1444  vec3 sum_dipoles = {0.0, 0.0, 0.0};
1445  for (size_t i = 0; i < positions.size(); i++) {
1446  sum_r_charges += positions[i] * charges[i];
1447  sum_dipoles += dipoles[i];
1448  }
1449  double sqDipoles =
1450  sum_r_charges.dot(sum_r_charges) + 2.0 * sum_r_charges.dot(sum_dipoles) + sum_dipoles.dot(sum_dipoles);
1451 
1452  return (2.0 * pi / (2.0 * eps_sur + 1.0) / volume * sqDipoles);
1453  }
1454 
1455 #ifdef NLOHMANN_JSON_HPP
1456  inline EwaldT(const nlohmann::json &j)
1457  : EwaldT(j.at("cutoff").get<double>(), j.at("alpha").get<double>(), j.value("epss", infinity),
1458  j.value("debyelength", infinity)) {}
1459 
1460  private:
1461  inline void _to_json(nlohmann::json &j) const override {
1462  j = {{ "alpha", eta/cutoff }};
1463  if (std::isinf(eps_sur))
1464  j["epss"] = "inf";
1465  else
1466  j["epss"] = eps_sur;
1467  }
1468 #endif
1469 };
1470 
1471 // -------------- Reaction-field ---------------
1472 
1476 class ReactionField : public EnergyImplementation<ReactionField> {
1477  private:
1478  double epsRF;
1479  double epsr;
1480  bool shifted;
1481  const double pi = 4.0 * std::atan(1.0);
1482 
1483  public:
1490  inline ReactionField(double cutoff, double epsRF, double epsr, bool shifted)
1491  : EnergyImplementation(Scheme::reactionfield, cutoff), epsRF(epsRF), epsr(epsr), shifted(shifted) {
1492  name = "Reaction-field";
1493  dipolar_selfenergy = true;
1494  doi = "10.1080/00268977300102101";
1495  //epsRF = epsRF;
1496  //epsr = epsr;
1497  //shifted = shifted;
1498  setSelfEnergyPrefactor({-3.0 * epsRF * double(shifted) / (4.0 * epsRF + 2.0 * epsr),
1499  -(2.0 * epsRF - 2.0 * epsr) / (2.0 * (2.0 * epsRF + epsr))});
1500  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1501  chi = -6.0 * cutoff * cutoff * pi * ((-10.0 * double(shifted) / 3.0 + 4.0) * epsRF + epsr) /
1502  ((5.0 * (2.0 * epsRF + epsr)));
1503  }
1504 
1505  inline double short_range_function(double q) const override {
1506  return (1.0 + (epsRF - epsr) * q * q * q / (2.0 * epsRF + epsr) -
1507  3.0 * epsRF * q / (2.0 * epsRF + epsr) * double(shifted));
1508  }
1509  inline double short_range_function_derivative(double q) const override {
1510  return (3.0 * (epsRF - epsr) * q * q / (2.0 * epsRF + epsr) -
1511  3.0 * epsRF * double(shifted) / (2.0 * epsRF + epsr));
1512  }
1513  inline double short_range_function_second_derivative(double q) const override {
1514  return 6.0 * (epsRF - epsr) * q / (2.0 * epsRF + epsr);
1515  }
1516  inline double short_range_function_third_derivative(double) const override {
1517  return 6.0 * (epsRF - epsr) / (2.0 * epsRF + epsr);
1518  }
1519 
1520 #ifdef NLOHMANN_JSON_HPP
1521 
1522  inline ReactionField(const nlohmann::json &j)
1523  : ReactionField(j.at("cutoff").get<double>(), j.at("epsRF").get<double>(), j.at("epsr").get<double>(),
1524  j.at("shifted").get<bool>()) {}
1525 
1526  private:
1527  inline void _to_json(nlohmann::json &j) const override {
1528  j = {{"epsr", epsr}, {"epsRF", epsRF}, { "shifted", shifted }};
1529  }
1530 #endif
1531 };
1532 
1533 // -------------- Zahn ---------------
1534 
1538 class Zahn : public EnergyImplementation<Zahn> {
1539  private:
1540  double alpha;
1541  double alphaRed, alphaRed2;
1542  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1543  const double pi = 4.0 * std::atan(1.0);
1544 
1545  public:
1551  inline Zahn(double cutoff, double alpha) : EnergyImplementation(Scheme::zahn, cutoff), alpha(alpha) {
1552  name = "Zahn";
1553  dipolar_selfenergy = false;
1554  doi = "10.1021/jp025949h";
1555  alphaRed = alpha * cutoff;
1556  alphaRed2 = alphaRed * alphaRed;
1557  setSelfEnergyPrefactor({-alphaRed * (1.0 - std::exp(-alphaRed2)) / pi_sqrt + 0.5 * std::erfc(alphaRed),
1558  0.0}); // Dipole self-energy undefined!
1559  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1560  chi = -(2.0 * (alphaRed * (alphaRed2 - 3.0) * std::exp(-alphaRed2) -
1561  0.5 * std::sqrt(pi) * ((7.0 - 3.0 / alphaRed2) * std::erf(alphaRed) - 7.0) * alphaRed2)) *
1562  cutoff * cutoff * std::sqrt(pi) / (3.0 * alphaRed2);
1563  }
1564 
1565  inline double short_range_function(double q) const override {
1566  return (std::erfc(alphaRed * q) -
1567  (q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1568  }
1569  inline double short_range_function_derivative(double q) const override {
1570  return (-(4.0 * (0.5 * std::exp(-alphaRed2 * q * q) * alphaRed +
1571  (alphaRed * std::exp(-alphaRed2) + 0.5 * pi_sqrt * std::erfc(alphaRed)) * (q - 0.5))) /
1572  pi_sqrt);
1573  }
1574  inline double short_range_function_second_derivative(double q) const override {
1575  return (4.0 * (alphaRed2 * alphaRed * q * std::exp(-alphaRed2 * q * q) - alphaRed * std::exp(-alphaRed2) -
1576  0.5 * pi_sqrt * std::erfc(alphaRed))) /
1577  pi_sqrt;
1578  }
1579  inline double short_range_function_third_derivative(double q) const override {
1580  return (-8.0 * std::exp(-alphaRed2 * q * q) * (alphaRed2 * q * q - 0.5) * alphaRed2 * alphaRed / pi_sqrt);
1581  }
1582 
1583 #ifdef NLOHMANN_JSON_HPP
1584 
1585  inline Zahn(const nlohmann::json &j) : Zahn(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1586 
1587  private:
1588  inline void _to_json(nlohmann::json &j) const override {
1589  j = {{ "alpha", alpha }};
1590  }
1591 #endif
1592 };
1593 
1594 // -------------- Fennell ---------------
1595 
1599 class Fennell : public EnergyImplementation<Fennell> {
1600  private:
1601  double alpha;
1602  double alphaRed, alphaRed2;
1603  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1604  const double pi = 4.0 * std::atan(1.0);
1605 
1606  public:
1611  inline Fennell(double cutoff, double alpha) : EnergyImplementation(Scheme::fennell, cutoff), alpha(alpha) {
1612  name = "Fennell";
1613  dipolar_selfenergy = false;
1614  doi = "10.1063/1.2206581";
1615  alphaRed = alpha * cutoff;
1616  alphaRed2 = alphaRed * alphaRed;
1617  setSelfEnergyPrefactor({-alphaRed * (1.0 + std::exp(-alphaRed2)) / pi_sqrt - std::erfc(alphaRed),
1618  0.0}); // Dipole self-energy undefined!
1619  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1620  chi = (2.0 * ((alphaRed2 + 3.0) * alphaRed * std::exp(-alphaRed2) +
1621  0.5 * (std::erf(alphaRed) * alphaRed2 - alphaRed2 - 3.0 * std::erf(alphaRed)) * std::sqrt(pi))) *
1622  std::sqrt(pi) * cutoff * cutoff / (3.0 * alphaRed2);
1623  }
1624 
1625  inline double short_range_function(double q) const override {
1626  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed) +
1627  (q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1628  }
1629  inline double short_range_function_derivative(double q) const override {
1630  return (2.0 * alphaRed * (2.0 * (q - 0.5) * std::exp(-alphaRed2) - std::exp(-alphaRed2 * q * q)) / pi_sqrt +
1631  2.0 * std::erfc(alphaRed) * (q - 1.0));
1632  }
1633  inline double short_range_function_second_derivative(double q) const override {
1634  return (4.0 * alphaRed * (alphaRed2 * q * std::exp(-alphaRed2 * q * q) + std::exp(-alphaRed2)) / pi_sqrt +
1635  2.0 * std::erfc(alphaRed));
1636  }
1637  inline double short_range_function_third_derivative(double q) const override {
1638  return 4.0 * alphaRed2 * alphaRed * (1.0 - 2.0 * alphaRed2 * q * q) * std::exp(-alphaRed2 * q * q) / pi_sqrt;
1639  }
1640 
1641 #ifdef NLOHMANN_JSON_HPP
1642 
1643  inline Fennell(const nlohmann::json &j) : Fennell(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1644 
1645  private:
1646  inline void _to_json(nlohmann::json &j) const override {
1647  j = {{ "alpha", alpha }};
1648  }
1649 #endif
1650 };
1651 
1652 // -------------- Zero-dipole ---------------
1653 
1657 class ZeroDipole : public EnergyImplementation<ZeroDipole> {
1658  private:
1659  double alpha;
1660  double alphaRed, alphaRed2;
1661  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1662  const double pi = 4.0 * std::atan(1.0);
1663 
1664  public:
1669  inline ZeroDipole(double cutoff, double alpha) : EnergyImplementation(Scheme::zerodipole, cutoff), alpha(alpha) {
1670  name = "ZeroDipole";
1671  dipolar_selfenergy = true;
1672  doi = "10.1063/1.3582791";
1673  alphaRed = alpha * cutoff;
1674  alphaRed2 = alphaRed * alphaRed;
1675  setSelfEnergyPrefactor({-alphaRed * (1.0 + 0.5 * std::exp(-alphaRed2)) / pi_sqrt - 0.75 * std::erfc(alphaRed),
1676  -alphaRed * (2.0 * alphaRed2 * (1.0 / 3.0) + std::exp(-alphaRed2)) / pi_sqrt -
1677  0.5 * std::erfc(alphaRed)});
1678  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1679  chi = cutoff * cutoff *
1680  ((6.0 * alphaRed2 - 15.0) * std::erf(alphaRed) * pi - 6.0 * pi * alphaRed2 +
1681  (8.0 * alphaRed2 + 30.0) * alphaRed * std::exp(-alphaRed2) * std::sqrt(pi)) /
1682  (15.0 * alphaRed2);
1683  }
1684 
1685  inline double short_range_function(double q) const override {
1686  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed) +
1687  0.5 * (q * q - 1.0) * q * (std::erfc(alphaRed) + 2.0 * alphaRed * std::exp(-alphaRed2) / pi_sqrt));
1688  }
1689  inline double short_range_function_derivative(double q) const override {
1690  return (alphaRed * ((3.0 * q * q - 1.0) * std::exp(-alphaRed2) - 2.0 * std::exp(-alphaRed2 * q * q)) / pi_sqrt +
1691  1.5 * std::erfc(alphaRed) * (q * q - 1.0));
1692  }
1693  inline double short_range_function_second_derivative(double q) const override {
1694  return (2.0 * alphaRed * q * (2.0 * alphaRed2 * std::exp(-alphaRed2 * q * q) + 3.0 * std::exp(-alphaRed2)) /
1695  pi_sqrt +
1696  3.0 * q * std::erfc(alphaRed));
1697  }
1698  inline double short_range_function_third_derivative(double q) const override {
1699  return (2.0 * alphaRed *
1700  (2.0 * alphaRed2 * (1.0 - 2.0 * alphaRed2 * q * q) * std::exp(-alphaRed2 * q * q) +
1701  3.0 * std::exp(-alphaRed2)) /
1702  pi_sqrt +
1703  3.0 * std::erfc(alphaRed));
1704  }
1705 
1706 #ifdef NLOHMANN_JSON_HPP
1707 
1708  inline ZeroDipole(const nlohmann::json &j)
1709  : ZeroDipole(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1710 
1711  private:
1712  inline void _to_json(nlohmann::json &j) const override {
1713  j = {{ "alpha", alpha }};
1714  }
1715 #endif
1716 };
1717 
1718 // -------------- Wolf ---------------
1719 
1723 class Wolf : public EnergyImplementation<Wolf> {
1724  private:
1725  double alpha;
1726  double alphaRed, alphaRed2;
1727  const double pi_sqrt = 2.0 * std::sqrt(std::atan(1.0));
1728  const double pi = 4.0 * std::atan(1.0);
1729 
1730  public:
1735  inline Wolf(double cutoff, double alpha) : EnergyImplementation(Scheme::wolf, cutoff), alpha(alpha) {
1736  name = "Wolf";
1737  dipolar_selfenergy = true;
1738  doi = "10.1063/1.478738";
1739  alphaRed = alpha * cutoff;
1740  alphaRed2 = alphaRed * alphaRed;
1741  setSelfEnergyPrefactor({-alphaRed / pi_sqrt - std::erfc(alphaRed) / 2.0,
1742  -powi(alphaRed, 3) * 2.0 / 3.0 / pi_sqrt});
1743  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1744  chi = 2.0 * std::sqrt(pi) * cutoff * cutoff *
1745  (3.0 * std::exp(-alphaRed2) * alphaRed -
1746  std::sqrt(pi) * (std::erfc(alphaRed) * alphaRed2 + 3.0 * std::erf(alphaRed) * 0.5)) /
1747  (3.0 * alphaRed2);
1748  }
1749 
1750  inline double short_range_function(double q) const override {
1751  return (std::erfc(alphaRed * q) - q * std::erfc(alphaRed));
1752  }
1753  inline double short_range_function_derivative(double q) const override {
1754  return (-2.0 * std::exp(-alphaRed2 * q * q) * alphaRed / pi_sqrt - std::erfc(alphaRed));
1755  }
1756  inline double short_range_function_second_derivative(double q) const override {
1757  return 4.0 * std::exp(-alphaRed2 * q * q) * alphaRed2 * alphaRed * q / pi_sqrt;
1758  }
1759  inline double short_range_function_third_derivative(double q) const override {
1760  return -8.0 * std::exp(-alphaRed2 * q * q) * alphaRed2 * alphaRed * (alphaRed2 * q * q - 0.5) / pi_sqrt;
1761  }
1762 
1763 #ifdef NLOHMANN_JSON_HPP
1764 
1765  inline Wolf(const nlohmann::json &j) : Wolf(j.at("cutoff").get<double>(), j.at("alpha").get<double>()) {}
1766 
1767  private:
1768  inline void _to_json(nlohmann::json &j) const override {
1769  j = {{ "alpha", alpha }};
1770  }
1771 #endif
1772 };
1773 
1774 // -------------- qPotential ---------------
1775 template <int order> class qPotentialFixedOrder : public EnergyImplementation<qPotentialFixedOrder<order>> {
1776  public:
1778  using base::chi;
1779  using base::name;
1780  using base::T0;
1785  inline qPotentialFixedOrder(double cutoff) : base(Scheme::qpotential, cutoff) {
1786  name = "qpotential";
1787  this->dipolar_selfenergy = true;
1788  this->doi = "10.1039/c9cp03875b";
1789  this->setSelfEnergyPrefactor({-0.5, -0.5});
1790  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1791  chi = -86459.0 * std::acos(-1.0) * cutoff * cutoff / 235620.0;
1792  }
1793 
1794  inline double short_range_function(double q) const override { return qPochhammerSymbol(q, 0, order); }
1795  inline double short_range_function_derivative(double q) const override {
1796  return qPochhammerSymbolDerivative(q, 0, order);
1797  }
1798  inline double short_range_function_second_derivative(double q) const override {
1799  return qPochhammerSymbolSecondDerivative(q, 0, order);
1800  }
1801  inline double short_range_function_third_derivative(double q) const override {
1802  return qPochhammerSymbolThirdDerivative(q, 0, order);
1803  }
1804 
1805 #ifdef NLOHMANN_JSON_HPP
1806  inline qPotentialFixedOrder(const nlohmann::json &j) : qPotentialFixedOrder(j.at("cutoff").get<double>()) {}
1807 
1808  private:
1809  inline void _to_json(nlohmann::json &j) const override {
1810  j = {{ "order", order }};
1811  }
1812 #endif
1813 };
1814 
1825 class qPotential : public EnergyImplementation<qPotential> {
1826  private:
1827  int order;
1828 
1829  public:
1834  inline qPotential(double cutoff, int order) : EnergyImplementation(Scheme::qpotential, cutoff), order(order) {
1835  name = "qpotential";
1836  dipolar_selfenergy = true;
1837  doi = "10.1039/c9cp03875b";
1838  setSelfEnergyPrefactor({-0.5, -0.5});
1839  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
1840  chi = 0.0; // -Pi*Rc^2 * [ 2/3 7/15 17/42 146/385 86459/235620 ]
1841  }
1842 
1843  inline double short_range_function(double q) const override { return qPochhammerSymbol(q, 0, order); }
1844  inline double short_range_function_derivative(double q) const override {
1845  return qPochhammerSymbolDerivative(q, 0, order);
1846  }
1847  inline double short_range_function_second_derivative(double q) const override {
1848  return qPochhammerSymbolSecondDerivative(q, 0, order);
1849  }
1850  inline double short_range_function_third_derivative(double q) const override {
1851  return qPochhammerSymbolThirdDerivative(q, 0, order);
1852  }
1853 
1854 #ifdef NLOHMANN_JSON_HPP
1855 
1856  inline qPotential(const nlohmann::json &j)
1857  : qPotential(j.at("cutoff").get<double>(), j.at("order").get<double>()) {}
1858 
1859  private:
1860  inline void _to_json(nlohmann::json &j) const override {
1861  j = {{ "order", order }};
1862  }
1863 #endif
1864 };
1865 
1898 class Poisson : public EnergyImplementation<Poisson> {
1899  private:
1900  signed int C, D;
1901  double kappaRed, kappaRed2;
1902  double yukawa_denom, binomCDC;
1903  bool yukawa;
1904 
1905  public:
1912  inline Poisson(double cutoff, signed int C, signed int D, double debye_length = infinity)
1913  : EnergyImplementation(Scheme::poisson, cutoff, debye_length), C(C), D(D) {
1914  if ( C < 1 )
1915  throw std::runtime_error("`C` must be larger than zero");
1916  if ( ( D < -1 ) && ( D != -C ) )
1917  throw std::runtime_error("If `D` is less than negative one, then it has to equal negative `C`");
1918  if ( ( D == 0 ) && ( C != 1 ) )
1919  throw std::runtime_error("If `D` is zero, then `C` has to equal one ");
1920  name = "poisson";
1921  dipolar_selfenergy = true;
1922  if(C < 2)
1923  dipolar_selfenergy = false;
1924  doi = "10/c5fr";
1925  double a1 = -double(C + D) / double(C);
1926  kappaRed = 0.0;
1927  yukawa = false;
1928  if( !std::isinf(debye_length) ) {
1929  kappaRed = cutoff / debye_length;
1930  if (std::fabs(kappaRed) > 1e-6) {
1931  yukawa = true;
1932  kappaRed2 = kappaRed * kappaRed;
1933  yukawa_denom = 1.0 / (1.0 - std::exp(2.0 * kappaRed));
1934  a1 *= -2.0 * kappaRed * yukawa_denom;
1935  }
1936  }
1937  binomCDC = 0.0;
1938  if( D != -C )
1939  binomCDC = double(binomial(C + D, C) * D);
1940  setSelfEnergyPrefactor({0.5 * a1, 0.0}); // Dipole self-energy seems to be 0 for C >= 2
1941  if (C == 2)
1942  setSelfEnergyPrefactor({0.5 * a1, -double(D) * (double(D * D) + 3.0 * double(D) + 2.0) / 12.0});
1943  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) +
1944  short_range_function(0.0); // Is this OK for Yukawa-interactions?
1945  chi = -2.0 * std::acos(-1.0) * cutoff * cutoff * (1.0 + double(C)) * (2.0 + double(C)) /
1946  (3.0 * double(D + 1 + C) *
1947  double(D + 2 + C)); // not confirmed, but have worked for all tested values of 'C' and 'D'
1948  }
1949 
1950  inline double short_range_function(double q) const override {
1951  if( D == -C )
1952  return 1.0;
1953  double tmp = 0;
1954  double qp = q;
1955  if (yukawa)
1956  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
1957  if( ( D == 0 ) && ( C == 1 ) )
1958  return ( 1.0 - qp );
1959  for (signed int c = 0; c < C; c++)
1960  tmp += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
1961  return powi(1.0 - qp, D + 1) * tmp;
1962  }
1963 
1964  inline double short_range_function_derivative(double q) const override {
1965  if( D == -C )
1966  return 0.0;
1967  if( ( D == 0 ) && ( C == 1 ) )
1968  return 0.0;
1969  double qp = q;
1970  double dqpdq = 1.0;
1971  if (yukawa) {
1972  double exp2kq = std::exp(2.0 * kappaRed * q);
1973  qp = (1.0 - exp2kq) * yukawa_denom;
1974  dqpdq = -2.0 * kappaRed * exp2kq * yukawa_denom;
1975  }
1976  double tmp1 = 1.0;
1977  double tmp2 = 0.0;
1978  for (int c = 1; c < C; c++) {
1979  double _fact = double(binomial(D - 1 + c, c)) * double(C - c) / double(C);
1980  tmp1 += _fact * powi(qp, c);
1981  tmp2 += _fact * double(c) * powi(qp, c - 1);
1982  }
1983  double dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
1984  return dSdqp * dqpdq;
1985  }
1986 
1987  inline double short_range_function_second_derivative(double q) const override {
1988  if( D == -C )
1989  return 0.0;
1990  if( ( D == 0 ) && ( C == 1 ) )
1991  return 0.0;
1992  double qp = q;
1993  double dqpdq = 1.0;
1994  double d2qpdq2 = 0.0;
1995  double dSdqp = 0.0;
1996  if (yukawa) {
1997  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
1998  dqpdq = -2.0 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
1999  d2qpdq2 = -4.0 * kappaRed2 * std::exp(2.0 * kappaRed * q) * yukawa_denom;
2000  double tmp1 = 1.0;
2001  double tmp2 = 0.0;
2002  for (signed int c = 1; c < C; c++) {
2003  tmp1 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
2004  tmp2 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * double(c) * powi(qp, c - 1);
2005  }
2006  dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
2007  }
2008  double d2Sdqp2 = binomCDC * powi(1.0 - qp, D - 1) * powi(qp, C - 1);
2009  return (d2Sdqp2 * dqpdq * dqpdq + dSdqp * d2qpdq2);
2010  };
2011 
2012  inline double short_range_function_third_derivative(double q) const override {
2013  if( D == -C )
2014  return 0.0;
2015  if( ( D == 0 ) && ( C == 1 ) )
2016  return 0.0;
2017  double qp = q;
2018  double dqpdq = 1.0;
2019  double d2qpdq2 = 0.0;
2020  double d3qpdq3 = 0.0;
2021  double d2Sdqp2 = 0.0;
2022  double dSdqp = 0.0;
2023  if (yukawa) {
2024  qp = (1.0 - std::exp(2.0 * kappaRed * q)) * yukawa_denom;
2025  dqpdq = -2.0 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
2026  d2qpdq2 = -4.0 * kappaRed2 * std::exp(2.0 * kappaRed * q) * yukawa_denom;
2027  d3qpdq3 = -8.0 * kappaRed2 * kappaRed * std::exp(2.0 * kappaRed * q) * yukawa_denom;
2028  d2Sdqp2 = binomCDC * powi(1.0 - qp, D - 1) * powi(qp, C - 1);
2029  double tmp1 = 1.0;
2030  double tmp2 = 0.0;
2031  for (signed int c = 1; c < C; c++) {
2032  tmp1 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * powi(qp, c);
2033  tmp2 += double(binomial(D - 1 + c, c)) * double(C - c) / double(C) * double(c) * powi(qp, c - 1);
2034  }
2035  dSdqp = (-double(D + 1) * powi(1.0 - qp, D) * tmp1 + powi(1.0 - qp, D + 1) * tmp2);
2036  }
2037  double d3Sdqp3 =
2038  binomCDC * powi(1.0 - qp, D - 2) * powi(qp, C - 2) * ((2.0 - double(C + D)) * qp + double(C) - 1.0);
2039  return (d3Sdqp3 * dqpdq * dqpdq * dqpdq + 3.0 * d2Sdqp2 * dqpdq * d2qpdq2 + dSdqp * d3qpdq3);
2040  };
2041 
2042 #ifdef NLOHMANN_JSON_HPP
2043 
2045  inline Poisson(const nlohmann::json &j)
2046  : Poisson(j.at("cutoff").get<double>(), j.at("C").get<int>(), j.at("D").get<int>(),
2047  j.value("debyelength", infinity)) {}
2048 
2049  private:
2050  inline void _to_json(nlohmann::json &j) const override {
2051  j = {{"C", C}, { "D", D }};
2052  }
2053 #endif
2054 };
2055 
2056 // -------------- Fanourgakis ---------------
2057 
2062 class Fanourgakis : public EnergyImplementation<Fanourgakis> {
2063  public:
2067  inline Fanourgakis(double cutoff) : EnergyImplementation(Scheme::fanourgakis, cutoff) {
2068  name = "fanourgakis";
2069  dipolar_selfenergy = true;
2070  doi = "10.1063/1.3216520";
2071  setSelfEnergyPrefactor({-0.875, 0.0});
2072  T0 = short_range_function_derivative(1.0) - short_range_function(1.0) + short_range_function(0.0);
2073  chi = -5.0 * std::acos(-1.0) * cutoff * cutoff / 18.0;
2074  }
2075 
2076  inline double short_range_function(double q) const override {
2077  double q2 = q * q;
2078  return powi(1.0 - q, 4) * (1.0 + 2.25 * q + 3.0 * q2 + 2.5 * q2 * q);
2079  }
2080  inline double short_range_function_derivative(double q) const override {
2081  return (-1.75 + 26.25 * powi(q, 4) - 42.0 * powi(q, 5) + 17.5 * powi(q, 6));
2082  }
2083  inline double short_range_function_second_derivative(double q) const override {
2084  return 105.0 * powi(q, 3) * powi(q - 1.0, 2);
2085  };
2086  inline double short_range_function_third_derivative(double q) const override {
2087  return 525.0 * powi(q, 2) * (q - 0.6) * (q - 1.0);
2088  };
2089 #ifdef NLOHMANN_JSON_HPP
2090 
2091  inline Fanourgakis(const nlohmann::json &j) : Fanourgakis(j.at("cutoff").get<double>()) {}
2092 
2093  private:
2094  inline void _to_json(nlohmann::json &) const override {}
2095 #endif
2096 };
2097 
2098 #ifdef NLOHMANN_JSON_HPP
2099 inline std::shared_ptr<SchemeBase> createScheme(const nlohmann::json &j) {
2100  const std::map<std::string, Scheme> m = {{"plain", Scheme::plain},
2101  {"qpotential", Scheme::qpotential},
2102  {"wolf", Scheme::wolf},
2103  {"poisson", Scheme::poisson},
2104  {"reactionfield", Scheme::reactionfield},
2105  {"spline", Scheme::spline},
2106  {"fanourgakis", Scheme::fanourgakis},
2107  {"fennell", Scheme::fennell},
2108  {"zahn", Scheme::zahn},
2109  {"zerodipole", Scheme::zerodipole},
2110  {"ewald", Scheme::ewald},
2111  {"ewaldt", Scheme::ewaldt}}; // map string keyword to scheme type
2112 
2113  std::string name = j.at("type").get<std::string>();
2114  auto it = m.find(name);
2115  if (it == m.end())
2116  throw std::runtime_error("unknown coulomb scheme " + name);
2117 
2118  std::shared_ptr<SchemeBase> scheme;
2119  switch (it->second) {
2120  case Scheme::plain:
2121  scheme = std::make_shared<Plain>(j);
2122  break;
2123  case Scheme::wolf:
2124  scheme = std::make_shared<Wolf>(j);
2125  break;
2126  case Scheme::zahn:
2127  scheme = std::make_shared<Zahn>(j);
2128  break;
2129  case Scheme::fennell:
2130  scheme = std::make_shared<Fennell>(j);
2131  break;
2132  case Scheme::zerodipole:
2133  scheme = std::make_shared<ZeroDipole>(j);
2134  break;
2135  case Scheme::fanourgakis:
2136  scheme = std::make_shared<Fanourgakis>(j);
2137  break;
2138  case Scheme::qpotential5:
2139  scheme = std::make_shared<qPotentialFixedOrder<5>>(j);
2140  break;
2141  case Scheme::qpotential:
2142  scheme = std::make_shared<qPotential>(j);
2143  break;
2144  case Scheme::ewald:
2145  scheme = std::make_shared<Ewald>(j);
2146  break;
2147  case Scheme::ewaldt:
2148  scheme = std::make_shared<EwaldT>(j);
2149  break;
2150  case Scheme::poisson:
2151  scheme = std::make_shared<Poisson>(j);
2152  break;
2153  case Scheme::reactionfield:
2154  scheme = std::make_shared<ReactionField>(j);
2155  break;
2156  default:
2157  break;
2158  }
2159  return scheme;
2160 }
2161 #endif
2162 
2163 // -------------- Splined ---------------
2164 
2181 class Splined : public EnergyImplementation<Splined> {
2182  private:
2183  std::shared_ptr<SchemeBase> pot;
2184  Tabulate::Andrea<double> splined_srf; // spline class
2185  std::array<Tabulate::TabulatorBase<double>::data, 4> splinedata; // 0=original, 1=first derivative, ...
2186 
2187  inline void generate_spline_data() {
2188  assert(pot);
2189  SchemeBase::operator=(*pot); // copy base data from pot -> Splined
2190  splinedata[0] = splined_srf.generate([pot = pot](double q) { return pot->short_range_function(q); }, 0, 1);
2191  splinedata[1] =
2192  splined_srf.generate([pot = pot](double q) { return pot->short_range_function_derivative(q); }, 0, 1);
2193  splinedata[2] = splined_srf.generate(
2194  [pot = pot](double q) { return pot->short_range_function_second_derivative(q); }, 0, 1);
2195  splinedata[3] =
2196  splined_srf.generate([pot = pot](double q) { return pot->short_range_function_third_derivative(q); }, 0, 1);
2197  }
2198 
2199  public:
2200  inline Splined() : EnergyImplementation<Splined>(Scheme::spline, infinity) {
2201  setTolerance(1e-3);
2202  }
2203 
2207  inline std::vector<size_t> numKnots() const {
2208  std::vector<size_t> n;
2209  for (auto &i : splinedata)
2210  n.push_back( i.numKnots() );
2211  return n;
2212  }
2213 
2217  inline void setTolerance(double tol) {
2218  splined_srf.setTolerance(tol);
2219  }
2220 
2227  template <class T, class... Args> void spline(Args &&... args) {
2228  pot = std::make_shared<T>(args...);
2229  generate_spline_data();
2230  }
2231  inline double short_range_function(double q) const override { return splined_srf.eval(splinedata[0], q); };
2232 
2233  inline double short_range_function_derivative(double q) const override {
2234  return splined_srf.eval(splinedata[1], q);
2235  }
2236  inline double short_range_function_second_derivative(double q) const override {
2237  return splined_srf.eval(splinedata[2], q);
2238  }
2239  inline double short_range_function_third_derivative(double q) const override {
2240  return splined_srf.eval(splinedata[3], q);
2241  }
2242 #ifdef NLOHMANN_JSON_HPP
2243  public:
2244  inline void to_json(nlohmann::json &j) const { pot->to_json(j); }
2245 
2246  private:
2247  inline void _to_json(nlohmann::json &) const override {}
2248 #endif
2249 };
2250 
2251 } // namespace CoulombGalore
constexpr double infinity
Numerical infinity.
Definition: coulombgalore.h:49
qPotential scheme
vec3 ion_field(double z, const vec3 &r) const override
electrostatic field from point charge
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)
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:85
Fanourgakis(const nlohmann::json &j)
ReactionField(const nlohmann::json &j)
Ewald real-space scheme using a truncated Gaussian screening-function.
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)
std::string name
Descriptive name.
Ewald real-space scheme using a Gaussian screening-function.
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
double quadrupole_potential(const mat33 &quad, const vec3 &r) const override
electrostatic potential from point dipole
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
vec3 ion_quadrupole_force(double z, const mat33 &quad, const vec3 &r) const override
interaction force between a point charge and a point quadrupole
double ion_quadrupole_energy(double z, const mat33 &quad, const vec3 &r) const override
interaction energy between a point charges and a point quadrupole
void spline(Args &&... args)
Spline given potential type.
vec3 dipole_field(const vec3 &mu, const vec3 &r) const override
electrostatic field from point dipole
Fennell(const nlohmann::json &j)
double qPochhammerSymbol(double q, int l=0, int P=300)
Help-function for the q-potential scheme.
Fanourgakis scheme.
double surface_energy(const std::vector< vec3 > &positions, const std::vector< double > &charges, const std::vector< vec3 > &dipoles, double volume) const
Surface-term.
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.
qPotential(double cutoff, int order)
double calc_dielectric(double M2V)
Calculate dielectric constant.
void setTolerance(double tol)
Set relative spline tolerance.
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.
double multipole_multipole_energy(double zA, double zB, const vec3 &muA, const vec3 &muB, const mat33 &quadA, const mat33 &quadB, const vec3 &r) const override
interaction energy between two multipoles with charges and dipole moments
Fennell scheme.
Eigen::Vector3d vec3
Definition: coulombgalore.h:46
double qPochhammerSymbolDerivative(double q, int l=0, int P=300)
Gives the derivative of the q-Pochhammer Symbol.
vec3 quadrupole_field(const mat33 &quad, const vec3 &r) const
electrostatic field from point quadrupole
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:73
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)
double self_energy(const std::array< double, 2 > &squared_moments) const
Self-energy for all type of interactions.
Ewald(double cutoff, double alpha, double eps_sur=infinity, double debye_length=infinity)
vec3 multipole_multipole_force(double zA, double zB, const vec3 &muA, const vec3 &muB, const mat33 &quadA, const mat33 &quadB, const vec3 &r) const override
interaction force between two point multipoles
EwaldT(double cutoff, double alpha, double eps_sur=infinity, double debye_length=infinity)
std::string doi
DOI for original citation.
vec3 multipole_field(double z, const vec3 &mu, const mat33 &quad, const vec3 &r) const
electrostatic field from point multipole
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)
std::vector< size_t > numKnots() const
Returns vector with number of spline knots the short-range-function and its derivatives.
Poisson(const nlohmann::json &j)
No truncation scheme, cutoff = infinity.
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.