diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 2723da8..2d3a01e 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -963,6 +963,13 @@ cdef class _ScatteringSystemAtOmegaK: cdef qpms_scatsys_at_omega_k_t *rawpointer(self): return &self.sswk + + property eta: + """Ewald parameter η""" + def __get__(self): + return self.sswk.eta + def __set__(self, double eta): + self.sswk.eta = eta cdef class _ScatteringSystemAtOmega: ''' @@ -1021,6 +1028,7 @@ cdef class _ScatteringSystemAtOmega: sswk.sswk.k[0] = k[0] sswk.sswk.k[1] = k[1] sswk.sswk.k[2] = k[2] + sswk.eta = qpms_ss_adjusted_eta(self.ssw[0].ss, self.ssw[0].wavenumber, sswk.sswk.k) return sswk property fecv_size: @@ -1033,20 +1041,6 @@ cdef class _ScatteringSystemAtOmega: def __get__(self): return self.ss_pyref.nirreps property wavenumber: def __get__(self): return self.ssw[0].wavenumber - property eta: - """Ewald parameter η (only relevant for periodic systems)""" - def __get__(self): - self.check_s() - if self.lattice_dimension: - return self.ssw[0].eta - else: - return None - def __set__(self, double eta): - self.check_s() - if self.lattice_dimension: - self.ssw[0].eta = eta - else: - raise AttributeError("Cannot set Ewald parameter for finite system") # different exception? def modeproblem_matrix_full(self, k=None): diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index 937b4a1..f9e44ee 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -630,7 +630,6 @@ cdef extern from "scatsystem.h": cdouble omega qpms_epsmu_t medium cdouble wavenumber - double eta qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym, cdouble omega, const qpms_tolerance_spec_t *tol) qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, cdouble omega) @@ -688,6 +687,7 @@ cdef extern from "scatsystem.h": struct qpms_scatsys_at_omega_k_t: const qpms_scatsys_at_omega_t *ssw double k[3] + double eta cdouble *qpms_scatsyswk_build_modeproblem_matrix_full(cdouble *target, const qpms_scatsys_at_omega_k_t *sswk) cdouble *qpms_scatsys_periodic_build_translation_matrix_full(cdouble *target, const qpms_scatsys_t *ss, cdouble wavenumber, const cart3_t *wavevector, double eta) qpms_ss_LU qpms_scatsyswk_build_modeproblem_matrix_full_LU(cdouble *target, int *target_piv, const qpms_scatsys_at_omega_k_t *sswk) @@ -703,6 +703,8 @@ cdef extern from "scatsystem.h": const cdouble *f_excitation_vector_full, cart3_t where) ccart3_t qpms_scatsysw_scattered_E__alt(const qpms_scatsys_at_omega_t *ssw, qpms_bessel_t btyp, const cdouble *f_excitation_vector_full, cart3_t where) + double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, cdouble wavenumber, const double *wavevector); + cdef extern from "ewald.h": struct qpms_csf_result: diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index db59b02..c9551f2 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -23,6 +23,7 @@ #include "kahansum.h" #include "tolerances.h" #include "beyn.h" +#include "tiny_inlines.h" #ifdef QPMS_SCATSYSTEM_USE_OWN_BLAS #include "qpmsblas.h" @@ -31,11 +32,15 @@ #define SERIAL_ZGEMM cblas_zgemm #endif -#define SQ(x) ((x)*(x)) #define QPMS_SCATSYS_LEN_RTOL 1e-13 #define QPMS_SCATSYS_TMATRIX_ATOL 1e-12 #define QPMS_SCATSYS_TMATRIX_RTOL 1e-12 +// This is used in adjustment of Ewald parameter to avoid high frequency breakdown. +// Very roughly, the value of 16 should lead to allowing terms containing incomplete Gamma +// functions with magnitudes around exp(16) == 8.9e6 +static const double QPMS_SCATSYS_EWALD_MAX_EXPONENT = 16.; + long qpms_scatsystem_nthreads_default = 4; long qpms_scatsystem_nthreads_override = 0; @@ -59,18 +64,23 @@ static inline void qpms_ss_ensure_nonperiodic_a(const qpms_scatsys_t *ss, const QPMS_ENSURE(ss->lattice_dimension == 0, "This method is applicable only to nonperiodic systems. Use %s instead.", s); } -// Adjust Ewald parameter to avoid high-frequency breakdown; -// TODO make this actually do something (other than just copying ss's eta.) -static inline double ss_adjusted_eta(const qpms_scatsys_t *ss, complex double omega) { +// Adjust Ewald parameter to avoid high-frequency breakdown +double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double k[3]) { qpms_ss_ensure_periodic(ss); - return ss->per.eta; + const double eta_default = ss->per.eta; + // FIXME here we silently assume that k lies in the first Brillioun zone, we should ensure that. + const double k2 = k[0]*k[0] + k[1]*k[1] + k[2] * k[2]; + const double kappa2 = SQ(cabs(wavenumber)); // maybe creal would be enough + if(kappa2 < k2) // This should happen only for pretty low frequencies + return eta_default; + const qpms_l_t maxj = ss->c->lMax; // Based on ewald.c:301, note that VSWF (c's) lMax is already half of corresponding translation matrix Ewald factors' (c->e3c's) lMax + const double mina = 0.5 * (ss->lattice_dimension - 1) - maxj; // minimum incomplete Gamma first argument, based on ewald.c:301; CHECKME whether this is fine also for 3D lattice + const double eta_min = sqrt(fabs((kappa2 - k2) * (mina - 1.) / QPMS_SCATSYS_EWALD_MAX_EXPONENT)); + return MAX(eta_default, eta_min); } // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() ------------- -#define MIN(x,y) (((x)<(y))?(x):(y)) -#define MAX(x,y) (((x)>(y))?(x):(y)) - // The following functions are just to make qpms_scatsys_apply_symmetry more readable. // They are not to be used elsewhere, as they do not perform any array capacity checks etc. @@ -549,8 +559,6 @@ qpms_scatsys_at_omega_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, ss->c = qpms_trans_calculator_init(lMax, normalisation); - ssw->eta = ss->lattice_dimension ? ss_adjusted_eta(ss, omega) : NAN; - return ssw; } @@ -604,7 +612,6 @@ qpms_scatsys_at_omega_t *qpms_scatsys_at_omega(const qpms_scatsys_t *ss, ssw->ss = ss; ssw->medium = qpms_epsmu_generator_eval(ss->medium, omega); ssw->wavenumber = qpms_wavenumber(omega, ssw->medium); - ssw->eta = ss->lattice_dimension ? ss_adjusted_eta(ss, omega) : NAN; QPMS_CRASHING_CALLOC(ssw->tm, ss->tm_count, sizeof(*ssw->tm)); qpms_tmatrix_t **tmatrices_preop; QPMS_CRASHING_CALLOC(tmatrices_preop, ss->tmg_count, sizeof(*tmatrices_preop)); @@ -1167,7 +1174,7 @@ complex double *qpms_scatsyswk_build_translation_matrix_full( const qpms_scatsys_t *ss = ssw->ss; qpms_ss_ensure_periodic(ss); const cart3_t k_cart3 = cart3_from_double_array(sswk->k); - return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3, ssw->eta); + return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3, sswk->eta); } complex double *qpms_scatsys_build_translation_matrix_e_full( @@ -1247,8 +1254,12 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full( complex double wavenumber, const cart3_t *wavevector, double eta) { QPMS_UNTESTED; qpms_ss_ensure_periodic(ss); - if (eta == 0 || isnan(eta)) - eta = ss->per.eta; + if (eta == 0 || isnan(eta)) { + double tmp[3]; + cart3_to_double_array(tmp, *wavevector); + eta = qpms_ss_adjusted_eta(ss, wavenumber, tmp); + } + const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); @@ -1273,7 +1284,8 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f /// Target memory with capacity for ss->fecv_size**2 elements. If NULL, new will be allocated. complex double *target, const qpms_scatsys_at_omega_t *ssw, - const double k[] // NULL if non-periodic + const double k[], // NULL if non-periodic + const double eta // ignored if non-periodic ) { const complex double wavenumber = ssw->wavenumber; @@ -1307,7 +1319,7 @@ static inline complex double *qpms_scatsysw_scatsyswk_build_modeproblem_matrix_f } else { // periodic case QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W(ss, piR, piC, wavenumber, k, tmp /*target*/, bspecC->n /*deststride*/, 1 /*srcstride*/, - QPMS_EWALD_FULL, ssw->eta)); + QPMS_EWALD_FULL, eta)); } cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, @@ -1334,14 +1346,14 @@ complex double *qpms_scatsysw_build_modeproblem_matrix_full( complex double *target, const qpms_scatsys_at_omega_t *ssw) { qpms_ss_ensure_nonperiodic_a(ssw->ss, "qpms_scatsyswk_build_modeproblem_matrix_full()"); return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full( - target, ssw, NULL); + target, ssw, NULL, NAN); } complex double *qpms_scatsyswk_build_modeproblem_matrix_full( complex double *target, const qpms_scatsys_at_omega_k_t *sswk) { qpms_ss_ensure_periodic_a(sswk->ssw->ss, "qpms_scatsysw_build_modeproblem_matrix_full()"); - return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(target, sswk->ssw, sswk->k); + return qpms_scatsysw_scatsyswk_build_modeproblem_matrix_full(target, sswk->ssw, sswk->k, sswk->eta); } @@ -2246,7 +2258,8 @@ static int qpms_scatsys_periodic_eval_Beyn_ImTW(complex double *target, QPMS_ENSURE(ssw != NULL, "qpms_scatsys_at_omega() returned NULL"); qpms_scatsys_at_omega_k_t sswk = { .ssw = ssw, - .k = {p->k[0], p->k[1], p->k[2]} + .k = {p->k[0], p->k[1], p->k[2]}, + .eta = qpms_ss_adjusted_eta(p->ss, ssw->wavenumber, p->k) }; QPMS_ASSERT(m == p->ss->fecv_size); QPMS_ENSURE(NULL != diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index dd96cac..4e99111 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -246,7 +246,6 @@ typedef struct qpms_scatsys_at_omega_t { complex double omega; ///< Angular frequency qpms_epsmu_t medium; ///< Background medium optical properties at the given frequency complex double wavenumber; ///< Background medium wave number - double eta; ///< Ewald parameter \f$ \eta \f$. } qpms_scatsys_at_omega_t; @@ -489,6 +488,7 @@ complex double *qpms_scatsys_scatter_solve( typedef struct qpms_scatsys_at_omega_k_t { const qpms_scatsys_at_omega_t *ssw; double k[3]; ///< The k-vector's cartesian coordinates. + double eta; ///< Ewald parameter η. } qpms_scatsys_at_omega_k_t; /// Creates the full \f$ (I - WS) \f$ matrix of the periodic scattering system. @@ -769,6 +769,10 @@ ccart3_t qpms_scatsysw_scattered_E__alt( cart3_t evalpoint ///< A point \f$ \vect r \f$, at which the field is evaluated. ); +/// Adjusted Ewadl parameter to avoid high-frequency breakdown. +// TODO DOC +double qpms_ss_adjusted_eta(const qpms_scatsys_t *ss, complex double wavenumber, const double wavevector[3]); + #if 0 /** Evaluates partial scattered fields (corresponding to a given irrep-reduced excitation vector) * at a given point.