SSWF evaluation functions for debugging Ewald sum.
This commit is contained in:
Marek Nečada
parent
6ccd785164
commit
3a8e5b99cb
|
|
@ -24,6 +24,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 version.c
|
||||
translations_dbg.c
|
||||
)
|
||||
use_c99()
|
||||
|
||||
|
|
|
|||
|
|
@ -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
|
||||
|
|
@ -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
|
||||
Loading…
Reference in New Issue