From 6ccd785164aefdd4dba9b62bd8c87091aa27e91f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 18 Nov 2020 18:08:57 +0200 Subject: [PATCH] Additional phase factors in Ewald sum for debugging --- qpms/cyewaldtest.pyx | 12 ++++++++++++ qpms/ewald.c | 14 ++++++++++++++ qpms/ewald.h | 5 +++++ 3 files changed, 31 insertions(+) diff --git a/qpms/cyewaldtest.pyx b/qpms/cyewaldtest.pyx index c88215b..1bb7105 100644 --- a/qpms/cyewaldtest.pyx +++ b/qpms/cyewaldtest.pyx @@ -7,6 +7,8 @@ cdef extern from "ewald.h": 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, bint bigimz) int complex_gamma_inc_e(double a, cdouble x, int xbranch, qpms_csf_result *result) + extern int ewald_factor_ipow_l + extern int ewald_factor_ipow_m 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 @@ -34,4 +36,14 @@ def gamma_inc(double a, cdouble x, int xbranch=0): complex_gamma_inc_e(a, x, xbranch, &res) return res.val +def ipow_l_factor(n = None): + global ewald_factor_ipow_l + if n is not None: + ewald_factor_ipow_l = n + return ewald_factor_ipow_l +def ipow_m_factor(n = None): + global ewald_factor_ipow_m + if n is not None: + ewald_factor_ipow_m = n + return ewald_factor_ipow_m diff --git a/qpms/ewald.c b/qpms/ewald.c index a6358b6..208b646 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -12,6 +12,17 @@ #include #include +/* + * Control additional phase factors + * that could possibly affect the resulting VSWF + * shape. + * Currently not used in ewald3_1_z_sigma_long() + */ + +int ewald_factor_ipow_l = 0; +int ewald_factor_ipow_m = 0; + + // parameters for the quadrature of integral in (4.6) #ifndef INTEGRATION_WORKSPACE_LIMIT #define INTEGRATION_WORKSPACE_LIMIT 30000 @@ -455,6 +466,7 @@ int ewald3_21_xy_sigma_long ( ckahanadd(&jsum, &jsum_c, summand); } jsum *= phasefac * factor1d; // PFC + jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m); ckahanadd(target + y, target_c + y, jsum); #ifdef EWALD_AUTO_CUTOFF kahanadd(&lsum, &lsum_c, cabs(jsum)); @@ -483,6 +495,7 @@ int ewald3_21_xy_sigma_long ( ckahanadd(&jsum, &jsum_c, jsummand); } jsum *= phasefac; // factor1d not here, off-axis sums not implemented/allowed. + jsum *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m); ckahanadd(target + y, target_c + y, jsum); #ifdef EWALD_AUTO_CUTOFF kahanadd(&lsum, &lsum_c, cabs(jsum)); @@ -903,6 +916,7 @@ int ewald3_sigma_short( kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown * interr[n])); // TODO include also other errors complex double thesummand = prefacn * R_pq_pown * leg * cintres[n] * e_beta_Rpq * e_imf * min1pow_m_neg(m); + thesummand *= ipow(n * ewald_factor_ipow_l) * ipow(m * ewald_factor_ipow_m); ckahanadd(target + y, target_c + y, thesummand); #ifdef EWALD_AUTO_CUTOFF kahanadd(&lsum, &lsum_c, cabs(thesummand)); diff --git a/qpms/ewald.h b/qpms/ewald.h index d7bad57..96121c8 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -331,4 +331,9 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep cart3_t particle_shift ); + +// If nonzero, adds an additional factor \f$ i^{nl} \f$ to the Ewald sum result (for debugging). +extern int ewald_factor_ipow_l; +// If nonzero, adds an additional factor \f$ i^{nm} \f$ to the Ewald sum result (for debubbing). +extern int ewald_factor_ipow_m; #endif //EWALD_H