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

Handle k->0 in Ewald sum. 

Former-commit-id: fb61b297ead4828541b763b2ecdb0c81de1d35f1
This commit is contained in:
Marek Nečada 2019-09-16 12:37:32 +03:00
parent ecbf5b10d4
commit e8a2426b55
3 changed files with 20 additions and 1 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
qpms/ewald.c
View File

@ -189,6 +189,17 @@ int ewald3_sigma0(complex double *result, double *err,
return 0; return 0;
} }
// Wrapper over cpow to correctly handle the k->0 behaviour
static inline complex double cpow_0lim_zi(complex double x, long z)
{
if (x == 0 && z == 0)
return 1;
else if (x == 0 && z > 0) // lol branch cut in cpow
return 0;
else
return cpow(x, z);
}
int ewald3_21_xy_sigma_long ( int ewald3_21_xy_sigma_long (
complex double *target, // must be c->nelem_sc long complex double *target, // must be c->nelem_sc long
double *err, double *err,
@ -293,7 +304,7 @@ int ewald3_21_xy_sigma_long (
double jsum_err, jsum_err_c; kahaninit(&jsum_err, &jsum_err_c); // TODO do I really need to kahan sum errors? 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]); 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 </<= ? for(qpms_l_t j = 0; j <= c->s1_jMaxes[y]/*(n-abs(m))/2*/; ++j) { // FIXME </<= ?
complex double summand = cpow(rbeta_pq/k, n-2*j) complex double summand = cpow_0lim_zi(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 * 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 * cpow(gamma_pq, 2*j-1) // * Gamma_pq[j] bellow (GGG) after error computation
* c->s1_constfacs[y][j]; * c->s1_constfacs[y][j];

5
tests/CMakeLists.txt
View File

@ -48,3 +48,8 @@ add_test( NAME scalar_ewald32_square_cplxk_negative_Im2 COMMAND test_scalar_ewal
# lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ... # lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ...
3 1 0 0 1 2.3 -0.1 2.7 1 0.0 0.0 -1 1e-8 1e-10 20 160 1.5 1.6 2.5 2.6 3 1 0 0 1 2.3 -0.1 2.7 1 0.0 0.0 -1 1e-8 1e-10 20 160 1.5 1.6 2.5 2.6
) )
add_test( NAME scalar_ewald32_square_zerok COMMAND test_scalar_ewald32
# lMax b1.x b1.y b2.x b2.y wavenum.real wavenum.imag k.x k.y particle_shift.x particle_shift.y csphase rtol atol maxR maxK eta1 ...
3 1 0 0 1 2.3 0 0 0 0.5 0.1325 -1 1e-8 1e-10 40 160 1.5 1.6 2.5 2.6
)

3
tests/test_scalar_ewald32.c
View File

@ -17,6 +17,8 @@
#include <string.h> #include <string.h>
#include <float.h> #include <float.h>
#include <gsl/gsl_sf_legendre.h> #include <gsl/gsl_sf_legendre.h>
#include <fenv.h>
typedef struct ewaldtest2d_params { typedef struct ewaldtest2d_params {
qpms_l_t lMax; qpms_l_t lMax;
point2d b1, b2; point2d b1, b2;
@ -66,6 +68,7 @@ int isclose_cmplx(complex double a, complex double b, double rtol, double atol)
ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p); ewaldtest2d_results *ewaldtest2d(const ewaldtest2d_params p);
int main(int argc, char **argv) { int main(int argc, char **argv) {
feenableexcept(FE_INVALID | FE_OVERFLOW);
bool verbose = !!getenv("QPMS_VERBOSE_TESTS"); bool verbose = !!getenv("QPMS_VERBOSE_TESTS");
gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler); gsl_set_error_handler(IgnoreUnderflowsGSLErrorHandler);
QPMS_ENSURE(argc >= 18, "At least 16 arguments expected, I see only %d.", argc-1); QPMS_ENSURE(argc >= 18, "At least 16 arguments expected, I see only %d.", argc-1);