Add section on pseudos and spherical harmonics to theory notes

This commit is contained in:
Matteo Giantomassi 2018-03-30 19:38:53 +02:00
parent 99938a8066
commit c0e8cbe9e6
12 changed files with 919 additions and 885 deletions

View File

@ -945,31 +945,31 @@ The bibtex file is available [here](../abiref.bib).
url = ""
if a.text is None: a.text = fragment
else:
if name.startswith("lesson_"):
self.warn("lesson_NAME is DEPRECATED, use lesson:name. %s in %s is deprecated" % (token, page_rpath))
# Handle [[lesson_gw1|text]]
url = "/tutorials/%s" % name.replace("lesson_" , " ", 1).strip()
if a.text is None: a.text = name
html_classes.append("lesson-wikilink")
#if name.startswith("lesson_"):
# self.warn("lesson_NAME is DEPRECATED, use lesson:name. %s in %s is deprecated" % (token, page_rpath))
# # Handle [[lesson_gw1|text]]
# url = "/tutorials/%s" % name.replace("lesson_" , " ", 1).strip()
# if a.text is None: a.text = name
# html_classes.append("lesson-wikilink")
elif name.startswith("topic_"):
# Handle [[topic_SelfEnergy|text]]
self.warn("topic_NAME is DEPRECATED, use topic:name. %s in %s is deprecated" % (token, page_rpath))
name = name.replace("topic_" , " ", 1).strip()
url = "/topics/%s" % name
if a.text is None: a.text = "%s topic" % name
html_classes.append("topic-wikilink")
add_popover(a, content=self.howto_topic[name])
#elif name.startswith("topic_"):
# # Handle [[topic_SelfEnergy|text]]
# self.warn("topic_NAME is DEPRECATED, use topic:name. %s in %s is deprecated" % (token, page_rpath))
# name = name.replace("topic_" , " ", 1).strip()
# url = "/topics/%s" % name
# if a.text is None: a.text = "%s topic" % name
# html_classes.append("topic-wikilink")
# add_popover(a, content=self.howto_topic[name])
elif name.startswith("help_"):
self.warn("help_NAME is DEPRECATED, use help:name. %s in %s is deprecated" % (token, page_rpath))
# Handle [[help_abinit|text]]
code = name.replace("help_" , " ", 1).strip()
url = "/guide/%s" % code
if a.text is None: a.text = "%s help file" % code
html_classes.append("user-guide-wikilink")
#elif name.startswith("help_"):
# self.warn("help_NAME is DEPRECATED, use help:name. %s in %s is deprecated" % (token, page_rpath))
# # Handle [[help_abinit|text]]
# code = name.replace("help_" , " ", 1).strip()
# url = "/guide/%s" % code
# if a.text is None: a.text = "%s help file" % code
# html_classes.append("user-guide-wikilink")
elif "@" in name:
if "@" in name:
# Handle [[dipdip@anaddb|text]]
vname, code = name.split("@")
var = self.codevars[code][vname]
@ -985,17 +985,17 @@ The bibtex file is available [here](../abiref.bib).
if a.text is None:
a.text = var.name if not var.is_internal else "%%%s" % var.name
elif name in self.bib_data.entries:
# Handle citation
self.warn("DEPRECATED citation without `cite` --> %s in %s is deprecated" % (token, page_rpath))
ref = self.bib_data.entries[name]
url = "/theory/bibliography#%s" % self.slugify(name)
content = ref.fields.get("title", "Unknown")
if content == "Unknown":
self.warn("Entry for %s does not provide title" % name)
add_popover(a, content=content) #+ "\n\n" + ref.authors
if a.text is None: a.text = "[%s]" % name
html_classes.append("citation-wikilink")
#elif name in self.bib_data.entries:
# # Handle citation
# self.warn("DEPRECATED citation without `cite` --> %s in %s is deprecated" % (token, page_rpath))
# ref = self.bib_data.entries[name]
# url = "/theory/bibliography#%s" % self.slugify(name)
# content = ref.fields.get("title", "Unknown")
# if content == "Unknown":
# self.warn("Entry for %s does not provide title" % name)
# add_popover(a, content=content) #+ "\n\n" + ref.authors
# if a.text is None: a.text = "[%s]" % name
# html_classes.append("citation-wikilink")
elif name.startswith("tests/") or name.startswith("~abinit/tests/"):
assert fragment is None
@ -1079,9 +1079,9 @@ The bibtex file is available [here](../abiref.bib).
if a.text is None: a.text = "%s_%s" % (namespace, name)
add_popover(a, content=self.howto_topic[name])
elif namespace in ("bib", "cite"):
if namespace == "bib":
self.warn("%s in %s is deprecated" % (token, page_rpath))
elif namespace == "cite":
#if namespace == "bib":
# self.warn("%s in %s is deprecated" % (token, page_rpath))
# Handle [[bib:biblio|bibliography]]
if name == "biblio":
url = "/theory/bibliography/"
@ -1101,9 +1101,9 @@ The bibtex file is available [here](../abiref.bib).
(exc.__class__, str(exc), token, page_rpath))
url, a.text = "FAKE_URL", "FAKE_URL"
elif namespace in ("theorydoc", "theory"):
if namespace == "theorydoc":
self.warn("%s in %s is deprecated" % (token, page_rpath))
elif namespace == "theory":
#if namespace == "theorydoc":
# self.warn("%s in %s is deprecated" % (token, page_rpath))
# TODO theorydoc --> theory
# Handle [[theorydoc:mbpt|text]]
url = "/theory/%s" % name
@ -1140,7 +1140,7 @@ The bibtex file is available [here](../abiref.bib).
html_classes.append("abifile-wikilink")
elif namespace == "ac":
# Handle [[ac:abiref_gnu_5.3_debug.ac]]
# Handle [[ac:abiref_gnu_5.3_debug.ac]]
# The following is incorrect : files in /build/config-examples are generated when */*/makemake is issued.
# url = "/build/config-examples/%s" % name
# By contrast, the following is a permanent reference

View File

@ -98,6 +98,7 @@
"ra": "{\\rangle}",
"la": "{\\langle}",
// matchcal version
mcG: "{\\mathcal{G}}",
mcR: "{\\mathcal{R}}",
mcS: "{\\mathcal{S}}",
mcC: "{\\mathcal{C}}",

View File

