diff --git a/qpms/cyewaldtest.pyx b/qpms/cyewaldtest.pyx index 5d286ff..4d545b1 100644 --- a/qpms/cyewaldtest.pyx +++ b/qpms/cyewaldtest.pyx @@ -3,12 +3,12 @@ from libc.stdlib cimport malloc, free, calloc import numpy as np cdef extern from "ewald.h": - void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, cdouble z) - void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, cdouble z) - void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, cdouble z) - int complex_gamma_inc_e(double a, cdouble x, int m, qpms_csf_result *result) + void ewald3_2_sigma_long_Delta(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) + void ewald3_2_sigma_long_Delta_series(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) + void ewald3_2_sigma_long_Delta_recurrent(cdouble *target, double *err, int maxn, cdouble x, int xbranch, cdouble z) + int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result) -def e32_Delta(int maxn, cdouble x, cdouble z, get_err=True, method='auto'): +def e32_Delta(int maxn, cdouble x, cdouble z, int xbranch = 0, get_err=True, method='auto'): cdef np.ndarray[double, ndim=1] err_np cdef double[::1] err_view cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((maxn+1,), dtype=complex, order='C') @@ -17,19 +17,19 @@ def e32_Delta(int maxn, cdouble x, cdouble z, get_err=True, method='auto'): err_np = np.empty((maxn+1,), order='C') err_view = err_np if method == 'recurrent': - ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + ewald3_2_sigma_long_Delta_recurrent(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z) elif method == 'series': - ewald3_2_sigma_long_Delta_series(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + ewald3_2_sigma_long_Delta_series(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z) else: - ewald3_2_sigma_long_Delta(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, z) + ewald3_2_sigma_long_Delta(&target_view[0], &err_view[0] if get_err else NULL, maxn, x, xbranch, z) if get_err: return target_np, err_np else: return target_np -def gamma_inc(double a, cdouble x, int m=0): +def gamma_inc(double a, cdouble x, int xbranch=0): cdef qpms_csf_result res - complex_gamma_inc_e(a, x, m, &res) + complex_gamma_inc_e(a, x, xbranch, &res) return res.val diff --git a/qpms/ewald.c b/qpms/ewald.c index c4e861b..cc36305 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -314,7 +314,7 @@ int ewald3_21_xy_sigma_long ( QPMS_UNTESTED; // TODO check the branches/phases! complex double z = I * k * gamma_pq * particle_shift.z; - ewald3_2_sigma_long_Delta(Gamma_pq, err ? Gamma_pq_err : NULL, lMax/2, x, z); + ewald3_2_sigma_long_Delta(Gamma_pq, err ? Gamma_pq_err : NULL, lMax/2, x, 0 /* FIXME */, z); } else { QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented"); } diff --git a/qpms/ewald.h b/qpms/ewald.h index 9485a2f..5c58024 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -217,13 +217,15 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result * unsuitable especially for big values of \a maxn. * */ -void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x, complex double z); +void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x, + int xbranch, complex double z); /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /** This function always uses Kambe's (corrected) recurrent formula. * For production, use ewald3_2_sigma_long_Delta() instead. */ -void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x, complex double z); +void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x, + int xbranch, complex double z); /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum. /** This function always uses Taylor expansion in \a z. @@ -233,7 +235,8 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_ * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth * of the error. */ -void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x, complex double z); +void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x, + int xbranch, 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 4889cae..8f50ee2 100644 --- a/qpms/ewaldsf.c +++ b/qpms/ewaldsf.c @@ -13,6 +13,7 @@ #include #include #include +#include "tiny_inlines.h" #ifndef COMPLEXPART_REL_ZERO_LIMIT #define COMPLEXPART_REL_ZERO_LIMIT 1e-14 @@ -290,11 +291,17 @@ int hyperg_2F2_series(const double a, const double b, const double c, const doub return GSL_SUCCESS; } +// Complex square root with branch selection +static inline complex double csqrt_branch(complex double x, int xbranch) { + return csqrt(x) * min1pow(xbranch); +} + // 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_recurrent(complex double *target, double *err, int maxn, complex double x, complex double z) { +void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *err, + int maxn, complex double x, int xbranch, 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 + complex double sqrtx = csqrt_branch(x, xbranch); // 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); @@ -305,7 +312,7 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *err, in target[1] = I / z * M_SQRTPI * expfac * (w_minus - w_plus); 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] = -(4 / (z*z)) * (-(0.5 - n) * target[n] + target[n-1] - cpow(x, 0.5 - n) * expfac); + target[n+1] = -(4 / (z*z)) * (-(0.5 - n) * target[n] + target[n-1] - sqrtx * cpow(x, -n) * expfac); } if (err) { // The error estimates for library math functions are based on @@ -334,7 +341,8 @@ void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *err, in } -void ewald3_2_sigma_long_Delta_series(complex double *target, double *err, int maxn, complex double x, complex double z) { +void ewald3_2_sigma_long_Delta_series(complex double *target, double *err, + int maxn, complex double x, int xbranch, complex double z) { complex double w = 0.25*z*z; double w_abs = cabs(w); int maxk; @@ -360,7 +368,7 @@ void ewald3_2_sigma_long_Delta_series(complex double *target, double *err, int m for(int j = 0; j <= maxn+maxk; ++j) { qpms_csf_result g; - QPMS_ENSURE_SUCCESS(complex_gamma_inc_e(0.5-j, x, 0 /* TODO branch choice */, &g)); + QPMS_ENSURE_SUCCESS(complex_gamma_inc_e(0.5-j, x, xbranch, &g)); Gammas[j] = g.val; if(err) { Gammas_abs[j] = cabs(g.val); @@ -393,11 +401,12 @@ void ewald3_2_sigma_long_Delta_series(complex double *target, double *err, int m } -void ewald3_2_sigma_long_Delta(complex double *target, double *err, int maxn, complex double x, complex double z) { +void ewald3_2_sigma_long_Delta(complex double *target, double *err, + int maxn, complex double x, int xbranch, complex double z) { double absz = cabs(z); if (absz < 2.) // TODO take into account also the other parameters - ewald3_2_sigma_long_Delta_series(target, err, maxn, x, z); + ewald3_2_sigma_long_Delta_series(target, err, maxn, x, xbranch, z); else - ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, z); + ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, xbranch, z); }