Additional phase factors in Ewald sum for debugging

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

qpms/ewald.c

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

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