@ -3627,3 +3627,84 @@ author = "Philip B. Allen and Božidar Mitrović"
year = {1975},
month = aug,
}
@article{Hobbs2000,
author = {Hobbs, D. and Kresse, G. and Hafner, J.},
doi = {10.1103/physrevb.62.11556},
number = {17},
pages = {11556-11570},
source = {Crossref},
url = {http://dx.doi.org/10.1103/physrevb.62.11556},
volume = {62},
journal = {Phys. Rev. B},
publisher = {American Physical Society (APS)},
title = {Fully unconstrained noncollinear magnetism within the projector augmented-wave method},
issn = {0163-1829, 1095-3795},
year = {2000},
month = nov,
}
@article{Perdew1992,
author = {Perdew, John P. and Chevary, J. A. and Vosko, S. H. and Jackson, Koblar A. and Pederson, Mark R. and Singh, D. J. and Fiolhais, Carlos},
doi = {10.1103/physrevb.46.6671},
number = {11},
pages = {6671-6687},
source = {Crossref},
url = {http://dx.doi.org/10.1103/physrevb.46.6671},
volume = {46},
journal = {Phys. Rev. B},
publisher = {American Physical Society (APS)},
title = {Atoms, molecules, solids, and surfaces: {Applications} of the generalized gradient approximation for exchange and correlation},
issn = {0163-1829, 1095-3795},
year = {1992},
month = sep,
}
@article{Blanco1997,
title = "Evaluation of the rotation matrices in the basis of real spherical harmonics",
journal = "Journal of Molecular Structure: THEOCHEM",
volume = "419",
number = "1",
pages = "19 - 27",
year = "1997",
issn = "0166-1280",
doi = "https://doi.org/10.1016/S0166-1280(97)00185-1",
url = "http://www.sciencedirect.com/science/article/pii/S0166128097001851",
author = "Miguel A. Blanco and M. Flórez and M. Bermejo",
keywords = "Wigner functions, Real spherical harmonics, Electronic structure calculations, Rotation matrices"
}
@book{bassani1975,
title={Electronic states and optical transitions in solids},
author={Bassani, G.F. and Parravicini, G.P.},
isbn={9780080168463},
lccn={gb75022672},
series={International series of monographs in the science of the solid state},
url={https://books.google.be/books?id=cGh5AAAAIAAJ},
year={1975},
publisher={Pergamon Press}
}
@Article{Bhagavantam1964,
author="Bhagavantam, S. and Pantulu, P. V.",
title="Crystal symmetry and physical properties: Application of group theory",
journal="Proceedings of the Indian Academy of Sciences - Section A",
year="1964",
month="Jul",
day="01",
volume="60",
number="1",
pages="1--10",
abstract="A group theoretical method given by one of us in earlier publications for studying the effect of crystal symmetry on physical properties is now extended to cover the galvanomagnetic, thermomagnetic and piezo-galvanomagnetic effects in single crystals. These effects have gained importance in semi-conductor physics and such studies are therefore of current interest and may be expected to yield valuable new information.",
issn="0370-0089",
doi="10.1007/BF03046362",
url="https://doi.org/10.1007/BF03046362"
}
@book{Landau1984,
title={Electrodynamics of continuous media: },
author={Landau, L.D. and Lifshit͡s, E.M. and Lifshit︠s︡, E.M. and Pitaevski{\u\i}, L.P.},
isbn={9780080302751},
lccn={83024997},
series={Pergamon international library of science, technology, engineering, and social studies},
url={https://books.google.be/books?id=j7nvAAAAMAAJ},
year={1984},
publisher={Pergamon}
}

Binary file not shown.

Before

Width:  |  Height:  |  Size: 1.3 KiB

View File

@ -29,8 +29,8 @@ related to the irreducible polarizability by the following relation
from which the inverse dielectric function is obtained via matrix inversion.
The _**macroscopic**_ dielectric function, _є M_LF(ω), can be directly related
to the inverse of a single element, the first ( **G** _1_ =0, **G** _2_ =0),
The _**macroscopic**_ dielectric function, $\ee_M^{\text{LF}}(\ww)$, can be directly related
to the inverse of a single element, the first ($\GG_1 = 0, \GG_2 =0$) element
of the inverse of the _**microscopic**_ dielectric matrix by means of:
\begin{equation}\label{eq:abs_LFE}
@ -39,8 +39,8 @@ of the inverse of the _**microscopic**_ dielectric matrix by means of:
The microscopic dielectric matrix is the one usually calculated within the RPA in
a GW calculation. The optical absorption spectrum is simply given by the
imaginary part of _є M_LF(ω). Sometimes, especially when comparing with
experimental spectra, the real part is simply called _є 1_ and the imaginary part _є 2_.
imaginary part of $\ee_M^{\text{LF}}(\ww)$ . Sometimes, especially when comparing with
experimental spectra, the real part is simply called $\ee_1$ and the imaginary part $\ee_2$.
Note that above equation differs from
@ -49,12 +49,15 @@ Note that above equation differs from
\end{equation}
due to the so called local-field effects introduced by the presence of the crystalline environment.
The former quantity _with_ local fields (LF) is more accurate than the latter one
with _no_ local field (NLF). The reason the two equations are different is
because in the first, the expression in the denominator is the first element
of the inverse of the _whole_ microscopic dielectric matrix. This element
depends on the _entire_ matrix and cannot simply be calculated by taking the
inverse of the first element of the non-inverted matrix.
The former quantity _with_ local fields (LF) is more accurate than the latter one with _no_ local field (NLF).
!!! important
The reason the two equations are different is
because in the first, the expression in the denominator is the first element
of the inverse of the _whole_ microscopic dielectric matrix. This element
depends on the _entire_ matrix and cannot simply be calculated by taking the
inverse of the first element of the non-inverted matrix.
In the [[theory:mbt#RPA_Fourier_space|GW_notes]], we have discussed how
to calculate the irreducible polarizability and thus the absorption spectrum
@ -75,8 +78,9 @@ The fundamental reason for the failure of the RPA can be traced back to the
neglect of the electron-hole interaction that, in the many-body language,
should be included through the vertex function. Replacing the vertex function
with a local and instantaneous function is a too crude and unrealistic
approximation for the many-body polarizability. In the next section we will
discuss how to obtain an improved approximation for the vertex and therefore
approximation for the many-body polarizability.
In the next section we will discuss how to obtain an improved approximation for the vertex and therefore
an improved approximation for the polarizability that takes into account many-
body effects due to the electron-hole interaction
@ -89,12 +93,12 @@ spectra available in ABINIT is given below.
The Bethe-Salpeter theory is formulated in terms of two-particle propagators.
These are four-point functions describing the motion through the system of two
particles at once. We will label them _L 0_ and _L_ , where the subscript zero
particles at once. We will label them $L^0$ and $L$ , where the subscript zero
indicates that the particles are non-interacting. By restricting the starting
and ending point of the particles, the two-point contraction of _L 0_ gives
and ending point of the particles, the two-point contraction of $L^0$ gives
the reducible independent-particle polarizability according to
$$ \chi^0(12) = L^0(11, 22) $$ while the two-point contraction of _L_ gives
$$ \chi^0(12) = L^0(11, 22) $$ while the two-point contraction of $L$ gives
the reducible many-body polarizability. $$ \hat\chi(12) = L(11, 22) $$ that
should not be confused with the _irreducible_ polarizability of the many-body
system (denoted with a tilde in the formulas in the [[theory:mbt#RPA_Fourier_space|GW_notes]]).
@ -104,21 +108,20 @@ _**with**_ local field effects is obtained directly from the reducible
polarizability using
\begin{equation}
\ee_M^{\text{LF}}(\ww) = 1 - \lim_{\qq
\rightarrow 0} v(\qq)\,\tchi_{00}(\qq;\ww)
\ee_M^{\text{LF}}(\ww) = 1 - \lim_{\qq \rightarrow 0} v(\qq)\,\tchi_{00}(\qq;\ww)
\end{equation}
Computing the reducible _L_ directly, if possible, thus allows one to avoid
Computing the reducible $L$ directly, if possible, thus allows one to avoid
the costly matrix inversion of the dielectric function that should otherwise
be performed for each frequency.
One can show that _L_ satisfies the Dyson-like equation:
One can show that $L$ satisfies the Dyson-like equation:
\begin{equation}
L = L^0 + L^0 K L \Longrightarrow L = \bigl [ 1 - L^0 K]^{-1} L^0
\end{equation}
that involves the Bethe-Salpeter kernel _K_ :
that involves the Bethe-Salpeter kernel $K$:
\begin{equation}\label{eq:BSE_kernel_LF}
K(1234) = \underbrace{\delta(12)\delta(34)\bar v(13)}_{Exchange} -
@ -126,7 +129,7 @@ K(1234) = \underbrace{\delta(12)\delta(34)\bar v(13)}_{Exchange} -
\end{equation}
where the bar symbol signifies that the Coulomb interaction has to be taken without its
long-range Fourier component at **G** =0:
long-range Fourier component at $\GG =0$:
\begin{equation} \begin{cases} {\bar
v(\qq)} = v(\qq) \quad {\text{if}}\; \qq \neq 0 \\\ \\\ {\bar v(\qq=0)} = 0
@ -140,11 +143,10 @@ rewriting of the kernel in this way is mostly an algebraic trick which allows
one to easily separate the contributions and calculate the optical spectrum
both with and without the LF.
The Dyson-like equation for _L_ becomes a matrix problem once a particular
The Dyson-like equation for $L$ becomes a matrix problem once a particular
basis set is chosen to expand the four-point functions involved in the
problem. In the standard approach, the polarisabilities and the BS kernel are
expressed in the so-called _transition_ space (i.e. as products of single-
particle orbitals) using:
expressed in the so-called _transition_ space (i.e. as products of single-particle orbitals) using:
\begin{equation} F(1234) =
\sum_{ \substack{(n_1 n_2) \\\ (n_3 n_4)} }
@ -158,7 +160,7 @@ and
} = \int F(1234) \psi_{n_1}(1) \psi_{n_2}^\*(2) \psi_{n_3}^\*(3)
\psi_{n_4}(4) \dd (1234) \end{equation}
where _n i_ is a shorthand notation to denote band, **k** -point and spin index.
where $n_i$ is a shorthand notation to denote band, **k**-point and spin index.
The expansion is exact provided that
the set of single-particle orbitals from a complete basis set in the Hilbert space.
@ -173,14 +175,14 @@ f_{n_1}) } { (\ee_{n_2} - \ee_{n_1} - \ww)} \delta_{n_1 n_3} \delta_{n_2 n_4}
thus leading to a considerable simplification in the mathematical formalism.
After some algebra (see [[cite:Onida2002]] for a more complete derivation) one
finds that, in a system with an energy gap, the Dyson-like equation for _L_
can be expressed in terms of an effective two-particle Hamiltonian, _H_ , according to
finds that, in a system with an energy gap, the Dyson-like equation for $L$
can be expressed in terms of an effective two-particle Hamiltonian, $H$ , according to
\begin{equation}
L = \bigl [ H-\ww \bigr]^{-1}\,F
\end{equation}
The explicit form for _H_ and _F_ in the transition space is given by
The explicit form for $H$ and $F$ in the transition space is given by
\begin{equation}
H =
@ -204,9 +206,9 @@ F =
\right)
\end{equation}
where valence states are indicated by the indices _v_ , _v'_ while _c_ , _c'_
are used for conduction bands. The _R_ sub-matrix is Hermitian and is usually
referred to as the resonant block while _C_ (the so-called coupling block) is
where valence states are indicated by the indices $v$ , $v'$ while $c$ , $c'$
are used for conduction bands. The $R$ sub-matrix is Hermitian and is usually
referred to as the resonant block while $C$ (the so-called coupling block) is
a complex symmetric sub-matrix. Note that, due to the presence of the coupling
term, the BS Hamiltonian is non-Hermitian although the operator possesses a real spectrum.
@ -222,7 +224,7 @@ W_{(vc\kk)(v'c'\kk')}
with single particle transition energies on the diagonal.
The matrix elements of v and _W_ are defined as:
The matrix elements of $v$ and $W$ are defined as:
\begin{equation}
{\bar v}_{(n_1 n_2) (n_3 n_4)} = \delta_{\sigma_1 \sigma_2}\,
@ -233,12 +235,12 @@ v(\rr-\rr')} \psi_{n_3}^*(\rr') \psi_{n_4}(\rr') \dd \rr \dd \rr'
{W(\rr,\rr',\ww=0)} \psi^*_{n_2}(\rr') \psi_{n_4}(\rr') \dd \rr \dd \rr'
\end{equation}
Note, in particular, that only the static limit ( ω = 0) of _W_ is involved in the last expression.
Note, in particular, that only the static limit ($\omega = 0$) of $W$ is involved in the last expression.
The coupling matrix elements are usually smaller than the resonant ones. This
is especially true in bulk systems due to the different spatial behavior of
conduction and valence states. In solid state calculations, it is therefore
common practice to ignore the _C_ block (the so-called Tamm-Dancoff
common practice to ignore the $C$ block (the so-called Tamm-Dancoff
approximation). Under this assumption the excitonic Hamiltonian is given by a Hermitian operator.
![](bse_assets/TDA.svg)
@ -246,8 +248,8 @@ approximation). Under this assumption the excitonic Hamiltonian is given by a He
The variable [[bs_coupling]] defines whether the coupling block should be included in the calculation.
The macroscopic dielectric function is obtained by contracting the many-body
_L_ and then taking the optical limit of the **G** = **G'** =0 component along
a particular direction in **q** -space. The final result reads:
$L$ and then taking the optical limit of the $\GG = \GG' = 0$ component along
a particular direction in $\qq$-space. The final result reads:
\begin{equation}
\ee_M(\ww) = 1 - \lim_{\qq \rightarrow 0} v(\qq)\,\langle
@ -273,7 +275,7 @@ change the default directions using the variables [[gw_nqlwl]] and [[gw_qlwl]].
At this point it is clear that the evaluation of the macroscopic dielectric
function within the BS formalism is a two-step process:
1. Construction of the _H_ matrix in the transition space.
1. Construction of the $H$ matrix in the transition space.
2. Computation of the macroscopic dielectric function using the two equations reported at the end of the previous section.
@ -291,7 +293,7 @@ For BS computations, it is common practice to simulate the self-energy
corrections by employing the scissors operator whose value can be obtained
either from experiments or from ab-initio calculations. The scissors operator
allows one to avoid a costly _GW_ calculation that should performed for all
the **k** -points and bands included in the transition space (the optional
the $\kk$-points and bands included in the transition space (the optional
path on the right indicated with yellow arrows that corresponds to [[optdriver]]=4).
The construction of the Hamiltonian matrix represents a significant portion of
@ -305,7 +307,7 @@ matrix inversion for each frequency. The variable [[bs_algorithm]] is used to
select among the different possibilities.
[[bs_algorithm]]=1 employs standard (Sca)LAPACK routines to obtain the
spectral representation of _H_ in terms of eigenvalues and right eigenvectors of _H_:
spectral representation of $H$ in terms of eigenvalues and right eigenvectors of $H$:
\begin{equation}
\begin{cases} H |\lambda\rangle = \ee_\lambda |\lambda\rangle
@ -314,7 +316,7 @@ spectral representation of _H_ in terms of eigenvalues and right eigenvectors of
\langle\lambda'| \end{cases}
\end{equation}
Then, as discussed in [2], the inverse of [ _H_ -ω] is obtained according to
Then, as discussed in [2], the inverse of $H - \ww$ is obtained according to
\begin{equation} \bigl[ H -\ww
\bigr]^{-1} = \sum_{\lambda \lambda'} |\lambda\rangle \dfrac{O_{\lambda
@ -371,7 +373,7 @@ We refer the reader to [[cite:Gruning2009]] for a generalization of the method t
the case in which the coupling term cannot be neglected. The main drawback of
the method is that it does not give direct access to the excitonic spectrum
hence it cannot be used to calculate binding energies or to plot the excitonic
wavefunctions. Moreover, for the time being, [[bs_algorithm]]=2 cannot be used
wavefunctions. Moreover, for the time being, [[bs_algorithm]] = 2 cannot be used
for calculations in which the coupling term is included.
[[bs_algorithm]]=3 employs the conjugate-gradient method to calculate the
@ -381,9 +383,8 @@ testing and does not support calculations with the coupling term.
## 4 Kernel matrix elements in reciprocal space
Our implementation employs planewaves to expand the periodic part of the Bloch
states, _u_ , and the two-point function _W(r,r')_ that becomes a **q**
-dependent matrix in reciprocal-space. The conventions used for the transforms
are documented in [this section](mbt.md#notations) of the _GW_ notes.
states, $u$ , and the two-point function $W(\rr,\rr')$ that becomes a $\qq$-dependent matrix in reciprocal-space.
The conventions used for the transforms are documented in [this section](mbt.md#notations) of the $GW$ notes.
The matrix elements of the exchange term are evaluated in reciprocal space using:
@ -414,21 +415,20 @@ calculations without local field effects. The variable [[bs_coulomb_term]]
governs the computation of the Coulomb term, the most CPU-demanding part due
to the presence of the double sum over **G** -vectors.
It is also important to stress that, in the two equations above, the **k**
-point index runs over the full Brillouin zone hence the size of the
It is also important to stress that, in the two equations above, the **k**-point
index runs over the full Brillouin zone hence the size of the
Hamiltonian is defined by the number of point in the full Brillouin zone and
not by the number of points in the irreducible wedge.
## 5 Matrix elements of the dipole operator
The accurate calculation of optical properties require the correct treatment
of the optical limit ( **G** =0, **q** -> 0) of the oscillator matrix
of the optical limit ($\GG = 0, |\qq| \rightarrow 0$) of the oscillator matrix
elements. The computation of these terms deserves some clarification, due to
the presence of the fully nonlocal pseudopotential term in the Kohn-Sham
Hamiltonian.
A linear expansion up to the first order in **q** of the oscillator matrix
element results in:
A linear expansion up to the first order in **q** of the oscillator matrix element results in:
\begin{equation}
\langle b_1,\kmq|e^{-i\qq\cdot\rr}|b_2,\kk \rangle
@ -436,9 +436,9 @@ element results in:
\rangle + \mcO(q^2)
\end{equation}
where we have assumed _b 1 ≠ b2 _and the difference between the two wave functions
at **k** - **q** and **k** has been
neglected because it only introduces terms that are quadratic in **q**.
where we have assumed $b_1 \neq b_2$ and the difference between the two wave functions
at $\kk - \qq$ and $\kk$ has been
neglected because it only introduces terms that are quadratic in $\qq$.
Unfortunately, the above expression cannot be directly used because the matrix
elements of the position operator are ill-defined when wavefunctions obey
@ -458,7 +458,7 @@ the pseudopotential is more involved and much more CPU-demanding.
The role played by this additional term is usually marginal in the case of GW
calculations: the QP corrections are obtained by performing an integration in
**q** -space and only the **q** -> 0 component of the inverse dielectric
$\qq$-space and only the $\qq \rightarrow 0$ component of the inverse dielectric
matrix is affected by the commutator of the non-local part of the pseudopotential.
For this reason it is common practice, especially during the _GW_ convergence

View File

@ -1,17 +1,17 @@
---
authors: X. Gonze, Y. Suzukawa, M. Mikami
authors: M. Giantomassi, X. Gonze, Y. Suzukawa, M. Mikami
---
# Geometric considerations
## Real space
The three primitive translation vectors are ${\bf R}_1$, ${\bf R}_2$, ${\bf R}_3$.
The three primitive translation vectors are $\RR_1$, $\RR_2$, $\RR_3$.
Their representation in Cartesian coordinates (atomic units) is:
$$ {\bf R}_1 \rightarrow {rprimd(:, 1)} $$
$$ {\bf R}_2 \rightarrow {rprimd(:, 2)} $$
$$ {\bf R}_3 \rightarrow {rprimd(:, 3)} $$
$$ \RR_1 \rightarrow {rprimd(:, 1)} $$
$$ \RR_2 \rightarrow {rprimd(:, 2)} $$
$$ \RR_3 \rightarrow {rprimd(:, 3)} $$
Related input variables: [[acell]], [[rprim]], [[angdeg]].
@ -33,7 +33,7 @@ Representation in reduced coordinates:
Related input variables: [[xangst]], [[xcart]], [[xred]]
The volume of the primitive unit cell (ucvol) is
The volume of the primitive unit cell (called *ucvol* in the code) is
\begin{eqnarray*}
\Omega &=& {\bf R}_1 \cdot ({\bf R}_2 \times {\bf R}_3)
@ -77,7 +77,7 @@ $$ {\bf R}^{met}_{ij} \rightarrow {rmet(i,j)} $$
## Reciprocal space
The three primitive translation vectors in reciprocal space are
${\bf G}_{1}$,${\bf G}_{2}$,${\bf G}_{3}$
$\GG_1$, $\GG_2$,$\GG_3$
\begin{eqnarray*}
{\bf G}_{1}&=&\frac{1}{\Omega}({\bf R}_{2}\times{\bf R}_{3}) \rightarrow {gprimd(:,1)} \\
@ -85,7 +85,7 @@ ${\bf G}_{1}$,${\bf G}_{2}$,${\bf G}_{3}$
{\bf G}_{3}&=&\frac{1}{\Omega}({\bf R}_{1}\times{\bf R}_{2}) \rightarrow {gprimd(:,3)}
\end{eqnarray*}
This definition is such that ${\bf G}_{i}\cdot{\bf R}_{j}=\delta_{ij}$
This definition is such that $\GG_i \cdot \RR_j = \delta_{ij}$
!!! warning
Often, a factor of $2\pi$ is present in definition of
@ -141,53 +141,33 @@ where ${\bf G}^{met}_{ij}$ is the metric tensor in reciprocal space:
$$ {\bf G}^{met}_{ij} \rightarrow {gmet(i,j)} $$
## Symmetries
## Fourier series for periodic lattice quantities
A symmetry operation in real space sends the point ${\bf r}$ to the point
${\bf r'}={\bf S_t}\{{\bf r}\}$ whose coordinates are
$({\bf r'})_{\alpha}=\sum_{\beta}S_{\alpha \beta}r_{\beta} + t_{\alpha}$ (Cartesian coordinates).
Any function with the periodicity of the lattice i.e any function fullfilling the property:
The symmetry operations that preserves the crystalline structure are
those that send every atom location on an atom location with the same atomic type.
$$ u(\rr + \RR) = u(\rr) $$
The application of a symmetry operation to a function of spatial coordinates $f$ is such that:
can be represented with the discrete Fourier series:
$$
({\bf S_t}f)({\bf r}) = f(({\bf S_t})^{-1}{\bf r})
$$
$$ u(\rr)= \sum_\GG u(\GG)e^{i\GG\cdot\rr} $$
$$
({\bf S_t})^{-1}({\bf r}) =
\sum_{\beta} S^{-1}_{\alpha \beta}(r_{\beta}-t_{\beta})
$$
where the Fourier coefficient, $u(\GG)$, is given by:
For each symmetry operation,$isym=1 \dots nsym$, the $3\times3$
${\bf S}^{red}$ matrix is stored in {symrel(:,:,isym)}.
$$ u(\GG) = \frac{1}{\Omega} \int_\Omega u(\rr)e^{-i\GG\cdot\rr}\dd\rr $$
[in reduced coordinates:
$r'^{red}_{\alpha}=\sum_{\beta}S^{red}_{\alpha \beta}
r^{red}_{\beta}+t^{red}_{\beta}$ ]
<!--
This appendix reports the conventions used in this work for the Fourier
representation in frequency- and momentum-space.
The volume of the unit cell is denoted with $\Omega$, while $V$ is the
total volume of the crystal simulated employing Born-von K\'arman periodic boundary condition~\cite{Ashcroft1976}.
and the vector ${\bf t}^{red}$ is stored in {tnons(:,isym)}.
\begin{equation}\label{eq:IFT_2points_convention}
f(\rr_1,\rr_2)= \frac{1}{V} \sum_{\substack{\qq \\ \GG_1 \GG_2}}
e^{i (\qq +\GG_1) \cdot \rr_1}\,f_{\GG_1 \GG_2}(\qq)\,e^{-i (\qq+\GG_2) \cdot \rr_2}
\end{equation}
The conversion between reduced coordinates and Cartesian coordinates is
$$
r'_{\gamma}=\sum_{\alpha \beta}(R_{\alpha p})_{\gamma}
[S^{red}_{\alpha \beta} r^{red}_{\beta}+t^{red}_{\alpha}]
$$
with $[{\rm as} \ G_{ip} \cdot R_{jp}=\delta_{ij}]$
$$
r_{\delta}=\sum_{\alpha}(R_{\alpha p})_\delta
r^{red}_{\alpha} \rightarrow
\sum_{\beta}(G_{\beta p})_{\delta} r_{\delta}=r^{red}_{\beta}
$$.
So
$$
S_{\gamma \delta}=\sum_{\alpha \beta}(R_{\alpha p})_{\gamma}
S^{red}_{\alpha \beta}(G_{\beta p})_{\gamma}
$$
\begin{equation}\label{eq:FT_2points_convention}
f_{\GG_1\GG_2}(\qq) = \frac{1}{V} \iint_V
e^{-i(\qq+\GG_1) \cdot \rr_1}\,f(\rr_1, \rr_2)\,e^{i (\qq+\GG_2) \cdot \rr_2}\dd\rr_1\dd\rr_2
\end{equation}
-->

View File

@ -7,7 +7,9 @@ authors: MG, MS
The aim of this section is to introduce the Green's function formalism, the
concept of self-energy and the set of coupled equations proposed by Hedin
whose self-consistent solution, in principle, gives the exact Green's function
of the interacting system. We mainly focus on the aspects of the theory that
of the interacting system.
We mainly focus on the aspects of the theory that
are important for understanding the different steps of the calculation and the
role played by the input variables used to control the run. A much more
consistent and rigorous introduction to the physical concept of Green's
@ -47,9 +49,9 @@ E_i^{(N+1)} - E_0^N & \text{if}\,\varepsilon_i > \mu \\\ E_0^N-E_i^{(N-1)} & \te
\end{cases}
\end{equation}
where _E 0N_ is the ground state energy of
the electron system with _N_ electrons, and _i_ is the set of quantum numbers
labeling the excited states with _N ±_1 electrons. Finally, _μ_ is the
where $E_0^N$ is the ground state energy of
the electron system with $N$ electrons, and $i$ is the set of quantum numbers
labeling the excited states with $N \pm 1$ electrons. Finally, $\mu$ is the
chemical potential of the system. Other important quantities that will be used
in the following are the so-called Lehmann amplitudes defined, in the Schrodinger representation, by
@ -67,8 +69,8 @@ G (\rr_1,\rr_2;\ww) = \sum_i \frac{\Psi_i(\rr_1)\Psi_i^\*(\rr_2)} {\ww-\ee_i
\end{equation}
makes it clear that, in the frequency domain, the time-ordered Green's function
contains the complete excitation spectrum corresponding to excitations of an (
_N_ -1)-particle and an ( _N_ +1)-particle system. Hence, locating the poles
contains the complete excitation spectrum corresponding to excitations of an
$(N-1)$-particle and an $(N+1)$-particle system. Hence, locating the poles
of the Green's function in the complex plane provides the information needed
to interpret those processes measured in experiments in which a single
electron is inserted to or removed from the system. The figure below gives a
@ -87,14 +89,14 @@ The Dyson equation
G(12) = \Go(12) + \int \Go(13)\,\Sigma(34)\,G(42)\dd34.
\end{equation}
establishes a connection between the fully interacting _G_ and the propagator, _G 0_, of an approximate
establishes a connection between the fully interacting $G$ and the propagator, $\Go$, of an approximate
non-interacting system through a (non-local, non-Hermitian and time dependent)
potential called the self-energy, Σ. Since _G 0_ is supposed to be known
exactly, the problem of calculating _G_ (12) has now been reduced to the calculation of the self-energy.
potential called the self-energy, $\Sigma$. Since $\Go$ is supposed to be known
exactly, the problem of calculating $G (12)$ has now been reduced to the calculation of the self-energy.
The self-energy is not a mere mathematical device used in a roundabout way to
obtain _G_ but is has a direct physical meaning. The knowledge of the self-
energy operator, allows one to describe the quantum-mechanical state of a
obtain $G$ but is has a direct physical meaning. The knowledge of the self-energy operator,
allows one to describe the quantum-mechanical state of a
renormalized electron in the many-body system by solving the quasiparticle
(QP) equation [[cite:Onida2002]]:
@ -104,7 +106,7 @@ renormalized electron in the many-body system by solving the quasiparticle
\end{equation}
The QP eigenstates so obtained can be used
to construct _G_ according to the Lehmann representation. Note that the QP
to construct $G$ according to the Lehmann representation. Note that the QP
equation departs from the Kohn Sham equation since the QP eigenvectors and
eigenvalues do have a direct physical meaning: they can be used to obtain both
the charge density of the interacting system and to describe the properties of charged excitations.
@ -113,7 +115,7 @@ the charge density of the interacting system and to describe the properties of c
In 1965 Hedin [[cite:Hedin1965]] showed how to derive a set of coupled integro-
differential equations whose self-consistent solution, in principle, gives the
exact self-energy of the system and therefore the exact _G_.
exact self-energy of the system and therefore the exact $G$.
The fundamental building blocks employed in the formalism are the irreducible polarizability:
@ -124,14 +126,14 @@ The fundamental building blocks employed in the formalism are the irreducible po
which describes the linear response of the density to changes in the total
effective potential (the superposition of the external potential plus the
internal classical Hartree potential) and the dynamically screened
interaction, _W_ , that is related to the bare Coulomb interaction, _v_ , and
interaction, $W$ , that is related to the bare Coulomb interaction, $v$ , and
to the inverse dielectric function through:
\begin{equation} \label{W_def}
W(12) \equiv \int \ee^{-1}(13)\,v(32)\dd3.
\end{equation}
The dielectric matrix є(12) is related to the irreducible polarizability _χ_ (12) by the
The dielectric matrix $\ee(12)$ is related to the irreducible polarizability $\tchi(12)$ by the
following relation:
\begin{equation}
@ -143,12 +145,12 @@ The pentagon sketched in the figure below shows how the various physical quantit
![](mbt_assets/hedin_pentagon.svg)
The polarization function renormalises the bare interaction resulting in the
screened interaction _W_ (12). The screened interaction, _W_ (12), the many-
body propagator _G_ (12), and the vertex function, Γ(12;3), which describe the
screened interaction $W (12)$. The screened interaction, $W (12)$, the many-
body propagator $G (12)$, and the vertex function, $\Gamma(12;3)$, which describe the
interactions between virtual hole and electron excitations [[cite:Onida2002]], are
the essential ingredients for the determination of Σ(12).
the essential ingredients for the determination of $\Sigma(12)$.
The iteration starts by setting _G_ = _G 0_. Then the set of equations should
The iteration starts by setting $G = \Go$. Then the set of equations should
in principle be iterated until self-consistency in all terms is reached.
## The GW approximation
@ -160,7 +162,7 @@ challenging. The set of equations can, however, be iterated assuming that only
a few iterations are actually needed to obtain physically meaningful results.
A widely used approach to the approximate solution of Hedin's equations is the
so-called _GW_ approximation [[cite:Hedin1965]], which consists in approximating
so-called $GW$ approximation [[cite:Hedin1965]], which consists in approximating
the vertex function with a local and instantaneous function:
\begin{equation}
@ -172,7 +174,7 @@ the full set of Hedin's equations, leads to a considerable simplification in the
![](mbt_assets/gw_pentagon.svg)
Thanks to the neglect of vertex corrections, the irreducible polarizability _χ_ (12) is now given by
Thanks to the neglect of vertex corrections, the irreducible polarizability $\tchi (12)$ is now given by
\begin{equation}\label{eq:RPA_with_G}
\tchi = -i\,G(12)G(21^+).
@ -182,7 +184,7 @@ which, once rewritten in terms of orbitals and energies, reduces to the RPA
expression proposed by Adler [[cite:Adler1962]] and Wiser [[cite:Wiser1963]].
In real space, the self-energy reduces to a simple direct product of the
dressed electron propagator, _G_ (12), and the dynamically screened interaction, _W_ (12):
dressed electron propagator, $G (12)$, and the dynamically screened interaction, $W (12)$:
\begin{equation}
\Sigma(12) = i G(12)\,W(1^+2).
@ -209,16 +211,16 @@ self-consistent approach would need the inclusion of some kind of vertex
correction during the solution of the equations.
For this reason, the most common approach employed in the _ab initio_
community consists of using the best available approximation for _G_ and _W_
community consists of using the best available approximation for $G$ and $W$
as a starting point and performing only a single-iteration of the
parallelogram (the so-called one-shot _GW_ method, or _G 0W0_). In this case
parallelogram (the so-called one-shot $GW$ method, or $\Go\Wo$). In this case
the self-energy is simply given by:
\begin{equation}
\Sigma(12) = i\Go^\KS(12)\Wo(1^+2)
\end{equation}
where _G0KS_(12) is the independent-particle propagator of the Kohn-Sham (KS)
where $\Go^\KS(12)$ is the independent-particle propagator of the Kohn-Sham (KS)
Hamiltonian, and the screened interaction is approximated with the RPA
calculated with KS energies and wave functions:
@ -236,11 +238,11 @@ are usually in qualitative agreement with experiment.
This observation suggests that a simple, albeit accurate, solution for the QP
energies can be obtained using first-order perturbation theory, treating the
exchange and correlation potential, _V xc_, as a zeroth-order approximation to
exchange and correlation potential, $V_{xc}$, as a zeroth-order approximation to
the non-local and energy dependent self-energy [[cite:Hybertsen1985]], [[cite:Hybertsen1986]]
Under the assumption that the QP wavefunctions equal the KS orbitals, we can
expand the self-energy operator around _ε KS_ obtaining a closed expression for _ε QP_:
expand the self-energy operator around $\ee^\KS$ obtaining a closed expression for $\ee^\QP$:
\begin{equation} \label{eq:implicit_QP_energy}
\ee^\QP = \ee^\KS + Z \langle\Psi^\KS|\Sigma(\ee^\KS)-\vxc|\Psi^\KS\rangle.
@ -267,9 +269,9 @@ particle polarizability assumes the form:
![](mbt_assets/chi0_sc_tr.svg)
where only the transitions between valence ( _v_ ) and conduction states (_c_) contribute
(for simplicity we have assumed a semiconductor with time-
reversal invariance, the conventions used for the Fourier transform are
where only the transitions between valence ($v$) and conduction states ($c$) contribute
(for simplicity we have assumed a semiconductor with time-reversal invariance,
the conventions used for the Fourier transform are
discussed in the [notation](#notations) section).
The number of bands used to compute the polarizability is specified by
@ -277,21 +279,22 @@ The number of bands used to compute the polarizability is specified by
divergences in the denominators. The frequency mesh is defined by the set of
variables [[nfreqre]], [[nfreqim]], [[freqremax]], and [[freqremin]] (a number
of more exotic grid choices are available through input variables beginning
with gw_... or cd_..., e.g. [[gw_frqim_inzgrid]]).
with `gw_...` or `cd_...`, e.g. [[gw_frqim_inzgrid]]).
_M_ is a shorthand notation to denote the matrix element of a plane wave
$M$ is a shorthand notation to denote the matrix element of a plane wave
sandwiched between two wavefunctions (i.e. oscillator matrix elements).
The number of planewaves (PW) used to describe the wavefunctions is determined by
[[ecutwfn]] while the number of _G_ -vectors used to describe the
[[ecutwfn]] while the number of $\GG$-vectors used to describe the
polarizability (i.e. the number of _G_ vectors in the oscillator matrix
elements) is determined by [[ecuteps]].
The oscillators are ubiquitous in the Many-Body part of ABINIT and their calculation
represents one of the most CPU intensive part of the execution. For this reason we
devote [[#6|section 6]] to the discussion of some important technical details concerning their computation.
In principle, the set of **q** -points in the screening matrix is given by all
In principle, the set of $\qq$-points in the screening matrix is given by all
the possible differences between two crystalline momenta of the wavefunctions
stored in the KSS file, so it is controlled by the chosen **k** -point grid.
stored in the KSS file, so it is controlled by the chosen $\kk$-point grid.
The code, however, exploits the invariance of the two-point function under the
action of any symmetry operation of the crystalline space group:
@ -300,7 +303,7 @@ action of any symmetry operation of the crystalline space group:
\bigl(\Ri(\rr_1-\tt),\Ri(\rr_2-\tt)\bigr)
\end{equation}
so that only the **q** -points in the irreducible Brillouin zone (IBZ) have to be calculated explicitly.
so that only the $\qq$-points in the irreducible Brillouin zone (IBZ) have to be calculated explicitly.
In frequency and reciprocal space, the microscopic dielectric function is
related to the irreducible polarizability by the following relation
@ -311,7 +314,7 @@ related to the irreducible polarizability by the following relation
\end{equation}
from which the inverse dielectric function is obtained via matrix inversion. Following Adler
[[cite:Adler1962]], the macroscopic dielectric function, єMLF(ω), can be directly
[[cite:Adler1962]], the macroscopic dielectric function, $\ee_M^{\text{LF}}(\ww)$, can be directly
related to the inverse of the microscopic dielectric matrix by means of:
\begin{equation} \label{eq:abs_LFE}
@ -319,7 +322,7 @@ related to the inverse of the microscopic dielectric matrix by means of:
\end{equation}
The optical absorption spectrum -- the quantity one can compare with experiments -- is
given by the imaginary part: Im[єMLF(ω)].
given by the imaginary part.
Note that the equation above differs from
@ -353,13 +356,13 @@ link and use the FFTW3 library in GW calculations (controlled by setting
[[fftalg]] 312). The performance of the various FFT libraries for a given type
of calculation can be benchmarked with the **fftprof** utility.
For a given set of indices ( _b 1_, _b 2_, **k**, **q** ), the calculation of
For a given set of indices ( $b_1$, $b_2$, $\kk$, $\qq$ ), the calculation of
the oscillator is done in four different steps:
1. The two wavefunctions in the irreducible wedge are FFT transformed from the **G** -space to the real space representation,
2. The orbitals are rotated in real space on the FFT mesh to obtain the points **k** and **k-q** in the full Brillouin zone.
1. The two wavefunctions in the irreducible wedge are FFT transformed from the $\GG$-space to the real space representation,
2. The orbitals are rotated in real space on the FFT mesh to obtain the points $\kk$ and $\kk-\qq$ in the full Brillouin zone.
3. Computation of the wavefunction product.
4. FFT transform of the product to obtain _M_
4. FFT transform of the product to obtain $M$
Each oscillator thus requires three different FFTs (two transforms to
construct the product, one FFT to get M). The number of FFTs can be
@ -371,10 +374,10 @@ time instead. This option is controlled by the second digit of the input variabl
The third term in the equation defining the oscillators makes it clear that
the product of the periodic part of the orbitals has non-zero Fourier
components in a sphere whose radius is 2 × _R wfn_ where _R wfn_ is the radius
of the **G** -sphere used for the wavefunctions (set by [[ecutwfn]]). To avoid
components in a sphere whose radius is $2 × R_{wfn}$ where $R_{wfn}$ is the radius
of the $\GG$-sphere used for the wavefunctions (set by [[ecutwfn]]). To avoid
aliasing errors in the FFT one should therefore use an FFT box that encloses
the sphere with radius 2× _R wfn_, but this leads to a significant increase in
the sphere with radius $2 × R_{wfn}$, but this leads to a significant increase in
the computing effort as well as in the memory requirements. The input variable
[[fftgw]] specifies how to setup the FFT box for the oscillators and should be
used to test how the aliasing errors affect the final results. The default
@ -428,7 +431,7 @@ expensive frequency integration (a Hilbert transform):
The price to be paid, however, is that a large table for the spectral function
has to be stored in memory and a Hilbert transform has to be performed for
each pair ( **G** 1, **G** 2). Since the computing time required for the
each pair ( $\GG_1, \GG_2$). Since the computing time required for the
transform scales quadratically with the number of vectors in the
polarizability (governed by [[ecuteps]]), the CPU time spent in this part will
overcome the computing time of the standard Adler-Wiser formalism for large
@ -450,8 +453,8 @@ W = v + (\ee^{-1} - 1) v
where matrix notation is used.
This particular decomposition of _W_ , once inserted in the convolution
defining Σ, leads to the split of the self-energy into two different
This particular decomposition of $W$ , once inserted in the convolution
defining $\Sigma$, leads to the split of the self-energy into two different
contributions (exchange and correlation):
\begin{equation}
@ -467,7 +470,7 @@ quasiparticle amplitudes
\sum_\nu^\text{occ} \Psi_{n\kk}(\rr_1){\Psi^\*_{n\kk}}(\rr_2)\,v(\rr_1,\rr_2)
\end{equation}
while the dynamic part Σc(ω) accounts for correlation effects beyond Σx.
while the dynamic part $\Sigma_c(\ww)$ accounts for correlation effects beyond $\Sigma_x$.
It is important to stress that, for computational efficiency, the code does
not compute the full self-energy operator by default. Only its matrix elements
@ -489,13 +492,14 @@ converged independently of others if necessary, and given a much larger value
in comparison to [[ecut]], [[ecutwfn]] and [[ecuteps]].
Another point worth noting is the presence in the expression of the Coulomb
singularity for | **q** + **G** | -> 0\. From a mathematical point of view,
the integral is well-defined since the singularity is integrable in three-
dimensional space once the thermodynamical limit, _N_ **q** -> ∞, is reached.
On the other hand, only a finite number of **q** -points can be used for
singularity for $| \qq + \GG | \rightarrow 0$ .
From a mathematical point of view, the integral is well-defined since the singularity is integrable in three-
dimensional space once the thermodynamical limit, $N_\qq \rightarrow \infty$, is reached.
On the other hand, only a finite number of $\qq$-points can be used for
practical applications, and a careful numerical treatment is needed to avoid
an exceedingly slow convergence with respect to the BZ sampling. To accelerate
the convergence in the number of **q** -points, the code implements several
the convergence in the number of $\qq$-points, the code implements several
techniques proposed in the literature. We refer to the documentation of
[[icutcoul]] for a more extensive discussion.
@ -512,7 +516,7 @@ approximate way, alternatively, it is possible to use the more sophisticated
frequency integration of the contour deformation method [[cite:Lebegue2003]] for
accurate QP calculations (see the related variables [[ppmodel]] and [[gwcalctyp]]).
The double sum over **G** -vectors is performed for all the plane waves
The double sum over $\GG$-vectors is performed for all the plane waves
contained within a sphere of energy [[ecuteps]] (it cannot be larger than the
value used to generate the SCR file). For each state, the correlated matrix
elements are evaluated on a linear frequency mesh centered around the initial
@ -520,7 +524,7 @@ KS energy and the derivative needed for the renormalization factor is obtained
numerically (see [[nomegasrd]] and [[omegasrdmax]]).
Note that here, in contrast to the exchange term, the sum over the band index
_n_ should extend up to infinity although in practice only a finite number of
$n$ should extend up to infinity although in practice only a finite number of
states can be used (specified by [[nband]]).
## Plasmon-pole models
@ -528,14 +532,14 @@ states can be used (specified by [[nband]]).
One of the major computational efforts in self-energy calculations is
represented by the calculation of the frequency dependence of the screened
interaction, which is needed for the evaluation of the convolution. Due to the
ragged behavior of _G_ (ω) and _W_ (ω) along the real axis, numerous real
ragged behavior of $G(\ww)$ and $W (\ww)$ along the real axis, numerous real
frequencies are in principle required to converge the results (note that the
maximum frequencies needed are now reported if [[prtvol]] > 9). On the other
hand, since the fine details of _W_ (ω) are integrated over, it is reasonable
hand, since the fine details of $W(\ww)$ are integrated over, it is reasonable
to expect that approximate models, able to capture the main physical features
of the screened interaction, should give sufficiently accurate results with a
considerably reduction of computational effort. This is the basic idea behind
the so-called plasmon-pole models in which the frequency dependence of _W_ (ω)
the so-called plasmon-pole models in which the frequency dependence of $W(\ww)$
is modelled in terms of analytic expressions depending on coefficients that
are derived from first principles, i.e. without any adjustable external parameters.
@ -558,7 +562,7 @@ function is modeled according to
The two models differ in the approach used to compute the
parameters. [[ppmodel]]=1 derives the parameters such that the inverse
dielectric matrix is correctly reproduced at two different explicitly
calculated frequencies: the static limit (ω = 0) and an additional imaginary
calculated frequencies: the static limit ($\ww = 0$) and an additional imaginary
point located at [[ppmfrq]]. Unless the user overrides this, the default value
is calculated from the average electronic density of the system. The plasmon-
pole parameters of [[ppmodel]]=2 are calculated so as to reproduce the static
@ -589,8 +593,8 @@ contributions coming from the poles of the integrand lying inside the contour:
\end{equation}
In the above equation, the first sum is restricted to the poles lying inside the path
_C_. _W_ c( _z_ ) represents the frequency dependent part of the screened
In the above equation, the first sum is restricted to the poles lying inside the path $\mcC$.
$W^c(z)$ represents the frequency dependent part of the screened
interaction, whose expression in reciprocal space is given by:
\begin{equation}
@ -619,7 +623,7 @@ v(\rr_1,\rr_2)\,\delta(t_1-t_2) $$
$$ 1^+ = (\rr_1,t_1 + \eta)_{\eta \rightarrow 0^+} $$
where _v_ ( **r** 1, **r** 2) represents the bare Coulomb interaction, and η is a positive infinitesimal.
where $v(\rr_1, \rr_2)$ represents the bare Coulomb interaction, and $\eta$ is a positive infinitesimal.
The Fourier transforms for periodic lattice quantities are defined as
@ -640,6 +644,6 @@ f_{\GG_1\GG_2}(\qq) =
(\qq+\GG_2) \cdot \rr_2}\dd\rr_1\dd\rr_2
\end{equation}
The volume of the unit cell is denoted with Ω, while _V_ is the total volume of the crystal simulated
The volume of the unit cell is denoted with $\Omega$, while $V$ is the total volume of the crystal simulated
employing Born-von Karman periodic boundary condition. Unless otherwise
specified, Hartree atomic units will be used throughout.

View File

@ -19,7 +19,7 @@ $$ \rho(\rr)=\sum_{\alpha} \rho^{\alpha\alpha}(\rr) $$
and the magnetization density $\vec m(\rr)$ (in units of $\hbar /2$) whose components are:
$$ m_i(\rr) = \sum_{\alpha\beta} \rho^{\alpha\beta}(\rr) \sigma_i^{\alpha\beta} $$,
$$ m_i(\rr) = \sum_{\alpha\beta} \rho^{\alpha\beta}(\rr) \sigma_i^{\alpha\beta}, $$
where the $\sigma_i$ are the Pauli matrices.
@ -99,7 +99,8 @@ Note that only the forurier transform are performed in *mkrho.f*, while the fina
$\rho(\rr)$, $\vec m(\rr)$ is performed in *symrhg.f*.
The computation of $V_{xc}$ is performed in *rhohxc.f*. The only transformation to this routine, is
to compute $|\vec m(\rr)|$ as input of the usual (i.e spin polarized) {\tt rhohxc.f} and yield back
the four component $V_{xc}$, from the expression of ${\delta E_{xc} \over \delta |m (\rr)| }$.
to compute $|\vec m(\rr)|$ and yield back the four component $V_{xc}$, from the expression
of ${\delta E_{xc} \over \delta |m (\rr)| }$.
For more information, see: Hobbs et al., PRB, 62, 11556; Perdew et al. PRB, 46, 6671 (for the xc functional)
For more information about noncollinear magnetism see [[cite:Hobbs2000]]
and [[cite:Perdew1992]] for the xc functional.

View File

@ -0,0 +1,137 @@
---
authors: MG
---
$
\newcommand{\aa}{\alpha}
\newcommand{\PS}{{\text{PS}}}
\newcommand{\Ylm}{{Y_m^l}}
\newcommand{\Vloc}{V_{\text{loc}}}
\newcommand{\vv}{\hat v}
\newcommand{\vloc}{\vv_{\text{loc}}}
\newcommand{\Vnl}{V_{\text{nl}}}
\newcommand{\lm}{{lm}}
\newcommand{\KK}{{\bf K}}
\newcommand{\KKp}{{\bf{K'}}}
\newcommand{\KKhat}{\widehat \KK}
\newcommand{\KKphat}{\widehat{\KK}'}
\newcommand{\jl}{j_l}
\newcommand{\rrhat}{{\widehat\rr}}
$
## Norm-conserving Pseudopotentials in the Kleynman Bylander form
This page reports basic definitions and equations needed
to evalute the matrix elements of the commutator between the position operator, $\rr$,
and the nonlocal part of the Hamiltonian in the particular case in which norm-conserving pseudopotentials are used.
Providing a consistent introduction to the theory of pseudopotentials is beyond the purposes of this section,
a more complete discussion of the topic can be found in specialized
articles available in the literature ~\cite{Hamann1979,Bachelet1982,Kleinman1982,Hamann1989,Troullier1991,Gonze1991,Fuchs1999}.
For our purposes, we can summarize by saying that a pseudopotential in constructed in order to
replace the atomic all-electron potential such that core states are eliminated
and valence electrons are described by nodeless pseudo wavefunctions whose
representation in the Fourier space decays rapidly.
Modern norm-conserving pseudopotentials are constructed following the procedure
described in~\cite{Hamann1979,Hamann1989,Troullier1991} such that the scattering properties of the
all-electron atom are reproduced around the reference energy configuration, up to first order in energy.
In modern ab initio codes, the fully separable form proposed by Kleynman and Bylander
is usually employed as it drastically reduces the number of operations
required for the application of the nonlocal part of the Hamiltonian
as well as the memory required to store the operator~\cite{ Kleinman1982}.
In the KB form, the interaction between a valence electron and an ion
of type $\aa$ is described by means of the operator:
\begin{equation}\label{eq:KBsingleatom}
\vv_\PS^\aa =
\vloc^\aa(r) + \sum_\lm |\chi_\lm^\aa\ra\,E_l^\aa\,\la\chi_\lm^\aa|,
% \vloc_\aa(\rr) + \sum_\lm |\Ylm\,u_{l\aa}^\KB\ra E_{l\aa}^\KB \lau_{l\aa}^\KB\,\Ylm|
\end{equation}
where $\vloc^\aa(r)$ is a purely local potential with a Coulomb tail $\gamma/r$.
!!! note
For simplicity, the discussion is limited to pseudopotentials with a single projector per angular channel.
The generalization to multi-projector pseudopotentials requirens introducing an additional $n$ index
in the KB energies i.e. $E_{nl}^\aa$ and in the projectors $\chi_{n\lm}^\aa(\rr)$.
The so-called Kleinman-Bylander energies, $E_l^\aa$,
measure the strength of the nonlocal component with respect to the local part.
The projectors $\chi_\lm^\aa(\rr)$ are
short-ranged functions, expressed in terms of a complex spherical harmonic $\Ylm(\theta,\phi)$
multiplied by an $l$- and atom-dependent radial function
\begin{equation}
\label{eq:KB_projectors}
\chi_\lm^\aa(\rr) = \dfrac{1}{r}\,\chi_l^\aa(r)\,\Ylm(\theta,\phi).
\end{equation}
<!--
where
\begin{equation}\label{eq:KBfunctionU}
u_{l\aa}^\KB (\rr) =
\frac{\Delta v_{l\aa}(\rr)\,u_{l\aa}^\PS (\rr)}
{\norm{u_{l\aa}(\rr)\,\Delta_{l\aa}}^{1/2}}
\end{equation}
are localized functions defined in terms of the short-ranged ...
and the pseudo eigenfuncions of the reference atom.
\begin{equation}\label{eq:KBenergy}
E_{l\aa}^\KB =
\frac{\la u_{l\aa}^\PS \Delta v_{l\aa} | \Delta v_{l\aa} u_{l\aa}^\PS \ra }
{\la u_{l\aa}^\PS | \Delta v_{l\aa} | u_{l\aa}^\PS \ra}
\end{equation}
-->
The total nonlocal part of the Hamiltonian is obtained
by summing the different atom-centered contributions.
The final expression in real space reads:
\begin{equation}
\Vnl(\rr_1,\rr_2)=
\sum_{\substack{\RR \\ \aa\tt_\aa}}
\sum_\lm \la\rr_1-\RR-\tt_\aa|\chi_\lm^\aa\ra E_l^\aa\la\chi_\lm^\aa|\rr_2-\RR-\tt_\aa\ra,
\end{equation}
where $\RR$ runs over the sites of the Bravais lattice, and $\tt_\aa$ over the
positions of the atoms with the same type $\aa$ located inside the unit cell.
Due to the nonlocality of the operator, its Fourier transform depends
on two separate indices, $\kk+\GG_1$ and $\kk+\GG_2$, instead of the simple difference $\GG_1-\GG_2$.
The shorthand notation $\KK = \kk + \GG$ is used in the following to simplify the derivation.
Using the Rayleigh expansion of planewaves in terms of spherical
harmonics $\Ylm(\theta,\phi)$ and spherical Bessel functions $\jl(Kr)$:
\begin{equation}\label{eq:PWinYlmPl}
e^{i\KK\cdot\rr} =
4\pi\,\sum_\lm i^l \jl(Kr) \Ylm^\*(\KKhat)\,\Ylm(\rrhat),
\end{equation}
the Fourier representation of the projector defined by \ref{eq:KB_projectors} can be expressed as:
\begin{equation}
\la\KK|\chi_\lm^\aa\ra =
\frac{4\pi}{\Omega^{1/2}}\, (-1)^l\,\Ylm(\KKhat)
\int_0^\infty\, r\jl(Kr)\,\chi_l^\aa(r)r\,\dd =
\frac{4\pi}{\Omega^{1/2}}\, (-1)^l\,\Ylm(\KKhat) F_l^\aa(K).
\end{equation}
<!--
%where the form factors $F_l^\aa(K)$ related to the atom of type $\aa$ is defined by
%\begin{equation}\label{wq:defformfactors}
% F_l^\aa(K) \df
% \frac{\int_0^\infty r\,j_l (Kr)\,u_{l\aa} \Delta v_{l\aa}\,\dd r}
% {\norm{u_{l\aa} \Delta v_{l\aa}}^{1/2} }
%\end{equation}
-->
Finally, the expression for the total nonlocal operator in reciprocal space is:
\begin{equation}
\label{eq:VnlKBmatrixelements}
\Vnl(\KK,\KKp)= \frac{(4\pi)^2}{\Omega} \sum_{\aa \tt_\aa} \sum_\lm
\,e^{-i(\KK-\KKp)\cdot\tt_\aa}\,
\Ylm(\KKhat)\Ylm^\*(\KKphat)\, E_l^\aa F_l^\aa(K) F_l^\aa(K').
\end{equation}

View File

@ -0,0 +1,222 @@
---
authors: MG
---
$$
\newcommand{\lm}{{lm}}
\newcommand{\Ylm}{{Y_m^l}}
\newcommand{\rY}{{\mathcal Y}}
\newcommand{\rYlm}{{\mathcal Y}_m^l}
\newcommand{\rYLM}{{\mathcal Y}_M^L}
\newcommand{\lp}{{l'}}
\newcommand{\rrhat}{{\widehat\rr}}
\newcommand{\rrphat}{{{\widehat\rr}'}}
\newcommand{\mcP}{{\mathcal{P}}}
\newcommand{\df}{\equiv} %this required mathtools :=
\newcommand{\li}{{l_i}}
\newcommand{\mi}{{m_i}}
\newcommand{\nj}{{n_j}}
\newcommand{\lj}{{l_j}}
\newcommand{\mj}{{m_j}}
\renewcommand{\mp}{{m'}}
\newcommand{\mcG}{{\mathcal G}}
\newcommand{\omcR}{{\hat\mcR}}
\newcommand{\mcR}{{\mathcal{R}}}
\newcommand{\ddO}{{\dd\Omega}}
\newcommand{\LM}{{LM}}
\newcommand{\mcRm}{\mcR^{-1}}
\newcommand{\kpGhat}{\widehat{\kpG}}
\newcommand{\rrp}{{\bf r'}}
\newcommand{\Gaunt}{{\mathcal G}}
$$
# Complex and Real Spherical Harmonics
## Complex Spherical Harmonics
Complex spherical harmonics, $\Ylm$, are defined as the eigenfunctions of the orbital
angular momentum operators, $\hat L^2$ and $\hat L_z$.
They can be labeled by two quantum numbers, $l$ and $m$, related to the eigenvalues of
$\hat L^2$ and $\hat L_z$:
\begin{equation}
\hat L^2 \Ylm(\theta,\phi) = l(l+1) \Ylm(\theta,\phi),
\end{equation}
and
\begin{equation}
\hat L_z \Ylm(\theta,\phi) = m \Ylm(\theta,\phi).
\end{equation}
Their explicit expression is:
\begin{equation}
\Ylm(\theta,\phi) = (-1)^m N_\lm P_m^l(\cos\theta) e^{im\phi},
\end{equation}
where $l$ takes non-negative integer values and the possible values for $m$ are
integers from $-l$ to $l$.
The two symbols $\theta$ and $\phi$ are used to denote angular spherical coordinates.
$N_\lm$ is the normalization factor:
\begin{equation}
N_\lm = \sqrt{\frac{2l+1}{4\pi} \dfrac{(l-m)!}{(l+m)!}}
\end{equation}
with $P_m^l$ being the associated Legendre polynomials.
Complex spherical harmonics are orthogonal functions on the sphere, i.e.:
\begin{equation}
\la\Ylm|Y_\mp^\lp\ra_\Omega = \delta_{l\lp}\delta_{m\mp},
\end{equation}
where $\la.|.\ra_\Omega$ denotes the integration over the unit sphere.
Furthermore, they form a complete basis set since any square-integrable
function $f(\theta,\phi)$ can be expanded in series according to:
\begin{equation}
f(\theta,\phi) = \sum_{l=0}^\infty \sum_{m=-l}^l a_{\lm}\Ylm(\theta,\phi).
\end{equation}
Complex spherical harmonics satisfy two useful properties commonly employed in practical applications.
The parity properties:
\begin{equation}\label{eq:parity_property_Ylm}
{Y_m^l}^\dagger(\rrhat) = (-1)^m Y_{-m}^l (\rrhat),
\end{equation}
\begin{equation}
\Ylm(-\rrhat) = (-1)^l \Ylm(\rrhat),
\end{equation}
and the addition theorem:
\begin{equation}\label{eq:addition_theorem_CSH}
\sum_m \Ylm(\rrhat)\, \Ylm^\*(\rrphat) =
\dfrac{2l+1}{4\pi}\,\mcP_l (\rrhat\cdot\rrphat),
\end{equation}
where $\mcP_l(q)$ are Legendre polynomials.
Equation \ref{eq:addition_theorem_CSH} is generally employed in the
application of the nonlocal part of the Hamiltonian in the case of
norm-conserving pseudopotentials when [[useylm]] = 0.
## Real spherical harmonics <a id="RSH"><a>
Real spherical harmonics (RSH) are obtained by combining complex conjugate
functions associated to opposite values of $m$.
RSH are the most adequate basis functions for calculations in which atomic
symmetry is important since they can be directly related to the irreducible
representations of the subgroups of $D_3$ [[cite:Blanco1997]].
Moreover, being real, they have half the memory requirement of complex spherical harmonics.
This is clearly an advantage if high angular momenta are needed or several RHS values have to be stored in memory.
A possible definition for RHS is [[cite:Blanco1997]]:
\begin{equation}\label{eq:Definition_real_harmonics}
\rYlm \df
\begin{cases}
\quad \dfrac{(-1)^m}{\sqrt{2}}\,(\Ylm + \Ylm^\*) \quad m > 0
\\\\
\quad Y_0^l \quad m = 0
\\\\
\quad \dfrac{(-1)^m}{i\sqrt{2}}\,(Y_{\lvert{m}\rvert}^l - {Y_{\lvert{m}\rvert}^l}^\*) \quad m < 0.
\end{cases}
\end{equation}
Equation \ref{eq:Definition_real_harmonics} can be rewritten in matrix notation as
\begin{equation}\label{eq:Unitaty_transformation_Y_Re2Cmplx}
\underline\rY^l = U^l \underline{Y}^l,
\end{equation}
where $\underline\rY^l$ and $\underline{Y}^l$ are column vectors containing real and
spherical harmonics ordered by increasing $m$ values.
$U^l$ is a unitary matrix whose explicit expression can be found in [[cite:Blanco1997]].
The basic properties of RSH can be easily derived from
the properties of complex spherical harmonics by means of \ref{eq:Definition_real_harmonics}.
Reality
:
\begin{equation}
{\rYlm}^\*(\rrhat) = \rYlm(\rrhat),
\end{equation}
Parity property
:
\begin{equation}
\rYlm(-\rrhat) = (-1)^l \rYlm (\rrhat),
\end{equation}
Orthonormality
:
\begin{equation}
\la\rYlm|\rY_{\mp}^\lp\ra_\Omega = \delta_{l\lp} \delta_{m\mp}.
\end{equation}
Certain on-site matrix elements include the integration of three real spherical harmonics, which
can be evaluated exactly:
\begin{equation}\label{eq:Gaunt_def}
\Gaunt_{\li\mi\lj\mj}^{\LM} \df \int \rY_\mi^\li\,\rY_M^L\,\rY_\mj^\lj\ddO.
\end{equation}
The above integral is known in the literature under the name of Gaunt coefficient.
Due to selection rules, the integral in equation \ref{eq:Gaunt_def} is nonzero only if $\mi+\mj = M$
and $\lvert{\li+\lj}\rvert \ge (L,L+2,L+4, \dots) \ge \lvert{\li-\lj}\rvert$.
!!! note
The expression for $\mcG$ follows the internal conventions
used inside the ABINIT code and slightly differs from the one reported in standard text-books.
## Symmetry transformation of RSH
The use of symmetry properties is of fundamental importance to simplify the
determination of electronic structure as it greatly
reduces the number of non-equivalent integrals that have to be computed,
as well as the size of any tables stored in memory.
Frequently one has to evaluate the value of a RSH at a rotated point after
the application of a symmetry operation belonging to the space group of the crystal.
This usually occurs when $\kk$-dependent matrix elements have to be symmetrized in the full
Brillouin zone starting from the knowledge of their symmetrical image in the irreducible wedge.
The effect of a proper or improper rotation on a RSH can be deduced from:
\begin{equation}
\label{eq:RSH_rotation}
\omcR\,\rYlm(\rrhat) = \rYlm(\mcRm \rrhat) =
\sum_\alpha D^l_{\alpha m}(\mcR)\,\rYlm(\rrhat).
\end{equation}
That is, spherical harmonics of given $l$ are transformed into a linear combination
of RHS of same $l$ where the coefficients are given by:
\begin{equation}
D^l_{\alpha m}(\mcR) \df \la\rY^l_\alpha|\hat\mcR|\rYlm\ra.
\end{equation}
## Useful expansions in terms of RSH
By means of the unitary transformation given in \ref{eq:Definition_real_harmonics},
it is possible to rewrite the Rayleigh expansion of a plane wave in terms of RSH as
\begin{equation}
\label{eq:Rayleigh_expansion_real_Ylm}
e^{i(\kpG)\cdot\rr} =
4\pi \sum_\lm i^l j_l(\lvert{\kpG}\rvert r)\,\rYlm(\kpGhat)\,\rYlm(\rrhat)
\end{equation}
with $j_l(x)$ the spherical Bessel function of order $l$.
In a similar way, the electrostatic potentials can be expanded in a real spherical harmonics basis set.
The final expression is very similar to the one obtained in the case of complex spherical harmonics:
\begin{equation}\label{eq:Expansion_Coulombian_RSH}
\dfrac{1}{\lvert{\rr-\rrp}\rvert} = \sum_{l=0}^\infty
\dfrac{4\pi}{2l+1}
\sum_{m=-l}^{m=l}
\dfrac{r_<^l}{r_>^{l+1}} \rYlm(\rrhat)\rYlm(\rrphat),
\end{equation}
with $r_< \df \min(r,r')$ and $r_> \df \max(r,r')$.

File diff suppressed because it is too large Load Diff

View File

@ -199,8 +199,10 @@ pages:
- tutorial/wannier90.md
- Theory:
- Geometry: theory/geometry.md
- Spherical Harmonics: theory/spherical_harmonics.md
- Wavefunctions: theory/wavefunctions.md
- Noncollinear: theory/noncollinear.md
- Pseudopotentials: theory/pseudopotentials.md
- Many-Body Perturbation Theory: theory/mbt.md
- Bethe-Salpeter: theory/bse.md
- Documents: theory/documents.md