diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 12d2c27..2723da8 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -568,7 +568,10 @@ cdef class ScatteringSystem: return None property eta: - """Ewald parameter η""" + """Ewald parameter η (only relevant for periodic systems) + + This is the default value that will be copied into + _ScatteringSystemAtOmega by __call__""" def __get__(self): self.check_s() if self.lattice_dimension: @@ -723,7 +726,7 @@ cdef class ScatteringSystem: self.s, iri, 0) return target_np - def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS): + def translation_matrix_full(self, cdouble wavenumber, blochvector = None, J = QPMS_HANKEL_PLUS, double eta = 0): """Constructs full translation matrix of a scattering system. This method enables to use any wave number for the background medium ignoring the @@ -741,6 +744,11 @@ cdef class ScatteringSystem: J : BesselType Optionally, one can replace Hankel functions of the first kind with different Bessel functions. + + eta : float + Ewald parameter η, only relevant for periodic lattices; + if set to 0 (default), the actual value is determined + automatically. See Also -------- @@ -763,7 +771,7 @@ cdef class ScatteringSystem: raise NotImplementedError("Translation operators based on other than Hankel+ functions not supperted in periodic systems") blochvector_c = {'x': blochvector[0], 'y': blochvector[1], 'z': blochvector[2]} with pgsl_ignore_error(15): - qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c) + qpms_scatsys_periodic_build_translation_matrix_full(&target_view[0][0], self.s, wavenumber, &blochvector_c, eta) return target def translation_matrix_packed(self, cdouble wavenumber, qpms_iri_t iri, J = QPMS_HANKEL_PLUS): @@ -1025,6 +1033,20 @@ 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 8422ab1..937b4a1 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -630,6 +630,7 @@ 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,7 +689,7 @@ cdef extern from "scatsystem.h": const qpms_scatsys_at_omega_t *ssw double k[3] 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) + 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) beyn_result_t *qpms_scatsys_periodic_find_eigenmodes(const qpms_scatsys_t *ss, const double *k, cdouble omega_centre, double omega_rr, double omega_ri, size_t contour_npoints, diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index a89ebb7..db59b02 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -59,6 +59,13 @@ 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) { + qpms_ss_ensure_periodic(ss); + return ss->per.eta; +} + // ------------ Stupid implementation of qpms_scatsys_apply_symmetry() ------------- #define MIN(x,y) (((x)<(y))?(x):(y)) @@ -541,6 +548,9 @@ 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; } @@ -594,6 +604,7 @@ 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)); @@ -1156,7 +1167,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); + return qpms_scatsys_periodic_build_translation_matrix_full(target, ss, wavenumber, &k_cart3, ssw->eta); } complex double *qpms_scatsys_build_translation_matrix_e_full( @@ -1201,17 +1212,17 @@ complex double *qpms_scatsys_build_translation_matrix_e_full( static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss, qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const cart2_t kvector, complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride, - qpms_ewald_part parts) { + qpms_ewald_part parts, double eta) { const qpms_vswf_set_spec_t *srcspec = qpms_ss_bspec_pi(ss, psrc); const qpms_vswf_set_spec_t *destspec = qpms_ss_bspec_pi(ss, pdest); // This might be a bit arbitrary, roughly "copied" from Unitcell constructor. TODO review. - const double maxR = sqrt(ss->per.unitcell_volume) * 32; - const double maxK = 1024 * 2 * M_PI / maxR; + const double maxR = sqrt(ss->per.unitcell_volume) * 64; + const double maxK = 2048 * 2 * M_PI / maxR; return qpms_trans_calculator_get_trans_array_e32_e(ss->c, target, NULL /*err*/, destspec, deststride, srcspec, srcstride, - ss->per.eta, wavenumber, + eta, wavenumber, cart3xy2cart2(ss->per.lattice_basis[0]), cart3xy2cart2(ss->per.lattice_basis[1]), kvector, cart2_substract(cart3xy2cart2(ss->p[pdest].pos), cart3xy2cart2(ss->p[psrc].pos)), @@ -1221,21 +1232,23 @@ static inline int qpms_ss_ppair_W32xy(const qpms_scatsys_t *ss, static inline int qpms_ss_ppair_W(const qpms_scatsys_t *ss, qpms_ss_pi_t pdest, qpms_ss_pi_t psrc, complex double wavenumber, const double wavevector[], complex double *target, const ptrdiff_t deststride, const ptrdiff_t srcstride, - qpms_ewald_part parts) { + qpms_ewald_part parts, double eta) { if(ss->lattice_dimension == 2 && // Currently, we can only the xy-plane !ss->per.lattice_basis[0].z && !ss->per.lattice_basis[1].z && !wavevector[2]) return qpms_ss_ppair_W32xy(ss, pdest, psrc, wavenumber, cart2_from_double_array(wavevector), - target, deststride, srcstride, parts); + target, deststride, srcstride, parts, eta); else QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported"); } complex double *qpms_scatsys_periodic_build_translation_matrix_full( complex double *target, const qpms_scatsys_t *ss, - complex double wavenumber, const cart3_t *wavevector) { + 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; const size_t full_len = ss->fecv_size; if(!target) QPMS_CRASHING_MALLOC(target, SQ(full_len) * sizeof(complex double)); @@ -1248,7 +1261,7 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full( for (qpms_ss_pi_t ps = 0; ps < ss->p_count; ++ps) { QPMS_ENSURE_SUCCESS(qpms_ss_ppair_W32xy(ss, pd, ps, wavenumber, cart3xy2cart2(*wavevector), target + deststride * ss->fecv_pstarts[pd] + srcstride * ss->fecv_pstarts[ps], - deststride, srcstride, QPMS_EWALD_FULL)); + deststride, srcstride, QPMS_EWALD_FULL, eta)); } } else QPMS_NOT_IMPLEMENTED("Only 2D xy-lattices currently supported"); @@ -1294,7 +1307,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)); + QPMS_EWALD_FULL, ssw->eta)); } cblas_zgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, diff --git a/qpms/scatsystem.h b/qpms/scatsystem.h index b40e347..dd96cac 100644 --- a/qpms/scatsystem.h +++ b/qpms/scatsystem.h @@ -147,7 +147,12 @@ typedef struct qpms_scatsys_periodic_info_t { */ double unitcell_volume; - /// Ewald parameter \f$ \eta \f$. + /// Default Ewald parameter \f$ \eta \f$. + /** Normally, this just gets copied into qpms_scatsys_at_omega_t, + * which is then used in the Ewald sums. + * However, for higher frequencies it must be adjusted to avoid + * numerical instability. + */ double eta; } qpms_scatsys_periodic_info_t; @@ -241,6 +246,7 @@ 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; @@ -501,7 +507,8 @@ complex double *qpms_scatsys_periodic_build_translation_matrix_full( complex double *target, const qpms_scatsys_t *ss, complex double wavenumber, ///< Wave number to use in the translation matrix. - const cart3_t *wavevector ///< Wavevector / pseudomomentum in cartesian coordinates. + const cart3_t *wavevector, ///< Wavevector / pseudomomentum in cartesian coordinates. + double eta ///< Ewald parameter eta. Pass 0 or NaN to use the default value in \a ss. ); /// Global translation matrix.