From 1f7ac7c1c1604e79182e235a68fe61805f4e3db7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 19 Sep 2018 09:16:38 +0300 Subject: [PATCH] New hexlattice_ewald.c with real ewald sum Former-commit-id: cb398312fcb1d78f39cc374e3171e216d7604f90 --- qpms/apps/hexlattice_ewald.c | 171 +++++++++++++++++++++++++++++++++++ qpms/translations.c | 14 +-- qpms/translations.h | 1 + 3 files changed, 179 insertions(+), 7 deletions(-) create mode 100644 qpms/apps/hexlattice_ewald.c diff --git a/qpms/apps/hexlattice_ewald.c b/qpms/apps/hexlattice_ewald.c new file mode 100644 index 0000000..9eab2b5 --- /dev/null +++ b/qpms/apps/hexlattice_ewald.c @@ -0,0 +1,171 @@ +// c99 -o ew -Wall -I ../.. -O2 -ggdb -DLATTICESUMS32 hexlattice_ewald.c ../translations.c ../ewald.c ../ewaldsf.c ../gaunt.c ../lattices2d.c -lgsl -lm -lblas +#include +#include +#include +#include +#include + + + +#define MAXOMEGACOUNT 1000 +#define MAXKCOUNT 100 +#define KLAYERS 20 +#define RLAYERS 20 + +const double s3 = 1.732050807568877293527446341505872366942805253810380628055; + +// IMPORTANT: lattice properties here +const qpms_y_t lMax = 2; +const double REFINDEX = 1.52; +const double LATTICE_H = 576e-9; +static const double SCUFF_OMEGAUNIT = 3e14; +static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR; +static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE; +static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT; +static const TriangularLatticeOrientation rs_orientation = TRIANGULAR_VERTICAL; + +int main (int argc, char **argv) { + const double LATTICE_A = s3*LATTICE_H; + const double INVLATTICE_A = 4*M_PI / s3 / LATTICE_A; + + char *omegafile = argv[1]; + // char *kfile = argv[2]; // not used + char *outfile = argv[3]; + char *errfile = NULL; + if (argc > 4) + errfile = argv[4]; + //char *outlongfile = argv[4]; + //char *outshortfile = argv[5]; + double scuffomegas[MAXOMEGACOUNT]; + cart2_t klist[MAXKCOUNT]; + FILE *f = fopen(omegafile, "r"); + size_t omegacount = 0; + while (fscanf(f, "%lf", scuffomegas + omegacount) == 1){ + assert(omegacount < MAXOMEGACOUNT); + ++omegacount; + } + fclose(f); + /*f = fopen(kfile, "r"); + int kcount = 100; + while (fscanf(f, "%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) { + assert(kcount < MAXKCOUNT); + ++kcount; + } + fclose(f); + */ + + int kcount = MAXKCOUNT; + for (int i = 0; i < kcount; ++i) { // TODO this should depend on orientation... + klist[i].x = 0; + klist[i].y = 2. * 4. * M_PI / 3. / LATTICE_A / kcount * i; + } + + const double refindex = REFINDEX; + const double h = LATTICE_H; + const double a = h * s3; + const double unitcell_area = s3*a*a/2.; + const double rec_a = 4*M_PI/s3/a; + + triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, rs_orientation, true, 0); + triangular_lattice_gen_extend_to_steps(Rlg, RLAYERS); + triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(rec_a, reverseTriangularLatticeOrientation(rs_orientation), true, 0); + triangular_lattice_gen_extend_to_steps(Klg, KLAYERS); + + const point2d *Rpoints = Rlg->ps.base; + size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs]; + const point2d *Kpoints = Klg->ps.base; + size_t nK = Klg->ps.r_offsets[Klg->ps.nrs]; + + const point2d pshift0 = {0, 0}; + point2d pshiftAB = {0, 0}, pshiftBA = {0,0}; + if(rs_orientation == TRIANGULAR_VERTICAL) { // CHECKSIGN + pshiftAB.x = h; + pshiftBA.x = -h; + } else { // CHECKSIGN + pshiftAB.y = -h; + pshiftBA.y = h; + } + + qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER); // vai POWER_CS? + + FILE *out = fopen(outfile, "w"); + FILE *err = NULL; + if (errfile) + err = fopen(errfile, "w"); + + for (size_t omegai = 0; omegai < omegacount; ++omegai) { + const double scuffomega = scuffomegas[omegai]; + const double omega = scuffomega * SCUFF_OMEGAUNIT; + const double EeV = omega * hbar / eV; + const double k0_vac = omega / c0; + const double k0_eff = k0_vac * refindex; + const double eta = 4*rec_a; // FIXME quite arbitrary + + // indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy + complex double W[2][2][2][c->nelem][c->nelem]; + double Werr[2][2][2][c->nelem][c->nelem]; + + for (size_t ki = 0; ki < kcount; ++ki) { + cart2_t beta = klist[ki]; + memset(W, 0, sizeof(W)); + if(err) + memset(Werr, 0, sizeof(Werr)); + + const ptrdiff_t deststride = &(W[0][0][0][1][0]) - &(W[0][0][0][0][0]); + const ptrdiff_t srcstride = &(W[0][0][0][0][1]) - &(W[0][0][0][0][0]); + assert (srcstride == 1 && deststride == c->nelem); + + // A<-A + qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c, + &(W[0][0][0][0][0]), err ? &(Werr[0][0][0][0][0]) : NULL, // Adest, Aerr, + &(W[0][0][1][0][0]), err ? &(Werr[0][0][1][0][0]) : NULL, // Bdest, Berr, + deststride, srcstride, + eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshift0 + ); + // B<-B + for(qpms_y_t desty = 0; desty < c->nelem; ++desty) + for (qpms_y_t srcy = 0; srcy < c->nelem; ++srcy) + for (int t = 0; t < 2; ++t) { + W[1][1][t][desty][srcy] = W[0][0][t][desty][srcy]; + if (err) + Werr[1][1][t][desty][srcy] = Werr[0][0][t][desty][srcy]; + } + + // A<-B + qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c, + &(W[0][1][0][0][0]), err ? &(Werr[0][1][0][0][0]) : NULL, // Adest, Aerr, + &(W[0][1][1][0][0]), err ? &(Werr[0][1][1][0][0]) : NULL, // Bdest, Berr, + deststride, srcstride, + eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshiftAB + ); + // B<-A + qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(c, + &(W[1][0][0][0][0]), err ? &(Werr[1][0][0][0][0]) : NULL, // Adest, Aerr, + &(W[1][0][1][0][0]), err ? &(Werr[1][0][1][0][0]) : NULL, // Bdest, Berr, + deststride, srcstride, + eta, k0_eff, unitcell_area, nR, Rpoints, nK, Kpoints, beta, pshiftBA + ); + // TODO CHECK B<-A vs. A<-B relation + + fprintf(out, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t", + scuffomega, EeV, k0_eff, beta.x, beta.y); + if(err) fprintf(err, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t", + scuffomega, EeV, k0_eff, beta.x, beta.y); + size_t totalelems = sizeof(W) / sizeof(complex double); + for (size_t i = 0; i < totalelems; ++i) { + complex double w = ((complex double *)W)[i]; + fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w)); + if (err) + fprintf(err, "%.3g\t", ((double *)Werr)[i]); + } + fputc('\n', out); + if(err) fputc('\n', err); + } + } + + fclose(out); + if(err) fclose(err); + + triangular_lattice_gen_free(Klg); + triangular_lattice_gen_free(Rlg); +} diff --git a/qpms/translations.c b/qpms/translations.c index fed2e9f..564f4ad 100644 --- a/qpms/translations.c +++ b/qpms/translations.c @@ -1162,22 +1162,22 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra complex double * const Bdest, double * const Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, /* qpms_bessel_t J*/ // assume QPMS_HANKEL_PLUS - const double eta, const double k, + const double eta, const double k, const double unitcell_area, const size_t nRpoints, const cart2_t *Rpoints, - const size_t nKpoints, const cart2_t *Kpoinst, + const size_t nKpoints, const cart2_t *Kpoints, const cart2_t beta, const cart2_t particle_shift ) { const qpms_y_t nelem2_sc = qpms_lMax2nelem_sc(c->e32c->lMax); - const qpms_y_t nelem = qpms_lMax2nelem(c->lMax); + //const qpms_y_t nelem = qpms_lMax2nelem(c->lMax); const bool doerr = Aerr || Berr; const bool do_sigma0 = ((particle_shift.x == 0) && (particle_shift.y == 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)*nelem_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); @@ -1187,7 +1187,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra int retval; retval = ewald32_sigma_long_points_and_shift(sigmas_long, serr_long, - c->e32c, eta, k, nKpoints, Kpoints, beta, particle_shift); + c->e32c, eta, k, unitcell_area, nKpoints, Kpoints, beta, particle_shift); if (retval) abort(); retval = ewald32_sigma_short_points_and_shift(sigmas_short, serr_short, @@ -1216,7 +1216,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra ptrdiff_t desti = 0, srci = 0; for (qpms_l_t n = 1; n <= c->lMax; ++n) for (qpms_m_t m = -n; m <= n; ++m) { for (qpms_l_t nu = 1; nu <= c->lMax; ++nu) for (qpms_m_t mu = -nu; mu <= nu; ++mu){ - const size_t i = qpms_trans_calculator_index_mnmu(c, m, n, mu, nu); + const size_t i = qpms_trans_calculator_index_mnmunu(c, m, n, mu, nu); const size_t qmax = c->A_multipliers[i+1] - c->A_multipliers[i] - 1; complex double Asum, Asumc; ckahaninit(&Asum, &Asumc); double Asumerr, Asumerrc; if(Aerr) kahaninit(&Asumerr, &Asumerrc); @@ -1258,7 +1258,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra } break; default: - abort() + abort(); } free(sigmas_short); diff --git a/qpms/translations.h b/qpms/translations.h index 44129e6..05dbb1d 100644 --- a/qpms/translations.h +++ b/qpms/translations.h @@ -186,6 +186,7 @@ int qpms_trans_calculator_get_AB_arrays_e32_both_points_and_shift(const qpms_tra complex double *Bdest, double *Berr, const ptrdiff_t deststride, const ptrdiff_t srcstride, const double eta, const double k, + const double unitcell_area, const size_t nRpoints, const cart2_t *Rpoints, const size_t nKpoints, const cart2_t *Kpoinst, const cart2_t beta,