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

Additional phase factors in Ewald sum for debugging

This commit is contained in:
Marek Nečada 2020-11-18 18:08:57 +02:00
parent e486e04e2f
commit 6ccd785164
3 changed files with 31 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
qpms/cyewaldtest.pyx
View File

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

14
qpms/ewald.c
View File

@ -12,6 +12,17 @@
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_expint.h>
/*
* 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));

5
qpms/ewald.h
View File

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