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

Branch selection for Δ_n in Ewald sum 

Former-commit-id: 6c3d6e6aa9010bb66975f65397b8a061fac1b5ef
This commit is contained in:
Marek Nečada 2020-05-28 13:13:08 +03:00
parent c50f40a747
commit 975d23b557
4 changed files with 34 additions and 22 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
qpms/cyewaldtest.pyx
View File

@ -3,12 +3,12 @@ from libc.stdlib cimport malloc, free, calloc
import numpy as np import numpy as np
cdef extern from "ewald.h": 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(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, 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, 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 m, qpms_csf_result *result) 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 np.ndarray[double, ndim=1] err_np
cdef double[::1] err_view cdef double[::1] err_view
cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((maxn+1,), dtype=complex, order='C') 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_np = np.empty((maxn+1,), order='C')
err_view = err_np err_view = err_np
if method == 'recurrent': 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': 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: 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: if get_err:
return target_np, err_np return target_np, err_np
else: else:
return target_np 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 cdef qpms_csf_result res
complex_gamma_inc_e(a, x, m, &res) complex_gamma_inc_e(a, x, xbranch, &res)
return res.val return res.val

2
qpms/ewald.c
View File

@ -314,7 +314,7 @@ int ewald3_21_xy_sigma_long (
QPMS_UNTESTED; QPMS_UNTESTED;
// TODO check the branches/phases! // TODO check the branches/phases!
complex double z = I * k * gamma_pq * particle_shift.z; 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 { } else {
QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented"); QPMS_NOT_IMPLEMENTED("1D lattices in 3D space outside of the line not implemented");
} }

9
qpms/ewald.h
View File

@ -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. * 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. /// 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. /** This function always uses Kambe's (corrected) recurrent formula.
* For production, use ewald3_2_sigma_long_Delta() instead. * 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. /// 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. /** 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 * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
* of the error. * 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 // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space

25
qpms/ewaldsf.c
View File

@ -13,6 +13,7 @@
#include <float.h> #include <float.h>
#include <stdbool.h> #include <stdbool.h>
#include <Faddeeva.h> #include <Faddeeva.h>
#include "tiny_inlines.h"
#ifndef COMPLEXPART_REL_ZERO_LIMIT #ifndef COMPLEXPART_REL_ZERO_LIMIT
#define COMPLEXPART_REL_ZERO_LIMIT 1e-14 #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; 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 // 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] // \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 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 // These are used in the first two recurrences
complex double w_plus = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0); complex double w_plus = Faddeeva_w(+z/(2*sqrtx) + I*sqrtx, 0);
complex double w_minus = 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); target[1] = I / z * M_SQRTPI * expfac * (w_minus - w_plus);
for(int n = 1; n < maxn; ++n) { // The rest via recurrence 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 // 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) { if (err) {
// The error estimates for library math functions are based on // 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; complex double w = 0.25*z*z;
double w_abs = cabs(w); double w_abs = cabs(w);
int maxk; 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) { for(int j = 0; j <= maxn+maxk; ++j) {
qpms_csf_result g; 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; Gammas[j] = g.val;
if(err) { if(err) {
Gammas_abs[j] = cabs(g.val); 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); double absz = cabs(z);
if (absz < 2.) // TODO take into account also the other parameters 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 else
ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, z); ewald3_2_sigma_long_Delta_recurrent(target, err, maxn, x, xbranch, z);
} }