ESPResSo
Extensible Simulation Package for Research on Soft Matter Systems
Loading...
Searching...
No Matches
mmm-modpsi.cpp
Go to the documentation of this file.
1/*
2 * Copyright (C) 2010-2022 The ESPResSo project
3 * Copyright (C) 2002,2003,2004,2005,2006,2007,2008,2009,2010
4 * Max-Planck-Institute for Polymer Research, Theory Group
5 *
6 * This file is part of ESPResSo.
7 *
8 * ESPResSo is free software: you can redistribute it and/or modify
9 * it under the terms of the GNU General Public License as published by
10 * the Free Software Foundation, either version 3 of the License, or
11 * (at your option) any later version.
12 *
13 * ESPResSo is distributed in the hope that it will be useful,
14 * but WITHOUT ANY WARRANTY; without even the implied warranty of
15 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 * GNU General Public License for more details.
17 *
18 * You should have received a copy of the GNU General Public License
19 * along with this program. If not, see <http://www.gnu.org/licenses/>.
20 */
21
22#include "config/config.hpp"
23
24#ifdef ELECTROSTATICS
25
26#include "mmm1d.hpp"
27#include "specfunc.hpp"
28
29#include <cmath>
30#include <numbers>
31#include <vector>
32
33static void preparePolygammaEven(int n, double binom,
34 std::vector<double> &series) {
35 /* (-0.5 n) psi^2n/2n! (-0.5 n) and psi^(2n+1)/(2n)! series expansions
36 note that BOTH carry 2n! */
37 auto const deriv = static_cast<double>(2 * n);
38 if (n == 0) {
39 // psi^0 has a slightly different series expansion
40 double maxx = 0.25;
41 series.resize(1);
42 series[0] = 2. * (1. - std::numbers::egamma);
43 for (int order = 1;; order += 1) {
44 auto const x_order = static_cast<double>(2 * order);
45 auto const coeff = -2 * hzeta(x_order + 1, 2);
46 if (fabs(maxx * coeff) * (4.0 / 3.0) < ROUND_ERROR_PREC)
47 break;
48 series.push_back(coeff);
49
50 maxx *= 0.25;
51 }
52 } else {
53 // even, n > 0
54 double maxx = 1;
55 double pref = 2;
56
57 for (int order = 0;; order++) {
58 // only even exponents of x
59 auto const x_order = static_cast<double>(2 * order);
60 auto const coeff = pref * hzeta(1 + deriv + x_order, 2);
61 if ((fabs(maxx * coeff) * (4.0 / 3.0) < ROUND_ERROR_PREC) &&
62 (x_order > deriv))
63 break;
64 series.push_back(-binom * coeff);
65
66 maxx *= 0.25;
67 pref *= (1.0 + deriv / (x_order + 1));
68 pref *= (1.0 + deriv / (x_order + 2));
69 }
70 }
71}
72
73static void preparePolygammaOdd(int n, double binom,
74 std::vector<double> &series) {
75 auto const deriv = static_cast<double>(2 * n + 1);
76 auto maxx = 0.5;
77 // to get 1/(2n)! instead of 1/(2n+1)!
78 auto pref = 2 * deriv * (1 + deriv);
79
80 for (int order = 0;; order++) {
81 // only odd exponents of x
82 auto const x_order = static_cast<double>(2 * order + 1);
83 auto const coeff = pref * hzeta(1 + deriv + x_order, 2);
84 if ((fabs(maxx * coeff) * (4.0 / 3.0) < ROUND_ERROR_PREC) &&
85 (x_order > deriv))
86 break;
87
88 series.push_back(-binom * coeff);
89 maxx *= 0.25;
90 pref *= (1.0 + deriv / (x_order + 1));
91 pref *= (1.0 + deriv / (x_order + 2));
92 }
93}
94
95void CoulombMMM1D::create_mod_psi_up_to(int new_n) {
96 auto const old_n = static_cast<int>(modPsi.size() >> 1);
97 if (new_n > old_n) {
98 modPsi.resize(2 * new_n);
99
100 double binom = 1.0;
101 for (int n = 0; n < old_n; n++)
102 binom *= (-0.5 - n) / static_cast<double>(n + 1);
103
104 for (int n = old_n; n < new_n; n++) {
105 preparePolygammaEven(n, binom, modPsi[2 * n]);
106 preparePolygammaOdd(n, binom, modPsi[2 * n + 1]);
107 binom *= (-0.5 - n) / static_cast<double>(n + 1);
108 }
109 }
110}
111
112#endif // ELECTROSTATICS
This file contains the defaults for ESPResSo.
#define ROUND_ERROR_PREC
Precision for capture of round off errors.
Definition config.hpp:66
static void preparePolygammaOdd(int n, double binom, std::vector< double > &series)
static void preparePolygammaEven(int n, double binom, std::vector< double > &series)
MMM1D algorithm for long-range Coulomb interactions on the CPU.
double hzeta(double s, double q)
Hurwitz zeta function.
Definition specfunc.cpp:215
This file contains implementations for some special functions which are needed by the MMM family of a...