Enable separation of long- and short-range of Ewald sum in python

qpms/cycommon.pyx

@@ -20,6 +20,12 @@ class BesselType(enum.IntEnum):
     HANKEL_PLUS = QPMS_HANKEL_PLUS
     HANKEL_MINUS = QPMS_HANKEL_MINUS
 
+class EwaldPart(enum.IntEnum):
+    LONG_RANGE = QPMS_EWALD_LONG_RANGE
+    SHORT_RANGE = QPMS_EWALD_SHORT_RANGE
+    FULL = QPMS_EWALD_FULL
+    ZEROTERM = QPMS_EWALD_0TERM
+
 class PointGroupClass(enum.IntEnum):
     CN    = QPMS_PGS_CN
     S2N   = QPMS_PGS_S2N

qpms/ewald.h

@@ -36,14 +36,6 @@
 #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_c.pyx

@@ -11,7 +11,7 @@ from .qpms_cdefs cimport *
 from .cyquaternions cimport IRot3, CQuat
 from .cybspec cimport BaseSpec
 from .cycommon cimport make_c_string
-from .cycommon import string_c2py, PointGroupClass, BesselType
+from .cycommon import string_c2py, PointGroupClass, BesselType, EwaldPart
 from .cytmatrices cimport CTMatrix, TMatrixFunction, TMatrixGenerator, TMatrixInterpolator
 from .cymaterials cimport EpsMuGenerator, EpsMu
 from libc.stdlib cimport malloc, free, calloc
@@ -1023,7 +1023,8 @@ cdef class _ScatteringSystemAtOmegaK:
             self.sswk.eta = eta
     
     @boundscheck(False)
-    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS):
+    def scattered_E(self, scatcoeffvector_full, evalpos, btyp=QPMS_HANKEL_PLUS,
+            ewaldparts = EwaldPart.FULL):
         """Evaluate electric field for a given excitation coefficient vector (periodic system)
 
         Parameters
@@ -1042,6 +1043,7 @@ cdef class _ScatteringSystemAtOmegaK:
         if(btyp != QPMS_HANKEL_PLUS):
             raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
+        cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
         evalpos = np.array(evalpos, dtype=float, copy=False)
         if evalpos.shape[-1] != 3:
             raise ValueError("Last dimension of evalpos has to be 3")
@@ -1057,14 +1059,14 @@ cdef class _ScatteringSystemAtOmegaK:
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                res = qpms_scatsyswk_scattered_E(&self.sswk, btyp_c, &scv_view[0], pos)
+                res = qpms_scatsyswk_scattered_E_e(&self.sswk, btyp_c, &scv_view[0], pos, ewaldparts_c)
                 results[i,0] = res.x
                 results[i,1] = res.y
                 results[i,2] = res.z
         return results.reshape(evalpos.shape)
 
     @boundscheck(False)
-    def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS):
+    def scattered_field_basis(self, evalpos, btyp=QPMS_HANKEL_PLUS, ewaldparts=EwaldPart.FULL):
         # TODO examples
         """Evaluate scattered field "basis" (periodic system)
 
@@ -1086,6 +1088,7 @@ cdef class _ScatteringSystemAtOmegaK:
         if(btyp != QPMS_HANKEL_PLUS):
             raise NotImplementedError("Only first kind Bessel function-based fields are supported")
         cdef qpms_bessel_t btyp_c = BesselType(btyp)
+        cdef qpms_ewald_part ewaldparts_c = EwaldPart(ewaldparts)
         cdef Py_ssize_t fecv_size = self.fecv_size
         evalpos = np.array(evalpos, dtype=float, copy=False)
         if evalpos.shape[-1] != 3:
@@ -1101,7 +1104,7 @@ cdef class _ScatteringSystemAtOmegaK:
                 pos.x = evalpos_a[i,0]
                 pos.y = evalpos_a[i,1]
                 pos.z = evalpos_a[i,2]
-                qpms_scatsyswk_scattered_field_basis(res, &self.sswk, btyp_c, pos)
+                qpms_scatsyswk_scattered_field_basis_e(res, &self.sswk, btyp_c, pos, ewaldparts_c)
                 for j in range(fecv_size):
                     results[i,j,0] = res[j].x
                     results[i,j,1] = res[j].y

qpms/qpms_cdefs.pxd

@@ -147,6 +147,11 @@ cdef extern from "qpms_types.h":
         cdouble mu
     ctypedef enum qpms_coord_system_t:
         pass
+    ctypedef enum qpms_ewald_part:
+        QPMS_EWALD_LONG_RANGE
+        QPMS_EWALD_SHORT_RANGE
+        QPMS_EWALD_FULL
+        QPMS_EWALD_0TERM
     # maybe more if needed
 
 cdef extern from "qpms_error.h":
@@ -705,6 +710,8 @@ cdef extern from "scatsystem.h":
             const cdouble *f_excitation_vector_full, cart3_t where) nogil
     ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
             const cdouble *f_excitation_vector_full, cart3_t where) nogil
+    ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk, qpms_bessel_t btyp,
+            const cdouble *f_excitation_vector_full, cart3_t where, qpms_ewald_part parts) nogil
     qpms_errno_t qpms_scatsys_scattered_field_basis(ccart3_t *target, const qpms_scatsys_t *ss,
             qpms_bessel_t btyp, cdouble wavenumber, cart3_t where) nogil
     qpms_errno_t qpms_scatsysw_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_t *ssw,
