QPMS
/
qpms
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity

hexlattice_ewald: už to něco robí 

Former-commit-id: c55716f7ac150413bd0c5891eba4995ea521fffe
This commit is contained in:
Marek Nečada 2018-05-15 11:03:33 +03:00
parent 81803d16a4
commit 65e8271feb
1 changed files with 31 additions and 14 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

45
qpms/apps/hexlattice_ewald.c
View File

@ -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 <stdlib.h>
#include <string.h>
#include <stddef.h>
#include <math.h>
#include <stdio.h>
#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;