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

Test of the new api – results disagree with the old one. 

Former-commit-id: ed09750246bca9a71e810745f443612ea0b989e8
This commit is contained in:
Marek Nečada 2018-12-11 11:36:12 +00:00
parent 3972f0a2e3
commit 4695792772
3 changed files with 81 additions and 46 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
qpms/ewald.c
View File

@ -902,7 +902,8 @@ int ewald3_sigma_short(
break; break;
case LAT_2D_IN_3D_XYONLY: case LAT_2D_IN_3D_XYONLY:
assert((pgen_retdata.flags &PGEN_AT_XY) == PGEN_AT_XY); assert((pgen_retdata.flags &PGEN_AT_XY) == PGEN_AT_XY);
assert(Rpq_shifted_theta == M_PI_2); assert(fabs(Rpq_shifted_theta - M_PI_2) < DBL_EPSILON * 1024);
// assert(Rpq_shifted_theta == M_PI_2); // FIXME this should work as well
legendre_array = c->legendre0; legendre_array = c->legendre0;
break; break;
default: default:

2
tests/ewaldshift2.c
View File

@ -1,4 +1,4 @@
// c99 -ggdb -Wall -I ../ ewaldshift2.c ../qpms/ewald.c ../qpms/ewaldsf.c ../qpms/lattices2d.c -lgsl -lm -lblas // c99 -o ewaldshift2 -ggdb -Wall -I ../ ewaldshift2.c ../qpms/ewald.c ../qpms/ewaldsf.c ../qpms/lattices2d.c -lgsl -lm -lblas
// implementation of the [LT(4.16)] test // implementation of the [LT(4.16)] test
#include <math.h> #include <math.h>

118
tests/ewaldshift_3g.c
View File

@ -1,4 +1,4 @@
// c99 -ggdb -Wall -I ../ ewaldshift2.c ../qpms/ewald.c ../qpms/ewaldsf.c ../qpms/lattices2d.c -lgsl -lm -lblas // c99 -o ewaldshift_3g -ggdb -Wall -I ../ ewaldshift_3g.c ../qpms/ewald.c ../qpms/ewaldsf.c ../qpms/lattices2d.c ../qpms/latticegens.c -lgsl -lm -lblas
// implementation of the [LT(4.16)] test // implementation of the [LT(4.16)] test
#include <math.h> #include <math.h>
@ -188,29 +188,17 @@ void ewaldtest_triang_results_free(ewaldtest_triang_results *r) {
free(r); free(r);
} }
void dump_points2d_rordered(const points2d_rordered_t *ps, char *filename) {
FILE *f = fopen(filename, "w");
for (size_t i = 0; i < ps->nrs; ++i) {
fprintf(f, "# r = %.16g\n", ps->rs[i]);
for (ptrdiff_t j = ps->r_offsets[i]; j < ps->r_offsets[i+1]; ++j)
fprintf(f, "%.16g %.16g\n", ps->base[j].x, ps->base[j].y);
}
fclose(f);
}
static inline double san(double x) { static inline double san(double x) {
return fabs(x) < 1e-13 ? 0 : x; return fabs(x) < 1e-13 ? 0 : x;
} }
ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p); ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p);
int main() { int main() {
gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler); gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) { for (size_t i = 0; i < sizeof(paramslist)/sizeof(ewaldtest_triang_params); ++i) {
ewaldtest_triang_params p = paramslist[i]; ewaldtest_triang_params p = paramslist[i];
ewaldtest_triang_results *r = ewaldtest_triang(p); ewaldtest_triang_results *r = ewaldtest_triang_3g(p);
// TODO print per-test header here // TODO print per-test header here
printf("===============================\n"); printf("===============================\n");
printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), ps = (%g,%g), csphase = %g\n", printf("a = %g, K = %g, Kmax = %g, Rmax = %g, lMax = %d, eta = %g, k = %g, beta = (%g,%g), ps = (%g,%g), csphase = %g\n",
@ -246,36 +234,51 @@ int main() {
int ewaldtest_counter = 0; int ewaldtest_counter = 0;
ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) {
const double a = p.a; //const double a = p.h * sqrt(3); const double a = p.a; //const double a = p.h * sqrt(3);
cart3_t beta3 = cart22cart3xy(p.beta);
cart3_t particle_shift3 = cart22cart3xy(p.particle_shift);
const double A = sqrt(3) * a * a / 2.; // unit cell size cart2_t b1, b2, rb1, rb2;
const double K_len = 4*M_PI/a/sqrt(3); // reciprocal vector length if (p.orientation == TRIANGULAR_VERTICAL) {
b1.x = 0;
b1.y = a;
b2.x = a * M_SQRT3 * .5;
b2.y = a * .5;
} else {
b1.x = a;
b1.y = 0;
b2.x = a * .5;
b2.y = a * M_SQRT3 * .5;
}
if (QPMS_SUCCESS != l2d_reciprocalBasis2pi(b1, b2, &rb1, &rb2))
abort();
const double A = l2d_unitcell_area(b1, b2); // sqrt(3) * a * a / 2.; // unit cell size
const double K_len = cart2norm(rb1); //4*M_PI/a/sqrt(3); // reciprocal vector length
ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results)); ewaldtest_triang_results *results = malloc(sizeof(ewaldtest_triang_results));
results->p = p; results->p = p;
triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, p.orientation, true, 0); // N.B. orig is included //triangular_lattice_gen_t *Rlg = triangular_lattice_gen_init(a, p.orientation, true, 0); // N.B. orig is included
triangular_lattice_gen_extend_to_r(Rlg, p.maxR + a); //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_t *Klg = triangular_lattice_gen_init(K_len, reverseTriangularLatticeOrientation(p.orientation), true, 0);
triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len); //triangular_lattice_gen_extend_to_r(Klg, p.maxK + K_len);
point2d *Rpoints = Rlg->ps.base; //point2d *Rpoints = Rlg->ps.base;
size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs]; //size_t nR = Rlg->ps.r_offsets[Rlg->ps.nrs];
if (fabs(p.particle_shift.x) ==0 && fabs(p.particle_shift.y) == 0) { /*if (fabs(p.particle_shift.x) ==0 && fabs(p.particle_shift.y) == 0) {
points2d_rordered_t Rpos = points2d_rordered_annulus(&(Rlg->ps), 0, false, INFINITY, false); points2d_rordered_t Rpos = points2d_rordered_annulus(&(Rlg->ps), 0, false, INFINITY, false);
Rpoints = Rpos.base + Rpos.r_offsets[0]; Rpoints = Rpos.base + Rpos.r_offsets[0];
nR = Rpos.r_offsets[Rpos.nrs] - Rpos.r_offsets[0]; nR = Rpos.r_offsets[Rpos.nrs] - Rpos.r_offsets[0];
} }*/
//point2d *Kpoints = Klg->ps.base;
//size_t nK = Klg->ps.r_offsets[Klg->ps.nrs];
/*
point2d *Kpoints = Klg->ps.base;
size_t nK = Klg->ps.r_offsets[Klg->ps.nrs];
point2d particle_shift = p.particle_shift; point2d particle_shift = p.particle_shift;
point2d minus_ps = {-particle_shift.x, -particle_shift.y}; point2d minus_ps = {-particle_shift.x, -particle_shift.y};
point2d Rpoints_plus_shift[nR]; point2d Rpoints_plus_shift[nR];
@ -283,8 +286,17 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
Rpoints_plus_shift[i].x = Rpoints[i].x - particle_shift.x; Rpoints_plus_shift[i].x = Rpoints[i].x - particle_shift.x;
Rpoints_plus_shift[i].y = Rpoints[i].y - particle_shift.y; Rpoints_plus_shift[i].y = Rpoints[i].y - particle_shift.y;
} }
*/
// skip zeroth point if it coincides with origin
bool include_origin = !(fabs(p.particle_shift.x) == 0
&& fabs(p.particle_shift.y) == 0);
PGen Rlgen = PGen_xyWeb_new(b1, b2, BASIS_RTOL, CART2_ZERO, 0, include_origin, p.maxR + a, false);
//PGen Rlgen_plus_shift = PGen_xyWeb_new(b1, b2, BASIS_RTOL, cart2_scale(-1 /* CHECKSIGN */, particle_shift2), 0, include_origin, p.maxR + a, false);
PGen Klgen = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, CART2_ZERO, 0, true, p.maxK + K_len, false);
//PGen Klgen_plus_beta = PGen_xyWeb_new(rb1, rb2, BASIS_RTOL, beta2, 0, true, p.maxK + K_len, false);
qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax); qpms_y_t nelem_sc = qpms_lMax2nelem_sc(p.lMax);
@ -297,26 +309,37 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
qpms_ewald32_constants_t *c = qpms_ewald32_constants_init(p.lMax, p.csphase); 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, //points2d_rordered_t *Kpoints_plus_beta = points2d_rordered_shift(&(Klg->ps), p.beta,
8*DBL_EPSILON, 8*DBL_EPSILON); // 8*DBL_EPSILON, 8*DBL_EPSILON);
char filename[BUFSIZ]; //char filename[BUFSIZ];
sprintf(filename, "betalattice_%d.out", ewaldtest_counter); //sprintf(filename, "betalattice_%d.out", ewaldtest_counter);
dump_points2d_rordered(Kpoints_plus_beta, filename); //dump_points2d_rordered(Kpoints_plus_beta, filename);
if (0!=ewald3_sigma_long(results->sigmas_long,
results->err_sigmas_long, c, p.eta, p.k, A,
LAT_2D_IN_3D_XYONLY, &Klgen, false, beta3, particle_shift3))
abort();
#if 0
if (0!=ewald32_sigma_long_points_and_shift(results->sigmas_long, if (0!=ewald32_sigma_long_points_and_shift(results->sigmas_long,
results->err_sigmas_long, c, p.eta, p.k, A, results->err_sigmas_long, c, p.eta, p.k, A,
nK, Kpoints, nK, Kpoints,
p.beta, p.beta,
particle_shift /*minus_ps*/ )) particle_shift /*minus_ps*/ ))
abort(); abort();
if (0!=ewald32_sigma_short_points_and_shift( #endif
if (0!=ewald3_sigma_short(
results->sigmas_short, results->err_sigmas_short, c,
p.eta, p.k, LAT_2D_IN_3D_XYONLY, &Rlgen, false, beta3, particle_shift3))
abort();
/*if (0!=ewald32_sigma_short_points_and_shift(
results->sigmas_short, results->err_sigmas_short, c, results->sigmas_short, results->err_sigmas_short, c,
p.eta, p.k, p.eta, p.k,
nR, Rpoints, p.beta, particle_shift)) nR, Rpoints, p.beta, particle_shift))
abort(); abort();*/
if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k)) //if (0!=ewald32_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
if (0!=ewald3_sigma0(&(results->sigma0), &(results->err_sigma0), c, p.eta, p.k))
abort(); abort();
for(qpms_y_t y = 0; y < nelem_sc; ++y) { for(qpms_y_t y = 0; y < nelem_sc; ++y) {
results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y]; results->sigmas_total[y] = results->sigmas_short[y] + results->sigmas_long[y];
@ -329,6 +352,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double)); results->regsigmas_416 = calloc(nelem_sc, sizeof(complex double));
results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)]; results->regsigmas_416[0] = -2 * c->legendre0[gsl_sf_legendre_array_index(0,0)];
#if 0 // not yet implemented for the new API
{ {
double legendres[gsl_sf_legendre_array_n(p.lMax)]; double legendres[gsl_sf_legendre_array_n(p.lMax)];
points2d_rordered_t sel = points2d_rordered_t sel =
@ -359,11 +383,21 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) {
} }
} }
} }
#else
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;
results->regsigmas_416[y] = NAN;
}
#endif
points2d_rordered_free(Kpoints_plus_beta);
//points2d_rordered_free(Kpoints_plus_beta);
qpms_ewald32_constants_free(c); qpms_ewald32_constants_free(c);
triangular_lattice_gen_free(Klg); //triangular_lattice_gen_free(Klg);
triangular_lattice_gen_free(Rlg); //triangular_lattice_gen_free(Rlg);
++ewaldtest_counter; ++ewaldtest_counter;
return results; return results;
} }