qpms/apps/hexlattice_wrongewald_altin.c

   1 // c99 -ggdb -O2 -DLATTICESUMS -I .. hexlattice_ewald.c ../translations.c ../bessels.c ../lrhankel_recspace_dirty.c ../gaunt.c -lm  -lgsl -lblas
   2 #include <stdlib.h>
   3 #include <string.h>
   4 #include <stddef.h>
   5 #include <math.h>
   6 #include <stdio.h>
   7 #include <string.h>
   8 #include "kahansum.h"
   9 #include "vectors.h"
  10 #include <gsl/gsl_const_mksa.h>
  11 #include <gsl/gsl_math.h>
  12 #include "qpms_types.h"
  13 #include "translations.h"
  14
  15 const double s3 = 1.732050807568877293527446341505872366942805253810380628055;
  16
  17 // IMPORTANT: lattice properties here
  18 const qpms_y_t lMax = 2;
  19 const double REFINDEX = 1.52;
  20 const double LATTICE_H = 576e-9;
  21 static const double SCUFF_OMEGAUNIT = 3e14;
  22 static const double hbar = GSL_CONST_MKSA_PLANCKS_CONSTANT_HBAR;
  23 static const double eV = GSL_CONST_MKSA_ELECTRON_CHARGE;
  24 static const double c0 = GSL_CONST_MKSA_SPEED_OF_LIGHT;
  25 static const double CC = 0.1;
  26
  27 // For sorting the points by distance from origin / radius
  28 int cart2_cmpr (const void *p1, const void *p2) {
  29   const  cart2_t *p1t = (const  cart2_t *)p1;
  30   const  cart2_t *p2t = (const  cart2_t *)p2;
  31   double r21 = cart2norm(*p1t);
  32   double r22 = cart2norm(*p2t);
  33   if (r21 < r22) return -1;
  34   else if (r21 > r22) return 1;
  35   else return 0;
  36 }
  37
  38 typedef struct {
  39   ptrdiff_t npoints; // number of lattice points.
  40   ptrdiff_t capacity; // for how much points memory is allocated
  41   double maxR; // circle radius, points <= R
  42    cart2_t *points;
  43 } latticepoints_circle_t;
  44
  45 void sort_cart2points_by_r(cart2_t *points, size_t nmemb) {
  46   qsort(points, nmemb, sizeof(cart2_t), cart2_cmpr);
  47 }
  48
  49 void latticepoints_circle_free(latticepoints_circle_t *c) {
  50   free(c->points);
  51   c->capacity = 0;
  52 }
  53
  54 // "horizontal" orientation of the adjacent A, B points
  55 latticepoints_circle_t generate_hexpoints_hor(double h, double R, cart2_t offset /* if zero, an A is in the origin */) {
  56   latticepoints_circle_t lat;
  57   lat.maxR = R;
  58   lat.npoints = 0;
  59   int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve)
  60   double unitcellS = s3 * 3 / 2 * h * h; // unit cell area
  61   double flcapacity =  5 + 2 * (R + 5*h) * (R + 5*h) * M_PI / unitcellS; // should be enough with some random reserve
  62   lat.capacity = flcapacity;
  63
  64   lat.points = malloc(lat.capacity *sizeof(cart2_t));
  65
  66   cart2_t BAoffset = {h, 0};
  67
  68   cart2_t a1 = {-1.5*h, s3/2 *h};
  69   cart2_t a2 = {1.5*h, s3/2 *h};
  70
  71   for (ptrdiff_t i1 = -nmax; i1 <= nmax; ++i1)
  72     for (ptrdiff_t i2 = -nmax; i2 <= nmax; ++i2) {
  73       cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2)));
  74          if (lat.npoints >= lat.capacity)
  75           printf("%zd %zd %g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, h, unitcellS);     if (cart2norm(Apoint) <= R) {
  76         assert(lat.npoints < lat.capacity);
  77         lat.points[lat.npoints] = Apoint;
  78         lat.npoints++;
  79       }
  80       cart2_t Bpoint = cart2_add(Apoint, BAoffset);
  81       if (cart2norm(Bpoint) <= R) {
  82         assert(lat.npoints < lat.capacity);
  83         lat.points[lat.npoints] = Bpoint;
  84         lat.npoints++;
  85       }
  86     }
  87   sort_cart2points_by_r(lat.points, lat.npoints);
  88   return lat;
  89 }
  90
  91 latticepoints_circle_t generate_tripoints_ver(double a, double R, cart2_t offset /* if zero, an A is in the origin */) {
  92   double h = a / s3;
  93   latticepoints_circle_t lat;
  94   lat.maxR = R;
  95   lat.npoints = 0;
  96   int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve)
  97   double unitcellS = (s3 * 3) / 2 * h * h; // unit cell area
  98   double flcapacity =  5 +  (R + 3*a) * (R + 3*a) * M_PI / unitcellS; // should be enough with some random reserve
  99   lat.capacity = flcapacity; // should be enough with some random reserve
 100   lat.points = malloc(lat.capacity *sizeof(cart2_t));
 101
 102   cart2_t a1 = {-1.5*h, s3/2 *h};
 103   cart2_t a2 = {1.5*h, s3/2 *h};
 104
 105   for (ptrdiff_t i1 = -nmax; i1 <= nmax; ++i1)
 106     for (ptrdiff_t i2 = -nmax; i2 <= nmax; ++i2) {
 107       cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2)));
 108       if (cart2norm(Apoint) <= R) {
 109         if (lat.npoints >= lat.capacity)
 110           printf("%zd %zd %g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, a, unitcellS);
 111         assert(lat.npoints < lat.capacity);
 112         lat.points[lat.npoints] = Apoint;
 113         lat.npoints++;
 114       }
 115     }
 116   sort_cart2points_by_r(lat.points, lat.npoints);
 117   return lat;
 118 }
 119
 120 latticepoints_circle_t generate_tripoints_hor(double a, double R, cart2_t offset /* if zero, an A is in the origin */) {
 121   double h = a / s3;
 122   latticepoints_circle_t lat;
 123   lat.maxR = R;
 124   lat.npoints = 0;
 125   int nmax = R / (1.5 * h) + 2; // max no of lattice shifts in each direction (with reserve)
 126   double unitcellS = s3 * 3 / 2 * h * h; // unit cell area
 127   double flcapacity =  5 +  (R + 3*a) * (R + 3*a) * M_PI / unitcellS; // should be enough with some random reserve
 128   lat.capacity = flcapacity; // should be enough with some random reserve
 129   lat.points = malloc(lat.capacity *sizeof(cart2_t));
 130
 131   cart2_t a1 = {s3/2 *h, -1.5*h};
 132   cart2_t a2 = {s3/2 *h, 1.5 * h};
 133
 134   for (int i1 = -nmax; i1 <= nmax; ++i1)
 135     for (int i2 = -nmax; i2 <= nmax; ++i2) {
 136       if (lat.npoints >= lat.capacity)
 137           printf("%zd %zd %.12g %g %g %g\n", lat.npoints, lat.capacity, flcapacity, R, a, unitcellS);
 138       cart2_t Apoint = cart2_add(offset, cart2_add(cart2_scale(i1, a1), cart2_scale(i2, a2)));
 139       if (cart2norm(Apoint) <= R) {
 140         assert(lat.npoints < lat.capacity);
 141         lat.points[lat.npoints] = Apoint;
 142         lat.npoints++;
 143       }
 144     }
 145   sort_cart2points_by_r(lat.points, lat.npoints);
 146   return lat;
 147 }
 148
 149 int main (int argc, char **argv) {
 150   const double LATTICE_A = s3*LATTICE_H;
 151   const double INVLATTICE_A = 4 * M_PI / s3 / LATTICE_A;
 152   const double MAXR_REAL = 100 * LATTICE_H;
 153   const double MAXR_K = 100 * INVLATTICE_A;
 154
 155   if (argc < 2) abort();
 156
 157   char *omegastr = argv[1];
 158   // char *kfile = argv[2]; // not used;, will be read from stdin
 159   char *outprefix = argv[2];
 160
 161   double scuffomega = strtod(omegastr, NULL);
 162   char *outfile = outprefix;
 163   char outlongfile[strlen(outprefix)+10];
 164   char outshortfile[strlen(outprefix)+10];
 165   sprintf(outlongfile, "%s.long", outprefix);
 166   sprintf(outshortfile, "%s.short", outprefix);
 167
 168
 169
 170   //cart2_t klist[MAXKCOUNT];
 171   /*f = fopen(kfile, "r");
 172   int kcount = 100;
 173   while (fscanf(f, "%lf %lf", &(klist[kcount].x), &(klist[kcount].y)) == 2) {
 174     assert(kcount < MAXKCOUNT);
 175     ++kcount;
 176   }
 177   fclose(f);
 178   */
 179
 180   const double refindex = REFINDEX;
 181   const double h = LATTICE_H;
 182   const double a = h * s3;
 183   const double rec_a = 4*M_PI/s3/a;
 184
 185   // generation of the real-space lattices
 186   const cart2_t cart2_0 = {0, 0};
 187   const cart2_t ABoffset = {h, 0};
 188   const cart2_t BAoffset = {-h, 0};
 189   //const cart2_t ab_particle_offsets[2][2] = {{ {0, 0}, {h, 0} }, {-h, 0}, {0, 0}};
 190
 191   // THIS IS THE LATTICE OF r_b
 192   latticepoints_circle_t lattice_0offset = generate_tripoints_ver(a, MAXR_REAL, cart2_0);
 193   // these have to have the same point order, therefore we must make the offset verision manually to avoid sorting;
 194   latticepoints_circle_t lattice_ABoffset, lattice_BAoffset;
 195   lattice_ABoffset.points = malloc(lattice_0offset.npoints * sizeof(cart2_t));
 196   lattice_ABoffset.capacity = lattice_0offset.npoints * sizeof(cart2_t);
 197   lattice_ABoffset.npoints = lattice_ABoffset.capacity;
 198   lattice_BAoffset.points = malloc(lattice_0offset.npoints * sizeof(cart2_t));
 199   lattice_BAoffset.capacity = lattice_0offset.npoints * sizeof(cart2_t);
 200   lattice_BAoffset.npoints = lattice_BAoffset.capacity;
 201   for (int i = 0; i < lattice_0offset.npoints; ++i) {
 202     lattice_ABoffset.points[i] = cart2_add(lattice_0offset.points[i], ABoffset);
 203     lattice_BAoffset.points[i] = cart2_add(lattice_0offset.points[i], BAoffset);
 204   }
 205
 206   // reciprocal lattice, without offset – DON'T I NEED REFINDEX HERE? (I DON'T THINK SO.)
 207   latticepoints_circle_t reclattice = generate_tripoints_hor(rec_a, MAXR_K, cart2_0);
 208
 209   qpms_trans_calculator *c = qpms_trans_calculator_init(lMax, QPMS_NORMALISATION_POWER_CS);
 210
 211   FILE *out = fopen(outfile, "w");
 212   FILE *outlong = fopen(outlongfile, "w");
 213   FILE *outshort = fopen(outshortfile, "w");
 214
 215   // as in eq. (5) in my notes
 216   double WL_prefactor = 4*M_PI/(a*a)/s3 / /*??*/ (4*M_PI*M_PI);
 217
 218     //double scuffomega = scuffomegas[omegai];
 219     double omega = scuffomega * SCUFF_OMEGAUNIT;
 220     double EeV = omega * hbar / eV;
 221     double k0_vac = omega / c0;
 222     double k0_eff = k0_vac * refindex; // this one will be used with the real x geometries
 223     double cv = CC * k0_eff;
 224
 225     complex double Abuf[c->nelem][c->nelem], Bbuf[c->nelem][c->nelem];
 226     // indices : destpart (A/B-particle), srcpart (A/B-particle), coeff type (A/B- type), desty, srcy
 227     complex double WS[2][2][2][c->nelem][c->nelem];
 228     complex double WS_comp[2][2][2][c->nelem][c->nelem];
 229     complex double WL[2][2][2][c->nelem][c->nelem];
 230     complex double WL_comp[2][2][2][c->nelem][c->nelem];
 231
 232     //for (int ki = 0; ki < kcount; ++ki) {
 233     //  cart2_t k = klist[ki];
 234     cart2_t k;
 235     while(scanf("%lf %lf", &(k.x), &(k.y)) == 2) {
 236       memset(WS, 0, sizeof(WS));
 237       memset(WS_comp, 0, sizeof(WS_comp));
 238       memset(WL, 0, sizeof(WL));
 239       memset(WL_comp, 0, sizeof(WL_comp));
 240
 241       for (int bi = 0; bi < lattice_0offset.npoints; ++bi) {
 242         cart2_t point0 = lattice_0offset.points[bi];
 243         double phase = cart2_dot(k,point0);
 244         complex double phasefac = cexp(I*phase);
 245
 246         if (point0.x || point0.y) { // skip the singular point
 247           qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
 248               cart22sph(cart2_scale(k0_eff,lattice_0offset.points[bi])), 3, 2, 5, CC);
 249           for (int desty = 0; desty < c->nelem; ++desty)
 250             for (int srcy = 0; srcy < c->nelem; ++srcy) {
 251               ckahanadd(&(WS[0][0][0][desty][srcy]),&(WS_comp[0][0][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
 252               ckahanadd(&(WS[0][0][1][desty][srcy]),&(WS_comp[0][0][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
 253             }
 254         }
 255         qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
 256             cart22sph(cart2_scale(k0_eff,lattice_ABoffset.points[bi])), 3, 2, 5, CC);
 257         for (int desty = 0; desty < c->nelem; ++desty)
 258             for (int srcy = 0; srcy < c->nelem; ++srcy) {
 259               ckahanadd(&(WS[0][1][0][desty][srcy]),&(WS_comp[0][1][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
 260               ckahanadd(&(WS[0][1][1][desty][srcy]),&(WS_comp[0][1][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
 261             }
 262
 263         qpms_trans_calculator_get_shortrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
 264             cart22sph(cart2_scale(k0_eff,lattice_BAoffset.points[bi])), 3, 2, 5, CC);
 265         for (int desty = 0; desty < c->nelem; ++desty)
 266             for (int srcy = 0; srcy < c->nelem; ++srcy) {
 267               ckahanadd(&(WS[1][0][0][desty][srcy]),&(WS_comp[1][0][0][desty][srcy]),Abuf[desty][srcy] * phasefac);
 268               ckahanadd(&(WS[1][0][1][desty][srcy]),&(WS_comp[1][0][1][desty][srcy]),Bbuf[desty][srcy] * phasefac);
 269             }
 270         // WS[1][1] is the same as WS[0][0], so copy in the end rather than double-summing
 271       }
 272       for (int desty = 0; desty < c->nelem; ++desty)
 273         for (int srcy = 0; srcy < c->nelem; ++srcy)
 274           for (int ctype = 0; ctype < 2; ctype++)
 275             WS[1][1][ctype][desty][srcy] = WS[0][0][ctype][desty][srcy];
 276       // WS DONE
 277       for (int Ki = 0; Ki < reclattice.npoints; ++Ki) {
 278         cart2_t K = reclattice.points[Ki];
 279         cart2_t k_K = cart2_substract(k, K);
 280         double phase_AB =
 281 #ifdef SWAPSIGN1
 282           -
 283 #endif
 284           cart2_dot(k_K, ABoffset); // And maybe the sign is excactly opposite!!! FIXME TODO CHECK
 285         complex double phasefacs[2][2];
 286         phasefacs[0][0] = phasefacs[1][1] = 1;
 287         phasefacs[1][0] = cexp(I * phase_AB); // sign???
 288         phasefacs[0][1] = cexp(- I * phase_AB); // sign???
 289
 290         // FIXME should I skip something (such as the origin?)
 291         qpms_trans_calculator_get_2DFT_longrange_AB_arrays(c, (complex double *) Abuf, (complex double *) Bbuf, c->nelem, 1,
 292             cart22sph(k_K), 3, 2, 5, cv, k0_eff);
 293         for (int dp = 0; dp < 2; dp++)
 294           for (int sp = 0; sp < 2; sp++)
 295             for (int dy = 0; dy < c->nelem; dy++)
 296               for (int sy = 0; sy < c->nelem; sy++) {
 297                 ckahanadd(&(WL[dp][sp][0][dy][sy]), &(WL_comp[dp][sp][0][dy][sy]), phasefacs[dp][sp] * Abuf[dy][sy] * WL_prefactor);
 298                 ckahanadd(&(WL[dp][sp][1][dy][sy]), &(WL_comp[dp][sp][1][dy][sy]), phasefacs[dp][sp] * Bbuf[dy][sy] * WL_prefactor);
 299               }
 300       }
 301       fprintf(outshort, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 302           scuffomega, EeV, k0_eff, k.x, k.y);
 303       fprintf(outlong, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 304           scuffomega, EeV, k0_eff, k.x, k.y);
 305       fprintf(out, "%.16g\t%.16g\t%16g\t%.16g\t%.16g\t",
 306           scuffomega, EeV, k0_eff, k.x, k.y);
 307       size_t totalelems = sizeof(WL) / sizeof(complex double);
 308       for (int i = 0; i < totalelems; ++i) {
 309         complex double ws = ((complex double *)WS)[i];
 310         complex double wl = ((complex double *)WL)[i];
 311         complex double w = ws+wl;
 312         fprintf(outshort, "%.16g\t%.16g\t", creal(ws), cimag(ws));
 313         fprintf(outlong, "%.16g\t%.16g\t", creal(wl), cimag(wl));
 314         fprintf(out, "%.16g\t%.16g\t", creal(w), cimag(w));
 315       }
 316       fputc('\n', outshort);
 317       fputc('\n', outlong);
 318       fputc('\n', out);
 319       fflush(outshort);
 320       fflush(outlong);
 321       fflush(out);
 322     }
 323   fclose(out);
 324   fclose(outlong);
 325   fclose(outshort);
 326
 327 }
 328
 329
 330
 331
 332 #if 0
 333 int main (int argc, char **argv) {
 334   cart2_t offset = {0,0};
 335
 336   latticepoints_circle_t lat = generate_tripoints_ver(1, 200, offset);
 337   for (int i = 0; i < lat.npoints; ++i)
 338     printf("%g %g %g\n", lat.points[i].x, lat.points[i].y, cart2norm(lat.points[i]));
 339   latticepoints_circle_free(&lat);
 340 }
 341 #endif