Hexlattice ewald summation compiles
Former-commit-id: 4bdbd08da64527e77926f66656ef1b0b81546cdb
This commit is contained in:
Marek Nečada
parent
2b00af241e
commit
81803d16a4
|
|
@ -1,10 +1,28 @@
|
||||||
|
// c99 -ggdb -O2 -DLATTICESUMS -I .. hexlattice_ewald.c ../translations.c ../bessels.c ../lrhankel_recspace_dirty.c ../gaunt.c -lm -lgsl -lblas
|
||||||
#include <stdlib.h>
|
#include <stdlib.h>
|
||||||
|
#include <string.h>
|
||||||
#include <math.h>
|
#include <math.h>
|
||||||
#include <stdio.h>
|
#include <stdio.h>
|
||||||
#include "kahansum.h"
|
#include "kahansum.h"
|
||||||
#include "vectors.h"
|
#include "vectors.h"
|
||||||
|
#include <gsl/gsl_const_mksa.h>
|
||||||
|
#include <gsl/gsl_math.h>
|
||||||
|
#include "qpms_types.h"
|
||||||
|
#include "translations.h"
|
||||||
|
|
||||||
static const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
|
#define MAXOMEGACOUNT 1000
|
||||||
|
#define MAXKCOUNT 10000
|
||||||
|
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 double CC = 0.1;
|
||||||
|
|
||||||
// For sorting the points by distance from origin / radius
|
// For sorting the points by distance from origin / radius
|
||||||
int cart2_cmpr (const void *p1, const void *p2) {
|
int cart2_cmpr (const void *p1, const void *p2) {
|
||||||
|
|
@ -120,11 +138,190 @@ latticepoints_circle_t generate_tripoints_hor(double a, double R, cart2_t offset
|
||||||
}
|
}
|
||||||
|
|
||||||
int main (int argc, char **argv) {
|
int main (int argc, char **argv) {
|
||||||
double h = 1.2;
|
const double LATTICE_A = s3*LATTICE_H;
|
||||||
|
const double INVLATTICE_A = 4 * M_PI / s3 / LATTICE_A;
|
||||||
|
const double MAXR_REAL = 100 * LATTICE_H;
|
||||||
|
const double MAXR_K = 100 * INVLATTICE_A;
|
||||||
|
|
||||||
|
|
||||||
|
char *omegafile = argv[1];
|
||||||
|
char *kfile = argv[2];
|
||||||
|
char *outfile = argv[3];
|
||||||
|
char *outlongfile = argv[4];
|
||||||
|
char *outshortfile = argv[5];
|
||||||
|
double scuffomegas[MAXOMEGACOUNT];
|
||||||
|
cart2_t klist[MAXKCOUNT];
|
||||||
|
FILE *f = fopen(omegafile, "r");
|
||||||
|
int omegacount = 0;
|
||||||
|
while (fscanf(f, "%lf", scuffomegas + omegacount) == 1){
|
||||||
|
assert(omegacount < MAXOMEGACOUNT);
|
||||||
|
++omegacount;
|
||||||
|
}
|
||||||
|
fclose(f);
|
||||||
|
f = fopen(kfile, "r");
|
||||||
|
int kcount = 0;
|
||||||
|
while (fscanf(f, "%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
|
||||||
|
assert(kcount < MAXKCOUNT);
|
||||||
|
++kcount;
|
||||||
|
}
|
||||||
|
fclose(f);
|
||||||
|
|
||||||
|
const double refindex = REFINDEX;
|
||||||
|
const double h = LATTICE_H;
|
||||||
|
const double a = h * s3;
|
||||||
|
const double rec_a = 4*M_PI/s3/a;
|
||||||
|
|
||||||
|
// generation of the real-space lattices
|
||||||
|
const cart2_t cart2_0 = {0, 0};
|
||||||
|
const cart2_t ABoffset = {h, 0};
|
||||||
|
const cart2_t BAoffset = {-h, 0};
|
||||||
|
//const cart2_t ab_particle_offsets[2][2] = {{ {0, 0}, {h, 0} }, {-h, 0}, {0, 0}};
|
||||||
|
|
||||||
|
// THIS IS THE LATTICE OF r_b
|
||||||
|
latticepoints_circle_t lattice_0offset = generate_tripoints_ver(a, MAXR_REAL, cart2_0);
|
||||||
|
// these have to have the same point order, therefore we must make the offset verision manually to avoid sorting;
|
||||||
|
latticepoints_circle_t lattice_ABoffset, lattice_BAoffset;
|
||||||
|
lattice_ABoffset.points = malloc(lattice_0offset.npoints * sizeof(cart2_t));
|
||||||
|
lattice_ABoffset.capacity = lattice_0offset.npoints * sizeof(cart2_t);
|
||||||
|
lattice_ABoffset.npoints = lattice_ABoffset.capacity;
|
||||||
|
lattice_BAoffset.points = malloc(lattice_0offset.npoints * sizeof(cart2_t));
|
||||||
|
lattice_BAoffset.capacity = lattice_0offset.npoints * sizeof(cart2_t);
|
||||||
|
lattice_BAoffset.npoints = lattice_BAoffset.capacity;
|
||||||
|
for (int i = 0; i < lattice_0offset.npoints; ++i) {
|
||||||
|
lattice_ABoffset.points[i] = cart2_add(lattice_0offset.points[i], ABoffset);
|
||||||
|
lattice_BAoffset.points[i] = cart2_add(lattice_0offset.points[i], BAoffset);
|
||||||
|
}
|
||||||
|
|
||||||
|
// reciprocal lattice, without offset – DON'T I NEED REFINDEX HERE? (I DON'T THINK SO.)
|
||||||
|
latticepoints_circle_t reclattice = generate_tripoints_hor(rec_a, MAXR_K, cart2_0);
|
||||||
|
|
||||||
|
qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS);
|
||||||
|
|
||||||
|
FILE *out = fopen(outfile, "w");
|
||||||
|
FILE *outlong = fopen(outlongfile, "w");
|
||||||
|
FILE *outshort = fopen(outshortfile, "w");
|
||||||
|
|
||||||
|
// as in eq. (5) in my notes
|
||||||
|
double WL_prefactor = 4*M_PI/(a*a)/s3 / /*??*/ (4*M_PI*M_PI);
|
||||||
|
|
||||||
|
for (int omegai = 0; omegai < omegacount; ++omegai) {
|
||||||
|
double scuffomega = scuffomegas[omegai];
|
||||||
|
double omega = scuffomega * SCUFF_OMEGAUNIT;
|
||||||
|
double EeV = omega * hbar / eV;
|
||||||
|
double k0_vac = omega / c0;
|
||||||
|
double k0_eff = k0_vac * refindex; // this one will be used with the real x geometries
|
||||||
|
double cv = CC * k0_eff;
|
||||||
|
|
||||||
|
complex double Abuf[c->nelem][c->nelem], Bbuf[c->nelem][c->nelem];
|
||||||
|
// indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
|
||||||
|
complex double WS[2][2][2][c->nelem][c->nelem];
|
||||||
|
complex double WS_comp[2][2][2][c->nelem][c->nelem];
|
||||||
|
complex double WL[2][2][2][c->nelem][c->nelem];
|
||||||
|
complex double WL_comp[2][2][2][c->nelem][c->nelem];
|
||||||
|
|
||||||
|
for (int ki = 0; ki < kcount; ++ki) {
|
||||||
|
cart2_t k = klist[ki];
|
||||||
|
memset(WS, 0, sizeof(WS));
|
||||||
|
memset(WS_comp, 0, sizeof(WS_comp));
|
||||||
|
memset(WL, 0, sizeof(WL));
|
||||||
|
memset(WL_comp, 0, sizeof(WL_comp));
|
||||||
|
|
||||||
|
for (int bi = 0; bi < lattice_0offset.npoints; ++bi) {
|
||||||
|
cart2_t point0 = lattice_0offset.points[bi];
|
||||||
|
double phase = cart2_dot(k,point0);
|
||||||
|
complex double phasefac = cexp(I*phase);
|
||||||
|
|
||||||
|
if (point0.x || point0.y) { // skip the singular point
|
||||||
|
qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
|
||||||
|
cart22sph(cart2_scale(k0_eff,lattice_0offset.points[bi])), 3, 2, 5, CC);
|
||||||
|
for (int desty = 0; desty < c->nelem; ++desty)
|
||||||
|
for (int srcy = 0; srcy < c->nelem; ++srcy) {
|
||||||
|
ckahanadd(&(WS[0][0][0][desty][srcy]),&(WS_comp[0][0][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
|
||||||
|
ckahanadd(&(WS[0][0][1][desty][srcy]),&(WS_comp[0][0][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
|
||||||
|
cart22sph(cart2_scale(k0_eff,lattice_ABoffset.points[bi])), 3, 2, 5, CC);
|
||||||
|
for (int desty = 0; desty < c->nelem; ++desty)
|
||||||
|
for (int srcy = 0; srcy < c->nelem; ++srcy) {
|
||||||
|
ckahanadd(&(WS[0][1][0][desty][srcy]),&(WS_comp[0][1][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
|
||||||
|
ckahanadd(&(WS[0][1][1][desty][srcy]),&(WS_comp[0][1][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
|
||||||
|
}
|
||||||
|
|
||||||
|
qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
|
||||||
|
cart22sph(cart2_scale(k0_eff,lattice_BAoffset.points[bi])), 3, 2, 5, CC);
|
||||||
|
for (int desty = 0; desty < c->nelem; ++desty)
|
||||||
|
for (int srcy = 0; srcy < c->nelem; ++srcy) {
|
||||||
|
ckahanadd(&(WS[1][0][0][desty][srcy]),&(WS_comp[1][0][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
|
||||||
|
ckahanadd(&(WS[1][0][1][desty][srcy]),&(WS_comp[1][0][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
|
||||||
|
}
|
||||||
|
// WS[1][1] is the same as WS[0][0], so copy in the end rather than double-summing
|
||||||
|
}
|
||||||
|
for (int desty = 0; desty < c->nelem; ++desty)
|
||||||
|
for (int srcy = 0; srcy < c->nelem; ++srcy)
|
||||||
|
for (int ctype = 0; ctype < 2; ctype++)
|
||||||
|
WS[1][1][ctype][desty][srcy] = WS[0][0][ctype][desty][srcy];
|
||||||
|
// WS DONE
|
||||||
|
for (int Ki = 0; Ki < reclattice.npoints; ++Ki) {
|
||||||
|
cart2_t K = reclattice.points[Ki];
|
||||||
|
cart2_t k_K = cart2_substract(k, K);
|
||||||
|
double phase_AB =
|
||||||
|
#ifdef SWAPSIGN1
|
||||||
|
-
|
||||||
|
#endif
|
||||||
|
cart2_dot(k_K, ABoffset); // And maybe the sign is excactly opposite!!! FIXME TODO CHECK
|
||||||
|
complex double phasefacs[2][2];
|
||||||
|
phasefacs[0][0] = phasefacs[1][1] = 1;
|
||||||
|
phasefacs[1][0] = cexp(I * phase_AB); // sign???
|
||||||
|
phasefacs[0][1] = cexp(- I * phase_AB); // sign???
|
||||||
|
|
||||||
|
// FIXME should I skip something (such as the origin?)
|
||||||
|
qpms_trans_calculator_get_2DFT_longrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
|
||||||
|
cart22sph(k_K), 3, 2, 5, cv, k0_eff);
|
||||||
|
for (int dp = 0; dp < 2; dp++)
|
||||||
|
for (int sp = 0; sp < 2; sp++)
|
||||||
|
for (int dy = 0; dy < c->nelem; dy++)
|
||||||
|
for (int sy = 0; sy < c->nelem; sy++) {
|
||||||
|
ckahanadd(&(WL[dp][sp][0][dy][sy]), &(WL_comp[dp][sp][0][dy][sy]), phasefacs[dp][sp] * Abuf[dy][sy] * WL_prefactor);
|
||||||
|
ckahanadd(&(WL[dp][sp][1][dy][sy]), &(WL_comp[dp][sp][1][dy][sy]), phasefacs[dp][sp] * Bbuf[dy][sy] * WL_prefactor);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fprintf(outshort, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
|
||||||
|
scuffomega, EeV, k0_eff, k.x, k.y);
|
||||||
|
fprintf(outlong, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
|
||||||
|
scuffomega, EeV, k0_eff, k.x, k.y);
|
||||||
|
fprintf(out, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
|
||||||
|
scuffomega, EeV, k0_eff, k.x, k.y);
|
||||||
|
size_t totalelems = sizeof(WL) / sizeof(complex double);
|
||||||
|
for (int i = 0; i < totalelems; ++i) {
|
||||||
|
complex double ws = ((complex double *)WS)[i];
|
||||||
|
complex double wl = ((complex double *)WL)[i];
|
||||||
|
complex double w = ws+wl;
|
||||||
|
fprintf(outshort, "%.16g\t%.16g\t", creal(ws), cimag(ws));
|
||||||
|
fprintf(outlong, "%.16g\t%.16g\t", creal(wl), cimag(wl));
|
||||||
|
fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
|
||||||
|
}
|
||||||
|
fputc('\n', outshort);
|
||||||
|
fputc('\n', outlong);
|
||||||
|
fputc('\n', out);
|
||||||
|
}
|
||||||
|
}
|
||||||
|
fclose(out);
|
||||||
|
fclose(outlong);
|
||||||
|
fclose(outshort);
|
||||||
|
|
||||||
|
}
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
|
||||||
|
#if 0
|
||||||
|
int main (int argc, char **argv) {
|
||||||
cart2_t offset = {0,0};
|
cart2_t offset = {0,0};
|
||||||
|
|
||||||
latticepoints_circle_t lat = generate_tripoints_ver(h, 30, offset);
|
latticepoints_circle_t lat = generate_tripoints_ver(1, 200, offset);
|
||||||
for (int i = 0; i < lat.npoints; ++i)
|
for (int i = 0; i < lat.npoints; ++i)
|
||||||
printf("%g %g %g\n", lat.points[i].x, lat.points[i].y, cart2norm(lat.points[i]));
|
printf("%g %g %g\n", lat.points[i].x, lat.points[i].y, cart2norm(lat.points[i]));
|
||||||
latticepoints_circle_free(&lat);
|
latticepoints_circle_free(&lat);
|
||||||
}
|
}
|
||||||
|
#endif
|
||||||
|
|
|
||||||
|
|
@ -1422,7 +1422,7 @@ int qpms_trans_calculator_get_2DFT_longrange_AB_arrays_buf(const qpms_trans_calc
|
||||||
}
|
}
|
||||||
|
|
||||||
|
|
||||||
int qpms_trans_calculator_get_Fourier_longrange_AB_arrays(const qpms_trans_calculator *c,
|
int qpms_trans_calculator_get_2DFT_longrange_AB_arrays(const qpms_trans_calculator *c,
|
||||||
complex double *Adest, complex double *Bdest,
|
complex double *Adest, complex double *Bdest,
|
||||||
size_t deststride, size_t srcstride,
|
size_t deststride, size_t srcstride,
|
||||||
sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */,
|
sph_t k_sph, qpms_bessel_t J /* Only J=3 valid for now */,
|
||||||
|
|
|
||||||
Loading…
Reference in New Issue