Allow to select only either of long/short range Ewald sum parts

(cherry picked from commit 97b7782291380193835c23a4f0ea04ff5f44273e [formerly 41b1a76d5ed571b507abc0515b3cba60e8c0ccca])


Former-commit-id: 4d51186e653babfdc84fc45a4ab70a3ffdc02eee

qpms/ewald.h

@@ -35,6 +35,15 @@
 #include "qpms_types.h"
 #include "lattices.h"
 
+
+typedef enum {
+	QPMS_EWALD_LONG_RANGE = 1,
+	QPMS_EWALD_SHORT_RANGE = 2,
+	QPMS_EWALD_0TERM = 4,
+	QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
+} qpms_ewald_part;
+
+
 /// Use this handler to ignore underflows of incomplete gamma.
 gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
 

qpms/qpms_cdefs.pxd

@@ -562,6 +562,12 @@ cdef extern from "ewald.h":
         cdouble val
         double err
 
+    ctypedef enum qpms_ewald_part:
+        QPMS_EWALD_LONG_RANGE
+        QPMS_EWALD_SHORT_RANGE
+        QPMS_EWALD_FULL
+        QPMS_EWALD_0TERM
+
     struct qpms_ewald3_constants_t:
         qpms_l_t lMax
         qpms_y_t nelem_sc

qpms/translations.c

@@ -729,7 +729,7 @@ qpms_errno_t qpms_trans_calculator_get_trans_array(const qpms_trans_calculator *
   return retval;
 }
 
-qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
+qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
     complex double *target, double *err,
     /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
     const qpms_vswf_set_spec_t *destspec, size_t deststride,
@@ -739,7 +739,8 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
     cart2_t b1, cart2_t b2,
     const cart2_t beta,
     const cart2_t particle_shift,
-    double maxR, double maxK
+    double maxR, double maxK,
+    const qpms_ewald_part parts
     )
 {
   TROPS_ONLY_AB_SYMMETRIC_NORMS_IMPLEMENTED(c->normalisation);
@@ -757,9 +758,9 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
     QPMS_CRASHING_MALLOC(Aerr, c->nelem*c->nelem*sizeof(double));
     QPMS_CRASHING_MALLOC(Berr, c->nelem*c->nelem*sizeof(double));
   }
-  qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays_e32(c,
+  qpms_errno_t retval = qpms_trans_calculator_get_AB_arrays_e32_e(c,
       A, Aerr, B, Berr, ldAB, 1, 
-      eta, k, b1, b2, beta, particle_shift, maxR, maxK);
+      eta, k, b1, b2, beta, particle_shift, maxR, maxK, parts);
   for (size_t desti = 0; desti < destspec->n; ++desti) {
     qpms_y_t desty; qpms_vswf_type_t destt;
     if(QPMS_SUCCESS != qpms_uvswfi2ty(destspec->ilist[desti], &destt, &desty))
@@ -781,6 +782,22 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
   return retval;
 }
 
+qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
+    complex double *target, double *err,
+    /// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm.
+    const qpms_vswf_set_spec_t *destspec, size_t deststride,
+    /// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm.
+    const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
+    const double eta, const complex double k,
+    cart2_t b1, cart2_t b2,
+    const cart2_t beta,
+    const cart2_t particle_shift,
+    double maxR, double maxK
+    )
+{
+  return qpms_trans_calculator_get_trans_array_e32_e(c, target, err, destspec, deststride,
+      srcspec, srcstride, eta, k, b1, b2, beta, particle_shift, maxR, maxK, QPMS_EWALD_FULL);
+}
 
 
 qpms_errno_t qpms_trans_calculator_get_trans_array_lc3p( 
@@ -913,7 +930,7 @@ int qpms_trans_calculator_get_AB_arrays_e31z_both_points_and_shift(const qpms_tr
 // and GEN_KSHIFTEDPOINTS.
 // The results should be the same. The performance can slightly differ (especially
 // if some optimizations in the point generators are implemented.)
-int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
+int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
     complex double * const Adest, double * const Aerr,
     complex double * const Bdest, double * const Berr,
     const ptrdiff_t deststride, const ptrdiff_t srcstride,
@@ -922,7 +939,8 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
     const cart2_t b1, const cart2_t b2,
     const cart2_t beta,
     const cart2_t particle_shift,
-    double maxR, double maxK
+    double maxR, double maxK,
+    const qpms_ewald_part parts
     )
 {
 
@@ -945,50 +963,58 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
   cart2_t rb1, rb2; // reciprocal basis
   QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2));
 
-  PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, 
+  if (parts & QPMS_EWALD_LONG_RANGE) {
-#ifdef GEN_RSHIFTEDPOINTS
+    PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL,
-      cart2_scale(-1 /*CHECKSIGN*/, particle_shift),
-#else
-      CART2_ZERO,
-#endif
-      0, !do_sigma0, maxR, false);
-  PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL,
 #ifdef GEN_KSHIFTEDPOINTS
-      beta,
+        beta,
 #else
-     CART2_ZERO,
+        CART2_ZERO,
 #endif
-     0, true, maxK, false);
+        0, true, maxK, false);
 
-  QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k, 
+    QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k, 
-      unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
+          unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen,
 #ifdef GEN_KSHIFTEDPOINTS
