QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

WIP 2D-in-3D Ewald sum for z != 0 

Not tested; error estimates not yet implemented.


Former-commit-id: d6886f64eb8b7e137abf6f187f8cd75f21a5f591
This commit is contained in:
Marek Nečada 2020-05-19 16:00:08 +03:00
parent 932f30e9f9
commit 85dbf79caa
3 changed files with 51 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
qpms/ewald.c
View File

@ -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

6
qpms/ewald.h
View File

@ -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

24
qpms/ewaldsf.c
View File

@ -12,6 +12,7 @@
#include <gsl/gsl_errno.h>
#include <float.h>
#include <stdbool.h>
#include <Faddeeva.h>
#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
}
}