Merge branch 'beyn_pureblas'
[qpms.git] / qpms / ewald.h
blob8eda80c39434eedd2bcb2cfc1fc1b5ccb1d26296
1 /*! \file ewald.h
2 * \brief Lattice sums of spherical waves.
4 * Implementation of two-dimensional lattice sum in three dimensions
5 * according to:
6 * - [1] C.M. Linton, I. Thompson
7 * Journal of Computational Physics 228 (2009) 1815–1829
8 * - [2] C.M.Linton
9 * SIAM Review Vol 52, No. 4, pp. 630–674
11 * N.B.!!! currently, the long-range parts are calculated
12 * not according to [1,(4.5)], but rather
13 * according to the spherical-harmonic-normalisation-independent
14 * formulation in my notes notes/ewald.lyx.
15 * Both parts of lattice sums are then calculated with
16 * the \f$ P_n^{|m|} e^{im\phi} \f$
17 * (N.B. or \f$ P_n^{|m|} e^{imf} (-1)^m \f$ for negative m)
18 * substituted in place of \f$ Y_n^m \f$
19 * (this is quite a weird normalisation especially
20 * for negative \f$ |m| \f$, but it is consistent
21 * with the current implementation of the translation coefficients in
22 * @ref translations.c;
23 * in the long run, it might make more sense to replace it everywhere with normalised
24 * Legendre polynomials).
27 #ifndef EWALD_H
28 #define EWALD_H
29 #include <gsl/gsl_sf_result.h>
30 #include <stdlib.h>
31 #include <gsl/gsl_sf_legendre.h>
32 #include <gsl/gsl_errno.h>
33 #include <math.h> // for inlined lilgamma
34 #include <complex.h>
35 #include "qpms_types.h"
36 #include "lattices.h"
39 typedef enum {
40 QPMS_EWALD_LONG_RANGE = 1,
41 QPMS_EWALD_SHORT_RANGE = 2,
42 QPMS_EWALD_0TERM = 4,
43 QPMS_EWALD_FULL = QPMS_EWALD_LONG_RANGE | QPMS_EWALD_SHORT_RANGE | QPMS_EWALD_0TERM,
44 } qpms_ewald_part;
47 /// Use this handler to ignore underflows of incomplete gamma.
48 gsl_error_handler_t IgnoreUnderflowsGSLErrorHandler;
51 /// Object holding the Ewald sum constant factors.
52 /**
53 * Used internally by qpms_translation_calculator_t.
54 * Initialised by qpms_ewald3_constants_init() and freed by qpms_ewald3_constants_free().
56 typedef struct qpms_ewald3_constants_t {
57 qpms_l_t lMax;
58 qpms_y_t nelem_sc;
59 /// The values of maximum \a j's in the long-range part summation, `[(l-|m|/2)]`.
60 qpms_l_t *s1_jMaxes;
61 /// The constant factors for the long range part of a 2D Ewald sum.
62 complex double **s1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
63 /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
64 * for m + n EVEN:
66 * s1_constfacs[y(m,n)][j] =
68 * -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j
69 * -----------------------------------------------------------
70 * j! * ((n-m)/2 - j)! * ((n+m)/2 + j)!
72 * for m + n ODD:
74 * s1_constfacs[y(m,n)][j] = 0
76 complex double *s1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
77 // similarly for the 1D z-axis aligned case; now the indices are [n][j] (as m == 0)
78 /// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
79 /** If the summation points lie along a different direction, use the formula for
80 * 2D sum with additional factor of
81 * \f$ \sqrt{pi} \kappa \gamma(\abs{\vect{k}+\vect{K}}/\kappa) \f$.
84 ///=============== NEW GENERATION GENERAL 2D-IN-3D, including z != 0 =========================
86 // TODO indexing mechanisms
88 /// The constant factors for the long range part of a 2D Ewald sum.
89 complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
90 /* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
91 * for m + n EVEN:
93 * S1_constfacs[y(m,n)][x(j,s)] =
95 * -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j / j \
96 * ----------------------------------------------------------- | |
97 * j! * ((n - m - s)/2)! * ((n + m - s)/2)! \ 2j - s /
99 * for m + n ODD:
101 * S1_constfacs[y(m,n)][j] = 0
103 complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
104 /// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
105 /** If the summation points lie along a different direction, use the formula for
106 * 2D sum with additional factor of
107 * \f$ \sqrt{pi} \kappa \gamma(\abs{\vect{k}+\vect{K}}/\kappa) \f$.
111 complex double **s1_constfacs_1Dz;
112 /* These are the actual numbers now:
113 * s1_constfacs_1Dz[n][j] =
115 * -I**(n+1) (-1)**j * n!
116 * --------------------------
117 * j! * 2**(2*j) * (n - 2*j)!
119 complex double *s1_constfacs_1Dz_base; ///<Internal pointer holding memory for the 1D Ewald sum constant factors.
121 double *legendre0; /* now with GSL_SF_LEGENDRE_NONE normalisation, because this is what is
122 * what the multipliers from translations.c count with.
124 double *legendre_plus1; // needed? TODO; in any case, nonzero only for m=0
125 double *legendre_minus1; // needed? TODO; in any case, nonzero only for m=0
126 gsl_sf_legendre_t legendre_normconv;
127 int legendre_csphase; /* 1 or -1; csphase of the Legendre polynomials saved in legendre0 etc.
128 This is because I dont't actually consider this fixed in
129 translations.c */
131 } qpms_ewald3_constants_t;
133 /// Constructor for qpms_ewald3_constants_t.
134 qpms_ewald3_constants_t *qpms_ewald3_constants_init(qpms_l_t lMax, int csphase);
135 /// Destructor for qpms_ewald3_constants_t.
136 void qpms_ewald3_constants_free(qpms_ewald3_constants_t *);
139 /// Structure for holding complex-valued result of computation and an error estimate.
140 /** Similar to gsl_sf_result, but with complex val. */
141 typedef struct qpms_csf_result {
142 complex double val; ///< Calculation result.
143 double err; ///< Error estimate.
144 } qpms_csf_result;
147 // [1, (A.9)]
148 static inline complex double lilgamma(double t) {
149 t = fabs(t);
150 if (t >= 1)
151 return sqrt(t*t - 1);
152 else
153 return -I * sqrt(1 - t*t);
156 // [1, (A.8)], complex version of lilgamma()
157 static inline complex double clilgamma(complex double z) {
158 complex double a1 = z - 1, a2 = z + 1;
159 // ensure -pi/2 < arg(z + 1) < 3*pi/2
160 if (creal(a2) < 0 && cimag(a2) <= 0)
161 a2 = -csqrt(a2);
162 else
163 a2 = csqrt(a2);
164 // ensure -3*pi/2 < arg(z - 1) < pi/2
165 if (creal(a1) < 0 && cimag(a1) >= 0)
166 a1 = -csqrt(a1);
167 else
168 a1 = csqrt(a1);
169 return a1 * a2;
172 /// Incomplete Gamma function as a series.
173 /** DLMF 8.7.3 (latter expression) for complex second argument.
175 * The principal value is calculated. On the negative real axis
176 * (where the function has branch cut), the sign of the imaginary
177 * part is what matters (even if it is zero). Therefore one
178 * can have
179 * `cx_gamma_inc_series_e(a, z1) != cx_gamma_inc_series_e(a, z2)`
180 * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
181 * The side of the branch cut can be determined using `signbit(creal(z))`.
183 int cx_gamma_inc_series_e(double a, complex double z, qpms_csf_result * result);
185 /// Incomplete Gamma function as continued fractions.
186 /**
187 * The principal value is calculated. On the negative real axis
188 * (where the function has branch cut), the sign of the imaginary
189 * part is what matters (even if it is zero). Therefore one
190 * can have
191 * `cx_gamma_inc_CF_e(a, z1) != cx_gamma_inc_CF_e(a, z2)`
192 * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
193 * The side of the branch cut can be determined using `signbit(creal(z))`.
195 int cx_gamma_inc_CF_e(double a, complex double z, qpms_csf_result * result);
197 /// Incomplete gamma for complex second argument.
198 /**
199 * If x is (almost) real, it just uses gsl_sf_gamma_inc_e().
201 * On the negative real axis
202 * (where the function has branch cut), the sign of the imaginary
203 * part is what matters (even if it is zero). Therefore one
204 * can have
205 * `complex_gamma_inc_e(a, z1, m) != complex_gamma_inc_e(a, z2, m)`
206 * even if `z1 == z2`, because `-0 == 0` according to IEEE 754.
207 * The side of the branch cut can be determined using `signbit(creal(z))`.
209 * Another than principal branch can be selected using non-zero \a m
210 * argument.
212 int complex_gamma_inc_e(double a, complex double x,
213 /// Branch index.
214 /** If zero, the principal value is calculated.
215 * Other branches might be chosen using non-zero \a m.
216 * In such case, the returned value corresponds to \f[
217 * \Gamma(a,ze^{2\pi mi})=e^{2\pi mia} \Gamma(a,z)
218 * + (1-e^{2\pi mia}) \Gamma(a).
219 * \f]
221 * If \a a is non-positive integer, the limiting value should
222 * be used, but this is not yet implemented!
224 int m,
225 qpms_csf_result *result);
227 /// Exponential integral for complex second argument.
228 /** If x is (almost) positive real, it just uses gsl_sf_expint_En_e(). */
229 int complex_expint_n_e(int n, complex double x, qpms_csf_result *result);
232 /// Hypergeometric 2F2, used to calculate some errors.
233 int hyperg_2F2_series(double a, double b, double c, double d, double x,
234 gsl_sf_result *result);
236 #if 0
237 // The integral from (4.6); maybe should be static and not here.
238 int ewald32_sr_integral(double r, double k, double n, double eta, double *result, double *err, gsl_integration_workspace *workspace);
239 #endif
241 /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
242 /** \f[ \Delta_n = \int_n^\infty t^{-1/2 - n} \exp(-t + z^2/(4t))\ud t \f]
244 * \bug The current choice of method, based purely on the value of \a z, might be
245 * unsuitable especially for big values of \a maxn.
248 void ewald3_2_sigma_long_Delta(complex double *target, double *target_err, int maxn, complex double x,
249 int xbranch, complex double z);
251 /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
252 /** This function always uses Kambe's (corrected) recurrent formula.
253 * For production, use ewald3_2_sigma_long_Delta() instead.
255 void ewald3_2_sigma_long_Delta_recurrent(complex double *target, double *target_err, int maxn, complex double x,
256 int xbranch, complex double z, _Bool bigimz);
258 /// The Delta_n factor from [Kambe II], Appendix 3, used in 2D-in-3D long range sum.
259 /** This function always uses Taylor expansion in \a z.
260 * For production, use ewald3_2_sigma_long_Delta() instead.
262 * \bug The error estimate seems to be wrong (too small) at least in some cases: try
263 * parameters maxn = 40, z = 0.5, x = -3. This might be related to the exponential growth
264 * of the error.
266 void ewald3_2_sigma_long_Delta_series(complex double *target, double *target_err, int maxn, complex double x,
267 int xbranch, complex double z);
269 // General functions acc. to [2], sec. 4.6 – currently valid for 2D and 1D lattices in 3D space
271 /// The Ewald sum "self-interaction" term that appears in the lattice sums with zero (direct-space) Bravais lattice displacement.
272 int ewald3_sigma0(complex double *result, ///< Pointer to save the result (single complex double).
273 double *err, ///< Pointer to save the result error estimate (single double).
274 const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
275 double eta, ///< Ewald parameter.
276 complex double wavenumber ///< Wavenumber of the background medium.
279 /// Short-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{S}(\vect k,\vect s)\f$.
280 int ewald3_sigma_short(
281 complex double *target_sigmasr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{S} \f$, must be `c->nelem_sc` long.
282 double *target_sigmasr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
283 const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
284 double eta, ///< Ewald parameter.
285 complex double wavenumber, ///< Wavenumber of the background medium.
286 /// Lattice dimensionality.
287 /** Ignored apart from asserts and possible optimisations, as the SR formula stays the same. */
288 LatticeDimensionality latdim,
289 /// Lattice point generator for the direct Bravais lattice.
290 /** There is a possibility that the whole PGen is not consumed
291 * (this might happen if the summand start to be consistently smaller
292 * than the (partial) sums * DBL_EPSILON.
293 * In such case, it is the responsibility of the caller to deallocate
294 * the generator.
296 PGen *pgen_R,
297 /// Indicates whether pgen_R already generates shifted points.
298 /** If false, the behaviour corresponds to the old ewald32_sigma_short_points_and_shift(),
299 * so the function assumes that the generated points correspond to the unshifted Bravais lattice,
300 * and adds particle_shift to the generated points before calculations.
301 * If true, it assumes that they are already shifted (if calculating interaction between
302 * different particles in the unit cell).
304 bool pgen_generates_shifted_points,
305 /// Wave vector \f$\vect k\f$.
306 cart3_t k,
307 /// Lattice offset \f$\vect s\f$ wrt. the Bravais lattice.
308 cart3_t particle_shift
311 /// Long-range part of outgoing scalar spherical wavefunctions' lattice sum \f$ \sigma_{l,m}^\mathrm{L}(\vect k,\vect s)\f$.
312 int ewald3_sigma_long( // calls ewald3_21_sigma_long or ewald3_3_sigma_long, depending on latdim
313 complex double *target_sigmalr_y, ///< Target array for \f$ \sigma_{l,m}^\mathrm{L} \f$, must be `c->nelem_sc` long.
314 double *target_sigmalr_y_err, ///< Target array for error estimates, must be `c->nelem_sc` long or `NULL`.
315 const qpms_ewald3_constants_t *c, ///< Constant factors structure initialised by qpms_ewald3_constants_init().
316 double eta, ///< Ewald parameter.
317 complex double wavenumber, ///< Wavenumber of the background medium.
318 double unitcell_volume, ///< Volume of the (direct lattice) unit cell (with dimension corresponding to the lattice dimensionality).
319 /// Lattice dimensionality.
320 LatticeDimensionality latdim,
321 /// Lattice point generator for the reciprocal lattice.
322 /** There is a possibility that the whole PGen is not consumed
323 * (this might happen if the summand start to be consistently smaller
324 * than the (partial) sums * DBL_EPSILON.
325 * In such case, it is the responsibility of the caller to deallocate
326 * the generator.
328 PGen *pgen_K,
329 /// Indicates whether pgen_K already generates shifted points.
330 /** If false, the behaviour corresponds to the old ewald32_sigma_long_points_and_shift(),
331 * so the function assumes that the generated points correspond to the unshifted reciprocal Bravais lattice,
332 * and adds beta to the generated points before calculations.
333 * If true, it assumes that they are already shifted.
335 bool pgen_generates_shifted_points,
336 /// Wave vector \f$\vect k\f$.
337 cart3_t k,
338 /// Lattice offset \f$\vect s\f$ wrt. the Bravais lattice.
339 cart3_t particle_shift
342 #endif //EWALD_H