From 85dbf79caa98a512592a290e50dde3fe0daa6ba9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 19 May 2020 16:00:08 +0300 Subject: [PATCH] WIP 2D-in-3D Ewald sum for z != 0 Not tested; error estimates not yet implemented. Former-commit-id: d6886f64eb8b7e137abf6f187f8cd75f21a5f591 --- qpms/ewald.c | 30 +++++++++++++++++++++--------- qpms/ewald.h | 6 ++++++ qpms/ewaldsf.c | 24 ++++++++++++++++++++++++ 3 files changed, 51 insertions(+), 9 deletions(-) diff --git a/qpms/ewald.c b/qpms/ewald.c index c992c50..7cf52ca 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -257,7 +257,8 @@ int ewald3_21_xy_sigma_long ( #endif // recycleable values if rbeta_pq stays the same: complex double gamma_pq; - complex double z; + complex double z; // I?*κ*γ*particle_shift.z + complex double x; // κ*γ/(2*η) complex double factor1d = 1; // the "additional" factor for the 1D case (then it is not 1) // space for Gamma_pq[j]'s qpms_csf_result Gamma_pq[lMax/2+1]; @@ -291,17 +292,28 @@ int ewald3_21_xy_sigma_long ( // R-DEPENDENT BEGIN + //void ewald3_2_sigma_long_Delta(complex double *target, int maxn, complex double x, complex double z) { if (new_rbeta_pq) { gamma_pq = clilgamma(rbeta_pq/k); - z = csq(gamma_pq*k/(2*eta)); - for(qpms_l_t j = 0; j <= lMax/2; ++j) { - int retval = complex_gamma_inc_e(0.5-j, z, - // we take the branch which is principal for the Re z < 0, Im z < 0 quadrant, cf. [Linton, p. 642 in the middle] - QPMS_LIKELY(creal(z) < 0) && !signbit(cimag(z)) ? -1 : 0, Gamma_pq+j); - QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW); + complex double x = gamma_pq*k/(2*eta); + complex double x2 = x*x; + if(particle_shift.z == 0) { + for(qpms_l_t j = 0; j <= lMax/2; ++j) { + int retval = complex_gamma_inc_e(0.5-j, x2, + // we take the branch which is principal for the Re x < 0, Im x < 0 quadrant, cf. [Linton, p. 642 in the middle] + QPMS_LIKELY(creal(x2) < 0) && !signbit(cimag(x2)) ? -1 : 0, Gamma_pq+j); + QPMS_ENSURE_SUCCESS_OR(retval, GSL_EUNDRFLW); + } + if (latdim & LAT1D) + factor1d = M_SQRT1_2 * .5 * k * gamma_pq; + } else if (latdim & LAT2D) { + QPMS_UNTESTED; + // TODO check the branches/phases! + complex double z = I * k * gamma_pq * particle_shift.z; + ewald3_2_sigma_long_Delta(Gamma_pq, lMax/2, x, z); + } else { + QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented"); } - if (latdim & LAT1D) - factor1d = M_SQRT1_2 * .5 * k * gamma_pq; } // R-DEPENDENT END diff --git a/qpms/ewald.h b/qpms/ewald.h index adde9d7..722e838 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -210,6 +210,12 @@ int hyperg_2F2_series(double a, double b, double c, double d, double x, int ewald32_sr_integral(double r, double k, double n, double eta, double *result, double *err, gsl_integration_workspace *workspace); #endif +/// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. +/** \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f] + * + * NOTE: The error part is not yet implemented, NANs are returned instead! + */ +void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z); // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space diff --git a/qpms/ewaldsf.c b/qpms/ewaldsf.c index ddf0bc1..69bce2e 100644 --- a/qpms/ewaldsf.c +++ b/qpms/ewaldsf.c @@ -12,6 +12,7 @@ #include #include #include +#include #ifndef COMPLEXPART_REL_ZERO_LIMIT #define COMPLEXPART_REL_ZERO_LIMIT 1e-14 @@ -289,5 +290,28 @@ int hyperg_2F2_series(const double a, const double b, const double c, const doub return GSL_SUCCESS; } +// The Delta_n factor from [Kambe II], Appendix 3 +// \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f] +void ewald3_2_sigma_long_Delta(qpms_csf_result *target, int maxn, complex double x, complex double z) { + complex double expfac = cexp(-x + 0.25 * z*z / x); + complex double sqrtx = csqrt(x); // TODO check carefully, which branch is needed + // These are used in the first two recurrences + complex double w_plus = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0); + complex double w_minus = Faddeeva_w(-z/(2*sqrtx) + I*sqrtx, 0); + QPMS_ASSERT(maxn >= 0); + if (maxn >= 0) { + target[0].val = 0.5 * M_SQRTPI * expfac * (w_minus + w_plus); + target[0].err = NAN; //TODO + } + if (maxn >= 1) { + target[1].val = -I / z * M_SQRTPI * expfac * (w_minus - w_plus); + target[1].err = NAN; //TODO + } + for(int n = 1; n < maxn; ++n) { // The rest via recurrence + // TODO The cpow(x, 0.5 - n) might perhaps better be replaced with a recurrently computed variant + target[n+1].val = (4 / (z*z)) * ((0.5 - n) * target[n].val + target[n-1].val - cpow(x, 0.5 - n) * expfac); + target[n+1].err = NAN; //TODO + } +}