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

Rename ewald3_constants, move legacy code 

Former-commit-id: e83dcfa532f7b8d7345103752aca924a56ad7138
This commit is contained in:
Marek Nečada 2018-12-21 19:31:07 +00:00
parent fd1aed02ca
commit 328d22de89
8 changed files with 472 additions and 375 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

357
qpms/ewald.c
View File

@ -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 </<= ?
complex double summand = pow(rbeta_pq/k, n-2*j)
* e_imalpha_pq * c->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 </<= ?
complex double summand = pow(rbeta_pq/k, n-2*j)
* e_imalpha_pq * c->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

36
qpms/ewald.h
View File

@ -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

434
qpms/ewald_legacy.c Normal file
View File

@ -0,0 +1,434 @@
#include "ewald.h"
#include <stdlib.h>
#include "indexing.h"
#include "kahansum.h"
#include <assert.h>
#include <string.h>
#include <complex.h>
#include "tiny_inlines.h"
#include <gsl/gsl_integration.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_sf_legendre.h>
#include <gsl/gsl_sf_expint.h>
// 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 </<= ?
complex double summand = pow(rbeta_pq/k, n-2*j)
* e_imalpha_pq * c->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 </<= ?
complex double summand = pow(rbeta_pq/k, n-2*j)
* e_imalpha_pq * c->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

4
tests/ewalds.c
View File

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

4
tests/ewaldshift.c
View File

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

4
tests/ewaldshift2.c
View File

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

4
tests/ewaldshift_3g.c
View File

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

4
tests/ewaldshift_ng.c
View File

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