The mixed product basis can still be quite large. In the calculation of the screened interaction, each matrix element, when represented in the basis of Coulomb eigenfunctions, is multiplied by :math:`{\sqrt{v_\mu v_\nu}}` with the Coulomb eigenvalues :math:`{\{v_\mu\}}`. This gives an opportunity for reducing the basis-set size further by introducing a Coulomb cutoff :math:`{v_\mathrm{min}}`. The reduced basis set is then used for the polarization function, the dielectric function, and the screened interaction. The parameter :math:`{v_\mathrm{min}}` can be specified after the keyword ``OPTIMIZE MB`` in three ways: first, as a "pseudo" reciprocal cutoff radius :math:`{\sqrt{4\pi/v_\mathrm{min}}}` (which derives from the plane-wave Coulomb eigenvalues :math:`{v_\mathbf{G}=4\pi/G^2}), second, directly as the parameter :math:`{v_\mathrm{min}}` by using a negative real number, and, finally, as the number of basis functions that should be retained when given as an integer. The so-defined basis functions are mathematically close to plane waves. For testing purposes, one can also enforce the usage of plane waves (or rather projections onto plane waves) with the keyword ``OPTIMIZE PW``, in which case the Coulomb matrix is known analytically. No optimization of the basis is applied, if ``OPTIMIZE`` is omitted.

The mixed product basis can still be quite large. In the calculation of the screened interaction, each matrix element, when represented in the basis of Coulomb eigenfunctions, is multiplied by :math:`{\sqrt{v_\mu v_\nu}}` with the Coulomb eigenvalues :math:`{\{v_\mu\}}`. This gives an opportunity for reducing the basis-set size further by introducing a Coulomb cutoff :math:`{v_\mathrm{min}}`. The reduced basis set is then used for the polarization function, the dielectric function, and the screened interaction. The parameter :math:`{v_\mathrm{min}}` can be specified after the keyword ``OPTIMIZE MB`` in three ways: first, as a "pseudo" reciprocal cutoff radius :math:`{\sqrt{4\pi/v_\mathrm{min}}}` (which derives from the plane-wave Coulomb eigenvalues :math:`{v_\mathbf{G}=4\pi/G^2}`), second, directly as the parameter :math:`{v_\mathrm{min}}` by using a negative real number, and, finally, as the number of basis functions that should be retained when given as an integer. The so-defined basis functions are mathematically close to plane waves. For testing purposes, one can also enforce the usage of plane waves (or rather projections onto plane waves) with the keyword ``OPTIMIZE PW``, in which case the Coulomb matrix is known analytically. No optimization of the basis is applied, if ``OPTIMIZE`` is omitted.

(*) The LAPW method relies on a partitioning of space into MT spheres and the interstitial region. The basis functions are defined differently in the two regions, interstitial plane waves in the latter and numerical functions in the spheres with radial parts :math:`{u(r)}, {\dot{u}(r)=\partial u(r)/\partial\epsilon}`, :math:`{u^\mathrm{LO}(r)}` and spherical harmonics :math:`{Y_{lm}(\hat{\mathbf{r}})}` The plane waves and the angular part of the MT functions can be converged straightforwardly with the reciprocal cutoff radius :math:`{g_\mathrm{max}}` and the maximal l quantum number :math:`{l_\mathrm{max}}`, respectively, whereas the radial part of the MT functions is not converged as easily. The standard LAPW basis is restricted to the functions :math:`{u}` and :math:`{\dot{u}}`. Local orbitals :math:`{u^\mathrm{LO}}` can be used to extend the basis set, to enable the description of semicore and high-lying conduction states. The accuracy of the radial MT basis can be analyzed with the keyword MTACCUR e1 e2 which gives the MT representation error [Phys. Rev. B 83, 081101] in the energy range between e1 and e2. (If unspecified, e1 and e2 are chosen automatically.) The results are written to the output files spex.mt.t where t is the atom type index, or spex.mt.s.t with the spin index s(=1 or 2) for spin-polarized calculations. The files contain sets of data for all l quantum numbers, which can be plotted separately with gnuplot (e.g., plot "spex.mt.1" i 3 for :math:`{l=3}`

(*) The LAPW method relies on a partitioning of space into MT spheres and the interstitial region. The basis functions are defined differently in the two regions, interstitial plane waves in the latter and numerical functions in the spheres with radial parts :math:`{u(r)}, {\dot{u}(r)=\partial u(r)/\partial\epsilon}`, :math:`{u^\mathrm{LO}(r)}` and spherical harmonics :math:`{Y_{lm}(\hat{\mathbf{r}})}` The plane waves and the angular part of the MT functions can be converged straightforwardly with the reciprocal cutoff radius :math:`{g_\mathrm{max}}` and the maximal l quantum number :math:`{l_\mathrm{max}}`, respectively, whereas the radial part of the MT functions is not converged as easily. The standard LAPW basis is restricted to the functions :math:`{u}` and :math:`{\dot{u}}`. Local orbitals :math:`{u^\mathrm{LO}}` can be used to extend the basis set, to enable the description of semicore and high-lying conduction states. The accuracy of the radial MT basis can be analyzed with the keyword MTACCUR e1 e2 which gives the MT representation error [Phys. Rev. B 83, 081101] in the energy range between e1 and e2. (If unspecified, e1 and e2 are chosen automatically.) The results are written to the output files spex.mt.t where t is the atom type index, or spex.mt.s.t with the spin index s(=1 or 2) for spin-polarized calculations. The files contain sets of data for all l quantum numbers, which can be plotted separately with gnuplot (e.g., plot "spex.mt.1" i 3 for :math:`{l=3}`)

@@ -63,9 +63,9 @@ The two GW values for the quasiparticle energy follow two common methods to appr

where :math:`{v^\mathrm{ext}}`, :math:`{v^\mathrm{H}}`, :math:`{\Sigma^\mathrm{xc}}`, :math:`{\psi_{\mathbf{k}n}}`, :math:`{E_{\mathbf{k}n})}` are the external and Hartree potential, the ``GW`` self-energy, and the quasiparticle wavefunction and energy, respectively. Taking the difference :math:`{\Sigma^\mathrm{xc}-v^\mathrm{xc}}` as a small perturbation, we can write the quasiparticle energy as a correction on the mean-field eigenvalue

with the single-particle wavefunction :math:`{\phi_{\mathbf{k}n}}` and the frequency-independent potential :math:`{v^{\mathrm{xc}}}`, which in the case of a KS solution would correspond to the local exchange-correlation potential; the nonlocal Hartree-Fock exchange potential and the ''hermitianized'' self-energy of QSGW (see below) are other examples. :math:`{Z_{\mathbf{k}n}=[1-\partial\Sigma^{\mathrm{xc}}/\partial\omega(\epsilon_{\mathbf{k}n})]^{-1}}` is called the renormalization factor. The two expressions on the right-hand side correspond to the "linearized" and "direct" (iterative) solutions given in the output. The direct solution takes into account the non-linearity of the quasiparticle equation and is thus considered the more accurate result. However, there is usually only little difference between the two values.

Up to this point, the job syntax for Hartree Fock (``JOB HF``), PBE0 (``JOB PBE0``), screened exchange (``JOB SX``), COHSEX (``JOB COSX``), and ''GT'' (``JOB GT`` and ``JOB GWT``) calculations is identical to the one of ``GW`` calculations, e.g., ``JOB HF FULL X:(1-10)``. Except the latter (``GT``), all of these methods are mean-field approaches, so one only gets one single-particle energy (instead of a ''linearized'' and a ''direct'' solution) for each band.

...

...

@@ -77,11 +77,13 @@ SPECTRAL (SENERGY)

(*) It should be pointed out that the quasiparticle energies given in the output rely on the quasiparticle approximation. The more fundamental equation, as it were, is the Dyson equation

which links the interacting Green function {$G$} to the non-interacting KS one {$G_0$} and which, in principle, requires the self-energy to be known on the complete :math:`{\omega}` axis. The spectral function measured in photoelectron spectroscopy is directly related to the Green function by

which links the interacting Green function {$G$} to the non-interacting KS one :math:`{G_0}` and which, in principle, requires the self-energy to be known on the complete :math:`{\omega}` axis. The spectral function measured in photoelectron spectroscopy is directly related to the Green function by

where the trace (tr) is over the eigenstates and :math:`{\epsilon_\mathrm{F}}` is the Fermi energy. The spectral function can be evaluated using the line

...

...

@@ -103,13 +105,13 @@ Alhough ``GW`` calculations can be readily performed with the default settings,

can be understood as a scattering potential that contains the exact exchange potential and correlation effects through the inclusion of ''W'', the dynamically screened interaction, which incorporates the screening of the many-electron system into an effective dynamical potential, obtained from the dielectric function :math:`{\varepsilon}` through

This integral equation turns into a matrix equation

...

...

@@ -117,13 +119,13 @@ This integral equation turns into a matrix equation

if the quantities :math:`{W}`, :math:`{\varepsilon}`, and :math:`{v}` ; are represented in the :ref:`mbp`, which thus has to be converged properly. The dielectric function, in turn, describes the change of the internal potential through screening processes and is related to the polarization matrix by :math:`{\varepsilon(\mathbf{k},\omega)=1-P(\mathbf{k},\omega)v(\mathbf{k})}` in matrix notation. The polarization function is one of the main quantities in the ``GW`` scheme. Its evaluation is described in the section :ref:`polar`.

The self-energy can be written as the sum of two terms, the first of which is the exact non-local exchange potential of Hartree-Fock theory, the remainder can be interpreted as a correlation self-energy and has the mathematical form of the [[#Eq:Selfene|self-energy]] with :math:{W(\omega)} replaced by :math:`{W^\mathrm{c}(\omega)=W(\omega)-v}`. The frequency integration is carried out analytically for the exchange part (by summing over the residues). The correlation part is more complex to evaluate because of the frequency dependence of the interaction. (Fast electrons experience a different potential than slow electrons.)

The self-energy can be written as the sum of two terms, the first of which is the exact non-local exchange potential of Hartree-Fock theory, the remainder can be interpreted as a correlation self-energy and has the mathematical form of the [[#Eq:Selfene|self-energy]] with :math:`{W(\omega)}` replaced by :math:`{W^\mathrm{c}(\omega)=W(\omega)-v}`. The frequency integration is carried out analytically for the exchange part (by summing over the residues). The correlation part is more complex to evaluate because of the frequency dependence of the interaction. (Fast electrons experience a different potential than slow electrons.)

There are several ways to represent the self-energy as a function of frequency. The default method is ``analytic continuation``, in which the screened interaction and the self-energy are evaluated on a mesh of purely imaginary frequencies. The self-energy is then analytically continued to the complete complex frequency plane (:math:`{\Sigma^\mathrm{xc}(i\omega)\rightarrow\Sigma^\mathrm{xc}(z)}`, :math:`{\omega\in\cal{R}}`, :math:`{z\in\cal{C}}`). This has several advantages over the usage of the real frequency axis. First, :math:`{W(\omega)}` is a hermitian (or real-symmetric) matrix if :math:`{\omega}` is purely imaginary. Second, ``W`` and :math:`{\Sigma^{\mathrm{xc}}}` show a lot of structure along the real axis, whereas they are much smoother on the imaginary axis, thereby making it easier to sample and interpolate these functions. Third, after the analytic continuation the self-energy is known, in principle, as an analytic function on the complete complex plane. And fourth, the method requires only few parameters and is, therefore, easy to handle. The main disadvantage lies is the badly controlled extrapolation of the Pade approximants, which can sometimes produce "outlier values", with a potential adverse effect on the accuracy and reliability of the method. Therefore, there is a more sophisticated but also more complex method called ``contour integration``, in which the frequency convolution is performed explicitly, yielding the self-energy directly for selected frequencies on the real axis. In this method, we also mostly integrate along the imaginary frequency axis.

MESH (SENERGY)

--------------

All methods employ a mesh of purely imaginary frequencies, which extends from 0 to some maximal :math:`{i\omega_\mathrm{max}}`. The number of mesh points :math:{N} and the maximal frequency must be provided as parameters, for example ``MESH 10 10.0``, which is the default. The mesh is defined by :math:`{i\omega_n=i\omega_\mathrm{max}f_n/f_N}` with :math:`{f_n=\{(N-1)/[0.9(n-1)]-1\}^{-1}}`, :math:`{n=1,2,...}`. It is fine for small :math:`{\omega}`, where the quantities have a lot of structure, and coarse for large :math:{\omega}. Sometimes it is helpful to make the mesh even finer for small :math:`{\omega}`. This is possible by specifying, for example, ``10+3``, which would yield three, two, and one extra equidistant frequencies in the ranges [:math:`{\omega_1,\omega_2}`], [:math:`{\omega_2,\omega_3}`], and [:math:`{\omega_3,\omega_4}`], respectively. If the second argument is defined negative (:math:`{-\omega_\mathrm{max}}`), then :math:`{f_n=\{N/(n-1)-1\}^{-1}}`. The latter definition is rarely used. One can also employ a user-defined mesh provided in a file (one frequency per line, comments ``#...`` are allowed).

All methods employ a mesh of purely imaginary frequencies, which extends from 0 to some maximal :math:`{i\omega_\mathrm{max}}`. The number of mesh points :math:`{N}` and the maximal frequency must be provided as parameters, for example ``MESH 10 10.0``, which is the default. The mesh is defined by :math:`{i\omega_n=i\omega_\mathrm{max}f_n/f_N}` with :math:`{f_n=\{(N-1)/[0.9(n-1)]-1\}^{-1}}`, :math:`{n=1,2,...}`. It is fine for small :math:`{\omega}`, where the quantities have a lot of structure, and coarse for large :math:`{\omega}`. Sometimes it is helpful to make the mesh even finer for small :math:`{\omega}`. This is possible by specifying, for example, ``10+3``, which would yield three, two, and one extra equidistant frequencies in the ranges [:math:`{\omega_1,\omega_2}`], [:math:`{\omega_2,\omega_3}`], and [:math:`{\omega_3,\omega_4}`], respectively. If the second argument is defined negative (:math:`{-\omega_\mathrm{max}}`), then :math:`{f_n=\{N/(n-1)-1\}^{-1}}`. The latter definition is rarely used. One can also employ a user-defined mesh provided in a file (one frequency per line, comments ``#...`` are allowed).

| Examples | ``MESH 12 15.0`` | Use a mesh containing twelve frequencies for [0,i15] htr |

...

...

@@ -344,11 +346,13 @@ One can also go beyond the perturbative solution of the [[#Eq:qpeq|quasiparticle

[[#QSGW]]If the job definition contains ``FULL`` and [[#IBZlabel|``IBZ``]], the full ``GW`` self-energy matrix is evaluated for the whole IBZ, which enables self-consistent calculations in the framework of the ``quasiparticle self-consistent GW`` (QSGW) approach. In this approach, one creates a mean-field system from the ``GW`` self-energy whose single-particle energies are as close as possible to the quasiparticle energies. This mean-field system is subsequently solved to self-consistency in a DFT code. The resulting solution can then serve as a starting point for a new one-shot ``GW`` calculation, which constitutes the second iteration and so on until the quasiparticle energies are converged. The construction of the mean-field system is, to some extent, arbitrary. We use the following definition, which is slightly modified from the original work [PRL 93, 126406]:

for off-diagonal elements. The ``hermitianized`` QSGW operator is then obtained from :math:`{\Sigma^\mathrm{xc,H}=(A+A^\dagger)/2}`. The difference to the original definition is the inclusion of the renormalization factor to better reproduce the ``GW`` quasiparticle energies. The ``hermitianized`` matrix, or rather the difference :math:`{\Sigma^{\mathrm{xc,H}}-v^\mathrm{xc}}`, is written to the file "spex.qsgw", which is later read by the DFT code. In Fleur, the following steps are required:

...

...

@@ -389,7 +393,9 @@ Polarization function

=====================

The polarization function gives the linear change in the electronic density of a non-interacting system with respect to changes in the effective potential. It is, thus, a fundamental quantity in the calculation of screening properties of a many-electron systems. For example, the dielectric function, instrumental in the calculation of spectroscopic quantities (e.g. EELS) and the screened interaction needed in ``GW``, is related to the polarization matrix through :math:`{\varepsilon(\mathbf{k},\omega)=1-P(\mathbf{k},\omega)v(\mathbf{k})}`, here given in matrix notation. The corresponding explicit formula for matrix elements of ``P`` in the mixed product basis is

We have implicitly defined the spectral function ``S`` in the last equation, an explicit expression for which is basically given by the formula in the middle with the :math:`{1/(\omega...)}` replaced by :math:`{\delta(\omega...)}`. (Technically, the :math:`{M_{\mathbf{k}\mu}(\mathbf{r})}` form the eigenbasis of the Coulomb matrix, so they are linear combinations of the mixed product basis functions.)

...

...

@@ -430,7 +436,7 @@ The most important keyword in the calculation of the polarization function is ``

MULTDIFF (SUSCEP)

-------------------

(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of [[#Eq:P|``P``]] approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in ``'k``' becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}-\epsilon_{\mathbf{q}n})^{-1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and non-divergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the [[#Eq:P|frequency integration]], thereby avoiding the numerical difficulties. (As an alternative, the energy differences can also be incorporated into the integration weights, which is arguably even more stable numerically, see option ``INT`` below.) By default, Spex does that only for k=0. The behavior can be changed with the keyword ``MULTDIFF`` in the section ``SUSCEP``.

(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of [[#Eq:P|``P``]] approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in :math:`k` becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}-\epsilon_{\mathbf{q}n})^{-1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and non-divergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the [[#Eq:P|frequency integration]], thereby avoiding the numerical difficulties. (As an alternative, the energy differences can also be incorporated into the integration weights, which is arguably even more stable numerically, see option ``INT`` below.) By default, Spex does that only for k=0. The behavior can be changed with the keyword ``MULTDIFF`` in the section ``SUSCEP``.