oldtests/ewalds.c

   1 // c99 -ggdb -Wall -I ../ ewalds.c ../qpms/ewald.c ../qpms/ewaldsf.c  ../qpms/lattices2d.c -lgsl -lm -lblas
   2
   3 // implementation of the [LT(4.16)] test
   4 #include <math.h>
   5 #define M_SQRTPI 1.7724538509055160272981674833411452
   6 #include <qpms/ewald.h>
   7 #include <qpms/tiny_inlines.h>
   8 #include <qpms/indexing.h>
   9 #include <stdlib.h>
  10 #include <stdio.h>
  11 #include <float.h>
  12 #include <gsl/gsl_sf_legendre.h>
  13 typedef struct ewaldtest_triang_params {
  14   qpms_l_t lMax;
  15   point2d beta;
  16   double k;
  17   double a;
  18   double eta;
  19   double maxR;
  20   double maxK;
  21   double csphase;
  22   TriangularLatticeOrientation orientation;
  23 } ewaldtest_triang_params;
  24
  25 typedef struct ewaldtest_triang_results {
  26   ewaldtest_triang_params p;
  27   complex double *sigmas_short,
  28           *sigmas_long,
  29           sigma0,
  30           *sigmas_total;
  31   double *err_sigmas_short,
  32     *err_sigmas_long,
  33     err_sigma0,
  34     *err_sigmas_total;
  35   complex double *regsigmas_416;
  36 } ewaldtest_triang_results;
  37
  38
  39 ewaldtest_triang_params paramslist[] = {
  40 // lMax, beta, k, a, eta, maxR, maxK, csphase, orientation
  41 #if 0
  42   { 3, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 20, 1., TRIANGULAR_VERTICAL},
  43   { 3, {1.1, 0.23}, 2.3, 0.97, 0.5,  30, 30, -1., TRIANGULAR_VERTICAL},
  44   { 3, {1.1, 0.23}, 2.3, 0.97, 0.9,  30, 30, 1., TRIANGULAR_VERTICAL},
  45   { 3, {1.1, 0.23}, 2.3, 0.97, 1.3,  30, 30, 1., TRIANGULAR_VERTICAL},
  46   { 6, {1.1, 0.23}, 2.3, 0.97, 1.9,  30, 30, 1., TRIANGULAR_VERTICAL},
  47   { 6, {1.1, 0.23}, 2.3, 0.97, 2.5,  30, 30, 1., TRIANGULAR_VERTICAL},
  48   { 6, {1.1, 0.23}, 2.3, 0.97, 3.2,  30, 30, 1., TRIANGULAR_VERTICAL},
  49   { 6, {1.1, 0.23}, 2.3, 0.97, 3.8,  30, 30, 1., TRIANGULAR_VERTICAL},
  50   { 6, {1.1, 0.23}, 2.3, 0.97, 4.5,  30, 30, 1., TRIANGULAR_VERTICAL},
  51   { 6, {1.1, 0.23}, 2.3, 0.97, 4.5,  40, 40, 1., TRIANGULAR_VERTICAL},
  52   { 6, {1.1, 0.23}, 2.3, 0.97, 2.3,  100, 100, 1., TRIANGULAR_VERTICAL},
  53   { 6, {1.1, 0.23}, 2.3, 0.97, 2.9,  100, 100, 1., TRIANGULAR_VERTICAL},
  54 #endif
  55   { 2, {1.1, 0.23}, 2.3, 0.97, 0.5,  20, 160, 1., TRIANGULAR_VERTICAL},
  56   { 2, {-1.1, -0.23}, 2.3, 0.97, 0.5,  20, 160, 1., TRIANGULAR_VERTICAL},
  57   { 2, {0, 1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  58   { 2, {0, -1.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  59   { 2, {-1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  60   { 2, {1.1, 0}, 2.3, 0.97, 5.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  61   { 2, {1.1, 0}, 2.3, 0.97, 5.5, 10, 80, 1., TRIANGULAR_VERTICAL},
  62   { 2, {1.1, 0}, 2.3, 0.97, 5.5, 5, 40, 1., TRIANGULAR_VERTICAL},
  63   { 2, {1.1, 0}, 2.3, 0.97, 5.5, 2, 16, 1., TRIANGULAR_VERTICAL},
  64   { 2, {1.1, 0}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  65   { 2, {1.1, 0}, 2.3, 0.97, 2.5, 10, 80, 1., TRIANGULAR_VERTICAL},
  66   { 2, {1.1, 0}, 2.3, 0.97, 2.5, 5, 40, 1., TRIANGULAR_VERTICAL},
  67   { 2, {1.1, 0}, 2.3, 0.97, 2.5, 2, 16, 1., TRIANGULAR_VERTICAL},
  68   { 2, {1.1, 0}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  69   { 2, {1.1, 0}, 2.3, 0.97, 0.5, 10, 80, 1., TRIANGULAR_VERTICAL},
  70   { 2, {1.1, 0}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL},
  71   { 2, {1.1, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL},
  72   { 2, {0.3, 0}, 2.3, 0.97, 0.5, 2, 16, 1., TRIANGULAR_VERTICAL},
  73   { 2, {2.7, 1}, 2.3, 0.97, 0.5, 5, 40, 1., TRIANGULAR_VERTICAL},
  74   { 2, {2.7, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  75   { 2, {2.7, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  76   { 2, {2.7, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  77   { 2, {2.7, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  78   { 2, {1.1, 1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  79   { 2, {1.1, 1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  80   { 2, {1.1, 1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  81   { 2, {1.1, 1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  82
  83   { 2, {0, 3.1}, 2.3, 0.97, 0.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  84   { 2, {0, 3.1}, 2.3, 0.97, 1.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  85   { 2, {0, 3.1}, 2.3, 0.97, 2.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  86   { 2, {0, 3.1}, 2.3, 0.97, 3.5, 20, 160, 1., TRIANGULAR_VERTICAL},
  87
  88
  89 // end:
  90 //  { 0,  {0, 0}, 0, 0, 0, 0, 0, 0, 0}
  91 };
  92
  93 void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
  94   free(r->sigmas_short);
  95   free(r->sigmas_long);
  96   free(r->sigmas_total);
  97   free(r->err_sigmas_long);
  98   free(r->err_sigmas_total);
  99   free(r->err_sigmas_short);
 100   free(r->regsigmas_416);
 101   free(r);
 102 }
 103
 104
 105 void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
 106   FILE *f = fopen(filename, "w");
 107   for (size_t i = 0; i < ps->nrs; ++i) {
 108     fprintf(f, "# r = %.16g\n", ps->rs[i]);
 109     for (ptrdiff_t j = ps->r_offsets[i]; j < ps->r_offsets[i+1]; ++j)
 110       fprintf(f, "%.16g %.16g\n", ps->base[j].x, ps->base[j].y);
 111   }
 112   fclose(f);
 113 }
 114
 115
 116 static inline double san(double x) {
 117   return fabs(x) < 1e-13 ? 0 : x;
 118 }
 119
 120 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p);
 121
 122 int main() {
 123   gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
 124   for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) {
 125     ewaldtest_triang_params p = paramslist[i];
 126     ewaldtest_triang_results *r = ewaldtest_triang(p);
 127     // TODO print per-test header here
 128     printf("===============================\n");
 129     printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), csphase = %g\n",
 130         p.a, 4*M_PI/sqrt(3)/p.a, p.maxK, p.maxR, p.lMax, p.eta, p.k, p.beta.x, p.beta.y, p.csphase);
 131     printf("sigma0: %.16g%+.16gj\n", creal(r->sigma0), cimag(r->sigma0));
 132     for (qpms_l_t n = 0; n <= p.lMax; ++n) {
 133       for (qpms_m_t m = -n; m <= n; ++m){
 134         if ((m+n)%2) continue;
 135         qpms_y_t y = qpms_mn2y_sc(m,n);
 136         qpms_y_t y_conj = qpms_mn2y_sc(-m,n);
 137         // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon
 138         printf("%zd %d %d: T:%.16g%+.16gj(%.3g) L:%.16g%+.16gj(%.3g) S:%.16g%+.16gj(%.3g) \n| predict %.16g%+.16gj \n| actual  %.16g%+.16gj\n",
 139             y, n, m, creal(san(r->sigmas_total[y])), san(cimag(r->sigmas_total[y])),
 140             r->err_sigmas_total[y],
 141             san(creal(r->sigmas_long[y])), san(cimag(r->sigmas_long[y])),
 142             r->err_sigmas_long[y],
 143             san(creal(r->sigmas_short[y])), san(cimag(r->sigmas_short[y])),
 144             r->err_sigmas_short[y],
 145             san(creal(r->regsigmas_416[y])), san(cimag(r->regsigmas_416[y])),
 146             san(creal(r->sigmas_total[y]) + creal(r->sigmas_total[y_conj])),
 147             san(cimag(r->sigmas_total[y]) - cimag(r->sigmas_total[y_conj]))
 148         );
 149       }
 150     }
 151     ewaldtest_triang_results_free(r);
 152   }
 153   return 0;
 154 }
 155
 156
 157 int ewaldtest_counter = 0;
 158
 159
 160 ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
 161   const double a = p.a; //const double a = p.h * sqrt(3);
 162
 163   const double A = sqrt(3) * a * a / 2.; // unit cell size
 164   const double K_len = 4*M_PI/a/sqrt(3); // reciprocal vector length
 165
 166
 167   ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results));
 168   results->p = p;
 169
 170   triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, p.orientation, false, 0); // N.B. orig is not included (not directly usable for the honeycomb lattice)
 171   triangular_lattice_gen_extend_to_r(Rlg, p.maxR + a);
 172   triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0);
 173   triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len);
 174
 175   point2d *Rpoints = Rlg->ps.base; //point2d  *Kpoints = Klg->ps.base;
 176   size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs],
 177          nK = Klg->ps.r_offsets[Klg->ps.nrs];
 178
 179   qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax);
 180
 181   results->sigmas_short = malloc(sizeof(complex double)*nelem_sc);
 182   results->sigmas_long = malloc(sizeof(complex double)*nelem_sc);
 183   results->sigmas_total = malloc(sizeof(complex double)*nelem_sc);
 184   results->err_sigmas_short = malloc(sizeof(double)*nelem_sc);
 185   results->err_sigmas_long = malloc(sizeof(double)*nelem_sc);
 186   results->err_sigmas_total = malloc(sizeof(double)*nelem_sc);
 187
 188   qpms_ewald3_constants_t *c = qpms_ewald3_constants_init(p.lMax, p.csphase);
 189
 190   points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
 191       8*DBL_EPSILON, 8*DBL_EPSILON);
 192
 193   char filename[BUFSIZ];
 194   sprintf(filename, "betalattice_%d.out", ewaldtest_counter);
 195   dump_points2d_rordered(Kpoints_plus_beta, filename);
 196
 197   point2d particle_shift = {0,0}; // TODO make this a parameter
 198
 199   if (0!=ewald32_sigma_long_shiftedpoints(results->sigmas_long,
 200         results->err_sigmas_long, c, p.eta, p.k, A,
 201         nK, Kpoints_plus_beta->base,
 202         p.beta,
 203         particle_shift))
 204     abort();
 205   if (0!=ewald32_sigma_short_shiftedpoints(
 206         results->sigmas_short, results->err_sigmas_short, c,
 207         p.eta, p.k,
 208         nR, Rpoints, p.beta, particle_shift))
 209     abort();
 210   if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
 211     abort();
 212   for(qpms_y_t y = 0; y < nelem_sc; ++y) {
 213     results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y];
 214     results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y];
 215   }
 216   results->sigmas_total[0] += results->sigma0;
 217   results->err_sigmas_total[0] += results->err_sigma0;
 218
 219   // Now calculate the reference values [LT(4.16)]
 220   results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double));
 221   results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)];
 222
 223   {
 224     double legendres[gsl_sf_legendre_array_n(p.lMax)];
 225     points2d_rordered_t sel =
 226       points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false);
 227     if (0 != sel.nrs)
 228     {
 229       point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0];
 230       size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0];
 231       for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) {
 232         point2d beta_pq = beta_pq_lessthan_k[i];
 233         double rbeta_pq = cart2norm(beta_pq);
 234         double arg_pq = atan2(beta_pq.y, beta_pq.x);
 235         double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq);
 236         if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE,
 237               p.lMax, denom/p.k, p.csphase, legendres) != 0)
 238           abort();
 239         for (qpms_y_t y = 0; y < nelem_sc; ++y) {
 240           qpms_l_t n; qpms_m_t m;
 241           qpms_y2mn_sc_p(y, &m, &n);
 242           if ((m+n)%2 != 0)
 243             continue;
 244           complex double eimf = cexp(I*m*arg_pq);
 245           results->regsigmas_416[y] +=
 246             4*M_PI*ipow(n)/p.k/A
 247             * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m)
 248             / denom;
 249         }
 250       }
 251     }
 252   }
 253
 254   points2d_rordered_free(Kpoints_plus_beta);
 255   qpms_ewald3_constants_free(c);
 256   triangular_lattice_gen_free(Klg);
 257   triangular_lattice_gen_free(Rlg);
 258   ++ewaldtest_counter;
 259   return results;
 260 }