-      true,
+          true,
 #else
-      false,
+          false,
 #endif
-      cart22cart3xy(beta), cart22cart3xy(particle_shift)));
+          cart22cart3xy(beta), cart22cart3xy(particle_shift)));
-  if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
+    if(Kgen.stateData) // PGen not consumed entirely (converged earlier)
-    PGen_destroy(&Kgen);
+      PGen_destroy(&Kgen);
+  }
 
-  QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
+  if (parts & QPMS_EWALD_SHORT_RANGE) {
-      LAT_2D_IN_3D_XYONLY, &Rgen, 
+    PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, 
 #ifdef GEN_RSHIFTEDPOINTS
-      true,
+        cart2_scale(-1 /*CHECKSIGN*/, particle_shift),
 #else
-      false,
+        CART2_ZERO,
 #endif
-      cart22cart3xy(beta), cart22cart3xy(particle_shift)));
+        0, !do_sigma0, maxR, false);
-  if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
+
-    PGen_destroy(&Rgen);
+    QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k,
+          LAT_2D_IN_3D_XYONLY, &Rgen, 
+#ifdef GEN_RSHIFTEDPOINTS
+          true,
+#else
+          false,
+#endif
+          cart22cart3xy(beta), cart22cart3xy(particle_shift)));
+
+    if(Rgen.stateData) // PGen not consumed entirely (converged earlier)
+      PGen_destroy(&Rgen);
+  }
 
   for(qpms_y_t y = 0; y < nelem2_sc; ++y)
-    sigmas_total[y] = sigmas_short[y] + sigmas_long[y];
+    sigmas_total[y] = ((parts & QPMS_EWALD_SHORT_RANGE) ? sigmas_short[y] : 0)
+                    + ((parts & QPMS_EWALD_LONG_RANGE)  ? sigmas_long[y] : 0);
   if (doerr) for(qpms_y_t y = 0; y < nelem2_sc; ++y)
-    serr_total[y]  = serr_short[y] + serr_long[y];
+    serr_total[y]  = ((parts & QPMS_EWALD_SHORT_RANGE) ? serr_short[y] : 0)
+                    + ((parts & QPMS_EWALD_LONG_RANGE)  ? serr_long[y] : 0);
 
   complex double sigma0 = 0; double sigma0_err = 0;
-  if (do_sigma0) {
+  if (do_sigma0 && (parts & QPMS_EWALD_0TERM)) {
     QPMS_ENSURE_SUCCESS(ewald3_sigma0(&sigma0, &sigma0_err, c->e3c, eta, k));
     const qpms_l_t y = qpms_mn2y_sc(0,0);
     sigmas_total[y] += sigma0;
@@ -1051,6 +1077,22 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
   return 0;
 }
 
+int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
+    complex double * const Adest, double * const Aerr,
+    complex double * const Bdest, double * const Berr,
+    const ptrdiff_t deststride, const ptrdiff_t srcstride,
+    /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS
+    const double eta, const complex double k,
+    const cart2_t b1, const cart2_t b2,
+    const cart2_t beta,
+    const cart2_t particle_shift,
+    double maxR, double maxK) 
+{
+  return qpms_trans_calculator_get_AB_arrays_e32_e(
+      c, Adest, Aerr, Bdest, Berr, deststride, srcstride, 
+      eta, k, b1, b2, beta, particle_shift, maxR, maxK, QPMS_EWALD_FULL);
+}
+
 #endif // LATTICESUMS32
 
 

qpms/translations.h

@@ -166,6 +166,18 @@ int qpms_trans_calculator_get_AB_arrays_e32(const qpms_trans_calculator *c,
 		double maxR, double maxK
 		);
 
+int qpms_trans_calculator_get_AB_arrays_e32_e(const qpms_trans_calculator *c,
+		complex double *Adest, double *Aerr,
+		complex double *Bdest, double *Berr,
+		const ptrdiff_t deststride, const ptrdiff_t srcstride,
+		const double eta, const complex double k,
+		cart2_t b1, cart2_t b2,
+		const cart2_t beta,
+		const cart2_t particle_shift,
+		double maxR, double maxK,
+		qpms_ewald_part parts
+		);
+
 // Convenience functions using VSWF base specs
 qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculator *c,
 		complex double *target, double *err,
@@ -180,6 +192,20 @@ qpms_errno_t qpms_trans_calculator_get_trans_array_e32(const qpms_trans_calculat
 		double maxR, double maxK
 		);
 
+qpms_errno_t qpms_trans_calculator_get_trans_array_e32_e(const qpms_trans_calculator *c,
+		complex double *target, double *err,
+		/// Must be destspec->lMax <= c-> lMax && destspec->norm == c->norm. 
+		const qpms_vswf_set_spec_t *destspec, size_t deststride,
+		/// Must be srcspec->lMax <= c-> lMax && srcspec->norm == c->norm. 
+		const qpms_vswf_set_spec_t *srcspec, size_t srcstride,
+		const double eta, const complex double k,
+		cart2_t b1, cart2_t b2,
+		const cart2_t beta,
+		const cart2_t particle_shift,
+		double maxR, double maxK,
+		qpms_ewald_part parts
+		);
+
 #endif //LATTICESUMS32
 
 #ifdef LATTICESUMS31