diff --git a/qpms/apps/hexlattice_ewald.c b/qpms/apps/hexlattice_ewald.c index ec0fae9..346bbdd 100644 --- a/qpms/apps/hexlattice_ewald.c +++ b/qpms/apps/hexlattice_ewald.c @@ -1,6 +1,7 @@ // c99 -ggdb -O2 -DLATTICESUMS -I .. hexlattice_ewald.c ../translations.c ../bessels.c ../lrhankel_recspace_dirty.c ../gaunt.c -lm -lgsl -lblas #include #include +#include #include #include #include "kahansum.h" @@ -36,8 +37,8 @@ int cart2_cmpr (const void *p1, const void *p2) { } typedef struct { - int npoints; // number of lattice points. - int capacity; // for how much points memory is allocated + ptrdiff_t npoints; // number of lattice points. + ptrdiff_t capacity; // for how much points memory is allocated double maxR; // circle radius, points <= R cart2_t *points; } latticepoints_circle_t; @@ -58,7 +59,9 @@ latticepoints_circle_t generate_hexpoints_hor(double h, double R, cart2_t offset lat.npoints = 0; int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve) double unitcellS = s3 * 3 / 2 * h * h; // unit cell area - lat.capacity = 5 + 2 * (R + 1) * (R + 1) * M_PI / unitcellS; // should be enough with some random reserve + double flcapacity = 5 + 2 * (R + 5*h) * (R + 5*h) * M_PI / unitcellS; // should be enough with some random reserve + lat.capacity = flcapacity; + lat.points = malloc(lat.capacity *sizeof(cart2_t)); cart2_t BAoffset = {h, 0}; @@ -66,10 +69,11 @@ latticepoints_circle_t generate_hexpoints_hor(double h, double R, cart2_t offset cart2_t a1 = {-1.5*h, s3/2 *h}; cart2_t a2 = {1.5*h, s3/2 *h}; - for (int i1 = -nmax; i1 <= nmax; ++i1) - for (int i2 = -nmax; i2 <= nmax; ++i2) { + for (ptrdiff_t i1 = -nmax; i1 <= nmax; ++i1) + for (ptrdiff_t i2 = -nmax; i2 <= nmax; ++i2) { cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2))); - if (cart2norm(Apoint) <= R) { + if (lat.npoints >= lat.capacity) + printf("%zd %zd %g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, h, unitcellS); if (cart2norm(Apoint) <= R) { assert(lat.npoints < lat.capacity); lat.points[lat.npoints] = Apoint; lat.npoints++; @@ -91,17 +95,20 @@ latticepoints_circle_t generate_tripoints_ver(double a, double R, cart2_t offset lat.maxR = R; lat.npoints = 0; int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve) - double unitcellS = s3 * 3 / 2 * h * h; // unit cell area - lat.capacity = 5 + (R + 1) * (R + 1) * M_PI / unitcellS; // should be enough with some random reserve + double unitcellS = (s3 * 3) / 2 * h * h; // unit cell area + double flcapacity = 5 + (R + 3*a) * (R + 3*a) * M_PI / unitcellS; // should be enough with some random reserve + lat.capacity = flcapacity; // should be enough with some random reserve lat.points = malloc(lat.capacity *sizeof(cart2_t)); cart2_t a1 = {-1.5*h, s3/2 *h}; cart2_t a2 = {1.5*h, s3/2 *h}; - for (int i1 = -nmax; i1 <= nmax; ++i1) - for (int i2 = -nmax; i2 <= nmax; ++i2) { + for (ptrdiff_t i1 = -nmax; i1 <= nmax; ++i1) + for (ptrdiff_t i2 = -nmax; i2 <= nmax; ++i2) { cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2))); if (cart2norm(Apoint) <= R) { + if (lat.npoints >= lat.capacity) + printf("%zd %zd %g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, a, unitcellS); assert(lat.npoints < lat.capacity); lat.points[lat.npoints] = Apoint; lat.npoints++; @@ -118,7 +125,8 @@ latticepoints_circle_t generate_tripoints_hor(double a, double R, cart2_t offset lat.npoints = 0; int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve) double unitcellS = s3 * 3 / 2 * h * h; // unit cell area - lat.capacity = 5 + (R + 1) * (R + 1) * M_PI / unitcellS; // should be enough with some random reserve + double flcapacity = 5 + (R + 3*a) * (R + 3*a) * M_PI / unitcellS; // should be enough with some random reserve + lat.capacity = flcapacity; // should be enough with some random reserve lat.points = malloc(lat.capacity *sizeof(cart2_t)); cart2_t a1 = {s3/2 *h, -1.5*h}; @@ -126,6 +134,8 @@ latticepoints_circle_t generate_tripoints_hor(double a, double R, cart2_t offset for (int i1 = -nmax; i1 <= nmax; ++i1) for (int i2 = -nmax; i2 <= nmax; ++i2) { + if (lat.npoints >= lat.capacity) + printf("%zd %zd %.12g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, a, unitcellS); cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2))); if (cart2norm(Apoint) <= R) { assert(lat.npoints < lat.capacity); @@ -145,7 +155,7 @@ int main (int argc, char **argv) { char *omegafile = argv[1]; - char *kfile = argv[2]; + // char *kfile = argv[2]; // not used char *outfile = argv[3]; char *outlongfile = argv[4]; char *outshortfile = argv[5]; @@ -158,13 +168,20 @@ int main (int argc, char **argv) { ++omegacount; } fclose(f); - f = fopen(kfile, "r"); - int kcount = 0; + /*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 = 1000; + for (int i = 0; i < kcount; ++i) { + 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;