From b707f65d83ea4a6b901af85805614d696052cfb6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 5 Sep 2018 18:50:02 +0000 Subject: [PATCH] Ewald sum first test compiles Former-commit-id: bb63ffdade407900478c26f2e10bc0ee8efb5154 --- qpms/ewald.c | 1 + qpms/ewald.h | 4 +- qpms/lattices.h | 15 ++++ qpms/lattices2d.c | 8 +++ tests/ewalds.c | 170 +++++++++++++++++++++++++++++++++++++++++++--- 5 files changed, 187 insertions(+), 11 deletions(-) diff --git a/qpms/ewald.c b/qpms/ewald.c index 17ec200..215e962 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -156,6 +156,7 @@ int ewald32_sigma_long_shiftedpoints ( const qpms_ewald32_constants_t *c, const double eta, const double k, const double unitcell_area, const size_t npoints, const point2d *Kpoints_plus_beta, + //const point2d beta, // not needed const point2d particle_shift // target - src ) { diff --git a/qpms/ewald.h b/qpms/ewald.h index d0ae271..bebebe9 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -87,13 +87,13 @@ int ewald32_sigma0(complex double *result, double *err, // TODO make "compressed versions" where the (m+n)-odd terms (which are zero) // are not included. -int ewald32_sigma_long_shiftedpoints_e ( +int ewald32_sigma_long_shiftedpoints ( complex double *target_sigmalr_y, // must be c->nelem_sc long double *target_sigmalr_y_err, // must be c->nelem_sc long or NULL const qpms_ewald32_constants_t *c, double eta, double k, double unitcell_area, size_t npoints, const point2d *Kpoints_plus_beta, - point2d beta, + //point2d beta, point2d particle_shift ); int ewald32_sigma_long_points_and_shift (//NI diff --git a/qpms/lattices.h b/qpms/lattices.h index d7012d9..7e2d1e3 100644 --- a/qpms/lattices.h +++ b/qpms/lattices.h @@ -16,6 +16,7 @@ #include #include #include +#include #define M_SQRT3 1.7320508075688772935274463415058724 #define M_SQRT3_2 (M_SQRT3/2) #define M_1_SQRT3 0.57735026918962576450914878050195746 @@ -190,6 +191,20 @@ typedef enum { TRIANGULAR_HORIZONTAL // there is a lattice base vector parallel to the x-axis } TriangularLatticeOrientation; +static inline TriangularLatticeOrientation reverseTriangularLatticeOrientation(TriangularLatticeOrientation o){ + switch(o) { + case TRIANGULAR_VERTICAL: + return TRIANGULAR_HORIZONTAL; + break; + case TRIANGULAR_HORIZONTAL: + return TRIANGULAR_VERTICAL; + break; + default: + abort(); + } + abort(); +} + // implementation data structures; not needed in the header file typedef struct triangular_lattice_gen_privstuff_t triangular_lattice_gen_privstuff_t; diff --git a/qpms/lattices2d.c b/qpms/lattices2d.c index 657843e..fd0f4bd 100644 --- a/qpms/lattices2d.c +++ b/qpms/lattices2d.c @@ -505,6 +505,10 @@ const points2d_rordered_t * triangular_lattice_gen_getpoints(const triangular_la return &(g->ps); } +int triangular_lattice_gen_extend_to_r(triangular_lattice_gen_t * g, const double maxr) { + return triangular_lattice_gen_extend_to_steps(g, maxr/g->a); +} + int triangular_lattice_gen_extend_to_steps(triangular_lattice_gen_t * g, int maxsteps) { if (maxsteps <= g->priv->maxs) // nothing needed @@ -605,6 +609,10 @@ void honeycomb_lattice_gen_free(honeycomb_lattice_gen_t *g) { free(g); } +int honeycomb_lattice_gen_extend_to_r(honeycomb_lattice_gen_t *g, double maxr) { + return honeycomb_lattice_gen_extend_to_steps(g, maxr/g->a); /*CHECKME whether g->a is the correct denom.*/ +} + int honeycomb_lattice_gen_extend_to_steps(honeycomb_lattice_gen_t *g, const int maxsteps) { if (maxsteps <= g->tg->priv->maxs) // nothing needed return 0; diff --git a/tests/ewalds.c b/tests/ewalds.c index 818bdf2..f67fc1e 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -1,21 +1,173 @@ // implementation of the [LT(4.16)] test - +#include +#define M_SQRTPI 1.7724538509055160272981674833411452 #include - -typedef struct ewaldtest_hex_params { +#include +#include +#include +#include +#include +#include +typedef struct ewaldtest_triang_params { qpms_l_t lMax; - poinnt2d beta; + point2d beta; double k; - double h; + double a; double eta; -} ewaldtest_hex_params; + double maxR; + double maxK; + double csphase; + TriangularLatticeOrientation orientation; +} ewaldtest_triang_params; + +typedef struct ewaldtest_triang_results { + ewaldtest_triang_params p; + complex double *sigmas_short, + *sigmas_long, + sigma0, + *sigmas_total; + double *err_sigmas_short, + *err_sigmas_long, + err_sigma0, + *err_sigmas_total; + complex double *regsigmas_416; +} ewaldtest_triang_results; -ewaldtest_hex_paraps paramslist = { - { 3, {1.1, 0.23}, 2.3, 0.97, 0.3}, +ewaldtest_triang_params paramslist[] = { + { 3, {1.1, 0.23}, 2.3, 0.97, 0.3, 1., TRIANGULAR_HORIZONTAL}//, // end: - { 0, {0, 0}, 0, 0, 0} +// { 0, {0, 0}, 0, 0, 0, 0, 0} +}; + +void ewaldtest_triang_results_free(ewaldtest_triang_results *r) { + free(r->sigmas_short); + free(r->sigmas_long); + free(r->sigmas_total); + free(r->err_sigmas_long); + free(r->err_sigmas_total); + free(r->err_sigmas_short); + free(r->regsigmas_416); + free(r); +} + +ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p); + +int main() { + for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) { + ewaldtest_triang_params p = paramslist[i]; + ewaldtest_triang_results *r = ewaldtest_triang(p); + // TODO print per-test header here + for (qpms_l_t n = 0; n <= p.lMax; ++n) { + for (qpms_m_t m = -n; m <= n; ++m){ + qpms_y_t y = qpms_mn2y_sc(m,n); + qpms_y_t y_conj = qpms_mn2y_sc(-m,n); + // y n m sigma_total (err), regsigmas_416 regsigmas_415_recon + printf("%zd %d %d: %.16g%+.16gj (%.3g) %.16g%+.16gj %.16g%+.16g\n", + y, n, m, creal(r->sigmas_total[y]), cimag(r->sigmas_total[y]), + r->err_sigmas_total[y], + creal(r->regsigmas_416[y]), cimag(r->regsigmas_416[y]), + creal(r->sigmas_total[y])+creal(r->sigmas_total[y_conj]), + cimag(r->sigmas_total[y])-cimag(r->sigmas_total[y_conj]) + ); + } + } + ewaldtest_triang_results_free(r); + } + return 0; } +ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { + const double a = p.a; //const double a = p.h * sqrt(3); + + const double A = sqrt(3) * a * a / 2.; // unit cell size + const double K_len = 4*M_PI/a/sqrt(3); // reciprocal vector length + + + ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results)); + results->p = p; + 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) + triangular_lattice_gen_extend_to_r(Rlg, p.maxR + a); + triangular_lattice_gen_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0); + triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len); + + point2d *Rpoints = Rlg->ps.base, *Kpoints = Klg->ps.base; + size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs], + nK = Klg->ps.r_offsets[Klg->ps.nrs]; + + qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax); + + results->sigmas_short = malloc(sizeof(complex double)*nelem_sc); + results->sigmas_long = malloc(sizeof(complex double)*nelem_sc); + results->sigmas_total = malloc(sizeof(complex double)*nelem_sc); + results->err_sigmas_short = malloc(sizeof(double)*nelem_sc); + results->err_sigmas_long = malloc(sizeof(double)*nelem_sc); + results->err_sigmas_total = malloc(sizeof(double)*nelem_sc); + + qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase); + + points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta, + 8*DBL_EPSILON, 8*DBL_EPSILON); + + point2d particle_shift = {0,0}; // TODO make this a parameter + + if (0!=ewald32_sigma_long_shiftedpoints(results->sigmas_long, + results->err_sigmas_long, c, p.eta, p.k, A, + nK, Kpoints_plus_beta->base, + //p.beta, + particle_shift)) + abort(); + if (0!=ewald32_sigma_short_shiftedpoints( + results->sigmas_short, results->err_sigmas_short, c, + p.eta, p.k, + nR, Rpoints, p.beta, particle_shift)) + abort(); + if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k)) + abort(); + for(qpms_y_t y = 0; y < nelem_sc; ++y) { + results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y]; + results->err_sigmas_total[y] = results->err_sigmas_short[y] + results->err_sigmas_long[y]; + } + results->sigmas_total[0] += results->sigma0; + results->err_sigmas_total[0] += results->err_sigma0; + + // Now calculate the reference values [LT(4.16)] + results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double)); + results->regsigmas_416[0] = -1/M_SQRTPI; + + { + double legendres[gsl_sf_legendre_array_n(p.lMax)]; + points2d_rordered_t sel = + points2d_rordered_annulus(Kpoints_plus_beta, 0, true, p.k, false); + point2d *beta_pq_lessthan_k = sel.base + sel.r_offsets[0]; + size_t beta_pq_lessthan_k_count = sel.r_offsets[sel.nrs] - sel.r_offsets[0]; + for(size_t i = 0; i < beta_pq_lessthan_k_count; ++i) { + point2d beta_pq = beta_pq_lessthan_k[i]; + double rbeta_pq = cart2norm(beta_pq); + double arg_pq = atan2(beta_pq.y, beta_pq.x); + double denom = sqrt(p.k*p.k - rbeta_pq*rbeta_pq); + if( gsl_sf_legendre_array_e(GSL_SF_LEGENDRE_NONE, + p.lMax, denom/p.k, p.csphase, legendres) != 0) + abort(); + for (qpms_y_t y = 0; y < nelem_sc; ++y) { + qpms_l_t n; qpms_m_t m; + qpms_y2mn_sc_p(y, &m, &n); + if ((m+n)%2 != 0) + continue; + complex double eimf = cexp(I*m*arg_pq); + results->regsigmas_416[y] += + 4*M_PI*ipow(n)/p.k/A + * eimf * legendres[gsl_sf_legendre_array_index(n,m)] + / denom; + } + } + } + + points2d_rordered_free(Kpoints_plus_beta); + qpms_ewald32_constants_free(c); + triangular_lattice_gen_free(Klg); + triangular_lattice_gen_free(Rlg); + return results; +}