diff --git a/qpms/CMakeLists.txt b/qpms/CMakeLists.txt index 01284bd..7bbd52c 100644 --- a/qpms/CMakeLists.txt +++ b/qpms/CMakeLists.txt @@ -16,6 +16,7 @@ add_library (qpms SHARED translations.c tmatrices.c vecprint.c vswf.c wigner.c e lattices2d.c gaunt.c error.c legendre.c symmetries.c vecprint.c bessel.c own_zgemm.c parsing.c scatsystem.c materials.c drudeparam_data.c lll.c beyn.c trivialgroup.c + translations_dbg.c ) use_c99() diff --git a/qpms/translations_dbg.c b/qpms/translations_dbg.c new file mode 100644 index 0000000..a2f998c --- /dev/null +++ b/qpms/translations_dbg.c @@ -0,0 +1,148 @@ +#include "qpms_types.h" +#include "qpms_specfunc.h" +#include "translations_dbg.h" +#include "indexing.h" +#include "tiny_inlines.h" +#include "qpms_error.h" +#include "normalisation.h" + + +int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c, + complex double *dest, const csph_t kdlj, const qpms_bessel_t J) { + if (0 == kdlj.r && J != QPMS_BESSEL_REGULAR) { + for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l) + for (qpms_m_t m = 0; m <= l; ++m) { + const qpms_y_t y_sc = qpms_mn2y_sc(m, l); + dest[y_sc] = NAN+I*NAN; + } + // TODO warn? different return value? + return 0; + } + double legendre_buf[gsl_sf_legendre_array_n(c->lMax+c->lMax+1)]; + complex double bessel_buf[2*c->lMax+2]; // maximum order is 2n+1 for B coeffs, plus the zeroth. + double costheta = cos(kdlj.theta); + QPMS_ENSURE_SUCCESS(gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,2*c->lMax+1, + costheta,-1,legendre_buf)); + QPMS_ENSURE_SUCCESS(qpms_sph_bessel_fill(J, 2*c->lMax+1, kdlj.r, bessel_buf)); + for (qpms_l_t l = 0; l <= 2*c->lMax+1; ++l) + for (qpms_m_t m = -l; m <= l; ++m) { + const qpms_y_t y_sc = qpms_mn2y_sc(m, l); + dest[y_sc] = bessel_buf[l] + * legendre_buf[gsl_sf_legendre_array_index(l, abs(m))] + * cexp(I * m * kdlj.phi); + } + return 0; +} + +int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c, + complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J) { + csph_t kdlj = cart2csph(r); + kdlj.r *= k; + return qpms_trans_calculator_test_sswf(c, dest, kdlj, J); +} + +#ifdef LATTICESUMS32 +int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c, + complex double * const sigmas_total, double * serr_total, + const double eta, const complex double k, + const cart2_t b1, const cart2_t b2, + const cart2_t beta, + const cart3_t particle_shift, + double maxR, double maxK, + const qpms_ewald_part parts) +{ + const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e3c->lMax); + //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax); + const bool doerr = !!serr_total; + const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 0) && (particle_shift.z == 0)); // FIXME ignoring the case where particle_shift equals to lattice vector + + complex double *sigmas_short = malloc(sizeof(complex double)*nelem2_sc); + complex double *sigmas_long = malloc(sizeof(complex double)*nelem2_sc); + //complex double *sigmas_total = malloc(sizeof(complex double)*nelem2_sc); + double *serr_short, *serr_long;//, *serr_total; + if(doerr) { + serr_short = malloc(sizeof(double)*nelem2_sc); + serr_long = malloc(sizeof(double)*nelem2_sc); + //serr_total = malloc(sizeof(double)*nelem2_sc); + } else serr_short = serr_long = serr_total = NULL; + + const double unitcell_area = l2d_unitcell_area(b1, b2); + cart2_t rb1, rb2; // reciprocal basis + QPMS_ENSURE_SUCCESS(l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2)); + + if (parts & QPMS_EWALD_LONG_RANGE) { + PGen Kgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, +#ifdef GEN_KSHIFTEDPOINTS + beta, +#else + CART2_ZERO, +#endif + 0, true, maxK, false); + + QPMS_ENSURE_SUCCESS(ewald3_sigma_long(sigmas_long, serr_long, c->e3c, eta, k, + unitcell_area, LAT_2D_IN_3D_XYONLY, &Kgen, +#ifdef GEN_KSHIFTEDPOINTS + true, +#else + false, +#endif + cart22cart3xy(beta), particle_shift)); + if(Kgen.stateData) // PGen not consumed entirely (converged earlier) + PGen_destroy(&Kgen); + } + + if (parts & QPMS_EWALD_SHORT_RANGE) { + PGen Rgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, +#ifdef GEN_RSHIFTEDPOINTS + cart2_scale(-1 /*CHECKSIGN*/, cart3xy2cart2(particle_shift)), +#else + CART2_ZERO, +#endif + 0, !do_sigma0, maxR, false); +#ifdef GEN_RSHIFTEDPOINTS // rather ugly hacks, LPTODO cleanup + if (particle_shift.z != 0) { + const cart3_t zshift = {0, 0, -particle_shift.z /*CHECKSIGN*/}; + Rgen = Pgen_shifted_new(Rgen, zshift); + } +#endif + + + QPMS_ENSURE_SUCCESS(ewald3_sigma_short(sigmas_short, serr_short, c->e3c, eta, k, + particle_shift.z ? LAT_2D_IN_3D : LAT_2D_IN_3D_XYONLY, &Rgen, +#ifdef GEN_RSHIFTEDPOINTS + true, +#else + false, +#endif + cart22cart3xy(beta), 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] = ((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] = ((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 && (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; + if(doerr) serr_total[y] += sigma0_err; + } + + free(sigmas_short); + free(sigmas_long); + //free(sigmas_total); + if(doerr) { + free(serr_short); + free(serr_long); + //free(serr_total); + } + return 0; +} +#endif // LATTICESUMS32 diff --git a/qpms/translations_dbg.h b/qpms/translations_dbg.h new file mode 100644 index 0000000..13aa0fc --- /dev/null +++ b/qpms/translations_dbg.h @@ -0,0 +1,28 @@ +#ifndef TRANSLATIONS_DBG_H + +// the following functions evaluate the SSWFs the same way as +// certain functions in translations.h, but they leave them +// as they are, without rearranging them into the translation +// operators. + +#include "translations.h" + +int qpms_trans_calculator_test_sswf(const qpms_trans_calculator *c, + complex double *dest, const csph_t kdlj, const qpms_bessel_t J); + +int qpms_trans_calculator_test_sswf_lc3p(const qpms_trans_calculator *c, + complex double *dest, const complex double k, const cart3_t r, const qpms_bessel_t J); + + +#ifdef LATTICESUMS32 +int qpms_trans_calculator_test_sswf_e32(const qpms_trans_calculator *c, + complex double * const sigmas_total, double * serr_total, + const double eta, const complex double k, + const cart2_t b1, const cart2_t b2, + const cart2_t beta, + const cart3_t particle_shift, + double maxR, double maxK, + const qpms_ewald_part parts); +#endif // LATTICESUMS32 + +#endif // TRANSLATIONS_DBG_H