WIP Ewald 2D in 3D general z != 0 constant factors.

Former-commit-id: 787689f357bd8670948ba8ce7d8dc1205ca77d0f
This commit is contained in:
Marek Nečada 2020-05-29 15:55:52 +03:00
parent 975d23b557
commit 97977dbb46
4 changed files with 688 additions and 1 deletions

View File

@ -0,0 +1,245 @@
#LyX 2.4 created this file. For more info see https://www.lyx.org/
\lyxformat 584
\begin_document
\begin_header
\save_transient_properties true
\origin unavailable
\textclass article
\use_default_options true
\maintain_unincluded_children false
\language finnish
\language_package default
\inputencoding utf8
\fontencoding auto
\font_roman "default" "default"
\font_sans "default" "default"
\font_typewriter "default" "default"
\font_math "auto" "auto"
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_roman_osf false
\font_sans_osf false
\font_typewriter_osf false
\font_sf_scale 100 100
\font_tt_scale 100 100
\use_microtype false
\use_dash_ligatures true
\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize default
\use_hyperref false
\papersize default
\use_geometry false
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine basic
\cite_engine_type default
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\use_minted 0
\use_lineno 0
\index Index
\shortcut idx
\color #008000
\end_index
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\is_math_indent 0
\math_numbering_side default
\quotes_style english
\dynamic_quotes 0
\papercolumns 1
\papersides 1
\paperpagestyle default
\tablestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header
\begin_body
\begin_layout Standard
\lang english
\begin_inset FormulaMacro
\newcommand{\uoft}[1]{\mathfrak{F}#1}
\end_inset
\begin_inset FormulaMacro
\newcommand{\uaft}[1]{\mathfrak{\mathbb{F}}#1}
\end_inset
\begin_inset FormulaMacro
\newcommand{\usht}[2]{\mathbb{S}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\bsht}[2]{\mathrm{S}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\sgn}{\operatorname{sgn}}
{\mathrm{sgn}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\pht}[2]{\mathfrak{\mathbb{H}}_{#1}#2}
\end_inset
\begin_inset FormulaMacro
\newcommand{\vect}[1]{\mathbf{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ud}{\mathrm{d}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\basis}[1]{\mathfrak{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\dc}[1]{Ш_{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\rec}[1]{#1^{-1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\recb}[1]{#1^{\widehat{-1}}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ints}{\mathbb{Z}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\nats}{\mathbb{N}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\reals}{\mathbb{R}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ush}[2]{Y_{#1,#2}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\hgfr}{\mathbf{F}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\hgf}{F}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ghgf}[2]{\mbox{}_{#1}F_{#2}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ghgfr}[2]{\mbox{}_{#1}\mathbf{F}_{#2}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\ph}{\mathrm{ph}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\kor}[1]{\underline{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\koru}[1]{\utilde{#1}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\swv}{\mathscr{H}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\expint}{\mathrm{E}}
\end_inset
\end_layout
\begin_layout Standard
\lang english
\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)!}\times\nonumber \\
& & \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}\nonumber \\
& = & -\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\nonumber \\
& & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\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}\nonumber \\
& = & -\frac{i^{n+1}}{k^{2}\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\times\nonumber \\
& & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\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}\label{eq:2D Ewald in 3D long-range part}
\end{eqnarray}
\end_inset
For
\begin_inset Formula $z\ne0$
\end_inset
\begin_inset Formula
\begin{align*}
& =-\frac{i^{n+1}}{k^{2}\mathscr{A}}\sqrt{\pi}2\left(\left(n-m\right)/2\right)!\left(\left(n+m\right)/2\right)!\\
& \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}Y_{n}^{m}\left(\frac{\pi}{2},\phi_{\vect{\beta}_{pq}}\right)\sum_{j=0}^{n-\left|m\right|}\frac{\Delta_{npq}}{j!}\left(-1\right)^{j}\left(\gamma_{pq}\right)^{2j-1}\sum_{s\overset{*}{=}j}^{\min(2j,n-\left|m\right|)}\binom{j}{2j-s}\frac{\left(-\kappa z\right)^{2j-s}\left(\beta_{pq}/k\right)^{n-s}}{\left(\frac{1}{2}\left(n-m-s\right)\right)!\left(\frac{1}{2}\left(n+m-s\right)\right)!}
\end{align*}
\end_inset
\end_layout
\end_body
\end_document

375
notes/kambe_linton_dict.lyx Normal file
View File

@ -0,0 +1,375 @@
#LyX 2.4 created this file. For more info see https://www.lyx.org/
\lyxformat 584
\begin_document
\begin_header
\save_transient_properties true
\origin unavailable
\textclass article
\use_default_options true
\maintain_unincluded_children false
\language finnish
\language_package default
\inputencoding utf8
\fontencoding auto
\font_roman "default" "default"
\font_sans "default" "default"
\font_typewriter "default" "default"
\font_math "auto" "auto"
\font_default_family default
\use_non_tex_fonts false
\font_sc false
\font_roman_osf false
\font_sans_osf false
\font_typewriter_osf false
\font_sf_scale 100 100
\font_tt_scale 100 100
\use_microtype false
\use_dash_ligatures true
\graphics default
\default_output_format default
\output_sync 0
\bibtex_command default
\index_command default
\paperfontsize default
\use_hyperref false
\papersize default
\use_geometry false
\use_package amsmath 1
\use_package amssymb 1
\use_package cancel 1
\use_package esint 1
\use_package mathdots 1
\use_package mathtools 1
\use_package mhchem 1
\use_package stackrel 1
\use_package stmaryrd 1
\use_package undertilde 1
\cite_engine basic
\cite_engine_type default
\use_bibtopic false
\use_indices false
\paperorientation portrait
\suppress_date false
\justification true
\use_refstyle 1
\use_minted 0
\use_lineno 0
\index Index
\shortcut idx
\color #008000
\end_index
\secnumdepth 3
\tocdepth 3
\paragraph_separation indent
\paragraph_indentation default
\is_math_indent 0
\math_numbering_side default
\quotes_style english
\dynamic_quotes 0
\papercolumns 1
\papersides 1
\paperpagestyle default
\tablestyle default
\tracking_changes false
\output_changes false
\html_math_output 0
\html_css_as_file 0
\html_be_strict false
\end_header
\begin_body
\begin_layout Standard
\lang english
\begin_inset FormulaMacro
\newcommand{\vect}[1]{\mathbf{#1}}
\end_inset
\lang finnish
\begin_inset FormulaMacro
\newcommand{\Kambe}[1]{#1^{\mathrm{K}}}
\end_inset
\begin_inset FormulaMacro
\newcommand{\Linton}[1]{#1^{\mathrm{L}}}
\end_inset
\end_layout
\begin_layout Standard
Here and in Kambe's papers,
\begin_inset Formula $\kappa$
\end_inset
is the wavenumber (
\begin_inset Formula $k$
\end_inset
in Linton).
Here
\begin_inset Formula $\vect K_{p}$
\end_inset
is a point of the reciprocal lattice (
\begin_inset Formula $\vect K_{p}=\Kambe{\vect K_{pt}}=\Linton{\vect{\beta}_{\mu}}$
\end_inset
)
\end_layout
\begin_layout Section
\begin_inset Quotes eld
\end_inset
Gammas
\begin_inset Quotes erd
\end_inset
\end_layout
\begin_layout Standard
For
\begin_inset Formula $\kappa$
\end_inset
positive,
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\Kambe{\Gamma_{p}}\equiv\begin{cases}
\sqrt{\kappa^{2}-\left|\vect K_{p}\right|^{2}} & \kappa^{2}-\left|\vect K_{p}\right|^{2}>0\\
i\sqrt{\left|\vect K_{p}\right|^{2}-\kappa^{2}} & \kappa^{2}-\left|\vect K_{p}\right|^{2}<0
\end{cases}
\]
\end_inset
\begin_inset Formula
\[
\Linton{\gamma_{\mu}}\equiv\begin{cases}
\sqrt{\left(\frac{\vect K_{p}}{\kappa}\right)^{2}-1} & \kappa-\left|\vect K_{p}\right|\le0\\
-i\sqrt{1-\left(\frac{\vect K_{p}}{\kappa}\right)^{2}} & \kappa-\left|\vect K_{p}\right|>0
\end{cases}
\]
\end_inset
hence
\begin_inset Formula
\[
\Kambe{\Gamma_{p}}=-i\kappa\Linton{\gamma_{\mu}},
\]
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\[
\Linton{\gamma_{\mu}}=i\frac{\Kambe{\Gamma_{p}}}{\kappa}.
\]
\end_inset
\end_layout
\begin_layout Section
D vs sigma
\end_layout
\begin_layout Standard
In-plane sums [Linton 2009, (4.5)], replacing
\begin_inset Formula $n,m\rightarrow L,M$
\end_inset
,
\begin_inset Formula $k\rightarrow\kappa$
\end_inset
\end_layout
\begin_layout Standard
\lang english
\begin_inset Formula
\begin{eqnarray*}
\sigma_{L}^{M(1)} & = & -\frac{i^{L+1}}{2\kappa^{2}\mathscr{A}}\left(-1\right)^{\left(L+M\right)/2}\sqrt{\left(2L+1\right)\left(L-M\right)!\left(L+M\right)!}\times\\
& & \times\sum_{\vect K_{pq}\in\Lambda^{*}}^{'}\sum_{j=0}^{\left[\left(L-\left|M\right|/2\right)\right]}\frac{\left(-1\right)^{j}\left(\beta_{pq}/2\kappa\right)^{L-2j}e^{iM\phi_{\vect{\beta}_{pq}}}\Gamma_{j,pq}}{j!\left(\frac{1}{2}\left(L-M\right)-j\right)!\left(\frac{1}{2}\left(L+M\right)-j\right)!}\left(\frac{\gamma_{pq}}{2}\right)^{2j-1}
\end{eqnarray*}
\end_inset
[Kambe II, (3.17)], replacing
\lang finnish
\begin_inset Formula $n\rightarrow j$
\end_inset
\lang english
,
\lang finnish
\begin_inset Formula $A\rightarrow\mathscr{A}$
\end_inset
,
\begin_inset Formula $\vect K_{pt}\to\vect K_{p}$
\end_inset
,
\begin_inset Formula $\Gamma\left(\frac{1}{2}-j,e^{-i\pi}\Gamma_{p}^{2}\omega/2\right)\to\Gamma_{j,p}$
\end_inset
and performing little typographic modifications
\lang english
\begin_inset Formula
\begin{align*}
D_{LM} & =-\frac{1}{\mathscr{A}\kappa}i^{\left|M\right|+1}2^{-L}\sqrt{\left(2L+1\right)\left(L+\left|M\right|\right)!\left(L-\left|M\right|\right)!}\times\\
& \quad\times\sum_{p}e^{i\vect K_{p}\cdot\vect c_{ijt}}e^{-iM\phi_{K_{p}}}\sum_{j=0}^{\left(L-\left|M\right|\right)/2}\frac{\left(\Gamma_{p}/\kappa\right)^{2j-1}\left(K_{p}/\kappa\right)^{L-2j}\Gamma_{j,p}}{j!\left(\frac{1}{2}\left(L-\left|M\right|\right)-j\right)!\left(\frac{1}{2}\left(L+\left|M\right|\right)-j\right)!}
\end{align*}
\end_inset
Using the relations between
\begin_inset Formula $\Kambe{\Gamma_{p}}=-i\kappa\Linton{\gamma_{\mu}}$
\end_inset
, we have (also, we replace the
\begin_inset Formula $\mu$
\end_inset
index with
\begin_inset Formula $p$
\end_inset
)
\begin_inset Formula
\begin{align*}
D_{LM} & =-\frac{1}{\mathscr{A}\kappa}i^{\left|M\right|+1}2^{-L}\sqrt{\left(2L+1\right)\left(L+\left|M\right|\right)!\left(L-\left|M\right|\right)!}\times\\
& \quad\times\sum_{p}e^{i\vect K_{p}\cdot\vect c_{ijt}}e^{-iM\phi_{K_{p}}}\sum_{j=0}^{\left(L-\left|M\right|\right)/2}\frac{\left(-i\gamma_{p}\right)^{2j-1}\left(K_{p}/\kappa\right)^{L-2j}\Gamma_{j,p}}{j!\left(\frac{1}{2}\left(L-\left|M\right|\right)-j\right)!\left(\frac{1}{2}\left(L+\left|M\right|\right)-j\right)!}
\end{align*}
\end_inset
and now, trying to make the exponents look the same as in Linton,
\begin_inset Formula $2^{-1}2^{2j-L}2^{1-2j}=2^{-L}$
\end_inset
(OK),
\begin_inset Formula $K_{p}^{L-2j}=K_{p}^{L-2j}$
\end_inset
(OK),
\begin_inset Formula
\begin{align*}
D_{LM} & =-\frac{1}{2\kappa\mathscr{A}}i^{\left|M\right|+1}\sqrt{\left(2L+1\right)\left(L+\left|M\right|\right)!\left(L-\left|M\right|\right)!}\times\\
& \quad\times\sum_{p}e^{i\vect K_{p}\cdot\vect c_{ij}}e^{-iM\phi_{K_{p}}}\sum_{j=0}^{\left(L-\left|M\right|\right)/2}\frac{\left(-i\right)^{2j-1}\left(K_{p}/2\kappa\right)^{L-2j}\Gamma_{j,p}}{j!\left(\frac{1}{2}\left(L-\left|M\right|\right)-j\right)!\left(\frac{1}{2}\left(L+\left|M\right|\right)-j\right)!}\left(\frac{\gamma_{p}}{2}\right)^{2j-1}
\end{align*}
\end_inset
There are now these differences left:
\end_layout
\begin_layout Itemize
\lang english
Additional
\begin_inset Formula $\kappa$
\end_inset
factor in
\begin_inset Formula $D_{LM}$
\end_inset
\end_layout
\begin_layout Itemize
\lang english
\begin_inset Formula $i^{L+1}\left(-1\right)^{\left(L+M\right)/2}\left(-1\right)^{j}$
\end_inset
vs.
\begin_inset Formula $i^{\left|M\right|+1}\left(-i\right)^{2j-1}$
\end_inset
\end_layout
\begin_layout Itemize
\lang english
Opposite phase in the angular part.
\end_layout
\begin_layout Itemize
\lang english
Plane wave factor in
\begin_inset Formula $D_{LM}$
\end_inset
\end_layout
\begin_layout Standard
\lang english
Let's look at the
\begin_inset Formula $i,-1$
\end_inset
factors (note that
\begin_inset Formula $L+M$
\end_inset
is odd):
\begin_inset Formula $\left(-i\right)^{2j}=\left(-1\right)^{j},$
\end_inset
leaving
\begin_inset Formula $i^{L+1}\left(-1\right)^{\left(L+M\right)/2}$
\end_inset
vs.
\begin_inset Formula $i^{\left|M\right|+1}i$
\end_inset
.
So there is might be a phase difference due to different conventions, but
it does not depend on
\begin_inset Formula $j$
\end_inset
, so one should be able to transplant the
\begin_inset Formula $z\ne0$
\end_inset
sum from Kambe without major problems.
\end_layout
\begin_layout Section
Ewald parameter (integration limits)
\end_layout
\end_body
\end_document

View File

@ -26,11 +26,12 @@
#endif
#ifndef M_SQRTPI
#define M_SQRTPI 1.7724538509055160272981674833411452
#define M_SQRTPI 1.7724538509055160272981674833411452L
#endif
// sloppy implementation of factorial
// We prefer to avoid tgamma/lgamma, as their errors are about 4+ bits
static inline double factorial(const int n) {
assert(n >= 0);
if (n < 0)
@ -45,6 +46,44 @@ static inline double factorial(const int n) {
return tgamma(n + 1); // hope it's precise and that overflow does not happen
}
// sloppy implementation of double factorial n!!
static inline double double_factorial(int n) {
assert(n >= 0);
if (n <= 25) {
double fac = 1;
while (n > 0) {
fac *= n;
n -= 2;
}
return fac;
} else {
if (n % 2) { // odd, (2*k - 1)!! = 2**k * Γ(k + 0.5) / sqrt(п)
const int k = n / 2 + 1;
return pow(2, k) * tgamma(k + 0.5) / M_SQRTPI;
} else { // even, n!! = 2**(n/2) * (n/2)!
const int k = n/2;
return pow(2, k) * factorial(k);
}
}
}
// sloppy implementation of (n/2)! = Γ(n/2 + 1)
// It is _slightly_ more precise than direct call of tgamma for small odd n
static inline double factorial_of_half(const int n2) {
assert(n2 >= 0);
if (n2 % 2 == 0) return factorial(n2/2);
else {
if (n2 <= 50) { // odd, use (k - 0.5)! = Γ(k + 0.5) = 2**(-k) (2*k - 1)!! sqrt(п) for small n2
const int k = n2 / 2 + 1;
double fac2 = 1;
for(int j = 2*k - 1; j > 0; j -= 2)
fac2 *= j;
return fac2 * pow(2, -k) * M_SQRTPI;
}
else return tgamma(1. + 0.5*n2);
}
}
static inline complex double csq(complex double x) { return x * x; }
static inline double sq(double x) { return x * x; }

View File

@ -80,6 +80,34 @@ typedef struct qpms_ewald3_constants_t {
* 2D sum with additional factor of
* \f$ \sqrt{pi} \kappa \gamma(\abs{\vect{k}+\vect{K}}/\kappa) \f$.
*/
///=============== NEW GENERATION GENERAL 2D-IN-3D, including z != 0 =========================
// TODO indexing mechanisms
/// The constant factors for the long range part of a 2D Ewald sum.
complex double **S1_constfacs; // indices [y][j] where j is same as in [1, (4.5)]
/* These are the actual numbers now: (in the EWALD32_CONSTANTS_AGNOSTIC version)
* for m + n EVEN:
*
* s1_constfacs[y(m,n)][x(j,s)] =
*
* -2 * I**(n+1) * sqrt(π) * ((n-m)/2)! * ((n+m)/2)! * (-1)**j / j \
* ----------------------------------------------------------- | |
* j! * ((n - m - s)/2)! * ((n + m - s)/2)! \ 2j - s /
*
* for m + n ODD:
*
* s1_constfacs[y(m,n)][j] = 0
*/
complex double *S1_constfacs_base; ///< Internal pointer holding memory for the 2D Ewald sum constant factors.
/// The constant factors for the long range part of a 1D Ewald sum along the \a z axis.
/** If the summation points lie along a different direction, use the formula for
* 2D sum with additional factor of
* \f$ \sqrt{pi} \kappa \gamma(\abs{\vect{k}+\vect{K}}/\kappa) \f$.
*/
complex double **s1_constfacs_1Dz;
/* These are the actual numbers now:
* s1_constfacs_1Dz[n][j] =