@@ -715,6 +722,8 @@ cdef extern from "scatsystem.h":
             qpms_ss_pi_t pi, qpms_bessel_t btyp, cart3_t where) nogil
     qpms_errno_t qpms_scatsyswk_scattered_field_basis(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
             qpms_bessel_t btyp, cart3_t where) nogil
+    qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
+            qpms_bessel_t btyp, cart3_t where, qpms_ewald_part parts) nogil
     double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector) nogil
 
 
@@ -723,12 +732,6 @@ 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/qpms_types.h

@@ -426,5 +426,12 @@ typedef struct qpms_epsmu_t {
 
 struct qpms_tolerance_spec_t; // See tolerances.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;
+
 #define lmcheck(l,m) assert((l) >= 1 && abs(m) <= (l))
 #endif // QPMS_TYPES

qpms/scatsystem.c

@@ -2173,10 +2173,11 @@ ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw,
 
 // For periodic lattices, we use directly the "alternative" implementation, 
 // using translation operator and regular dipole waves at zero
-ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
+ccart3_t qpms_scatsyswk_scattered_E_e(const qpms_scatsys_at_omega_k_t *sswk,
     qpms_bessel_t btyp,
     const complex double *cvf, 
-    const cart3_t where
+    const cart3_t where,
+    const qpms_ewald_part parts
     ) {
   QPMS_UNTESTED;
   if (btyp != QPMS_HANKEL_PLUS) 
@@ -2208,13 +2209,13 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
     const double maxR = sqrt(ss->per.unitcell_volume) * 64;
     const double maxK = 2048 * 2 * M_PI / maxR;
 
-    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
+    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
           ss->c, s, NULL, 
           &dipspec, 1, bspec, dipspec.n,
           sswk->eta, sswk->ssw->wavenumber,
           cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
           cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, 
-          maxR, maxK));
+          maxR, maxK, parts));
 
     for(size_t i = 0; i < bspec->n; ++i)
       for(size_t j = 0; j < dipspec.n; ++j){
@@ -2228,11 +2229,17 @@ ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
   return res;
 }
 
-qpms_errno_t qpms_scatsyswk_scattered_field_basis(
+ccart3_t qpms_scatsyswk_scattered_E(const qpms_scatsys_at_omega_k_t *sswk,
+    qpms_bessel_t btyp, const complex double *cvf, const cart3_t where) {
+  return qpms_scatsyswk_scattered_E_e(sswk, btyp, cvf, where, QPMS_EWALD_FULL);
+}
+
+qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
     ccart3_t *target,
     const qpms_scatsys_at_omega_k_t *sswk,
     const qpms_bessel_t btyp,
-    const cart3_t where
+    const cart3_t where,
+    const qpms_ewald_part parts
     ) {
   QPMS_UNTESTED;
   if (btyp != QPMS_HANKEL_PLUS) 
@@ -2266,13 +2273,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
     const double maxR = sqrt(ss->per.unitcell_volume) * 64;
     const double maxK = 2048 * 2 * M_PI / maxR;
 
-    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32(
+    QPMS_ENSURE_SUCCESS(qpms_trans_calculator_get_trans_array_e32_e(
           ss->c, s, NULL, 
           &dipspec, 1, bspec, dipspec.n,
           sswk->eta, sswk->ssw->wavenumber,
           cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]),
           cart2_from_double_array(sswk->k), cart3_substract(where, particle_pos) /*CHECKSIGN*/, 
-          maxR, maxK));
+          maxR, maxK, parts));
 
     for(size_t i = 0; i < bspec->n; ++i)
       for(size_t j = 0; j < dipspec.n; ++j){
@@ -2288,6 +2295,13 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
   return QPMS_SUCCESS;
 }
 
+qpms_errno_t qpms_scatsyswk_scattered_field_basis(
+    ccart3_t *target, const qpms_scatsys_at_omega_k_t *sswk,
+    const qpms_bessel_t btyp, const cart3_t where) {
+  return qpms_scatsyswk_scattered_field_basis_e(
+      target, sswk, btyp, where, QPMS_EWALD_FULL);
+}
+
 #if 0
 ccart3_t qpms_scatsys_scattered_E_irrep(const qpms_scatsys_t *ss,
    qpms_iri_t iri, const complex double *cvr, cart3_t where) {

qpms/scatsystem.h

@@ -892,6 +892,14 @@ ccart3_t qpms_scatsyswk_scattered_E(
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated.
 		);
 
+ccart3_t qpms_scatsyswk_scattered_E_e(
+		const qpms_scatsys_at_omega_k_t *sswk,
+		qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supported).
+		const complex double *scatcoeff_full, ///< Full vector of the scattered field coefficients \f$ \wckcout \f$.
+		cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the field is evaluated.
+		qpms_ewald_part parts
+		);
+
 /// Evaluates "scattered" field basis functions in a periodic system.
 /**
  *
@@ -904,6 +912,14 @@ qpms_errno_t qpms_scatsyswk_scattered_field_basis(
 		cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the basis is evaluated.
 		);
 
+qpms_errno_t qpms_scatsyswk_scattered_field_basis_e(
+		ccart3_t *target, ///< Target array of size sswk->ssw->ss->fecv_size
+		const qpms_scatsys_at_omega_k_t *sswk,
+		qpms_bessel_t typ, ///< Bessel function kind to use (only QPMS_HANKEL_PLUS is currently supponted).
+		cart3_t evalpoint, ///< A point \f$ \vect r \f$, at which the basis is evaluated.
+		qpms_ewald_part parts
+		);
+
 
 /// Adjusted Ewadl parameter to avoid high-frequency breakdown.
 // TODO DOC