From 328d22de89d9bd17aea2e855fd771cdd5b5915ce Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 21 Dec 2018 19:31:07 +0000 Subject: [PATCH] Rename ewald3_constants, move legacy code Former-commit-id: e83dcfa532f7b8d7345103752aca924a56ad7138 --- qpms/ewald.c | 357 +--------------------------------- qpms/ewald.h | 36 ++-- qpms/ewald_legacy.c | 434 ++++++++++++++++++++++++++++++++++++++++++ tests/ewalds.c | 4 +- tests/ewaldshift.c | 4 +- tests/ewaldshift2.c | 4 +- tests/ewaldshift_3g.c | 4 +- tests/ewaldshift_ng.c | 4 +- 8 files changed, 472 insertions(+), 375 deletions(-) create mode 100644 qpms/ewald_legacy.c diff --git a/qpms/ewald.c b/qpms/ewald.c index d4d6a29..46e1c60 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -57,14 +57,14 @@ typedef enum { * spherical harmonic. See notes/ewald.lyx. */ -} ewald32_constants_option; +} ewald3_constants_option; -static const ewald32_constants_option type = EWALD32_CONSTANTS_AGNOSTIC; +static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC; -qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, const ewald32_constants_option type */, +qpms_ewald3_constants_t *qpms_ewald3_constants_init(const qpms_l_t lMax /*, const ewald3_constants_option type */, const int csphase) { - qpms_ewald32_constants_t *c = malloc(sizeof(qpms_ewald32_constants_t)); + qpms_ewald3_constants_t *c = malloc(sizeof(qpms_ewald3_constants_t)); //if (c == NULL) return NULL; // Do I really want to do this? c->lMax = lMax; c->nelem_sc = qpms_lMax2nelem_sc(lMax); @@ -160,7 +160,7 @@ qpms_ewald32_constants_t *qpms_ewald32_constants_init(const qpms_l_t lMax /*, co return c; } -void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) { +void qpms_ewald3_constants_free(qpms_ewald3_constants_t *c) { free(c->legendre0); free(c->legendre_plus1); free(c->legendre_minus1); @@ -175,7 +175,7 @@ void qpms_ewald32_constants_free(qpms_ewald32_constants_t *c) { int ewald3_sigma0(complex double *result, double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k) { qpms_csf_result gam; @@ -189,198 +189,10 @@ int ewald3_sigma0(complex double *result, double *err, return 0; } -int ewald32_sigma0(complex double *result, double *err, - const qpms_ewald32_constants_t *c, - const double eta, const double k) -{ - return ewald3_sigma0(result, err, c, eta, k); -} - - - -int ewald32_sigma_long_shiftedpoints ( - complex double *target, // must be c->nelem_sc long - double *err, - 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 - ) -{ - const qpms_y_t nelem_sc = c->nelem_sc; - const qpms_l_t lMax = c->lMax; - - // Manual init of the ewald summation targets - complex double *target_c = calloc(nelem_sc, sizeof(complex double)); - memset(target, 0, nelem_sc * sizeof(complex double)); - double *err_c = NULL; - if (err) { - err_c = calloc(nelem_sc, sizeof(double)); - memset(err, 0, nelem_sc * sizeof(double)); - } - - const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC) - assert(commonfac > 0); - - // space for Gamma_pq[j]'s - qpms_csf_result Gamma_pq[lMax/2+1]; - - // CHOOSE POINT BEGIN - for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP - const point2d beta_pq = Kpoints_plus_beta[i]; - const point2d K_pq = {beta_pq.x - beta.x, beta_pq.y - beta.y}; - const double rbeta_pq = cart2norm(beta_pq); - // CHOOSE POINT END - - const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!! - const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT - - // R-DEPENDENT BEGIN - const complex double gamma_pq = lilgamma(rbeta_pq/k); - const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné - for(qpms_l_t j = 0; j <= lMax/2; ++j) { - int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j); - // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above - if(creal(z) < 0) - Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above - if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort(); - } - // R-DEPENDENT END - - // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array - // and just fetched for each n, m pair - for(qpms_l_t n = 0; n <= lMax; ++n) - for(qpms_m_t m = -n; m <= n; ++m) { - if((m+n) % 2 != 0) // odd coefficients are zero. - continue; - const qpms_y_t y = qpms_mn2y_sc(m, n); - const complex double e_imalpha_pq = cexp(I*m*arg_pq); - complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c); - double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? - assert((n-abs(m))/2 == c->s1_jMaxes[y]); - for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop - * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation - * c->s1_constfacs[y][j]; - if(err) { - // FIXME include also other errors than Gamma_pq's relative error - kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand)); - } - summand *= Gamma_pq[j].val; // GGG - ckahanadd(&jsum, &jsum_c, summand); - } - jsum *= phasefac; // PFC - ckahanadd(target + y, target_c + y, jsum); - if(err) kahanadd(err + y, err_c + y, jsum_err); - } - } // END POINT LOOP - - free(err_c); - free(target_c); - for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above - target[y] *= commonfac; - if(err) - for(qpms_y_t y = 0; y < nelem_sc; ++y) - err[y] *= commonfac; - return 0; -} - -int ewald32_sigma_long_points_and_shift ( - complex double *target, // must be c->nelem_sc long - double *err, - const qpms_ewald32_constants_t *c, - const double eta, const double k, const double unitcell_area, - const size_t npoints, const point2d *Kpoints, - const point2d beta, - const point2d particle_shift // target - src - ) -{ - const qpms_y_t nelem_sc = c->nelem_sc; - const qpms_l_t lMax = c->lMax; - - // Manual init of the ewald summation targets - complex double *target_c = calloc(nelem_sc, sizeof(complex double)); - memset(target, 0, nelem_sc * sizeof(complex double)); - double *err_c = NULL; - if (err) { - err_c = calloc(nelem_sc, sizeof(double)); - memset(err, 0, nelem_sc * sizeof(double)); - } - - const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC) - assert(commonfac > 0); - - // space for Gamma_pq[j]'s - qpms_csf_result Gamma_pq[lMax/2+1]; - - // CHOOSE POINT BEGIN - for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP - // Only these following two lines differ from ewald32_sigma_long_points_and_shift()!!! WTFCOMMENT?! - const point2d K_pq = Kpoints[i]; - const point2d beta_pq = {K_pq.x + beta.x, K_pq.y + beta.y}; - const double rbeta_pq = cart2norm(beta_pq); - // CHOOSE POINT END - - const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!! - const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT - - // R-DEPENDENT BEGIN - const complex double gamma_pq = lilgamma(rbeta_pq/k); - const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné - for(qpms_l_t j = 0; j <= lMax/2; ++j) { - int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j); - // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above - if(creal(z) < 0) - Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above - if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort(); - } - // R-DEPENDENT END - - // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array - // and just fetched for each n, m pair - for(qpms_l_t n = 0; n <= lMax; ++n) - for(qpms_m_t m = -n; m <= n; ++m) { - if((m+n) % 2 != 0) // odd coefficients are zero. - continue; - const qpms_y_t y = qpms_mn2y_sc(m, n); - const complex double e_imalpha_pq = cexp(I*m*arg_pq); - complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c); - double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? - assert((n-abs(m))/2 == c->s1_jMaxes[y]); - for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop - * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation - * c->s1_constfacs[y][j]; - if(err) { - // FIXME include also other errors than Gamma_pq's relative error - kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand)); - } - summand *= Gamma_pq[j].val; // GGG - ckahanadd(&jsum, &jsum_c, summand); - } - jsum *= phasefac; // PFC - ckahanadd(target + y, target_c + y, jsum); - if(err) kahanadd(err + y, err_c + y, jsum_err); - } - } // END POINT LOOP - - free(err_c); - free(target_c); - for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above - target[y] *= commonfac; - if(err) - for(qpms_y_t y = 0; y < nelem_sc; ++y) - err[y] *= commonfac; - return 0; -} - int ewald3_21_xy_sigma_long ( complex double *target, // must be c->nelem_sc long double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k, const double unitcell_volume /* with the corresponding lattice dimensionality */, const LatticeDimensionality latdim, @@ -517,7 +329,7 @@ int ewald3_21_xy_sigma_long ( int ewald3_1_z_sigma_long ( complex double *target, // must be c->nelem_sc long double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k, const double unitcell_volume /* length (periodicity) in this case */, const LatticeDimensionality latdim, @@ -633,7 +445,7 @@ int ewald3_1_z_sigma_long ( int ewald3_sigma_long ( complex double *target, // must be c->nelem_sc long double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k, const double unitcell_volume /* with the corresponding lattice dimensionality */, const LatticeDimensionality latdim, @@ -749,139 +561,10 @@ static int ewald32_sr_integral_ck(double r, complex double k, int n, double eta, return retval; } -int ewald32_sigma_short_shiftedpoints( - complex double *target, // must be c->nelem_sc long - double *err, - const qpms_ewald32_constants_t *c, // N.B. not too useful here - const double eta, const double k, - const size_t npoints, const point2d *Rpoints_plus_particle_shift, - const point2d beta, - const point2d particle_shift // used only in the very end to multiply it by the phase - ) -{ - const qpms_y_t nelem_sc = c->nelem_sc; - const qpms_l_t lMax = c->lMax; - gsl_integration_workspace *workspace = - gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT); - - // Manual init of the ewald summation targets - complex double * const target_c = calloc(nelem_sc, sizeof(complex double)); - memset(target, 0, nelem_sc * sizeof(complex double)); - double *err_c = NULL; - if (err) { - err_c = calloc(nelem_sc, sizeof(double)); - memset(err, 0, nelem_sc * sizeof(double)); - } - - -// CHOOSE POINT BEGIN - for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP - const point2d Rpq_shifted = Rpoints_plus_particle_shift[i]; - const double r_pq_shifted = cart2norm(Rpq_shifted); -// CHOOSE POINT END - - const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT - const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT - - for(qpms_l_t n = 0; n <= lMax; ++n) { - const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end - const double R_pq_pown = pow(r_pq_shifted, n); - // TODO the integral here - double intres, interr; - int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta, - &intres, &interr, workspace); - if (retval) abort(); - for (qpms_m_t m = -n; m <= n; ++m){ - if((m+n) % 2 != 0) // odd coefficients are zero. - continue; // nothing needed, already done by memset - const complex double e_imf = cexp(I*m*Rpq_shifted_arg); - const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))]; - const qpms_y_t y = qpms_mn2y_sc(m,n); - if(err) - kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown - * interr)); // TODO include also other errors - ckahanadd(target + y, target_c + y, - prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m)); - } - - } - } - - gsl_integration_workspace_free(workspace); - if(err) free(err_c); - free(target_c); - return 0; -} - -int ewald32_sigma_short_points_and_shift( - complex double *target, // must be c->nelem_sc long - double *err, - const qpms_ewald32_constants_t *c, // N.B. not too useful here - const double eta, const double k, - const size_t npoints, const point2d *Rpoints, - const point2d beta, - const point2d particle_shift // used only in the very end to multiply it by the phase - ) -{ - const qpms_y_t nelem_sc = c->nelem_sc; - const qpms_l_t lMax = c->lMax; - gsl_integration_workspace *workspace = - gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT); - - // Manual init of the ewald summation targets - complex double * const target_c = calloc(nelem_sc, sizeof(complex double)); - memset(target, 0, nelem_sc * sizeof(complex double)); - double *err_c = NULL; - if (err) { - err_c = calloc(nelem_sc, sizeof(double)); - memset(err, 0, nelem_sc * sizeof(double)); - } - - -// CHOOSE POINT BEGIN - for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP - //const point2d Rpq_shifted = Rpoints_plus_particle_shift[i]; - const point2d Rpq_shifted = cart2_add(Rpoints[i], cart2_scale(-1,particle_shift)); // CHECKSIGN!!!! - const double r_pq_shifted = cart2norm(Rpq_shifted); -// CHOOSE POINT END - - const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT - const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT - - for(qpms_l_t n = 0; n <= lMax; ++n) { - const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end - const double R_pq_pown = pow(r_pq_shifted, n); - // TODO the integral here - double intres, interr; - int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta, - &intres, &interr, workspace); - if (retval) abort(); - for (qpms_m_t m = -n; m <= n; ++m){ - if((m+n) % 2 != 0) // odd coefficients are zero. - continue; // nothing needed, already done by memset - const complex double e_imf = cexp(I*m*Rpq_shifted_arg); - const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))]; - const qpms_y_t y = qpms_mn2y_sc(m,n); - if(err) - kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown - * interr)); // TODO include also other errors - ckahanadd(target + y, target_c + y, - prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m)); - } - - } - } - - gsl_integration_workspace_free(workspace); - if(err) free(err_c); - free(target_c); - return 0; -} - int ewald3_sigma_short( complex double *target, // must be c->nelem_sc long double *err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k /* TODO COMPLEX */, const LatticeDimensionality latdim, // apart from asserts and possible optimisations ignored, as the SR formula stays the same PGen *pgen_R, const bool pgen_generates_shifted_points @@ -1027,23 +710,3 @@ int ewald3_sigma_short( return 0; } - - -#if 0 - -int ewald32_sigma_long_shiftedpoints_rordered( - complex double *target_sigmalr_y, // must be c->nelem_sc long - const qpms_ewald32_constants_t *c, - double eta, double k, double unitcell_area, - const points2d_rordered_t *Kpoints_plus_beta_rordered, - point2d particle_shift - ); -int ewald32_sigma_short_points_rordered( - complex double *target_sigmasr_y, // must be c->nelem_sc long - const qpms_ewald32_constants_t *c, // N.B. not too useful here - double eta, double k, - const points2d_rordered_t *Rpoints_plus_particle_shift_rordered, - point2d particle_shift // used only in the very end to multiply it by the phase - ); - -#endif diff --git a/qpms/ewald.h b/qpms/ewald.h index a4a5a56..ca794b3 100644 --- a/qpms/ewald.h +++ b/qpms/ewald.h @@ -76,10 +76,10 @@ typedef struct { This is because I dont't actually consider this fixed in translations.c */ -} qpms_ewald32_constants_t; +} qpms_ewald3_constants_t; -qpms_ewald32_constants_t *qpms_ewald32_constants_init(qpms_l_t lMax, int csphase); -void qpms_ewald32_constants_free(qpms_ewald32_constants_t *); +qpms_ewald3_constants_t *qpms_ewald3_constants_init(qpms_l_t lMax, int csphase); +void qpms_ewald3_constants_free(qpms_ewald3_constants_t *); typedef struct { // as gsl_sf_result, but with complex val @@ -144,14 +144,14 @@ int ewald32_sr_integral(double r, double k, double n, double eta, double *result // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space int ewald3_sigma0(complex double *result, double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, double eta, complex double k ); int ewald3_sigma_short( complex double *target_sigmasr_y, // must be c->nelem_sc long double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, const double eta, const complex double k, const LatticeDimensionality latdim, // apart from asserts and possible optimisations ignored, as the SR formula stays the same PGen *pgen_R, const bool pgen_generates_shifted_points @@ -168,7 +168,7 @@ int ewald3_sigma_short( int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim 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, + const qpms_ewald3_constants_t *c, const double eta, const complex double k, const double unitcell_volume /* with the corresponding lattice dimensionality */, const LatticeDimensionality latdim, @@ -182,11 +182,10 @@ int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, dep const cart3_t particle_shift ); -/// !!!!!!!!!!!!!!! ZDE JSEM SKONČIL !!!!!!!!!!!!!!!!!!!!!!. - +#ifdef EWALD_LEGACY // moved to ewald_legacy.c, not even everything implemented int ewald32_sigma0(complex double *result, double *err, // actually, this should be only alias for ewald3_sigma0 - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, double eta, double k ); @@ -196,7 +195,7 @@ int ewald32_sigma0(complex double *result, double *err, // actually, this should 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, + const qpms_ewald3_constants_t *c, double eta, double k, double unitcell_area, size_t npoints, const point2d *Kpoints_plus_beta, point2d beta, @@ -205,7 +204,7 @@ int ewald32_sigma_long_shiftedpoints ( int ewald32_sigma_long_points_and_shift ( 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, + const qpms_ewald3_constants_t *c, double eta, double k, double unitcell_area, size_t npoints, const point2d *Kpoints, point2d beta, @@ -214,7 +213,7 @@ int ewald32_sigma_long_points_and_shift ( int ewald32_sigma_long_shiftedpoints_rordered(//NI 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, + const qpms_ewald3_constants_t *c, double eta, double k, double unitcell_area, const points2d_rordered_t *Kpoints_plus_beta_rordered, point2d particle_shift @@ -223,7 +222,7 @@ int ewald32_sigma_long_shiftedpoints_rordered(//NI int ewald32_sigma_short_shiftedpoints( complex double *target_sigmasr_y, // must be c->nelem_sc long double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, // N.B. not too useful here + const qpms_ewald3_constants_t *c, // N.B. not too useful here double eta, double k, size_t npoints, const point2d *Rpoints_plus_particle_shift, point2d beta, @@ -232,7 +231,7 @@ int ewald32_sigma_short_shiftedpoints( int ewald32_sigma_short_points_and_shift( complex double *target_sigmasr_y, // must be c->nelem_sc long double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, // N.B. not too useful here + const qpms_ewald3_constants_t *c, // N.B. not too useful here double eta, double k, size_t npoints, const point2d *Rpoints, point2d beta, @@ -241,7 +240,7 @@ int ewald32_sigma_short_points_and_shift( int ewald32_sigma_short_points_rordered(//NI complex double *target_sigmasr_y, // must be c->nelem_sc long double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, // N.B. not too useful here + const qpms_ewald3_constants_t *c, // N.B. not too useful here double eta, double k, const points2d_rordered_t *Rpoints_plus_particle_shift_rordered, point2d particle_shift // used only in the very end to multiply it by the phase @@ -252,7 +251,7 @@ int ewald32_sigma_short_points_rordered(//NI int ewald31z_sigma_long_points_and_shift ( 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, + const qpms_ewald3_constants_t *c, double eta, double k, double unitcell_area, size_t npoints, const double *Kpoints, double beta, @@ -261,16 +260,17 @@ int ewald31z_sigma_long_points_and_shift ( int ewald31z_sigma_short_points_and_shift( complex double *target_sigmasr_y, // must be c->nelem_sc long double *target_sigmasr_y_err, // must be c->nelem_sc long or NULL - const qpms_ewald32_constants_t *c, // N.B. not too useful here + const qpms_ewald3_constants_t *c, // N.B. not too useful here double eta, double k, size_t npoints, const double *Rpoints, double beta, double particle_shift ); int ewald31z_sigma0(complex double *result, double *err, - const qpms_ewald32_constants_t *c, + const qpms_ewald3_constants_t *c, double eta, double k ); // exactly the same as ewald32_sigma0 +#endif // EWALD_LEGACY #endif //EWALD_H diff --git a/qpms/ewald_legacy.c b/qpms/ewald_legacy.c new file mode 100644 index 0000000..c056456 --- /dev/null +++ b/qpms/ewald_legacy.c @@ -0,0 +1,434 @@ +#include "ewald.h" +#include +#include "indexing.h" +#include "kahansum.h" +#include +#include +#include +#include "tiny_inlines.h" +#include +#include +#include +#include + +// parameters for the quadrature of integral in (4.6) +#ifndef INTEGRATION_WORKSPACE_LIMIT +#define INTEGRATION_WORKSPACE_LIMIT 30000 +#endif + +#ifndef INTEGRATION_EPSABS +#define INTEGRATION_EPSABS 1e-13 +#endif + +#ifndef INTEGRATION_EPSREL +#define INTEGRATION_EPSREL 1e-13 +#endif + +#ifndef M_SQRTPI +#define M_SQRTPI 1.7724538509055160272981674833411452 +#endif + +// sloppy implementation of factorial +static inline double factorial(const int n) { + assert(n >= 0); + if (n < 0) + return 0; // should not happen in the functions below. (Therefore the assert above) + else if (n <= 20) { + double fac = 1; + for (int i = 1; i <= n; ++i) + fac *= i; + return fac; + } + else + return tgamma(n + 1); // hope it's precise and that overflow does not happen +} + +static inline complex double csq(complex double x) { return x * x; } +static inline double sq(double x) { return x * x; } + + +typedef enum { + EWALD32_CONSTANTS_ORIG, // As in [1, (4,5)], NOT USED right now. + EWALD32_CONSTANTS_AGNOSTIC /* Not depending on the spherical harmonic sign/normalisation + * convention – the $e^{im\alpha_pq}$ term in [1,(4.5)] being + * replaced by the respective $Y_n^m(\pi/2,\alpha)$ + * spherical harmonic. See notes/ewald.lyx. + */ + +} ewald3_constants_option; + +static const ewald3_constants_option type = EWALD32_CONSTANTS_AGNOSTIC; + + + +int ewald32_sigma0(complex double *result, double *err, + const qpms_ewald3_constants_t *c, + const double eta, const double k) +{ + return ewald3_sigma0(result, err, c, eta, k); +} + +int ewald32_sigma_long_shiftedpoints ( + complex double *target, // must be c->nelem_sc long + double *err, + const qpms_ewald3_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 + ) +{ + const qpms_y_t nelem_sc = c->nelem_sc; + const qpms_l_t lMax = c->lMax; + + // Manual init of the ewald summation targets + complex double *target_c = calloc(nelem_sc, sizeof(complex double)); + memset(target, 0, nelem_sc * sizeof(complex double)); + double *err_c = NULL; + if (err) { + err_c = calloc(nelem_sc, sizeof(double)); + memset(err, 0, nelem_sc * sizeof(double)); + } + + const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC) + assert(commonfac > 0); + + // space for Gamma_pq[j]'s + qpms_csf_result Gamma_pq[lMax/2+1]; + + // CHOOSE POINT BEGIN + for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP + const point2d beta_pq = Kpoints_plus_beta[i]; + const point2d K_pq = {beta_pq.x - beta.x, beta_pq.y - beta.y}; + const double rbeta_pq = cart2norm(beta_pq); + // CHOOSE POINT END + + const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!! + const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT + + // R-DEPENDENT BEGIN + const complex double gamma_pq = lilgamma(rbeta_pq/k); + const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné + for(qpms_l_t j = 0; j <= lMax/2; ++j) { + int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j); + // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above + if(creal(z) < 0) + Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above + if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort(); + } + // R-DEPENDENT END + + // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array + // and just fetched for each n, m pair + for(qpms_l_t n = 0; n <= lMax; ++n) + for(qpms_m_t m = -n; m <= n; ++m) { + if((m+n) % 2 != 0) // odd coefficients are zero. + continue; + const qpms_y_t y = qpms_mn2y_sc(m, n); + const complex double e_imalpha_pq = cexp(I*m*arg_pq); + complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c); + double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? + assert((n-abs(m))/2 == c->s1_jMaxes[y]); + for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop + * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation + * c->s1_constfacs[y][j]; + if(err) { + // FIXME include also other errors than Gamma_pq's relative error + kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand)); + } + summand *= Gamma_pq[j].val; // GGG + ckahanadd(&jsum, &jsum_c, summand); + } + jsum *= phasefac; // PFC + ckahanadd(target + y, target_c + y, jsum); + if(err) kahanadd(err + y, err_c + y, jsum_err); + } + } // END POINT LOOP + + free(err_c); + free(target_c); + for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above + target[y] *= commonfac; + if(err) + for(qpms_y_t y = 0; y < nelem_sc; ++y) + err[y] *= commonfac; + return 0; +} + +int ewald32_sigma_long_points_and_shift ( + complex double *target, // must be c->nelem_sc long + double *err, + const qpms_ewald3_constants_t *c, + const double eta, const double k, const double unitcell_area, + const size_t npoints, const point2d *Kpoints, + const point2d beta, + const point2d particle_shift // target - src + ) +{ + const qpms_y_t nelem_sc = c->nelem_sc; + const qpms_l_t lMax = c->lMax; + + // Manual init of the ewald summation targets + complex double *target_c = calloc(nelem_sc, sizeof(complex double)); + memset(target, 0, nelem_sc * sizeof(complex double)); + double *err_c = NULL; + if (err) { + err_c = calloc(nelem_sc, sizeof(double)); + memset(err, 0, nelem_sc * sizeof(double)); + } + + const double commonfac = 1/(k*k*unitcell_area); // used in the very end (CFC) + assert(commonfac > 0); + + // space for Gamma_pq[j]'s + qpms_csf_result Gamma_pq[lMax/2+1]; + + // CHOOSE POINT BEGIN + for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP + // Only these following two lines differ from ewald32_sigma_long_points_and_shift()!!! WTFCOMMENT?! + const point2d K_pq = Kpoints[i]; + const point2d beta_pq = {K_pq.x + beta.x, K_pq.y + beta.y}; + const double rbeta_pq = cart2norm(beta_pq); + // CHOOSE POINT END + + const complex double phasefac = cexp(I*cart2_dot(K_pq,particle_shift)); // POINT-DEPENDENT (PFC) // !!!CHECKSIGN!!! + const double arg_pq = atan2(beta_pq.y, beta_pq.x); // POINT-DEPENDENT + + // R-DEPENDENT BEGIN + const complex double gamma_pq = lilgamma(rbeta_pq/k); + const complex double z = csq(gamma_pq*k/(2*eta)); // Když o tom tak přemýšlím, tak tohle je vlastně vždy reálné + for(qpms_l_t j = 0; j <= lMax/2; ++j) { + int retval = complex_gamma_inc_e(0.5-j, z, Gamma_pq+j); + // we take the other branch, cf. [Linton, p. 642 in the middle]: FIXME instead use the C11 CMPLX macros and fill in -O*I part to z in the line above + if(creal(z) < 0) + Gamma_pq[j].val = conj(Gamma_pq[j].val); //FIXME as noted above + if(!(retval==0 ||retval==GSL_EUNDRFLW)) abort(); + } + // R-DEPENDENT END + + // TODO optimisations: all the j-dependent powers can be done for each j only once, stored in array + // and just fetched for each n, m pair + for(qpms_l_t n = 0; n <= lMax; ++n) + for(qpms_m_t m = -n; m <= n; ++m) { + if((m+n) % 2 != 0) // odd coefficients are zero. + continue; + const qpms_y_t y = qpms_mn2y_sc(m, n); + const complex double e_imalpha_pq = cexp(I*m*arg_pq); + complex double jsum, jsum_c; ckahaninit(&jsum, &jsum_c); + double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? + assert((n-abs(m))/2 == c->s1_jMaxes[y]); + for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME legendre0[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) // This line can actually go outside j-loop + * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation + * c->s1_constfacs[y][j]; + if(err) { + // FIXME include also other errors than Gamma_pq's relative error + kahanadd(&jsum_err, &jsum_err_c, Gamma_pq[j].err * cabs(summand)); + } + summand *= Gamma_pq[j].val; // GGG + ckahanadd(&jsum, &jsum_c, summand); + } + jsum *= phasefac; // PFC + ckahanadd(target + y, target_c + y, jsum); + if(err) kahanadd(err + y, err_c + y, jsum_err); + } + } // END POINT LOOP + + free(err_c); + free(target_c); + for(qpms_y_t y = 0; y < nelem_sc; ++y) // CFC common factor from above + target[y] *= commonfac; + if(err) + for(qpms_y_t y = 0; y < nelem_sc; ++y) + err[y] *= commonfac; + return 0; +} + +struct sigma2_integrand_params { + int n; + double k, R; +}; + +static double sigma2_integrand(double ksi, void *params) { + struct sigma2_integrand_params *p = (struct sigma2_integrand_params *) params; + return exp(-sq(p->R * ksi) + sq(p->k / ksi / 2)) * pow(ksi, 2*p->n); +} + +static int ewald32_sr_integral(double r, double k, int n, double eta, + double *result, double *err, gsl_integration_workspace *workspace) +{ + struct sigma2_integrand_params p; + + + const double R0 = r; // Maybe could be chosen otherwise, but fuck it for now. + p.n = n; + eta *= R0; + p.k = k * R0; + p.R = r / R0; // i.e. p.R = 1 + + gsl_function F; + F.function = sigma2_integrand; + F.params = &p; + + int retval = gsl_integration_qagiu(&F, eta, INTEGRATION_EPSABS, + INTEGRATION_EPSREL, INTEGRATION_WORKSPACE_LIMIT, + workspace, result, err); + double normfac = pow(R0, -2*p.n - 1); + *result *= normfac; + *err *= normfac; + return retval; +} + + +int ewald32_sigma_short_shiftedpoints( + complex double *target, // must be c->nelem_sc long + double *err, + const qpms_ewald3_constants_t *c, // N.B. not too useful here + const double eta, const double k, + const size_t npoints, const point2d *Rpoints_plus_particle_shift, + const point2d beta, + const point2d particle_shift // used only in the very end to multiply it by the phase + ) +{ + const qpms_y_t nelem_sc = c->nelem_sc; + const qpms_l_t lMax = c->lMax; + gsl_integration_workspace *workspace = + gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT); + + // Manual init of the ewald summation targets + complex double * const target_c = calloc(nelem_sc, sizeof(complex double)); + memset(target, 0, nelem_sc * sizeof(complex double)); + double *err_c = NULL; + if (err) { + err_c = calloc(nelem_sc, sizeof(double)); + memset(err, 0, nelem_sc * sizeof(double)); + } + + +// CHOOSE POINT BEGIN + for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP + const point2d Rpq_shifted = Rpoints_plus_particle_shift[i]; + const double r_pq_shifted = cart2norm(Rpq_shifted); +// CHOOSE POINT END + + const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT + const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT + + for(qpms_l_t n = 0; n <= lMax; ++n) { + const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end + const double R_pq_pown = pow(r_pq_shifted, n); + // TODO the integral here + double intres, interr; + int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta, + &intres, &interr, workspace); + if (retval) abort(); + for (qpms_m_t m = -n; m <= n; ++m){ + if((m+n) % 2 != 0) // odd coefficients are zero. + continue; // nothing needed, already done by memset + const complex double e_imf = cexp(I*m*Rpq_shifted_arg); + const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))]; + const qpms_y_t y = qpms_mn2y_sc(m,n); + if(err) + kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown + * interr)); // TODO include also other errors + ckahanadd(target + y, target_c + y, + prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m)); + } + + } + } + + gsl_integration_workspace_free(workspace); + if(err) free(err_c); + free(target_c); + return 0; +} + +int ewald32_sigma_short_points_and_shift( + complex double *target, // must be c->nelem_sc long + double *err, + const qpms_ewald3_constants_t *c, // N.B. not too useful here + const double eta, const double k, + const size_t npoints, const point2d *Rpoints, + const point2d beta, + const point2d particle_shift // used only in the very end to multiply it by the phase + ) +{ + const qpms_y_t nelem_sc = c->nelem_sc; + const qpms_l_t lMax = c->lMax; + gsl_integration_workspace *workspace = + gsl_integration_workspace_alloc(INTEGRATION_WORKSPACE_LIMIT); + + // Manual init of the ewald summation targets + complex double * const target_c = calloc(nelem_sc, sizeof(complex double)); + memset(target, 0, nelem_sc * sizeof(complex double)); + double *err_c = NULL; + if (err) { + err_c = calloc(nelem_sc, sizeof(double)); + memset(err, 0, nelem_sc * sizeof(double)); + } + + +// CHOOSE POINT BEGIN + for (size_t i = 0; i < npoints; ++i) { // BEGIN POINT LOOP + //const point2d Rpq_shifted = Rpoints_plus_particle_shift[i]; + const point2d Rpq_shifted = cart2_add(Rpoints[i], cart2_scale(-1,particle_shift)); // CHECKSIGN!!!! + const double r_pq_shifted = cart2norm(Rpq_shifted); +// CHOOSE POINT END + + const double Rpq_shifted_arg = atan2(Rpq_shifted.y, Rpq_shifted.x); // POINT-DEPENDENT + const complex double e_beta_Rpq = cexp(I*cart2_dot(beta, Rpq_shifted)); // POINT-DEPENDENT + + for(qpms_l_t n = 0; n <= lMax; ++n) { + const double complex prefacn = - I * pow(2./k, n+1) * M_2_SQRTPI / 2; // TODO put outside the R-loop and multiply in the end + const double R_pq_pown = pow(r_pq_shifted, n); + // TODO the integral here + double intres, interr; + int retval = ewald32_sr_integral(r_pq_shifted, k, n, eta, + &intres, &interr, workspace); + if (retval) abort(); + for (qpms_m_t m = -n; m <= n; ++m){ + if((m+n) % 2 != 0) // odd coefficients are zero. + continue; // nothing needed, already done by memset + const complex double e_imf = cexp(I*m*Rpq_shifted_arg); + const double leg = c->legendre0[gsl_sf_legendre_array_index(n, abs(m))]; + const qpms_y_t y = qpms_mn2y_sc(m,n); + if(err) + kahanadd(err + y, err_c + y, cabs(leg * (prefacn / I) * R_pq_pown + * interr)); // TODO include also other errors + ckahanadd(target + y, target_c + y, + prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m)); + } + + } + } + + gsl_integration_workspace_free(workspace); + if(err) free(err_c); + free(target_c); + return 0; +} + + + +#if 0 + +int ewald32_sigma_long_shiftedpoints_rordered( + complex double *target_sigmalr_y, // must be c->nelem_sc long + const qpms_ewald3_constants_t *c, + double eta, double k, double unitcell_area, + const points2d_rordered_t *Kpoints_plus_beta_rordered, + point2d particle_shift + ); +int ewald32_sigma_short_points_rordered( + complex double *target_sigmasr_y, // must be c->nelem_sc long + const qpms_ewald3_constants_t *c, // N.B. not too useful here + double eta, double k, + const points2d_rordered_t *Rpoints_plus_particle_shift_rordered, + point2d particle_shift // used only in the very end to multiply it by the phase + ); + +#endif diff --git a/tests/ewalds.c b/tests/ewalds.c index 5623403..32af265 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -185,7 +185,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { 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); + qpms_ewald3_constants_t *c = qpms_ewald3_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); @@ -252,7 +252,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { } points2d_rordered_free(Kpoints_plus_beta); - qpms_ewald32_constants_free(c); + qpms_ewald3_constants_free(c); triangular_lattice_gen_free(Klg); triangular_lattice_gen_free(Rlg); ++ewaldtest_counter; diff --git a/tests/ewaldshift.c b/tests/ewaldshift.c index 0d00e87..3bb2984 100644 --- a/tests/ewaldshift.c +++ b/tests/ewaldshift.c @@ -266,7 +266,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { 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); + qpms_ewald3_constants_t *c = qpms_ewald3_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); @@ -332,7 +332,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { } points2d_rordered_free(Kpoints_plus_beta); - qpms_ewald32_constants_free(c); + qpms_ewald3_constants_free(c); triangular_lattice_gen_free(Klg); triangular_lattice_gen_free(Rlg); ++ewaldtest_counter; diff --git a/tests/ewaldshift2.c b/tests/ewaldshift2.c index 77bd250..f9460b7 100644 --- a/tests/ewaldshift2.c +++ b/tests/ewaldshift2.c @@ -295,7 +295,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { 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); + qpms_ewald3_constants_t *c = qpms_ewald3_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); @@ -361,7 +361,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { } points2d_rordered_free(Kpoints_plus_beta); - qpms_ewald32_constants_free(c); + qpms_ewald3_constants_free(c); triangular_lattice_gen_free(Klg); triangular_lattice_gen_free(Rlg); ++ewaldtest_counter; diff --git a/tests/ewaldshift_3g.c b/tests/ewaldshift_3g.c index 83b51e9..7b493b9 100644 --- a/tests/ewaldshift_3g.c +++ b/tests/ewaldshift_3g.c @@ -307,7 +307,7 @@ ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) { 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); + qpms_ewald3_constants_t *c = qpms_ewald3_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); @@ -395,7 +395,7 @@ ewaldtest_triang_results *ewaldtest_triang_3g(const ewaldtest_triang_params p) { //points2d_rordered_free(Kpoints_plus_beta); - qpms_ewald32_constants_free(c); + qpms_ewald3_constants_free(c); //triangular_lattice_gen_free(Klg); //triangular_lattice_gen_free(Rlg); ++ewaldtest_counter; diff --git a/tests/ewaldshift_ng.c b/tests/ewaldshift_ng.c index e3d0991..c5172ad 100644 --- a/tests/ewaldshift_ng.c +++ b/tests/ewaldshift_ng.c @@ -295,7 +295,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { 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); + qpms_ewald3_constants_t *c = qpms_ewald3_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); @@ -361,7 +361,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { } points2d_rordered_free(Kpoints_plus_beta); - qpms_ewald32_constants_free(c); + qpms_ewald3_constants_free(c); triangular_lattice_gen_free(Klg); triangular_lattice_gen_free(Rlg); ++ewaldtest_counter;