From 84a2b7f64d6386d0c0278b72920bc41f1235479e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Mon, 10 Sep 2018 15:28:41 +0300 Subject: [PATCH] Fix spherical harmonics for negative m in ewald.c. Former-commit-id: 960dd721a009e2b2f82c986677c3faa65d36394c --- notes/ewald.lyx | 9 ++++++--- qpms/ewald.c | 4 ++-- qpms/tiny_inlines.h | 7 +++++++ tests/ewalds.c | 2 +- 4 files changed, 16 insertions(+), 6 deletions(-) diff --git a/notes/ewald.lyx b/notes/ewald.lyx index 43e0c5d..126cd3c 100644 --- a/notes/ewald.lyx +++ b/notes/ewald.lyx @@ -3104,9 +3104,12 @@ Y_{n}^{m}\left(\frac{\pi}{2},\phi\right) & = & \left(-1\right)^{m}\sqrt{\frac{\l Let us substitute this into [LT(4.5)] \begin_inset Formula \begin{eqnarray*} -\sigma_{n}^{m(1)} & = & -\frac{i^{n+1}}{2k^{2}\mathscr{A}}\left(-1\right)^{\left(n+m\right)/2}\sqrt{\left(2n+1\right)\left(n-m\right)!\left(n+m\right)!}\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/2k\right)^{n-2j}e^{im\phi_{\vect{\beta}_{pq}}}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\frac{\gamma_{pq}}{2}\right)^{2j-1}\\ - & = & -\frac{i^{n+1}}{2k^{2}\mathscr{A}}\sqrt{\pi}2^{n+1}\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!Y_{n}^{m}\left(0,\phi_{\vect{\beta}_{pq}}\right)\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/2k\right)^{n-2j}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\frac{\gamma_{pq}}{2}\right)^{2j-1}\\ - & = & -\frac{i^{n+1}}{k\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!Y_{n}^{m}\left(0,\phi_{\vect{\beta}_{pq}}\right)\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/k\right)^{n-2j}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\gamma_{pq}\right)^{2j-1}, +\sigma_{n}^{m(1)} & = & -\frac{i^{n+1}}{2k^{2}\mathscr{A}}\left(-1\right)^{\left(n+m\right)/2}\sqrt{\left(2n+1\right)\left(n-m\right)!\left(n+m\right)!}\times\\ + & & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/2k\right)^{n-2j}e^{im\phi_{\vect{\beta}_{pq}}}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\frac{\gamma_{pq}}{2}\right)^{2j-1}\\ + & = & -\frac{i^{n+1}}{2k^{2}\mathscr{A}}\sqrt{\pi}2^{n+1}\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\times\\ + & & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(0,\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/2k\right)^{n-2j}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\frac{\gamma_{pq}}{2}\right)^{2j-1}\\ + & = & -\frac{i^{n+1}}{k\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\times\\ + & & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(0,\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{\left[\left(n-\left|m\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/k\right)^{n-2j}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(n-m\right)-j\right)!\left(\frac{1}{2}\left(n+m\right)-j\right)!}\left(\gamma_{pq}\right)^{2j-1} \end{eqnarray*} \end_inset diff --git a/qpms/ewald.c b/qpms/ewald.c index e35942c..d8915e6 100644 --- a/qpms/ewald.c +++ b/qpms/ewald.c @@ -208,7 +208,7 @@ int ewald32_sigma_long_shiftedpoints ( 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))] // 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 * c->s1_constfacs[y][j]; if(err) { @@ -322,7 +322,7 @@ int ewald32_sigma_short_shiftedpoints( 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); + prefacn * R_pq_pown * leg * intres * e_beta_Rpq * e_imf * min1pow_m_neg(m)); } } diff --git a/qpms/tiny_inlines.h b/qpms/tiny_inlines.h index 0ac7848..3b01f0f 100644 --- a/qpms/tiny_inlines.h +++ b/qpms/tiny_inlines.h @@ -4,6 +4,13 @@ static inline int min1pow(int pow) { return (pow % 2) ? -1 : 1; } + +// This is useful for calculating spherical harmonics with negative m +// if spharm-normalised legendre functions for positive m are available. +static inline int min1pow_m_neg(int m) { + return (m < 0) ? min1pow(m) : 1; +} + // this has shitty precision: // static inline complex double ipow(int x) { return cpow(I, x); } diff --git a/tests/ewalds.c b/tests/ewalds.c index 3774785..8cc24a1 100644 --- a/tests/ewalds.c +++ b/tests/ewalds.c @@ -175,7 +175,7 @@ ewaldtest_triang_results *ewaldtest_triang(const ewaldtest_triang_params p) { complex double eimf = cexp(I*m*arg_pq); results->regsigmas_416[y] += 4*M_PI*ipow(n)/p.k/A - * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] + * eimf * legendres[gsl_sf_legendre_array_index(n,abs(m))] * min1pow_m_neg(m) / denom; } }