From 8abe7e335725be1a3d78da2426de9147fb66256d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Tue, 28 Jul 2020 09:39:22 +0300 Subject: [PATCH] Enable separation of long- and short-range of Ewald sum in python --- qpms/cycommon.pyx | 6 ++++++ qpms/ewald.h | 8 -------- qpms/qpms_c.pyx | 13 ++++++++----- qpms/qpms_cdefs.pxd | 15 +++++++++------ qpms/qpms_types.h | 7 +++++++ qpms/scatsystem.c | 30 ++++++++++++++++++++++-------- qpms/scatsystem.h | 16 ++++++++++++++++ 7 files changed, 68 insertions(+), 27 deletions(-) diff --git a/qpms/cycommon.pyx b/qpms/cycommon.pyx index 88c557d..31abcfa 100644 --- a/qpms/cycommon.pyx +++ b/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 diff --git a/qpms/ewald.h b/qpms/ewald.h index 8eda80c..d7bad57 100644 --- a/qpms/ewald.h +++ b/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; diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 7f41b30..f6f2362 100644 --- a/qpms/qpms_c.pyx +++ b/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 @@ -1026,7 +1026,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 @@ -1045,6 +1046,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") @@ -1060,14 +1062,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) @@ -1089,6 +1091,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: @@ -1104,7 +1107,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 diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index daa2f29..45a18e0 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -150,6 +150,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": @@ -708,6 +713,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, @@ -718,6 +725,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 @@ -726,12 +735,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 diff --git a/qpms/qpms_types.h b/qpms/qpms_types.h index e0b0b4f..4d465ba 100644 --- a/qpms/qpms_types.h +++ b/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 diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index 95743ba..ba84d07 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -2179,10 +2179,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) @@ -2214,13 +2215,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){ @@ -2234,11 +2235,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) @@ -2272,13 +2279,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){ @@ -2294,6 +2301,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) { diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index 3704ac9..78e90fd 100644 --- a/qpms/scatsystem.h +++ b/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