New hexlattice_ewald.c with real ewald sum

Former-commit-id: cb398312fcb1d78f39cc374e3171e216d7604f90

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 <stdio.h>
+#include <math.h>
+#include <qpms/translations.h>
+#include <qpms/lattices.h>
+#include <gsl/gsl_const_mksa.h>
+
+
+
+#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);
+}

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);

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,