diff git a/.gitignore b/.gitignore
index 01c3d69278d24be8005b2a4c2c42b10744717739..c5235b8431051393749ae0023582c81d17a9a879 100644
 a/.gitignore
+++ b/.gitignore
@@ 1,2 +1,4 @@
build
.gitlabci.yml
+sec
+*~
diff git a/source/_static/theme_overrides.css b/source/_static/theme_overrides.css
new file mode 100644
index 0000000000000000000000000000000000000000..727df0963c0850a0b2d7e1912b55a327fb5b9807
 /dev/null
+++ b/source/_static/theme_overrides.css
@@ 0,0 +1,13 @@
+/* override table width restrictions */
+@media screen and (minwidth: 767px) {
+
+ .wytableresponsive table td {
+ /* !important prevents the common CSS stylesheets from overriding
+ * this as on RTD they are loaded after this stylesheet */
+ whitespace: normal !important;
+ }
+
+ .wytableresponsive {
+ overflow: visible !important;
+ }
+}
diff git a/source/basis_sets.rst b/source/basis_sets.rst
index 072ef08c8758bd8c7402bf41bc31f6196951a632..c2dc18535d9eb215fdc4f89f6ac1ab2b4db79b7d 100644
 a/source/basis_sets.rst
+++ b/source/basis_sets.rst
@@ 9,95 +9,121 @@ LAPW basis
The reader is supposed to be familiar with the LAPW basis set. This section is intended as a short recapitulation and to define the notation used in the manual.
The LAPW method relies on a partitioning of space into nonoverlapping MT spheres centered at the atomic nuclei and the interstitial region. The basis functions are defined differently in the two regions, plane waves in the interstitial and numerical functions in the spheres, consisting of a radial part :math:`{u_{lp}^a(r)}` (:math:`{a}` is the atomic index) and the spherical harmonics :math:`{Y_{lm}(\hat{\mathbf{r}})}`. :math:`{p}` is an index to distinguish different types of radial functions: :math:`{u_{l0}^a=u_l^a}`, :math:`{\,u_{l1}^a=\dot{u}_l^a=\partial u_l^a/\partial\epsilon}`, :math:`{\,u_{lp}^a=u_l^{a,\mathrm{LO}},\,p\ge 2}`, where we have used the usual notation of LAPW. Augmented plane waves are formed by matching linear combinations of the :math:`{u}` and :math:`{\dot{u}}` to the interstitial plane waves, forming the standard LAPW basis set. Local orbitals :math:`{u^\mathrm{LO}}` can be used to extend the basis set, to enable the description of semicore and highlying conduction states. 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 more difficult to converge.
+The LAPW method relies on a partitioning of space into nonoverlapping MT spheres centered at the atomic nuclei and the interstitial region.
+The basis functions are defined differently in the two regions, plane waves in the interstitial and numerical functions in the spheres,
+consisting of a radial part :math:`{u_{lp}^a(r)}` (:math:`{a}` is the atomic index) and the spherical harmonics :math:`{Y_{lm}(\hat{\mathbf{r}})}`.
+:math:`{p}` is an index to distinguish different types of radial functions: :math:`{u_{l0}^a=u_l^a}`,
+:math:`{\,u_{l1}^a=\dot{u}_l^a=\partial u_l^a/\partial\epsilon}`, :math:`{\,u_{lp}^a=u_l^{a,\mathrm{LO}},\,p\ge 2}`,
+where we have used the usual notation of LAPW. Augmented plane waves are formed by matching linear combinations of the
+:math:`{u}` and :math:`{\dot{u}}` to the interstitial plane waves, forming the standard LAPW basis set.
+Local orbitals :math:`{u^\mathrm{LO}}` can be used to extend the basis set, to enable the description of semicore and highlying conduction states.
+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 more difficult to improve systematically, but it is possible, see [Comput. Phys. Commun. 184, 2670 (2013)].
MTACCUR

(*) 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 spinpolarized 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}`).

++++
 Examples  ``MTACCUR 1 2``   Calculate MT representation error 
    between 1 and 2 Hartree. 
++++
  ``MTACCUR``  Automatic energy range. 
++++

.. _mbp:

Mixed product basis
+(*) The accuracy of the radial MT basis can be analyzed with the keyword ``MTACCUR`` followed by two numbers :math:`e_1` and :math:`e_2`.
+Spex then calculates the MT representation error
+[Phys. Rev. B 83, 081101] in the energy range between :math:`e_1` and :math:`e_2`.
+(If unspecified, upper and lower bounds 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 spinpolarized 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}`).
+
+.. listtable:: Examples
+ :widths: 24 100
+
+ *  ``MTACCUR 1 2``
+  Calculate MT representation error between 1 and 2 Hartree.
+ *  ``MTACCUR``
+  Automatic energy range.
+
+.. _mpb:
+
+Mixed product basis
===================
The methods implemented in Spex distinguish themselves from conventional DFT with local or semilocal functionals by the fact that individual twoparticle scattering processes are described explicitly. In such processes, the state of an incoming particle and the state that it scatters into form products of wave functions. Many of the physical quantities, such as the Coulomb interaction, the polarization function, and the dielectric function, become matrices when represented in a suitable product basis. In the context of the FLAPW method, the mixed product basis is an optimal choice. It consists of two separate sets of functions that are defined only in one of the spatial regions, while being zero in the other, the first is defined by products of the LAPW MT functions and the second by interstitial plane waves. (The products of plane waves are plane waves again.) The corresponding parameters are defined in the section ``MBASIS`` of the input file.
+The methods implemented in Spex distinguish themselves from conventional DFT with local or semilocal functionals by the fact that individual twoparticle scattering processes are described explicitly. In such processes, the state of an incoming particle and the state that it scatters into form products of wavefunctions. Many of the physical quantities, such as the Coulomb interaction, the polarization function, and the dielectric function, become matrices when represented in a suitable product basis. In the context of the FLAPW method, the mixed product basis is an optimal choice. It consists of two separate sets of functions that are defined only in one of the spatial regions, while being zero in the other, the first is defined by products of the LAPW MT functions and the second by interstitial plane waves. (The products of plane waves are plane waves again.) The corresponding parameters are defined in the section ``MBASIS`` of the input file.
GCUT (MBASIS)

If :math:`{N}` is the number of LAPW basis functions, one would naively expect the number of product functions to be roughly :math:`{N^2}`. In the case of the interstitial plane waves, this is not so, since, with a cutoff :math:`{g_\mathrm{max}}`, the maximum momentum of the product would be :math:`{2g_\mathrm{max}}`, leading to :math:`{8N}` as the number of product plane waves. Fortunately, it turns out that the basis size can be much smaller in practice. Therefore, we introduce a reciprocal cutoff radius :math:`{G_\mathrm{max}}` for the interstitial plane waves and find that, instead of :math:`{G_\mathrm{max}=2g_\mathrm{max}}`, good convergence is achieved already with :math:`{G_\mathrm{max}=0.75g_\mathrm{max}}`, the default value. The parameter :math:`{G_\mathrm{max}}` can be set to a different value with the keyword ``GCUT`` in the section ``MBASIS`` of the input file.
++++
 Example  ``GCUT 2.9``   Use a reciprocal cutoff radius of 2.9 
    :math:`Bohr^{1}` for the product plane waves. 
++++

+.. listtable:: Example
+ :widths: 16 100
+
+ *  ``GCUT 2.9``
+  Use a reciprocal cutoff radius of 2.9 Bohr\ :sup:`1` for the product plane waves.
+
.. _lcut:
LCUT (MBASIS)
+LCUT (MBASIS)

The MT functions are constructed from the products :math:`{u_{lp}^a(r)Y_{lm}(\hat{\mathbf{r}})u_{l'p'}^a(r)Y_{l'm'}(\hat{\mathbf{r}})}`. (Obviously, the atom index :math:`{a}` must be indentical.) The product of spherical harmonics exands into a linear combination of :math:`{Y_{LM}(\hat{\mathbf{r}})}` with :math:`{L=ll',...,l+l'}` and :math:`{M\le L}`. In principle, this leads to a cutoff of :math:`{2l_\mathrm{max}}` for the products, but, as with the plane waves, one can afford to use a much smaller cutoff parameter in practice. We use :math:`{L_\mathrm{max}=l_\mathrm{max}/2}` as the default. The corresponding keyword in the input file is called ``LCUT``.

++++
 Example  ``LCUT 5 4``  Use l cutoffs of 5 and 4 for two atom types in the system. 
++++

+The MT functions are constructed from the products :math:`{u_{lp}^a(r)Y_{lm}(\hat{\mathbf{r}})u_{l'p'}^a(r)Y_{l'm'}(\hat{\mathbf{r}})}`. (Obviously, the atom index :math:`{a}` must be identical.) The product of spherical harmonics expands into a linear combination of :math:`{Y_{LM}(\hat{\mathbf{r}})}` with :math:`{L=ll',...,l+l'}` and :math:`{M\le L}`. In principle, this leads to a cutoff of :math:`{2l_\mathrm{max}}` for the products, but, as with the plane waves, one can afford to use a much smaller cutoff parameter in practice. We use :math:`{L_\mathrm{max}=l_\mathrm{max}/2}` as the default. The corresponding keyword in the input file is called ``LCUT``.
+
+.. listtable:: Example
+ :widths: 16 100
+
+ *  ``LCUT 5 4``
+  Use l cutoffs of 5 and 4 for two atom types in the system.
+
SELECT (MBASIS)

In the LAPW basis, the matching conditions at the MT boundaries require large :math:`{l}` cutoff radii, typically around :math:`{l=10}`, giving a correspondingly large number of radial functions. Not all of these functions must be considered in the product basis set. The keyword ``SELECT`` enables a specific choice of the radial functions. When specified, an entry for each atom type is required. The default is ``l;l+1`` with ``l`` being the ``LCUT``. The best way to describe the syntax of ``SELECT`` is with examples:

++++
 Examples  ``SELECT 2;3``   Use products of :math:`{u_{lp} u_{l'p'}}` with :math:`{l\le 2}` and :math:`{l'\le 3}` 
    for :math:`{p=0}` (socalled :math:`{u}`) and no function 
    of :math:`{p=1}` (socalled :math:`{\dot{u}}`). 
++++
  ``SELECT 2,1;3,2``   Same as before but also include the :math:`{\dot{u}}` functions 
    with :math:`{l\le 2}` and :math:`{l'\le 3}`. 
++++
  ``SELECT 2,,1100;3,,1111``   Same as first line but also include the local orbitals 
    :math:`{p\ge 2}`, which are selected (deselected) by "1" ("0"): 
    here, the first two and all four LOs, respectively. 
    The default behavior is to include semicore LOs but to 
    exclude the ones at higher energies. 
++++
  ``SELECT 2,1,1100;3,2,1111``   Same as second line with the LOs. 
++++

The LOs are selected by series of 1s and 0s. If there are many LOs, a definition like ``0000111`` can be abbreviated to ``4031`` (four times ``0``, three times ``1``). It is important to note that the wavefunction products that must be represented consist mostly of occupiedunoccupied pairs. So, the first set of parameters before the semicolon refers to the occupied state, the set after the semicolon refers to the unoccupied state.
+In the LAPW basis, the matching conditions at the MT boundaries require large l cutoff values, typically around :math:`{l=10}`, giving a correspondingly large number of radial functions. Not all of these functions must be considered in the product basis set. The keyword ``SELECT`` enables a specific choice of the radial functions. When specified, an entry for each atom type is required.
+The best way to describe the syntax of ``SELECT`` is with examples:
+
+.. listtable:: Examples
+ :widths: 50 100
+
+ *  ``SELECT 2;3``
+  Use products of :math:`{u_{lp} u_{l'p'}}` with :math:`{l\le 2}` and :math:`{l'\le 3}` for :math:`{p=0}` (socalled :math:`{u}`) and no function of :math:`{p=1}` (socalled :math:`{\dot{u}}`).
+ *  ``SELECT 2,1;3,2``
+  Same as before but also include the :math:`{\dot{u}}` functions with :math:`{l\le 1}` and :math:`{l'\le 2}`.
+ *  ``SELECT 2,,1100;3,,1111``
+  Same as first line but also include the local orbitals :math:`{p\ge 2}`, which are selected (deselected) by "1" ("0"): here, the first two and all four LOs, respectively. The default behavior is to include semicore LOs but to exclude the ones at higher energies.
+ *  ``SELECT 2,1,1100;3,2,1111``
+  Same as second line with the LOs.
+
+The LOs are selected by series of 1s and 0s. If there are many LOs, a definition like ``0000111``
+can be abbreviated by ``4031`` (four times ``0``, three times ``1``).
+
+One of the most important quantities to be calculated in *GW* calculations (and related methods)
+are the "projections"
+:math:`{\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle}`,
+which describe the coupling of a propagating particle to an interaction line. Often, the two states with
+:math:`n` and :math:`n'` are occupied and unoccupied, respectively.
+In this sense, the first set of ``SELECT`` parameters before the semicolon refers to the occupied states,
+the set after the semicolon refers to the unoccupied states.
+For :math:`l=2` with :math:`2l` or :math:`2l+1` being the l cutoff specified by ``LCUT``, the default
+would be ``2;3``.
TOL (MBASIS)

(*) The set of MT products selected by ``SELECT`` can still be highly linearly dependent. Therefore, in a subsequent optimization step one diagonalizes the MT overlap matrix and retains only those eigenfunctions whose eigenvalues exceed a predefined tolerance value. This tolerance is 0.0001 by default and can be changed with the keyword ``TOL`` in the input file.
++++
 Example  ``TOL 0.00001``   Remove linear dependencies that fall below a 
    tolerance of 0.00001 (see text). 
++++

+.. listtable:: Example
+ :widths: 22 100
+
+ *  ``TOL 0.00001``
+  Remove linear dependencies that fall below a tolerance of 0.00001 (see text).
+
OPTIMIZE (MBASIS)

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 basisset 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 planewave 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 sodefined 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.

++++
 Examples  ``OPTIMIZE MB 4.0``   Optimize the mixed product basis by removing eigenfunctions 
    with eigenvalues below ``4\pi/4.5^2``. 
++++
  ``OPTIMIZE MB 0.05``   Optimize the mixed product basis by removing eigenfunctions 
    with eigenvalues below ``0.05``. 
++++
  ``OPTIMIZE MB 80``  Retain only the eigenfunctions with the 80 largest eigenvalues. 
++++
  ``OPTIMIZE PW 4.5``   Use projected plane waves with the cutoff :math:`{4.5 \mathrm{Bohr}^{1}}` 
    (for testing only, can be quite slow). 
++++

In summary, there are a number of parameters that influence the accuracy of the basis set. Whenever a new physical system is investigated, it is recommendable to converge the basis set for that system. The parameters to consider in this respect are ``GCUT``, ``LCUT``, ``SELECT``, and ``OPTIMIZE``.
+Despite the favorable convergence behavior of the mixed product basis,
+it 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 basisset size further by introducing a Coulomb eigenvalue threshold :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 planewave 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 sodefined 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.
+
+.. listtable:: Examples
+ :widths: 34 100
+
+ *  ``OPTIMIZE MB 4.0``
+  Optimize the mixed product basis by removing eigenfunctions with eigenvalues below :math:`4\pi/4.0^2`.
+ *  ``OPTIMIZE MB 0.05``
+  Optimize the mixed product basis by removing eigenfunctions with eigenvalues below 0.05.
+ *  ``OPTIMIZE MB 80``
+  Retain only the eigenfunctions with the 80 largest eigenvalues.
+ *  ``OPTIMIZE PW 4.5``
+  Use projected plane waves with the cutoff 4.5 Bohr\ :sup:`1` (for testing only, can be quite slow).
+
+In summary, there are a number of parameters that influence the accuracy of the basis set. Whenever a new physical system is investigated, it is recommendable to converge the basis set for that system. The parameters to consider in this respect are ``GCUT``, ``LCUT``, ``SELECT``, and ``OPTIMIZE``.
Finally, we show a section ``MBASIS`` for a system with three atom types as an example
@@ 117,7 +143,7 @@ NOAPW (MBASIS)
ADDBAS (MBASIS)

(*) This keyword augments the mixed product basis by the radial :math:`{u_{lp}}` functions, the :math:`{l}` and :math:`{p}` parameters are according to first part of the ``SELECT`` definition. This improves the basis convergence for COHSEX calculations.
+(*) This keyword augments the mixed product basis by the radial :math:`{u_{lp}}` functions, the :math:`{l}` and :math:`{p}` parameters are according to the first part of the ``SELECT`` definition. This improves the basis convergence for COHSEX calculations.
CHKPROD (MBASIS)

@@ 125,38 +151,199 @@ CHKPROD (MBASIS)
WFADJUST (MBASIS)

(**) Removes wavefunction coefficients belonging to radial functions that are not used to construct the MT product basis. (Currently the coefficients are actually not removed from the memory.)
+(**) Removes wavefunction coefficients belonging to radial functions that are not used to construct the MT product basis. (Currently the coefficients are actually not removed from memory.)
FFT (WFPROD)

When the interaction potential is represented in the mixed product basis, the coupling to the singleparticle states involve projections of the form
:math:`{\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle\,.}`
The calculation of these projections can be quite expensive. Therefore, there are a number of keywords that can be used for acceleration. Most of them are, by now, somewhat obsolete. An important keyword, though, is ``FFT`` in the section ``WFPROD`` of the input file. When used, the interstitial terms are evaluated using Fast Fourier Transformations (FFTs), i.e., by transforming into real space (where the convolutions turn into products), instead of by explicit convolutions in reciprocal space. For small systems the latter is faster, but for large systems it is recommendable to use FFTs because of a better scaling with system size. A run with FFTs can be made to yield results identical to the explicit summation. This requires an FFT reciprocal cutoff radius of :math:`{2G_\mathrm{max}+g_\mathrm{max}}`, which can be achieved by setting ``FFT EXACT``, but such a calculation is quite costly. It is, therefore, advisable to use smaller cutoff radii, thereby sacrificing a bit of accuracy but speeding up the computations a lot. If given without an argument, Spex will use 2/3 of the above *exact* cutoff. One can also specify a cutoff by a realvalued argument explicitly, good compromises between accuracy and speed are values between 6 and 8 Bohr'^1^'.

++++
 Examples  ``FFT 6``  Use FFTs with the cutoff 6 Bohr'^1^'. 
++++
  ``FFT``  Use the default cutoff :math:`{\frac{2}{3}(2G_\mathrm{max}+g_\mathrm{max})}`. 
++++
  ``FFT EXACT``  Use the exact cutoff :math:`{2G_\mathrm{max}+g_\mathrm{max}}`. 
++++

+When the interaction potential is represented in the mixed product basis, the coupling to the singleparticle states involve projections of the form
+:math:`{\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle\,.}`
+The calculation of these projections can be quite expensive. Therefore, there are a number of keywords that can be used for acceleration.
+Most of them are, by now, somewhat obsolete. An important keyword, though, is ``FFT`` in the section ``WFPROD`` of the input file.
+When used, the interstitial terms are evaluated using Fast Fourier Transformations (FFTs), i.e., by transforming into real space
+(where the convolutions turn into products), instead of by explicit convolutions in reciprocal space.
+For small systems the latter is faster, but for large systems it is recommendable to use FFTs because of a better scaling with system size.
+A run with FFTs can be made to yield results identical to the explicit summation. This requires an FFT reciprocal cutoff radius of :math:`{2g_\mathrm{max}+G_\mathrm{max}}`, which can be achieved by setting ``FFT EXACT``, but such a calculation is quite costly. It is, therefore, advisable to use smaller cutoff radii, thereby sacrificing a bit of accuracy but speeding up the computations a lot. If given without an argument, Spex will use 2/3 of the above *exact* cutoff. One can also specify a cutoff by a realvalued argument explicitly, good compromises between accuracy and speed are values between 6 and 8 Bohr\ :sup:`1`.
+
+.. listtable:: Examples
+ :widths: 18 100
+
+ *  ``FFT 6``
+  Use FFTs with the cutoff 6 Bohr\ :sup:`1`.
+ *  ``FFT``
+  Use the default cutoff :math:`{\frac{2}{3}(2g_\mathrm{max}+G_\mathrm{max})}`.
+ *  ``FFT EXACT``
+  Use the exact cutoff :math:`{2g_\mathrm{max}+G_\mathrm{max}}`.
+
APPROXPW (WFPROD)

(**) Since version 02.03 the interstitial part of the projections is calculated exactly, leading to favorable convergence with respect to ``GCUT``. However, the exact evaluation takes considerably more time. If ``APPROXPW`` is set, the old and approximate calculation of the products is used instead. (Only included for testing.)
+(**) Since version 02.03 the interstitial part of the projections is calculated exactly, leading to favorable convergence with respect to ``GCUT``. However, the exact evaluation takes considerably more time. If ``APPROXPW`` is set, the old and approximate way of calculating the products is used instead. (Only included for testing.)
LCUT MTTHR (MBASIS)

(**) The MT part of the projections can be accelerated as well with the keywords ``LCUT`` and ``MTTHR`` in the section ``WFPROD``. (Note the other unrelated :ref:`lcut` keyword in ``MBASIS``.) The first can be used to restrict the l cutoff of the wave functions (only for the projections), one l cutoff for each atom type. The second is a threshold value below which wavefunction coefficients are neglected. By default, the two options are unused.
+LCUT and MTTHR (WFPROD)
+
+(**) The MT part of the projections can be accelerated as well using the the keywords ``LCUT`` and ``MTTHR``
+in the section ``WFPROD``. (The keyword ``LCUT`` should not be confused with the keyword of the same name
+in the section ``MBASIS``, :numref:`lcut`.) The first can be used to restrict the l cutoff of the wavefunctions (only for the projections), one l cutoff for each atom type. The second is a threshold value below which wavefunction coefficients are neglected. By default, the two options are unused.
MINCPW (MBASIS)
+MINCPW (WFPROD)

(**) The keyword ``MINCPW`` modifies the planewave wavefunction coefficients such that their absolute values become as small as possible. This reduces the error in the evaluation of the wavefunction products arising from the finite G cutoff and leads to better convergence with respect to ``GCUT`` and smoother band structures, when ``APPROXPW`` is used. The coefficients can be modified, because the set of IPWs is overcomplete. If set, the eigenvectors of the overlap matrix with the ``n`` smallest eigenvalues are used for this modification. Note that it only makes sense to use this together with ``APPROXPW``.
+(**) The keyword ``MINCPW`` modifies the planewave wavefunction coefficients such that their absolute values become as small as possible. This reduces the error in the evaluation of the wavefunction products arising from the finite G cutoff and leads to better convergence with respect to ``GCUT`` and smoother band structures, when ``APPROXPW`` is used. The coefficients can be modified, because the set of IPWs is overcomplete. If set, the eigenvectors of the overlap matrix with the ``n`` smallest eigenvalues are used for this modification. Note that it only makes sense to use this together with ``APPROXPW``.
.. _wannier:
Wannier orbitals
=============================
...
+.. note:: For Wannier functions to be available, the Spex configure script for compilation has to be run with the option ``withwan``. See :numref:`installation`.
+
+For many features of the Spex code, Wannier functions are required: Wannier interpolation, spinwave calculations,
+*GT* selfenergy, Hubbard *U* parameter.
+Wannier functions are constructed by a linear combination of singleparticle states.
+They can be viewed as Fourier
+transforms of the Bloch states. While Bloch states have a definite k vector and are delocalized over real space,
+the Wannier functions are localized at specific atomic sites and delocalized in k space.
+The mathematical definition is
+
+.. math::
+ \displaystyle w_{\mathbf{R}n}(\mathbf{r}) = \frac{1}{N} \sum_\mathbf{k} e^{i\mathbf{kR}}
+ \sum_m U_{\mathbf{k}m,n} \varphi_{\mathbf{k}m}(\mathbf{r})
+ :label: wannier
+with the lattice vector :math:`\mathbf{R}` (the Wannier orbital is localized in the respective unit cell),
+the orbital index :math:`n`, the number of k vectors :math:`N`, and the transformation matrix :math:`U`.
+(The latter is unitary if there are as many Wannier orbitals as there are Bloch bands included in the
+construction. There can be more Bloch bands but not less.)
+
+Wannier functions are defined in the section ``WANNIER``:
+
+.. codeblock:: bash
+
+ SECTION WANNIER
+ ORBITALS 1 4 (sp3)
+ MAXIMIZE
+ END
+
+ORBITALS (WANNIER)
+
+This line is the main definition of Wannier orbitals giving the orbital characters for the *first guess* orbitals and
+the energy window. The energy window is defined by the first two arguments, the indices of the lower and upper band
+[summation index :math:`m` in Eq. :eq:`wannier`].
+Note that if, in the above example, the fourth state is degenerate with the fifth at some k point, the fifth state will
+automatically be included as well.
+(This is the default behavior. It can be changed to excluding the fourth band or to ignoring
+degeneracies altogether by setting preprocessor macros, see the source file "wannier.f".) After the energy window follows the
+definition of the orbital characters of the projection functions used to calculate the *firstguess* Wannier functions.
+Although strictly not guaranteed, the orbital character is usually preserved not only after projection in the
+firstguess functions but even after the maximal localization
+procedure (see :numref:`MAXIMIZE`). The orbital character is defined for each atom (round brackets) or each atom type
+(square brackets). Empty brackets, such as ``()`` or ``[]``, are allowed. For example, ``(p) () (d)``
+gives three p orbitals in the first atom, five d orbitals in the third, and none in the second. The number of bands
+in the energy window must not be smaller than the number of Wannier orbitals (but can be larger).
+
+.. listtable:: Examples
+ :widths: 42 100
+
+ *  ``ORBITALS 1 18 (s,p,d)``
+  Define s, p, and d Wannier orbitals in the first (and perhaps only) atom.
+ *  ``ORBITALS 6 11 [t2g]``
+  Define t\ :sub:`2g` Wannier orbitals in the first atom type (here two atoms) from the sixth through the 11th band.
+ *  ``ORBITALS 1 9 () (d)``
+  Define d Wannier orbitals in the second atom.
+
+The above example for the ``WANNIER`` section is for bulk silicon, which has two atoms in the unit cell. Four sp\ :sup:`3` hybrid orbitals are defined
+for the first and none for the second. (Trailing empty brackets as in ``(sp3) ()`` can be omitted.)
+In this case, the four bands of the energy window (``1 4``) are all formed by bonding sp\ :sup:`3` orbitals.
+Therefore, by projection, the resulting Wannier orbitals are evenly distributed over the two atoms.
+This is an effect of the system's symmetry and not explicitly caused by the particular Wannier definition.
+If one wants to describe the antibonding orbitals as well, one would have to define, for example, ``(sp3) (sp3)``
+or just ``[sp3]``, and increase the energy window. Since the empty states of silicon are entangled, simply
+doubling the energy window (``1 8``) does not suffice. We will discuss the case of entangled bands below.
+
+The orbital characters are defined in the same way as in Wannier90. The following orbital labels are available:
+``s``, ``p``, ``d``, ``f``, ``px``, ``py``, ``pz``, ``dxy``, ``dyz``, ``dz2``, ``dxz``, ``dx2y2``, ``t2g``, ``eg``,
+``fz3``, ``fxz2``, ``fyz2``, ``fzxy``, ``fxyz``, ``fxxy``, ``fyxy``, ``sp``, ``sp2``, ``sp3``, ``sp3d``, and ``sp3d2``.
+Consult the Wannier90 manual for details.
+In addition, there are capitalized labels (``S``, ``P``, ``D``), in which case complex spherical harmonics are used
+instead of real spherical harmonics. (Obviously, ``S`` is equivalent to ``s`` and included only for completeness.
+In practice, real and complex spherical harmonics behave identically in the present context.)
+Wannier functions (except for s orbitals) have an angular dependence.
+Sometimes it is necessary to specify a particular spatial orientation of Wannier orbitals. For this purpose, one can
+give Euler angles (in degrees), for example, ``(s,px,30,60,90,d)``, through which the Wannier orbital (of p\ :sub:`x`
+symmetry in this case) is rotated.
+
+The sogenerated Wannier function ares *firstguess* Wannier functions. To be more precise,
+the expansion coefficients :math:`U_{\mathbf{k}n,m}` are obtained from projecting the Bloch functions
+:math:`\varphi_{\mathbf{k}n}` onto localized functions with the chosen orbital character. Then, an
+orthonormalization procedure is applied to make the Wannier set orthonormal. For many purposes, the
+resulting set of Wannier functions is already usable. However, it is not yet maximally localized.
+The following keyword enforces maximal localization with Wannier90 (using the Wannier90 library),
+yielding the set of maximally localized Wannier functions (MLWFs).
+
+.. _MAXIMIZE:
+
+MAXIMIZE (WANNIER)
+
+The maximal localization procedure of Wannier90 is applied starting from the set of firstguess
+Wannier functions. The input file for Wannier90 ("wannier.win") is generated automatically by
+Spex if it is not present yet. The user can still provide his own "wannier.win" and, in this way,
+finetune the localization procedure. However, some of the parameters (``num_wan`` and
+``num_bands``) might still be changed by Spex if they are inconsistent with the Spex input
+(in ``ORBITALS``). The output of Wannier90 is written to the file "wannier.wout". This file should
+be checked for convergence of the maximal localization. For further details, consult the Wannier90
+manual.
+
+FROZEN (WANNIER)
+
+In the case of entangled bands, it is often necessary to add an extra step in the construction
+of MLWFs, the disentanglement, which is also an iterative optimization procedure.
+To this end, one defines a *frozen energy window*, which lies inside the *outer* energy window
+defined by ``ORBITALS``. The frozen energy window is defined by two energies, the lower and the
+upper bound. The lower bound is assumed to coincide with the minimum energy of the lowest band
+of the outer window, i.e., the first argument of ``ORBITALS``. If ``FROZEN`` is given without an
+argument, the upper bound is *guessed* by Spex. This guess can be overridden by giving the upper
+energy as an argument to ``FROZEN``. Alternatively, one can specify the Wannier90 parameters
+``dis_froz_min`` and ``dis_froz_max`` in the file "wannier.win". Again, you should check the
+convergence of the disentanglement procedure in the file "wannier.wout". For further details,
+consult the Wannier90 manual.
+
+.. listtable:: Examples
+ :widths: 20 100
+
+ *  ``FROZEN``
+  Run disentanglement procedure. Frozen window boundaries are guessed by Spex.
+ *  ``FROZEN 5eV``
+  Run disentanglement procedure and use 5 eV as upper boundary.
+
+
+IRREP (WANNIER)
+
+(*) The set of Wannier functions should reflect the symmetry of the system. This can be used
+to define matrix representations that determine how the Wannier functions transform when a
+given symmetry operation is applied. These matrix representations can be used to speed
+up the evaluation of certain Wannierexpanded quantities, for example, in spinwave and
+*GT* calculations. (The name ``IRREP`` is a bit misleading because the matrices are
+not irreducible in general.) Unfortunately, Wannier90 tends to break the symmetry of the
+Wannier set, so ``IRREP`` can, in most cases, not be applied to MLWFs. However, unless the
+``ORBITALS`` definition itself does not break the symmetry, the firstguess Wannier functions
+reflect the full symmetry, and ``IRREP`` can be used.
+
+RESTART
+
+If ``RESTART`` is set, Spex writes the Wannier functions (the *U* matrix) to the file "spex.uwan"
+if it does not exist. Otherwise, the Wannier functions are read from "spex.uwan". In the same way
+the file "spex.kolap" is written or read containing the overlap of wavefunctions at neighboring
+k points. These quantities are needed for the maximalization procedure. (The file "spex.kolap" can be reused
+if the orbitals in ``ORBITALS`` are changed but not the outer energy window, i.e., the first two
+arguments.)
+
+UREAD (WANNIER)
+
+(**) This keyword makes Spex read Wannier functions generated by Fleur. (Works only with old version of Fleur. Currently
+unmaintained.)
+
+INTERPOL (WANNIER)
+
+The keyword INTERPOL has already been discussed in :numref:`INTERPOL_GW`. It is not related to the
+Wannier construction, but we explain it here briefly for completeness.
+When given, Spex performs a Wannier interpolation for the reference meanfield system (output to
+the file "bands0") and of the calculated *GW* (HF, COSX, et cetera) quasiparticle energies
+(output to the file "bands1").
diff git a/source/common_keywords.rst b/source/common_keywords.rst
index 77655f82512444ebdbbe5bb0da09d088ae260e75..a8e4bda3c84e59e37cdb46f26f5f69856440a34f 100644
 a/source/common_keywords.rst
+++ b/source/common_keywords.rst
@@ 4,32 +4,35 @@
Common Keywords
****************
.. warning:: This section discusses several settings that are common for all calculations. Some of them are explained in more detail in later sections.
+This section discusses several settings that are common for all calculations. Some of them are explained in more detail in later chapters.
General Settings
=================
MEM

Calculations based on manybody perturbation theory are costly, not only in terms of computation time but also in terms of memory. Since the problem sizes (e.g., number of atoms) are thus effectively restricted by the memory hardware of the computer, the latter constitutes in some ways a more serious limit than the computation time the user can afford to invest. Therefore, an upper limit for the memory that Spex is allowed to use can be defined with the keyword MEM x with x in MB. The default is MEM 1000. In parallel Spex, x is the upper limit per node. If possible, Spex will divide the work load into packets in such a way that the memory usage stays below the limit. (Actually, Spex might use a bit more memory than specified.)

++++
 Example  MEM 2000  Restrict memory usage to 2GB (Default: 1GB) 
++++

+Calculations based on manybody perturbation theory are costly, not only in terms of computation time but also in terms of memory. Since the problem sizes (e.g., number of atoms) are thus effectively restricted by the memory hardware of the computer, the latter constitutes in some ways a more serious limit than the computation time the user can afford to invest. Therefore, an upper limit for the memory that Spex is allowed to use can be defined with the keyword ``MEM x`` with x in MB. The default is ``MEM 1000``.
+In parallel Spex, the value would correspond to the upper limit per node. If possible, Spex will divide the work load into packets in such a way that the memory usage stays below the limit. (Actually, Spex might use a bit more memory than specified.)
+
+.. listtable:: Example
+ :widths: 16 100
+
+ *  ``MEM 2000``
+  Restrict memory usage to 2GB (Default: 1GB).
+
.. _job:
JOB

Each Spex run needs a job definition, which defines what Spex should do, e.g., ``JOB GW 1:(15)`` for a ``GW`` calculation of the specified set of bands. The available jobs are
+Each Spex run needs a job definition, which defines what Spex should do, e.g., ``JOB GW 1:(15)`` for a *GW* calculation of the specified set of bands. The available jobs are
* ``GW``  GW method
+* ``GW``  *GW* method
* ``HF``  HartreeFock method
* ``PBE0``  PBE0 hybrid functional
* ``SX``  Screened exchange
* ``COSX``  Coulombhole screenedexchange (COHSEX)
* ``GT``  GT method
* ``GWT``  GW and GT combined
+* ``GT``  *GT* method
+* ``GWT``  *GW* and *GT* combined
* ``HFENE``  Total HartreeFock exchange energy
* ``RPAENE``  RPA correlation energy
* ``SUSCEP``  Susceptibility (Polarization function)
@@ 41,43 +44,47 @@ Each Spex run needs a job definition, which defines what Spex should do, e.g., `
Details of these jobs are explained in subsequent sections. The job definition must not be omitted but may be empty: ``JOB``, in which case Spex will just read the wavefunctions and energies, perform some checks and some elemental calculations (e.g., Wannier interpolation), and stop. In principle, Spex supports multiple jobs such as ``JOB GW 1:(15) DIELEC 1:{0:1,0.01}``. This feature is, however, seldom used and is not guaranteed to work correctly in all versions.
++++
 Examples  ``JOB COSX 1:(15)``  Perform COHSEX calculation. 
++++
  ``JOB COSX 1:(15) SCREEN 2:{0:1,0.01}``   Subsequently, calculate the screened 
    interaction on a frequency mesh 
++++
  ``JOB``  Just perform some checks and stop 
++++

BZ


.. _wrtkpt:

WRTKPT

There are only two keywords that may not be omitted from the input file. The first is ``JOB``, the second ``BZ``, which defines the kpoint set for the sampling of the Brillouin zone. The keyword requires four integer arguments, BZ l m n, defining an lxmxn kpoint set. If the keyword ``WRTKPT`` is set (or, alternatively, with the command line option ``spex w``), Spex uses the ``BZ`` definition to construct the irreducible wedge of the kpoint set and writes the respective list of k points to a file, which is then read by the DFT code.

++++
 Example  ``BZ 8 4 2``  Define a 8x4x2 kpoint sampling 
++++

+.. listtable:: Examples
+ :widths: 95 100
+
+ *  ``JOB COSX 1:(15)``
+  Perform COHSEX calculation.
+ *  ``JOB COSX 1:(15) SCREEN 2:{0:1,0.01}``
+  Subsequently, calculate the screened interaction on a frequency mesh.
+ *  ``JOB``
+  Just perform some checks and stop.
+
+.. _WRTKPT:
+
+BZ and WRTKPT
+
+There are only two keywords that must not be omitted from the input file. The first is ``JOB``, the second ``BZ``, which defines the kpoint set for the sampling of the Brillouin zone. The keyword requires four integer arguments, ``BZ l m n``, defining an :math:`l\times m\times n` kpoint set. If the keyword ``WRTKPT`` is set (or, alternatively, with the command line option ``spex w``), Spex uses the ``BZ`` definition to construct the irreducible wedge of the kpoint set and writes the respective list of k points to a file, which is then read by the DFT code.
+
+.. listtable:: Example
+ :widths: 16 100
+
+ *  ``BZ 8 4 2``.
+  Define a :math:`8\times 4\times 2` kpoint sampling.
+
RESTART

This is an important keyword. It enables (a) continuing a calculation that has, for some reason, stopped or crashed and (b) reusing data of a previous calculation. If specified, Spex will store intermediate results in a file or read them if the particular data is already present in the file. Even if RESTART was not specified in the first run, Spex might still be able to reuse data from a previous run with ``RESTART 2``. (This option usually gives less freedom for a change of parameters. ``RESTART 1`` and RESTART are equivalent.) It depends on the particular job definition whether ``RESTART`` and/or ``RESTART 2`` are available and how they work. See below for details.

++++
 Examples  ``RESTART``  Read/write restart file 
++++
  ``RESTART 2``  Try to reuse data from standard output files 
++++

+This is an important keyword. It enables (a) continuing a calculation that has, for some reason, stopped or crashed and (b) reusing data of a previous calculation. If specified, Spex will store intermediate results in a file or read them if the particular data is already present in the file. Even if ``RESTART`` was not specified in the first run, Spex might still be able to reuse data from a previous run with ``RESTART 2``. (This option usually gives less freedom for a change of parameters. ``RESTART 1`` and ``RESTART`` are equivalent.) It depends on the particular job definition whether ``RESTART`` and/or ``RESTART 2`` are available and how they work. See below for details.
+
+.. listtable:: Examples
+ :widths: 18 100
+
+ *  ``RESTART``
+  Read/write restart file.
+ *  ``RESTART 2``
+  Try to reuse data from standard output files.
+
STOREIBZ

Spex reads wavefunctions and energies as a first step. In general, these have to be stored for all k vectors in the set (e.g., for all 64 points of a 4x4x4 set), filling up a lot of memory. This can be restricted to the irreducible zone (only eight points in the example) with the keyword ``STOREIBZ``. Spex then regenerates the wavefunctions on equivalent k points on the fly. This might slow down the calculation but is more economical with computer memory.
+Spex reads wavefunctions and energies as a first step. In general, these have to be stored for all k vectors in the set (e.g., for all 64 points of a :math:`4\times 4\times 4` set), filling up a lot of memory. This can be restricted to the irreducible zone (only eight points in the example) with the keyword ``STOREIBZ``. Spex then regenerates the wavefunctions on equivalent k points on the fly. This might slow down the calculation but is more economical with computer memory.
.. note:: ``enableload`` (*) Another possibility to reduce the memory consumption is to reconfigure (and recompile) the code with the option enableload (or recompile with the flag DLOAD in the Makefile). With this option, Spex stores the wavefunctions in memory only temporarily and rereads them from harddisc whenever they are needed. Obviously, this requires a fast harddisc and network connection. The HDF5 library is used for the data transfer. There is no possibility to set the option in the input file because it strongly affects the way the wavefunction data is handled in the code.
+enableload
+
+(*) Another possibility to reduce the memory consumption is to reconfigure (and recompile) the code with the option ``enableload`` (or recompile with the flag ``DLOAD`` in the Makefile). With this option, Spex stores the wavefunctions in memory only temporarily and rereads them from harddisc whenever they are needed. Obviously, this requires a fast harddisc and network connection. The HDF5 library is used for the data transfer. There is no possibility to set the option in the input file because it strongly affects the way the wavefunction data is handled in the code.
Input data
===========
@@ 87,78 +94,129 @@ Spex needs data from a previous selfconsistent solution of a meanfield system
NBAND

The Green function involves a summation over electronic bands, as do the polarization function and the selfenergy. The convergence with respect to the number of bands is a very important aspect in calculations based on manybody perturbation theory. The keyword NBAND n sets the number of bands to n. Obviously, n must not be larger than the number of bands provided by the input data. It should be noted that the actual number of bands can differ from n and also from k point to k point because Spex makes sure that degenerate subspaces are not cut in order to avoid symmetry breaking. If n is a real number (i.e., if it contains a decimal point: 3.0), it is interpreted as the maximal band energy. Then, all bands up to this energy are included. This is helpful if one wants to use a single parameter for calculations with differently sized unit cells. For example, when the unit cell doubles, one would have to include twice as many bands to reach the same accuracy, while the maximal band energy practically remains unchanged. If NBAND is unset, all available bands are used.
+The Green function involves a summation over electronic bands, as do the polarization function [Eq. :eq:`eqP`] and the selfenergy. The convergence with respect to the number of bands is a very important aspect in calculations based on manybody perturbation theory. The keyword ``NBAND n`` sets the number of bands to n. Obviously, n must not be larger than the number of bands provided by the input data. It should be noted that the actual number of bands can differ from n and also from k point to k point because Spex makes sure that degenerate subspaces are not cut in order to avoid symmetry breaking. If n is a real number (i.e., if it contains a decimal point: ``3.0``), it is interpreted as the maximal band energy. Then, all bands up to this energy are included. This is helpful if one wants to use a single parameter for calculations with differently sized unit cells. For example, when the unit cell doubles, one would have to include twice as many bands to reach the same accuracy, while the maximal band energy (practically) remains unchanged. If ``NBAND`` is unset, all available bands are used.
++++
 Examples  NBAND 100  Use 100 bands per k point 
++++
  NBAND 3.0  Use all bands up to 3 Hartree 
++++
  NBAND 50.0eV  Use all bands up to 50 eV 
++++
+.. listtable:: Examples
+ :widths: 24 100
+
+ *  ``NBAND 100``
+  Use 100 bands per k point.
+ *  ``NBAND 3.0``
+  Use all bands up to 3 Hartree.
+ *  ``NBAND 50.0eV``
+  Use all bands up to 50 eV.
CHKOLAP

+.. _CORESOC:
CHKMISM
+CORESOC

(*) The set of wavefunctions read from harddisc can be checked in terms of their orthonormality and continuity (mismatch) at the muffintin (MT) boundary with the keywords CHKOLAP and CHKMISM, respectively. The former test gives a value for the deviation from exact orthonormality. The smaller the value, the better the condition of orthonormality is fulfilled. If the deviation is too large, the program stops. The second keyword yields for each k point the average and maximal MT mismatch of the wavefunctions. If any of these tests yield too large errors, there could be a problem with the readin process or with the data written by the DFT code. The developers should be informed in this case.
+By default, Spex averages over the SOCsplitted states. For example, the two states 2p\ :sup:`1/2` and the four states 2p\ :sup:`3/2` are
+averaged to give the six quasinonrelativistic states of the 2p shell (or the three spinup 2p states and the three spindown 2p states
+in the case of spin polarization). Then, all core states of a given shell have the same energy and the same radial form of the wavefunction. The averaging procedure
+can be avoided by specifying ``CORESOC``. The SOCsplitted states and the different energies and wavefunctions are then retained and taken
+into account in the calculations of the selfenergy, polarization function, et cetera. In the case of spin polarization, the Weiss field gives
+rise to further lifting of degeneracies. Spex calculates the respective new splittings, energies, and wavefunctions.
+For reference, the corestate energies are written to the output file. Usually, the effect of ``CORESOC`` is numerically small, but it
+still might play an important role for *GW* calculations in the case of very small band gaps.
+
+
+CHKMISM and CHKOLAP
+
+(*) The set of wavefunctions read from harddisc can be checked in terms of their orthonormality and continuity (mismatch) at the muffintin (MT) boundary with the keywords ``CHKOLAP`` and ``CHKMISM``, respectively. The former test gives a value for the deviation from exact orthonormality. The smaller the value, the better the condition of orthonormality is fulfilled. If the deviation is too large, the program stops.
+The second keyword yields, for each k point, the average and maximal MT mismatch of the wavefunctions. If any of these tests yield too large errors, there could be a problem with the readin process or with the data written by the DFT code. The developers should be informed in this case.
MTACCUR

(*) 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 highlying 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 spinpolarized 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}`)

++++
 Examples  ``MTACCUR 1 2``  Calculate MT representation error between 1 and 2 Hartree 
++++
  ``MTACCUR``  Automatic energy range 
++++
+(*) 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 highlying 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 spinpolarized 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}`)
+
+.. listtable:: Examples
+ :widths: 24 100
+
+ *  ``MTACCUR 1 2``
+  Calculate MT representation error between 1 and 2 Hartree.
+ *  ``MTACCUR``
+  Automatic energy range.
BANDINFO
+
+(*) This keyword lets Spex analyze the orbital character of the wavefunctions. The results are given, first,
+in the form of a table in the output file containing the orbital decomposition for each eigenstate and,
+second, as projected densityofstates (DOS) plots in the files "spex.dos.NNN" (where NNN is a threedigit running index).
+As the table can be very large, there is a possibility to choose the atom type index at whose spheres the analysis is to be performed.
+
+.. listtable:: Examples
+ :widths: 20 100
+
+ *  ``BANDINFO``
+  Calculate orbital decomposition and DOS plots.
+ *  ``BANDINFO 2``
+  Restrict the orbitaldecomposition table to atom type 2.
+
+ENERGY

(*) In some cases, it may be necessary to replace the energy eigenvalues, provided by the meanfield (DFT) code, by energies (e.g., GW quasiparticle energies) obtained in a previous Spex calculation, for example, to determine the GW Fermi energy or to perform energyonly selfconsistent calculations. This can be achieved with the keyword ``ENERGY file``, where ``file`` contains the new energies in eV. The format of ``file`` corresponds to the output of the ``spex.extr`` utility: ``spex.extr g spex.out > file``. It must be made sure that ``file`` contains energy values for the whole irreducible Brillouin zone. Band energies not contained in ``file`` will be adjusted so that the energies are in accending order (provided that there is at least one energy value for the particular k point).

++++
 Example  ``ENERGY energy.inp``   Replace the meanfield energy eigenvalues by the 
    energies provided in the file ``energy.inp`` 
++++

+(*) In some cases, it may be necessary to replace the energy eigenvalues, provided by the meanfield (DFT) code, by energies (e.g., *GW* quasiparticle energies) obtained in a previous Spex calculation, for example, to determine the *GW* Fermi energy or to perform energyonly selfconsistent calculations. This can be achieved with the keyword ``ENERGY file``, where "file" contains the new energies in eV. The format of "file" corresponds to the output of the ``spex.extr`` utility: ``spex.extr g spex.out > file``. It must be made sure that "file" contains energy values for the whole irreducible Brillouin zone. Band energies not contained in "file" will be adjusted so that the energies are in accending order (provided that there is at least one energy value for the particular k point).
+
+.. listtable:: Example
+ :widths: 34 100
+
+ *  ``ENERGY energy.inp``
+  Replace the meanfield energy eigenvalues by the energies provided in the file "energy.inp".
+
DELTAEX

(*) This keyword modifies the exchange splitting of a collinear magnetic system, i.e., it shifts spinup and spindown energies relative to to each other so as to increase or decrease the exchange splitting. With ``DELTAEX x``, the spinup (spindown) energies are lowered (elevated) by x/2. The parameter x can be used to enforce the Goldstone condition in spinwave calculations [Phys. Rev. B 94, 064433 (2016)]
++++
 Example  ``DELTAEX 0.2eV``   Increase the exchange splitting by 0.2eV 
    (spinup/down energies are decreased/increased by 0.1eV) 
++++

+.. listtable:: Example
+ :widths: 26 100
+
+ *  ``DELTAEX 0.2eV``
+  Increase the exchange splitting by 0.2 eV (spinup/down energies are decreased/increased by 0.1 eV)
+
PLUSSOC

(*) With this keyword, Spex adds spinorbit coupling (SOC) to the meanfield solution using a oneshot secondvariation approach and stops, so it is noniterative but updates the wavefunctions and energies. This can be used, for example, to include SOC to a previous ``GW`` calculation that has been carried out without SOC. The new energies are written to standard output and, together with the new wavefunctions, to new input files on the harddisc, thus overwriting the old input data. In the case of Fleur, the modified files are ``gwa``, ``KS.hdf``, and ``spex.sym``. The latter is written only if the system is magnetic (noncollinear magnetism). It has the same form as ``sym.out``. Spex reads the file "spex.sym" if it exists and "sym.out" otherwise. After the ``PLUSSOC`` calculation, a subsequent Spex run will read the new input data, now including SOC. (The keyword PLUSSOC must be removed for this new run, otherwise Spex stops with an error, since it cannot add SOC to a solution that contains it already.)
+(*) With this keyword, Spex adds spinorbit coupling (SOC) to the meanfield solution using a oneshot secondvariation approach and stops,
+so it is noniterative but updates the wavefunctions and energies.
+This can be used, for example, to include SOC to a previous *GW* calculation that has been carried out without SOC.
+The new energies are written to standard output and, together with the new wavefunctions, to new input files on the harddisc,
+thus overwriting the old input data. In the case of Fleur, the modified files are ``gwa``, ``KS.hdf``, and ``spex.sym``.
+The latter is written only if the system is spin polarized (noncollinear magnetism). It has the same form as ``sym.out``.
+Spex reads the file "spex.sym" if it exists and "sym.out" otherwise. After the ``PLUSSOC`` calculation, a subsequent
+Spex run will read the new input data, now including SOC. (The keyword ``PLUSSOC`` must be removed for this new run, otherwise
+Spex stops with an error, since it cannot add SOC to a solution that contains it already.)
ITERATE

(*) If specified, Spex only reads the LAPW basis set from the input data, provided by the meanfield (DFT) code, but performs the diagonalization of the Hamiltonian at the k points itself. This calculation effectively replaces the second run of the DFT code. In this sense, the name of the keyword is a bit misleading, as the calculation is noniterative. The keyword ``ITERATE`` is mostly intended for testing and debugging. It is not available for executables compiled with ``DLOAD`` (configured with ``enableload``).

++++
 Examples  ``ITERATE NR``  Diagonalize a nonrelativistic Hamiltonian 
++++
  ``ITERATE SR``  Use scalarrelativity 
++++
  ``ITERATE FR``  Also include SOC 
++++
  ``ITERATE SR 1``   Diagonalize scalarrelativistic Hamiltonian and 
    neglect eigenvalues below 1 htr 
++++
  ``ITERATE FR STOP``  Diagonalize relativistic Hamiltonian (including SOC), then stop 
++++

Parallelized version
====================
The parallelized version of Spex uses the MPI standard 3.1. It can be run on several CPUs on the same node or on several nodes with the command mpirun or mpiexec or whatever MPI launcher your computer system uses. In principle, there are no restrictions with respect to the number of processes. A better performance is expected, though, if the number of processes is not a prime number but has a long prime factorization because this gives the code more freedom to distribute the work among the processes.

MPISPLIT

(*) The sharedmemory functionality of MPI 3.1 is used for several big arrays, which allows the same memory region to be accessed by several MPI processes. By default, all processes running on one node share the memory. It can be reasonable to change this behavior to, e.g., having processes on the same socket to share memory. This is possible with ``MPISPLIT SHRD=NODE`` (default), ``MPISPLIT SHRD=SOCKET`` (only works with ``OpenMPI``), or ``MPISPLIT SHRD=n``, where n is the number of the sharedmemory processes. Using this option increases the memory consumption but might be advantageous in terms of computation time. (There is not much experience yet with this option.)
+(*) If ``ITERATE`` is specified, Spex only reads the LAPW basis set from the input data, provided by the meanfield (DFT) code,
+but performs the diagonalization of the Hamiltonian at the k points itself. This calculation effectively replaces
+the second run of the DFT code. In this sense, the name of the keyword is a bit misleading, as the calculation is noniterative.
+The keyword ``ITERATE`` is mostly intended for testing and debugging. It is not available for executables compiled with ``DLOAD``
+(configured with ``enableload``).
+
+.. listtable:: Examples
+ :widths: 30 100
+
+ *  ``ITERATE NR``
+  Diagonalize a nonrelativistic Hamiltonian.
+ *  ``ITERATE SR``
+  Use scalarrelativity.
+ *  ``ITERATE FR``
+  Also include SOC.
+ *  ``ITERATE SR 1``
+  Diagonalize scalarrelativistic Hamiltonian and neglect eigenvalues below 1 htr.
+ *  ``ITERATE FR STOP``
+  Diagonalize relativistic Hamiltonian (including SOC), then stop.
+
+
diff git a/source/conf.py b/source/conf.py
index c292b38e6d9e1cd75c478405bb8d45ae1409e5e1..e566917131ae275ae1bc60d303653b7a980fecc7 100644
 a/source/conf.py
+++ b/source/conf.py
@@ 21,15 +21,14 @@ def setup(app):
#  Project information 
project = 'spex'
copyright = '2018, PGI, Forschungszentrum Juelich, Germany'
author = 'spex team'
+project = 'Spex'
+copyright = u'2019, PGI, Forschungszentrum J\u00fclich, Germany'
+author = ''
# The short X.Y version
version = '0.0.1'
# The full version, including alpha/beta/rc tags
release = '0.0.1'

+release = open('version').read()
#  General configuration 
@@ 106,6 +105,13 @@ html_theme_options = {
# so a file named "default.css" will overwrite the builtin "default.css".
html_static_path = ['_static']
+html_context = {
+ 'css_files': [
+ '_static/theme_overrides.css', # override wide tables in RTD theme
+ ],
+ }
+
+
# Custom sidebar templates, must be a dictionary that maps document names
# to template names.
#
@@ 147,8 +153,8 @@ latex_elements = {
# (source start file, target name, title,
# author, documentclass [howto, manual, or own class]).
latex_documents = [
 (master_doc, 'spex.tex', 'spex Documentation',
 'Spex Team', 'manual'),
+ (master_doc, 'spex.tex', 'Spex Manual',
+ '', 'manual'),
]
@@ 157,7 +163,7 @@ latex_documents = [
# One entry per manual page. List of tuples
# (source start file, name, description, authors, manual section).
man_pages = [
 (master_doc, 'spex', 'spex Documentation',
+ (master_doc, 'Spex', 'Spex Manual',
[author], 1)
]
@@ 168,8 +174,8 @@ man_pages = [
# (source start file, target name, title, author,
# dir menu entry, description, category)
texinfo_documents = [
 (master_doc, 'spex', 'spex Documentation',
 author, 'spex', 'One line description of project.',
+ (master_doc, 'Spex', 'Spex Manual',
+ author, 'Spex', 'One line description of project.',
'Miscellaneous'),
]
@@ 195,5 +201,5 @@ epub_copyright = copyright
epub_exclude_files = ['search.html']
latex_elements = {
 'releasename' : 'Version'
 }
+ 'releasename' : 'Version'
+ }
diff git a/source/general_information.rst b/source/general_information.rst
index ec8ce729816c7b6119a324cf2b8df9cef1d4f241..43ee7ce9157f0722934aaa20a19fa966b119763e 100644
 a/source/general_information.rst
+++ b/source/general_information.rst
@@ 4,9 +4,13 @@
General information
===================
Although the Spex executable can have different names (e.g., ``spex.inv``, ``spex.noinv``), we assume that spex is the name of the executable in the following. (This can be achieved by using the caller utility, see the file ``INSTALL``.)
+Although the Spex executable can have different names (e.g., "spex.inv", "spex.noinv"),
+we assume that spex is the name of the executable in the following.
+(This can be achieved by using the caller utility, see the file "INSTALL".)
.. note:: Paragraphs discussing advanced options are preceded with (*), and the ones about obsolete, unmaintained, or experimental options are marked with (**). You can safely skip the paragraphs marked with (*) at first reading.
+.. note:: Paragraphs discussing advanced options are preceded with (*), and the ones about obsolete, unmaintained, or experimental options are marked with (**). You can safely skip the paragraphs marked with (*) and (**) at first reading.
+
+.. _cmdlineopt:
Commandline options
+++++++++++++++++++++++++++++
@@ 14,24 +18,24 @@ The Spex code recognizes several commandline options (``spex [OPTIONS]``):
* ``version``: output version information and exit
* ``help``: display help and exit
* ``inp=file``: uSie wollen den Besitz und Betrieb von Medienunternehmen doch nicht ernsthaft mit der Entgegennahme von unversteuerten Geschenken von auslndischen Politikern vergleichen.se file as input file instead of the default "spex.inp"
+* ``inp=file``: use file as input file instead of the default "spex.inp"
* ``w``: set ``WRTKPT`` automatically
* ``o=opt``: use special option ``opt`` (for debugging and testing).
+* ``o=opt``: use special integer option parameter ``opt`` (for debugging and testing).
Input file
++++++++++++
By default, Spex reads input parameters from the file ``spex.inp``. The syntax is as follows.
+By default, Spex reads input parameters from the file "spex.inp". The syntax is as follows.
* The file ``spex.inp`` does not use a fixed format, and the keywords may be given in any order.
+* The file "spex.inp" does not use a fixed format, and the keywords may be given in any order.
* Each keyword is given on one line together with its parameters following it.
* Keywords must not be specified more than once.
* Some keywords are grouped in sections. A section starts with ``SECTION`` sectionname and ends with END.
+* Some keywords are grouped in sections. A section starts with ``SECTION sectionname`` and ends with ``END``.
* All keywords (and section names) are in capital letters.
* Everything after a ``#`` sign is interpreted as a comment. Empty lines are allowed.
* Lines are only interpreted up to the 80th column. Line continuation is possible using a single backslash ``\`` at the end of the previous line.
* Most keyword parameters have default values that are used if the keyword is not specified.
* By default, parameters are given in atomic units. In the case of energies, it is possible to use eV instead, e.g., 1.0eV.
* The input file can be truncated with the keyword EXIT. Everything after EXIT is ignored.
+* By default, parameters are given in atomic units. In the case of energies, it is possible to use eV instead, e.g., ``1.0eV``.
+* The input file can be truncated with the keyword ``EXIT``. Everything after ``EXIT`` is ignored.
* Currently, there are six section keywords: ``MBASIS``, ``WFPROD``, ``COULOMB``, ``SUSCEP``, ``SENERGY``, and ``WANNIER``.
An example for a section in the input file is
@@ 51,31 +55,38 @@ An example for a section in the input file is
Output files
++++++++++++
The main output file of a Spex calculation is written directly to standard output (stdout). Errors, warnings, and (additional) info(rmation) are written to standard error (stderr). Errors stop the program run. Errors, warnings, and infos are of the form::
+The main output file of a Spex calculation is written directly to standard output (stdout).
+Errors, warnings, and (additional) info(rmation) are written to standard error (stderr).
+Errors stop the program run. Errors, warnings, and infos are of the form::
 routine: Error message.
 routine: Warning! Warning message.
 routine: Info! Info message.
+ SPEXERROR (source.f:0123) Error message.
+ SPEXBUG (source.f:0123) Error message.
+ SPEXWARNING (source.f:0123) Warning message.
+ SPEXINFO (source.f:0123) Info message.
respectively, where routine is the name of the subroutine where the error (etc.) has occurred. An error message that ends on ... (bug?) indicates a possible bug in the code (please inform the developers).
+respectively, where "source.f" is the name of the source file and 0123 is the respective line where the error (etc.)
+has occurred.
+An error message starting with ``SPEXBUG`` is likely caused by a bug in the code (please inform the developers).
+A warning message indicates a possible problem in the calculation, but the calculation continues anyway.
+An info message informs about a less important issue, for example, about the usage of obsolete syntax in the input file.
It is recommendable to pipe the output to a file (in Bash)
Pipes stdout to spex.out and stderr to the screen::
+Pipes stdout to "spex.out" and stderr to the screen::
$ spex > spex.out
Pipes stdout and stderr to spex.out::
+Pipes stdout and stderr to "spex.out"::
$ spex &> spex.out
Pipes stdout to spex.out and stderr to spex.err::
+Pipes stdout to "spex.out" and stderr to "spex.err"::
$ spex > spex.out 2> spex.err
+For the parallelized version, Spex has to be run through an MPI launcher, for example,
+``mpiexec np 8 spex``. See :numref:`parallel` for details.
Depending on the job type and keywords, there are additional output files containing, for example, the resulting spectra. These output files can be in ASCII and binary format. Further details are given in the :ref:`tutorials`.



+Depending on the job type and keywords, there are additional output files containing, for example, the resulting spectra.
+These output files can be in ASCII and binary format. Further details are given in :numref:`tutorials`.
diff git a/source/getting_started.rst b/source/getting_started.rst
index 243ba47020d5a6e41ce70c2174f924008e0afb31..7f9fc387ceb762fc171c84734aa736243f9a8b6d 100644
 a/source/getting_started.rst
+++ b/source/getting_started.rst
@@ 6,21 +6,28 @@ Getting Started
This section is intended to give a first impression on how a typical calculation is carried out. Details are given in later sections.
Manybody perturbation theory is defined with respect to a formally unperturbed system, i.e., a system without electronelectron interaction. Any meanfield theory, e.g., HartreeFock or KohnSham (KS) DFT, provides such a system. Thus before running Spex we must find a selfconsistent meanfield solution. Normally this first step is fast and we can use a dense kpoint set. For the calculation with Spex we have to generate a new (and usually coarser) kpoint set and calculate the corresponding wave functions and energies with the meanfield program in a second run. Thus a typical procedure is
+Manybody perturbation theory is defined with respect to a formally unperturbed system, i.e., a system without electronelectron interaction.
+Any meanfield theory, e.g., HartreeFock or KohnSham (KS) DFT, provides such a system. Thus, before running Spex, we must find a selfconsistent
+meanfield solution. Normally this first step is fast and we can use a dense kpoint set. For the calculation with Spex, we have to generate a new
+(and usually coarser) kpoint set and calculate the corresponding wavefunctions and energies with the meanfield program in a second run.
+Thus, a typical procedure is
1. selfconsistent DFT calculation (e.g. with Fleur),
2. generate new kpoint set with Spex,
3. second run of the DFT program to generate the eigenstates at the new k points,
4. finally run Spex.
For the first and third step consult the manual of your DFT program. (Click `here`_ for the Fleur manual and also read the section `Spex and Fleur`_. You must set ``gw=1`` and ``gw=2`` for the first and third step, respectively.) For the second and fourth step you need to create an input file ``spex.inp`` for Spex. A very simple input file for a GW calculation of Si is this:
+Of course, for subsequent Spex runs, steps 13 need not be repeated, unless the reference meanfield system is to be
+recalculated with different parameters (e.g., different kpoint set or different crystal structure).
+For the first and third step consult the manual of your DFT program. (Click `here`_ for the Fleur manual and also read
+:numref:`SpexandFleur`. You must set ``gw=1`` and ``gw=2`` for the first and third step, respectively.)
+For the second and fourth step you need to create an input file "spex.inp" for Spex. A very simple input file for a *GW* calculation of Si is this:
.. _here: http://www.flapw.de/pm/index.php?n=UserDocumentation.FLEUR
.. _Spex and Fleur: http://www.flapw.de/pm/index.php?n=UserDocumentation.SPEXFLEUR#Fleur
.. codeblock:: bash
 # Quasiparticle energies of bands at the highsymmetry points Gamma, X and L.
+ # Quasiparticle energies of bands at the highsymmetry points Gamma, X, and L.
BZ 4 4 4
WRTKPT
@@ 28,10 +35,11 @@ For the first and third step consult the manual of your DFT program. (Click `her
KPT X=[0,0,1] L=1/2*[1,1,1]
JOB GW 1:(1,2,5) X:(1,3,5) L:(13,5)
.. note:: Refer to the tutorials and other chapter for a full list of available keywords
+.. note:: Refer to the tutorials and other chapters for a full list of available keywords
For the second step only the lines beginning with ``BZ`` and ``WRTKPT`` are necessary. The line ``BZ 4 4 4`` defines a 4x4x4 kpoint set, and the keyword ``WRTKPT`` tells Spex to write the (irreducible) kpoints to a file (``kpts`` for Fleur) and stop.
+For the second step, only the lines beginning with ``BZ`` and ``WRTKPT`` are necessary. The line ``BZ 4 4 4`` defines a :math:`4\times 4\times 4` kpoint set, and the keyword ``WRTKPT`` tells Spex to write the (irreducible) kpoints to a file ("kpts" for Fleur) and stop.
For the fourth step you must remove (or comment out) the keyword ``WRTKPT``. Now all other parameters are needed. The results are written to standard output and should be piped to a :ref:`spex.out` file. As long as you do not change the kpoint set, the first three steps need not be repeated for a new run of Spex.
+For the fourth step, you must remove (or comment out) the keyword ``WRTKPT``. Now all other parameters are needed.
+The results are written to standard output and should be piped to an output file (:numref:`spex.out`).
.. note:: As an alternative to ``WRTKPT``, the command line option w can be used (``spex w``). It sets the keyword ``WRTKPT`` automatically.
+.. note:: As an alternative to ``WRTKPT``, the command line option ``w`` can be used (``spex w``). It sets the keyword ``WRTKPT`` automatically (:numref:`cmdlineopt`).
diff git a/source/index.rst b/source/index.rst
index eb2793c1c1978069aa48bd243e2b15fdfcd247b8..51a16edd458fa394031e8fe3fe6afc7ec7f98c5e 100644
 a/source/index.rst
+++ b/source/index.rst
@@ 3,27 +3,37 @@
You can adapt this file completely to your liking, but it should at least
contain the root `toctree` directive.
Welcome to spex's documentation!
+Spex Manual
================================
Spex is a computer code based on manybody perturbation theory. It uses the allelectron fullpotential linearized augmented planewave method (FLAPW), which provides an accurate basis set for all kinds of materials including transition metals, oxides and even felectron systems.
Currently Spex can calculate quasiparticle properties (oneshot and selfconsistent) using the GW approximation, EELS and optical spectra as well as total energies in the RPA approximation, and spinwave spectra from the BetheSalpeter equation. (TDDFT is implemented but unmaintained at the moment.)
+Spex is a computer code based on manybody perturbation theory.
+It uses the allelectron fullpotential linearized augmented planewave method (FLAPW),
+which provides an accurate basis set for all kinds of materials including transition metals,
+oxides, and even felectron systems.
+
+Currently Spex can calculate quasiparticle properties (oneshot and selfconsistent)
+using the *GW* approximation, EELS and optical spectra as well as total energies
+in the RPA approximation, and spinwave spectra from the BetheSalpeter equation.
+(TDDFT is implemented but unmaintained at the moment.)
It needs input from a converged DFT calculation, which can be generated by Fleur but also by any other FLAPWDFT code.
If you use SPEX for your research, please cite the following work:
+If you use Spex for your research, please cite the following work:
.. highlights:: Christoph Friedrich, Stefan BlÃ¼gel, Arno Schindlmayr, "Efficient implementation of the GW approximation within the allelectron FLAPW method", *Phys. Rev. B 81, 125102 (2010)*.
+.. highlights:: Christoph Friedrich, Stefan BlÃ¼gel, Arno Schindlmayr, "Efficient implementation of the `GW` approximation within the allelectron FLAPW method", *Phys. Rev. B 81, 125102 (2010)*.
.. toctree::
:maxdepth: 2

+ :numbered:
+
installation
getting_started
general_information
common_keywords
tutorials
basis_sets
+ parallelized_version
+ spex_and_fleur
spex_faq
diff git a/source/installation.rst b/source/installation.rst
index bc0420b87ace2132af6ed8a98fe72619183cae01..d9a3896724db44c1ca72983a2e7638426579499f 100644
 a/source/installation.rst
+++ b/source/installation.rst
@@ 6,23 +6,29 @@ Installation
The Spex source code cannot be downloaded directly. Please contact c.friedrich@fzjuelich.de
Install prerequisite software
+++++++++++++++++++++++++++++
+Prerequisite software
+++++++++++++++++++++++
For compiling the source you need:
* a modern FORTRAN compiler (the code uses some FORTRAN 2008 functionality),
+* a modern FORTRAN compiler,
* the BLAS, LAPACK, FFTW, and HDF5 libraries as well as the Wannier90 library if you need Wannier support,
* an implementation of the MPI3 standard if you intend to compile the parallel version.
Configure the spex installation
++++++++++++++++++++++++++++++++
Once you have obtained the code (``spexVERSION.tgz``) unpack it with your local version of tar and change to the ``spexVERSION`` directory, where ``VERSION`` is the version number, e.g. "05.00". There you should find the source code and a configure shell script that helps to generate a Makefile for your system. Consult the file ``INSTALL`` for a description of the compilation process. In short, if all necessary libraries are in standard directories, ``./configure`` and make should suffice to create a working executable. Special options can be given to the configure script:
+Compilation
+++++++++++++
+Once you have obtained the code ("spexVERSION.tgz") unpack it and change
+to the "spexVERSION" directory, where "VERSION" is the version number, e.g. "05.00". There, you find
+a configure shell script that helps to generate a Makefile for your system. Consult
+the file "INSTALL" for a description of the compilation process. In short, if all necessary libraries are
+in standard directories, ``./configure`` and ``make`` should suffice to create working executables.
+Special options can be given to the configure script:
* ``enableparallel``: Build parallel version (requires an MPI3 compiler).
* ``enableload``: Load wave functions from harddisc when needed (otherwise stored in memory).
* ``withwan``: Build with Wannierfunction support(Wannierfunction support requires the Wannier90 library, whose location can be specified with ``â€“withwan=DIR``)
* ``prefix=PREF``: Executables will be installed in directory PREF (`/usr/local` is default).
+* ``enableparallel``: Build parallel version (requires an MPI3 compiler).
+* ``enableload``: Load wavefunctions from harddisc when needed (otherwise stored in memory).
+* ``withwan``: Build with Wannierfunction support (Wannierfunction support requires the Wannier90 library, whose location can be specified with ``withwan=DIR``).
+* ``prefix=PREF``: Executables will be installed in directory "PREF" (/usr/local is default).
.. note:: Enter ``./configure h`` for a full list of configuration options
+Enter ``./configure h`` for a full list of configuration options.
+.. note:: If the configuration fails, the file "config.log" should be consulted for more details about the failure. In most cases, the compiler does not find a required library or an include file, in which case the corresponding directories must be specified with the environment variables LDFLAGS and FCFLAGS, respectively.
diff git a/source/parallelized_version.rst b/source/parallelized_version.rst
new file mode 100644
index 0000000000000000000000000000000000000000..c4a6062a353e611907218a4aec64134e740b2a43
 /dev/null
+++ b/source/parallelized_version.rst
@@ 0,0 +1,89 @@
+.. _parallel:
+
+=====================
+Parallelized version
+=====================
+
+General remarks
+++++++++++++++++
+
+The parallelized version of Spex uses the MPI standard 3.1.
+It can be run on several CPUs on the same node or on several nodes with the command ``mpirun`` or ``mpiexec`` or
+whatever MPI launcher your computer system uses. In principle, there are no restrictions with respect to
+the number of processes. A better performance is expected, though, if the number of processes is not a
+prime number but has a long prime factorization because this gives the code more freedom to distribute
+the work among the processes.
+
+The default parallelization strategy of Spex is conservative in the sense that memory demand and
+load imbalances are minimized. Often, the parallelized run can be sped up substantially by using ``MPIKPT``
+or ``MPIBLK``, see below.
+
+Special MPI keywords
++++++++++++++++++++++
+
+MPISPLIT
+
+(*) The sharedmemory functionality of MPI 3.1 is used for several big arrays,
+which allows the same memory region to be accessed by several MPI processes.
+By default, all processes running on one node share the memory.
+It can be reasonable to change this behavior to, e.g., having processes on the same socket to share memory.
+This is possible with ``MPISPLIT SHRD=NODE`` (default), ``MPISPLIT SHRD=SOCKET`` (only works with ``OpenMPI``), or ``MPISPLIT SHRD=n``,
+where n is the number of the sharedmemory processes.
+Using this option increases the memory consumption but might be advantageous
+in terms of computation time. (There is not much experience yet with this option.)
+
+MPIKPT
+
+In many calculation types (*GW*, Hubbard *U* calculations, COHSEX, ...) there is an outer loop over the kpoint set.
+By default, Spex does not parallelize over this loop because different k points need different computation times depending on their symmetry,
+making the work distribution nontrivial. However, if there are many k points, it is recommendable to additionally
+parallelize over this loop. This can be enabled with the keyword ``MPIKPT``. The kloop parallelization is over nodes, not over processes.
+The computation for each individual k point runs in parallel over the processes on the respective node, in the same way as all processes would without ``MPIKPT``.
+If you only have a single node (or very few nodes) available, you can still use ``MPIKPT`` in conjunction with ``MPISPLIT``, which
+allows processes to be grouped into virtual nodes.
+
+MPIBLK (SENERGY)
+
+Another special parallelization layer is the parallelization over blocks of the selfenergy matrix (or over the diagonal elements).
+This may speed up the calculation if there are many blocks but may also result in work inbalance. (Different blocks need different
+computation times.) Parallelization over blocks is enabled with the keyword ``MPIBLK``. (An optional argument, e.g., ``MPIBLK 5``, can
+be used to finetune the work distribution. It gives the "relative computational overhead" of each block that does not
+scale with the number of bands.
+The default value is 10.) ``MPIBLK`` is enabled automatically for the jobs ``HF``, ``PBE0``, ``SX``, and ``COSX``.
+
+.. listtable:: Examples
+ :widths: 18 100
+
+ *  ``MPIBLK``
+  Enable parallelization over selfenergy blocks (or diagonal elements).
+ *  ``MPIBLK 50``
+  Enable parallelization with assumed large "computational overhead".
+ *  ``MPIBLK 0``
+  Disable parallellization over blocks.
+
+MPISYM (SENERGY)
+
+(*) Using PadÃ© approximants in the evaluation of the *GW* selfenergy (``CONTINUE`` or ``CONTOUR`` with PadÃ© approximant for *W* or
+``FREQINT PADE``)
+might lead to a slight symmetry breaking in the quasiparticle energies, leading to unphysical lifting of degeneracies.
+(This is caused by the fact that Thiele's continuedfraction PadÃ© formula is numerically unstable, especially for a large number
+of imaginary frequencies.)
+Since these errors are usually very small, this is not a big problem. Furthermore, when the full selfenergy matrix is calculated
+(e.g., ``GW FULL``), Spex performs a symmetrization of the selfenergy matrix, which enforces the correct degeneracies again.
+However, for testing purposes, it is possible to enforce the correct symmetries already in the evaluation of the selfenergy by
+using the keyword ``MPISYM``. This requires additional communication among the processes, potentially slowing down
+the calculation due to the necessary blocking synchronization.
+
+RESTART
+
+In the parallelized version, the ``RESTART`` option works in exactly the same way as for the serial version
+(see :numref:`RESTART_GW` and :numref:`RESTART_U`). However, the restart
+data might be written to separate files when ``MPIKPT`` is used. (The underlying reason for this is that binary or HDF5 files can be written
+in parallel, i.e., by all processes at the same time, only if the dataset sizes are known in advance. This is not the case for the restart data.)
+Instead of a single file "spex.cor", Spex writes the files "spex.cor.1", "spex.cor.2",
+et cetera, and a directory "spex.cor.map", which contains, for each k point, *links* to the respective *cor file* that contains the data.
+Furthermore, in addition to "spex.sigc" ("spex.sigx", "spex.wcou", "spex.ccou", "spex.core"),
+the files "spex.sigc.2", "spex.sigc.3", et cetera, might be
+written (and analogously for the other file names). These multiple files should be taken into account, when restart files are transferred.
+Switching between different numbers of processes, different numbers of nodes,
+or between serial and parallel runs should not lead to problems. Spex should be able to always read the correct data.
diff git a/source/spex_and_fleur.rst b/source/spex_and_fleur.rst
new file mode 100644
index 0000000000000000000000000000000000000000..a7c61c2da728e33c14b9362b7ed89f75f6b7d338
 /dev/null
+++ b/source/spex_and_fleur.rst
@@ 0,0 +1,99 @@
+.. _SpexandFleur:
+
+===============
+Spex and Fleur
+===============
+
+Spex needs input data about a reference meanfield system, for example, the KS system of DFT. A selfconsistent meanfield
+calculation must thus precede any Spex calculation. The following three steps prepare the input data for Spex.
+
+:Step 1: Selfconsistent meanfield calculation.
+:Step 2: Generate special kpoint set with Spex.
+:Step 3: Diagonalize meanfield Hamiltonian with the potential from step 1, which produces the eigensolutions at the new k points.
+
+Fleur
+++++++
+
+The Fleur code is an implementation of KSDFT based on the FLAPW method. Consult the `Fleur manual`_ for details.
+
+.. _Fleur manual: http://www.flapw.de/pm/index.php?n=UserDocumentation.FLEUR
+
+We distinguish between an old and a new generation of Fleur versions. In the long run, the old versions will become
+obsolete. However, as of now, the new Fleur versions have not been completely integrated.
+Therefore, both are documented here.
+
+Old Fleur
+==========
+In the Fleur input file there is a flag for writing out the necessary data for a Spex calculation. It is called ``gw`` and is appended to line 23 of the Fleur "inp" file:
+
+``23vchk=F,cdinf=F,pot8=F,gw=n,numbands=N``
+
+*n* can take four values:
+
+* ``gw=0``: No output files are generated (default; the same as omitting ``gw`` and ``numbands`` altogether).
+* ``gw=1``: Some basic parameters are written to several files. Otherwise Fleur runs as usual.
+* ``gw=2``: All output files necessary for Spex are created and Fleur stops after one iteration (without updating the potential). In this case you should also set ``pot8=T`` and the maximal number of bands ``N``. The latter overrides the default ``neigd`` value in "fl7para".
+* ``gw=3``: Selfconsistent cycle for QSGW. Same as ``gw=1`` but adds :math:`\Sigma^{\mathrm{xc,QSGW}}(\mathbf{r},\mathbf{r}')v^{\mathrm{xc}}(\mathbf{r})` to the Hamiltonian in each iteration.
+
+The output files are
+
+* "gwa": Basic parameters including the atomic numbers and positions in the unit cell, lattice parameters and basis vectors, FLAPW l cutoff, and localorbital parameters,
+* "LATTC": FLAPW G cutoff,
+* "radfun": radial basis functions,
+* "ecore": coreelectron functions,
+* "eig": k points, wavefunction coeffients (interstitial), energies,
+* "abcoeff": wavefunction coefficients (muffin tins),
+* "KS.hdf": alternative to "eig" and "abcoeff" (now the default),
+* "vxc": expectation values of the xc potential (diagonal elements, obsolete),
+* "vxcfull": matrix of the xc potential,
+* "qsgw": matrix of the QSGW selfenergy (only for QSGW calculations).
+
+Furthermore, Spex needs the Fleur file "sym.out".
+
+New Fleur
+==========
+Spex has to be compiled with ``DFleurMax`` to make it compatible with the new Fleur "MaX" versions.
+
+...(undocumented)...
+
+Remarks about MT basis
++++++++++++++++++++++++
+
+Conduction states
+==================
+As in a *GW* calculation the conduction states enter both the Green function and the screened interaction,
+they should be accurately described already in the KS system which is our starting point for the perturbation treatment.
+The LAPW basis, however, only guarantees well described occupied states. In order to obtain wellconverged *GW* results,
+one must make the basis set more flexible, especially in the MT regions, by introducing local orbitals either at high energy
+parameters or as higherorder energy derivatives [Phys. Rev. B 74, 045104 (2006), Comput. Phys. Commun., 184, 2670 (2013)].
+
+Semicore states
+=================
+Highlying semicore states can appear as badly described "ghost bands" in the valence band region leading to a wrong DFT ground state.
+A possible solution is to use local orbitals and extend the energy window, such that the semicore states are treated as valence states.
+On the other hand if the ghost bands are above the Fermi energy, they usually pose no problem in a DFT calculation.
+In a *GW* calculation, however, the conduction states are also relevant, and one must make sure that there are no ghost bands at all
+in the energy spectrum. Therefore, it might be necessary to treat more (and deeper) semicore states as valence states than are needed
+in the DFT calculation.
+As an indication for such a case the overlap of core states with basis and wavefunctions calculated in the routine "checkinput"
+should be checked. As an example, here is the output for Strontium Titanate with a Titanium 3s ghost band:
+
+.. codeblock:: bash
+
+ Overlap
+ Atom type 1
+ u(s) udot(s)
+ 1s 0.001431 0.000644
+ 2s 0.000830 0.000414
+ 3s 0.287789 0.878671
+ u(p) udot(p) ulo(p) ...
+ 2p 0.000665 0.000112 0.000494
+ ...
+ Maximum overlap at (band/kpoint)
+ Atom type 1
+ 1s 0.000611 (031/ 001)
+ 2s 0.000343 (031/ 001)
+ 3s 0.664265 (047/ 001)
+ 2p 0.000340 (003/ 036) 0.000481 (003/ 042) 0.000340 (003/ 036)
+ ...
+
diff git a/source/spex_faq.rst b/source/spex_faq.rst
index 5e4a607c9d56cb03c62bef1f076284463d4d7262..a1ef5bae830261eb2cfe48d276735855039fe602 100644
 a/source/spex_faq.rst
+++ b/source/spex_faq.rst
@@ 47,7 +47,7 @@ Usage
* In the calculation of the susceptibility groups of degenerate states are split into different packets. This can only be avoided with a larger ``MEM`` if possible.
* The tetrahedron method breaks the symmetry. You can use ``GAUSS`` to test this.
**My ``GW`` calculations do not converge wrt the kpoint set.**
+**My GW calculations do not converge wrt the kpoint set.**
.. highlights:: In systems with small band gaps (e.g. GaAs) the default treatment of the point k=0 might lead to bad kpoint convergence. In this case try the option ``NOZERO`` which guarantees a smooth (but slower) convergence.
@@ 65,27 +65,71 @@ Usage
Error Messages

.. highlights:: This is not a complete list of errors. In many error messages a possible solution is already indicated. These are not listed. The program might also stop at places where this is not supposed to happen. Then the message ends with ``...(bug?)``. In this case you should contact the programmers.
+.. highlights:: This is not a complete list of errors. In many error messages a possible solution is already indicated. These are not listed.
+ The program might also stop at places where this is not supposed to happen. Then, the error message starts with ``SPEXBUG`` (see :numref:`spex.out`).
+ In this case you should contact the programmers. Error messages have the form ``SPEXERROR (source.f:0123) Error message``, where
+ "source.f" is the name of the source file and 0123 is the respective line where the error occurred. We shorten this to ``(source.f) Error message`` in this list.
**error while loading shared libraries: lib...: cannot open shared object file: No such file or directory**
.. highlights:: The library "lib..." is not found in standard directories. In this case, the location of the library has to be defined in the environment variable ``LD_LIBRARY_PATH``.
+.. highlights:: The library "lib..." is not found in standard directories. In this case, the location of the library has to be defined in the
+ environment variable ``LD_LIBRARY_PATH`` of the shell. (Consult your Unix or Linux manual.)
**correlation: Solution of quasiparticle equation did not converge after 1000 iterations.**
+**(read_write.f) kpoint sets inconsistent (check ...)**
.. highlights:: You probably have poles very close to the real frequency axis (normally with a small weight). Usually a slight change of parameters solves the problem.
+.. highlights:: This error usually occurs if the kpoint set defined in "spex.inp" (keyword ``BZ``) conflicts with the one the DFT code has used
+ for writing the input data. Make sure that the kpoint sets are identical by generating the kpoint set with Spex and running the oneshot
+ DFT calculation in the proper order (steps 2 and 3 in :numref:`getting_started`). (Only for older versions:)
+ The error can also be caused by the internal structure of the directaccess binary file
+ containing the wavefunction data. Make sure that Spex and the DFT code are compiled in such a way that these files are treated identically.
**diagonalize_coulomb: Negative/Residual eigenvalue.**
+**(coulombmatrix.f) Negative eigenvalue.**
.. highlights:: Your Coulomb matrix seems to be underconverged. Try to increase the parameter ``LEXP`` or adjust the scaling with ``SCALE``. If this does not solve the problem, the reason might be that the LAPACK routine fails to solve the general eigenvalue problem because of a nearly singular overlap matrix. Then you might try the option ``CUTZERO`` which removes linear dependencies from the overlap matrix and reattempts the diagonalization.
+.. highlights:: Your Coulomb matrix seems to be underconverged. Try to increase the parameter ``LEXP``.
+ If this does not solve the problem, the reason might be that the LAPACK routine fails to solve the general eigenvalue problem
+ because of a nearly singular overlap matrix. Then you might try the option ``CUTZERO`` which removes linear dependencies from
+ the overlap matrix and reattempts the diagonalization. The accuracy of the Coulomb matrix can be tested with the (undocumented) keywords
+ ``CHKCOUL``, ``TSTCOUL``, and ``STEPRAD`` (all in section ``COULOMB``).
**intgrf: Negative exponent x in extrapolation a+c*r**x.**
+.. **intgrf: Negative exponent x in extrapolation a+c*r**x.**
+ The first (nonzero) point of your radial MT mesh might be too far away from zero. Modify the mesh accordingly. You can also try to increase ``TOL`` in section ``MBASIS``.
.. highlights:: The first (nonzero) point of your radial MT mesh might be too far away from zero. Modify the mesh accordingly. You can also try to increase ``TOL`` in section ``MBASIS``.
+**(global.f) argument not an element of kpoint set.**
**modulo1: argument not an element of kpoint set.**
+.. highlights:: This message is caused by an inconsistency in the kpoint set. Either the kpoint set defined by ``BZ``
+ and the one stored in the DFT output file are different or the k points defined by ``KPT`` are not elements of the set.
+ The latter happens, for example, if you use a definition with square brackets (such as ``KPT A=[1,1,0]``), but
+ the lattice constant is not properly defined (such definitions are interpreted with respect to the lattice constant).
+ In any case, it is recommended to use internal coordinates (such as ``KPT A=(0,0,1)``). See :numref:`KPT`.
.. highlights:: This message is caused by an error in the kpoint set. Either the kpoint sets defined by ``BZ`` and in the DFT output file are different or the k points defined by ``KPT`` are not elements of the set. The latter happens, if you use ``KPT ?=[?,?,?]`` and the lattice constant is not properly defined since ``[?,?,?]`` is interpreted with respect to the lattice constant.
+**(getinput.f) Defective symmetry transformation matrix.**
+
+.. highlights:: A calculated symmetry transformation matrix turned out to be nonunitary. This is usually caused by insufficient precision of
+ the lattice basis vectors provided by the DFT code. The obvious solution is to make sure that the lattice basis vectors are defined with
+ sufficient precision in the DFT code. In particular, they have to conform to the crystal symmetries with a precision of at least ten
+ significant digits.
+
+**(getinput.f) Symmetry operation ... rotates spin quantization axis.**
+
+.. highlights:: Similarly to the previous one, this error is caused by insufficient precision: Spex expects the spinquantization axis to
+ be invariant with respect to the crystal symmetries in the case of noncollinear magnetism. Obvious solution: Make sure that the spinquantization axis
+ is defined with sufficient precision in the DFT code.
+
+ **(irreps.f) Broken unitarity ...**
+ **(irreps.f) Deviation in characters; too large irrep threshold increase: ...**
+
+.. highlights:: Both errors are caused by eigenstates (read from the DFT code) that are symmetry broken, i.e., they do not reflect sufficiently
+ the crystal symmetries. This could be a problem of the DFT code. At the moment, there is no general solution to this problem. The second
+ error can be circumvented, though, by setting the parameter ``irrep_threshold`` in "irreps.f" to 0, which might lead to slower calculations.
+
+.. Found 2 states in subspace of eigenvalue 2.12487 while expecting 4.
+ Currently no solution to this error.
+
+**(quasiparticle.f) Solution of quasiparticle equation did not converge after 5000 iterations.**
+
+.. highlights:: If you use analytic continuation (``CONTINUE``), you probably have poles very close to the real frequency axis
+ (normally with a small weight). Usually a slight change of parameters solves the problem. You can also try to use the keyword ``SMOOTH`` (:numref:`SMOOTH`).
+ In the case of contour integration (``CONTOUR``), strong spectral structure in the selfenergy can cause this problem. For example in metallic systems, a small
+ Drude frequency can give rise to a deltalike peak in the selfenergy. A possible solution in this case is to use ``PLASMA METAL`` (:numref:`PLASMA`).
To be continued ...
diff git a/source/tutorials.rst b/source/tutorials.rst
index 97c6f3c3cac64be74ba33e6ab0293a5d97a72ef5..dcc7056d7a459c8cd5be08796b2b0d6f36051276 100644
 a/source/tutorials.rst
+++ b/source/tutorials.rst
@@ 5,34 +5,63 @@ Tutorials
===================
In this section, we provide detailed tutorials for the main features of the Spex code.
GW calculations
===============
The Spex code was started as a pure GW code. The GW approach remains the main computational method implemented in the code.
+*GW* calculations
+==================
+The Spex code was started as a pure *GW* code. The *GW* approach remains the main computational method implemented in the code.
Oneshot approach

A simple input file for a oneshot GW calculation for bulk silicon was already presented in :ref:`getting_started`. The following, equivalent input file is even more minimalistic.
+Oneshot *GW* approach
+
+A simple input file for a oneshot *GW* calculation for bulk silicon was already presented in :numref:`getting_started`.
+The following, equivalent input file is even more minimalistic.
.. codeblock:: bash
BZ 4 4 4
JOB GW 1:(1,2,5) 7:(1,3,5) 3:(13,5)
It only contains the two essential keywords that any Spex input file must contain: ``BZ`` and ``JOB``, specifying, respectively, the k mesh for the Brillouinzone sampling and the requested type of calculation. Here, we calculate quasiparticle corrections for the bands (1,2,5) of the first, (1,3,5) of the seventh, and (1,2,3,5) of the third k point. The band indices are defined according to the ordering in the meanfield input data. In the example, we have left out some band indices because of energy degeneracy. For example, bands 3 and 4 (of the first k point) are degenerate with band 2.
+It only contains the two essential keywords that any Spex input file must contain: ``BZ`` and ``JOB``, specifying, respectively,
+the k mesh for the Brillouinzone sampling and the requested type of calculation. Here, we calculate quasiparticle corrections
+for the bands (1,2,5) of the first, (1,3,5) of the seventh, and (1,2,3,5) of the third k point. The band indices are defined
+according to the ordering in the meanfield input data. In the example, we have left out some band indices because of energy
+degeneracy. For example, bands 3 and 4 (of the first k point) are degenerate with band 2.
In the case of magnetic systems, the quasiparticle energies for both spins are calculated by default. Alternatively, one can choose the spin index by using u (spin up) or d (spin down), e.g., ``JOB GW 1:u(1,2,5) 1:d(1,2,5)`` is identical to ``JOB GW 1:(1,2,5)``.
+In the case of a spinpolarized systems, the quasiparticle energies for both spins are calculated by default.
+Alternatively, one can choose the spin index by using u (spin up) or d (spin down), e.g., ``JOB GW 1:u(1,2,5) 1:d(1,2,5)``
+is identical to ``JOB GW 1:(1,2,5)``.
.. _kpt:
KPT
+KPT

The kpoints follow an internal ordering: first the irreducible parent k points, then the remaining equivalent k points. A list of the irreducible kvectors is written to the output file. In the present example, there are eight irreducible k points and 64 k points in total. All 64 kpoint indices can be used in the job definition. We have chosen the ones from the irreducible set. They correspond to the :math:`{\Gamma}` (index 1), X (index 7), and L point (index 3). Clearly, for a different kpoint set, e.g., 8x8x8, the kpoint indices would change (except for the index 1, which always corresponds to the :math:`{\Gamma}` point). Therefore, there is the possibility of defining kpoint labels with
+The kpoints follow an internal ordering: first the irreducible (parent) k points, then the remaining equivalent k points.
+A list of the irreducible kvectors is written to the output file. In the present example, there are eight irreducible
+k points and 64 k points in total. All 64 kpoint indices can be used in the job definition. We have chosen the ones
+from the irreducible set. They correspond to the :math:`{\Gamma}` (index 1), X (index 7), and L point (index 3).
+Clearly, for a different kpoint set, e.g., :math:`8\times 8\times 8`, the kpoint indices would change (except for the index 1, which
+always corresponds to the :math:`{\Gamma}` point). Therefore, there is the possibility of defining kpoint labels by
``KPT X=1/2*(0,1,1) L=1/2*(0,0,1)``
or
+or
``KPT X=[1,0,0] L=1/2*[1,1,1]``
The labels must be single (upper or lowercase) characters. The two definitions are equivalent. The first gives the k vectors in internal coordinates, i.e., in the basis of the reciprocal basis vector, the second in cartesian coordinates in units of :math:`{2\pi/a}` with the lattice constant a . In the case of the fcc lattice of silicon, the lattice basis vectors are :math:`{a_1=a(011)/2}` , :math:`a_2=a(101)/2` , and :math:`a_3=a(110)/2` . The reciprocal basis vectors are thus :math:`{b_1=2\pi/a(111)}` et cetera. For the X point, e.g., one then gets :math:`{(a_2+a_3)/2=2\pi/a(100)}` . In order for the cartesian coordinates (square brackets) to be interpreted correctly, Spex must obviously know the lattice constant a. (In Fleur, the lattice constant is taken to be the global scaling factor for the lattice vectors.) The sodefined k points must be elements of the k mesh defined by ``BZ``. We will later discuss how points outside the k mesh can be considered using the special label +. With the k labels, the above job definition can be written as ``JOB GW 1:(1,2,5) X:(1,3,5) L:(13,5)`` as in the input file described in :ref:`getting_started` There are two more special kpoint labels: ``IBZ`` and ``BZ`` (e.g., ``JOB GW IBZ:(1,2,5)``) stand for all k points in the irreducible and the full k set, respectively. (The label ``BZ`` is included for completeness but is not needed in practice.) The ``IBZ`` label is helpful when a selfconsistent GW calculation is to be performed, which requires the selfenergy to be calculated for the whole irreducible Brillouin zone.
+The labels must be single (upper or lowercase) characters. The two definitions are equivalent.
+The first gives the k vectors in internal coordinates, i.e., in the basis of the reciprocal basis vector,
+the second in cartesian coordinates in units of :math:`{2\pi/a}` with the lattice constant :math:`a`. In the case
+of the fcc lattice of silicon, the lattice basis vectors are :math:`{a_1=a(0,1,1)/2}` , :math:`a_2=a(1,0,1)/2` ,
+and :math:`a_3=a(1,1,0)/2` . The reciprocal basis vectors are thus :math:`{b_1=2\pi/a(1,1,1)}` et cetera.
+For the X point, e.g., one then gets :math:`{(a_2+a_3)/2=2\pi/a(1,0,0)}` . In order for the cartesian
+coordinates (square brackets) to be interpreted correctly, Spex obviously needs to know the lattice constant
+:math:`a`.
+(In Fleur, the lattice constant is taken to be the global scaling factor for the lattice vectors.)
+The sodefined k points must be elements of the k mesh defined by ``BZ``. We will later discuss how
+points outside the k mesh can be considered using the special label ``+``. With the k labels, the above
+job definition can be written as ``JOB GW 1:(1,2,5) X:(1,3,5) L:(13,5)`` as in the input file described
+in :numref:`getting_started`. There are two more special kpoint labels: ``IBZ`` and ``BZ`` (e.g., ``JOB GW IBZ:(1,2,5)``)
+stand for all k points in the irreducible and the full k set, respectively. (The label ``BZ`` is included
+for completeness but is not needed in practice.) The ``IBZ`` label is helpful when a Wannier interpolation
+of *GW* quasiparticle energies or
+a selfconsistent *GW* calculation is to be performed, which require the selfenergy
+to be calculated for the whole irreducible Brillouin zone.
The quasiparticle energies are written to standard output in tabular form::
@@ 47,279 +76,504 @@ The quasiparticle energies are written to standard output in tabular form::
The columns are
* ``Bd`` band index,
* ``vxc`` expectation value of the exchangecorrelation potential :math:`v^{xc}`, e.g., 12.15897 eV,
* ``sigmax`` expectation value of the exchange selfenergy :math:`\sum^x`, e.g., 19.20415 eV,
* ``sigmac`` expectation value of the correlation selfenergy :math:`\sum^c`, e.g., (6.60649+1.05993i) eV (note that the selfenergy is complex),
* ``Z`` renormalization factor Z, e.g., (0.651460.10556i) eV,
* ``KS`` meanfield energy eigenvalue (e.g., from the previous KSDFT calculation), e.g., 6.19011 eV,
* ``HF`` HartreeFock value (oneshot), e.g., 13.23529 eV,
* ``GW`` GW quasiparticle values, e.g., (6.36401+0.73681i ) eV (linearized solution) and (6.36806+0.72704i ) eV (direct solution).
+* ``vxc`` expectation value of the exchangecorrelation potential :math:`v^\mathrm{xc}`, e.g., 12.15897 eV,
+* ``sigmax`` expectation value of the exchange selfenergy :math:`\Sigma^x`, e.g., 19.20415 eV,
+* ``sigmac`` expectation value of the correlation selfenergy :math:`\Sigma^c`, e.g., (6.60649+1.05993\ *i*\ ) eV (note that the selfenergy is complex),
+* ``Z`` renormalization factor :math:`Z`, e.g., (0.651460.10556\ *i*\ ) eV,
+* ``KS`` meanfield energy eigenvalue (e.g., from the previous KSDFT calculation), e.g., 6.19011 eV,
+* ``HF`` HartreeFock value (oneshot), e.g., 13.23529 eV,
+* ``GW`` *GW* quasiparticle values, e.g., (6.36401+0.73681\ *i*\ ) eV (linearized solution) and (6.36806+0.72704\ *i*\ ) eV (direct solution).
The two GW values for the quasiparticle energy follow two common methods to approximately solve the quasiparticle equation
+The two *GW* values for the quasiparticle energy follow two common methods to approximately solve the quasiparticle equation
.. math::
 \displaystyle\left\{\frac{\nabla}{2}+v^\mathrm{ext}(\mathbf{r})+v^\mathrm{H}(\mathbf{r})\right\}\psi_{\mathbf{k}n}(\mathbf{r})+\int \Sigma^\mathrm{xc}(\mathbf{r},\mathbf{r}';E_{\mathbf{k}n})\psi_{\mathbf{k}n}(\mathbf{r}')d^3 r=E_{\mathbf{k}n}\psi_{\mathbf{k}n}(\mathbf{r})
+ \left\{\frac{\nabla}{2}+v^\mathrm{ext}(\mathbf{r})+v^\mathrm{H}(\mathbf{r})\right\}\psi_{\mathbf{k}n}(\mathbf{r})+\int \Sigma^\mathrm{xc}(\mathbf{r},\mathbf{r}';E_{\mathbf{k}n})\psi_{\mathbf{k}n}(\mathbf{r}')d^3 r=E_{\mathbf{k}n}\psi_{\mathbf{k}n}(\mathbf{r})
:label: qpeq
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`` selfenergy, 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 meanfield eigenvalue
+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* selfenergy, 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 meanfield eigenvalue
.. math::
 \displaystyle E_{\mathbf{k}n}=\epsilon_{\mathbf{k}n}+\langle\phi_{\mathbf{k}n}\Sigma^\mathrm{xc}(E_{\mathbf{k}n})v^\mathrm{xc}\phi_{\mathbf{k}n}\rangle\approx\epsilon_{\mathbf{k}n}+Z_{\mathbf{k}n}\langle\phi_{\mathbf{k}n}\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})v^\mathrm{xc}\phi_{\mathbf{k}n}\rangle
+ E_{\mathbf{k}n}=\epsilon_{\mathbf{k}n}+\langle\phi_{\mathbf{k}n}\Sigma^\mathrm{xc}(E_{\mathbf{k}n})v^\mathrm{xc}\phi_{\mathbf{k}n}\rangle\approx\epsilon_{\mathbf{k}n}+Z_{\mathbf{k}n}\langle\phi_{\mathbf{k}n}\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})v^\mathrm{xc}\phi_{\mathbf{k}n}\rangle
:label: qppert
with the singleparticle wavefunction :math:`{\phi_{\mathbf{k}n}}` and the frequencyindependent potential :math:`{v^{\mathrm{xc}}}`, which in the case of a KS solution would correspond to the local exchangecorrelation potential; the nonlocal HartreeFock exchange potential and the *hermitianized* selfenergy 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 righthand side correspond to the "linearized" and "direct" (iterative) solutions given in the output. The direct solution takes into account the nonlinearity 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:(110)``. Except the latter (``GT``), all of these methods are meanfield approaches, so one only gets one singleparticle energy (instead of a ''linearized'' and a ''direct'' solution) for each band.

.. _spectral:

SPECTRAL (SENERGY)
+with the singleparticle wavefunction :math:`{\phi_{\mathbf{k}n}}` and the frequencyindependent potential
+:math:`{v^{\mathrm{xc}}}`, which in the case of a KS solution would correspond to the
+local exchangecorrelation potential; the nonlocal HartreeFock exchange potential and
+the *hermitianized* selfenergy 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 righthand side correspond to the "linearized" and "direct" (iterative) solutions given in the output. The direct solution takes into account the nonlinearity 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:(110)``.
+(The keyword ``FULL`` is explained in :numref:`beyondpert`.)
+Except the latter (*GT*), all of these methods are meanfield approaches,
+so one only gets one singleparticle energy (instead of a ''linearized'' and a ''direct'' solution) for each band.
+
+.. _SPECTRAL:
+
+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
+(*) 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
.. math::
 \displaystyle G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)v^{\mathrm{xc}}]G(\omega)

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

.. math::
 A(\mathbf{k},\omega)=\pi^{1}\,\text{sgn}(\omega\epsilon_\mathrm{F})\,\,\text{tr}\left\{\text{Im}[\omega IH^\mathrm{KS}(\mathbf{k})\Sigma^{\mathrm{xc}}(\mathbf{k},\omega)]^{1}\right\}\,,

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

``SPECTRAL``

in the section ``SENERGY`` of the input file. This option does not require the ``FULL`` keyword. Optionally, one can define a frequency mesh by
the keyword ``SPECTRAL`` in the section ``SENERGY``, in the examples below from 10 eV to 1 eV in steps of 0.01 eV. If there is no explicit mesh, Spex chooses one automatically. The spectral function is then written to the file "spectral", one block of data for each k point given in the job definition. There is another optional parameter given at the end of the line (second example below), which can be used to bound the imaginary part of the selfenergy and, thus, the quasiparticle peak widths from below by this value (unset if not given). This can be helpful in the case of very sharp quasiparticle peaks that are otherwise hard to catch with the frequency mesh.


++++
 Examples  ``SPECTRAL {10eV:1eV,0.01eV}``   Write spectral function :math:`{\text{Im}G(\omega)}` on the 
    specified frequency mesh to the file "spectral". 
++++
  ``SPECTRAL {10eV:1eV,0.01eV} 0.002``   Bound imaginary part from below by 0.02, preventing 
    peak widths to become too small to be plotted. 
++++
+ G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)v^{\mathrm{xc}}]G(\omega)
Alhough ``GW`` calculations can be readily performed with the default settings, the user should be familiar to some degree with the details of the computation and how he/she can influence each step of the program run. Also note that the default settings might work well for one physical system but be unsuitable for another.
+which links the interacting Green function :math:`G` to the noninteracting KS one :math:`{G_0}` and which,
+in principle, requires the selfenergy to be known on the complete :math:`{\omega}` axis.
+The spectral function measured in photoelectron spectroscopy is directly related to the Green function by
The ``GW`` selfenergy

.. math:: \displaystyle\Sigma^{\mathrm{xc}}(\mathbf{r},\mathbf{r}';\omega)=\frac{i}{2\pi}\int_{\infty}^\infty G(\mathbf{r},\mathbf{r}';\omega+\omega')\,W(\mathbf{r},\mathbf{r}';\omega')e^{i\eta\omega'}\,d\omega'\,.
 :label: selfene

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 manyelectron system into an effective dynamical potential, obtained from the dielectric function :math:`{\varepsilon}` through

.. math::
 \displaystyle W(\mathbf{r},\mathbf{r}';\omega)=\int \varepsilon^{1}(\mathbf{r},\mathbf{r}'';\omega) v(\mathbf{r}'',\mathbf{r}')\,d^3 r''\,.

This integral equation turns into a matrix equation

:math:`{\displaystyle W_{\mu\nu}(\mathbf{k},\omega)=\sum_\eta \varepsilon^{1}_{\mu\eta}(\mathbf{k},\omega)v_{\eta\nu}(\mathbf{k})}`

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)=1P(\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 selfenergy can be written as the sum of two terms, the first of which is the exact nonlocal exchange potential of HartreeFock theory, the remainder can be interpreted as a correlation selfenergy and has the mathematical form of the equation :eq:`selfene` 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.)
+.. math::
+ A(\mathbf{k},\omega)=\pi^{1}\,\text{sgn}(\omega\epsilon_\mathrm{F})\,\,\text{tr}\left\{\text{Im}[\omega IH^\mathrm{KS}(\mathbf{k})\Sigma^{\mathrm{xc}}(\mathbf{k},\omega)]^{1}\right\}\,,
+
+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 keyword ``SPECTRAL``
+in the section ``SENERGY``.
+Optionally, one can define a frequency mesh by
+the keyword ``SPECTRAL`` in the section ``SENERGY``, in the examples below from 10 eV to 1 eV in steps of 0.01 eV.
+If there is no explicit mesh, Spex chooses one automatically.
+The spectral function is then written to the file "spectral",
+one block of data for each k point given in the job definition.
+There is another optional parameter given at the end of the line (third example below),
+which can be used to bound the imaginary part of the selfenergy and, thus,
+the quasiparticle peak widths from below by this value (unset if not given).
+This can be helpful in the case of very sharp quasiparticle peaks that are otherwise hard to catch with the frequency mesh.
+The ``SPECTRAL`` keyword can be used in conjunction with bandstructure calculations, see :numref:`GWbandstr` for details.
+
+.. listtable:: Examples
+ :widths: 95 100
+
+ *  ``SPECTRAL``
+  Write spectral function :math:`{\text{Im}G(\omega)}` to the file "spectral"; frequency mesh automatically defined.
+ *  ``SPECTRAL {10eV:1eV,0.01eV}``
+  Write spectral function on the specified frequency mesh.
+ *  ``SPECTRAL {10eV:1eV,0.01eV} 0.002eV``
+  Bound imaginary part from below by 0.002 eV, preventing peak widths to become too small to be plotted.
+
+*GW* basics
+
+Alhough *GW* calculations can be readily performed with the default settings,
+the user should be familiar to some degree with the details of the computation
+and how he/she can influence each step of the program run.
+Also note that the default settings might work well for one physical system but be unsuitable for another.
+To lay the basis for the subsequent sections, we briefly recapitulate the basics of the *GW* method.
+
+The *GW* selfenergy
+
+.. math:: \Sigma^{\mathrm{xc}}(\mathbf{r},\mathbf{r}';\omega)=\frac{i}{2\pi}\int_{\infty}^\infty G(\mathbf{r},\mathbf{r}';\omega+\omega')\,W(\mathbf{r},\mathbf{r}';\omega')e^{i\eta\omega'}\,d\omega'\,.
+ :label: selfene
+
+can be understood as a scattering potential that contains the exact exchange potential and correlation effects
+through the inclusion of :math:`W`, the dynamically screened interaction, which incorporates the screening of
+the manyelectron system into an effective dynamical potential, obtained from the dielectric function :math:`\varepsilon` through
+
+.. math:: W(\mathbf{r},\mathbf{r}';\omega)=\int \varepsilon^{1}(\mathbf{r},\mathbf{r}'';\omega) v(\mathbf{r}'',\mathbf{r}')\,d^3 r''\,.
+
+This integral equation turns into a matrix equation
+
+.. math:: W_{\mu\nu}(\mathbf{k},\omega)=\sum_\eta \varepsilon^{1}_{\mu\eta}(\mathbf{k},\omega)v_{\eta\nu}(\mathbf{k})
+
+if the quantities :math:`{W}`, :math:`{\varepsilon}`, and :math:`{v}`
+are represented in the mixed product basis (:numref:`mpb`), 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)=1P(\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 :numref:`polar`.
+
+The selfenergy can be written as the sum of two terms,
+the first of which is the exact nonlocal exchange potential of HartreeFock theory,
+the remainder can be interpreted as a correlation selfenergy and has the mathematical
+form of Eq. :eq:`selfene` 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 selfenergy as a function of frequency.
+The default method is *analytic continuation*, in which the screened interaction and the selfenergy
+are evaluated on a mesh of purely imaginary frequencies.
+The selfenergy 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 realsymmetric) matrix if :math:`{\omega}` is purely imaginary. Second, :math:`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 selfenergy 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 in the badly controlled extrapolation of the PadÃ© 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 selfenergy directly for selected
+frequencies on the real axis. This method requires the screened interaction :math:`W` to be calculated in a small
+range on the real frequency axis, another contribution comes from an integration along the imaginary frequency axis.
+As a third alternative, there is a *hybrid* method that is nearly as accurate as the "explicit" contour integration
+method, but it is faster and avoids the calculation of :math:`W` for real frequencies altogether.
There are several ways to represent the selfenergy as a function of frequency. The default method is ``analytic continuation``, in which the screened interaction and the selfenergy are evaluated on a mesh of purely imaginary frequencies. The selfenergy 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 realsymmetric) 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 selfenergy 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 selfenergy 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=\{(N1)/[0.9(n1)]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/(n1)1\}^{1}}`. The latter definition is rarely used. One can also employ a userdefined 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 
++++
  ``MESH 12+2 15.0``   Add points at low frequencies, two (one) between :math:`{\omega_1}` and :math:`{\omega_2}` 
    (:math:`{\omega_2}` and :math:`{\omega_3}` ) 
++++
  ``MESH 12 15.0``  Use alternative mesh definition. 
++++
  ``MESH mesh.dat``  Read frequency mesh from file "mesh.dat" 
++++

+All methods to calculate the selfenergy 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=\{(N1)/[0.9(n1)]1\}^{1}}`, :math:`{n=1,2,...}`.
+It is dense 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/(n1)1\}^{1}}`. The latter definition is rarely used. One can also employ a userdefined mesh provided
+in a file (one frequency per line, comments ``#...`` are allowed).
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``MESH 12 15.0``
+  Use a mesh containing twelve frequencies for the range :math:`[0,i15]` htr.
+ *  ``MESH 12+2 15.0``
+  Add points at low frequencies, two (one) between :math:`{\omega_1}` and :math:`{\omega_2}` (:math:`{\omega_2}` and :math:`{\omega_3}`).
+ *  ``MESH 12 15.0``
+  Use alternative mesh definition.
+ *  ``MESH mesh.dat``
+  Read frequency mesh from file "mesh.dat".
+
CONTINUE (SENERGY)

The keyword ``CONTINUE`` in the section ``SENERGY`` chooses the analytic continuation method. It can have optional parameters. If given without parameters (or if not specified at all), Pade approximants are used. An integer number, such as ``CONTINUE 2``, lets Spex perform a fit to an npole function (here, n=2) with the Newton method. The latter approach is somewhat obsolete by now and recommended only for test calculations. The argument can have additional flags +c*. Any combination is possible. They are explained in the examples.

++++
 Examples  ``CONTINUE``  Use Pade approximants (Default). 
++++
  ``CONTINUE 2``  Fit to the twopole function :math:`{a_1/(\omegab_1)+a_2/(\omegab_2)}` 
++++
  ``CONTINUE 2+``  Include a constant in the fit function :math:`{a_1/(\omegab_1)+a_2/(\omegab_2)+c}` 
++++
  ``CONTINUE 2c``   Take constraints (continuity of value and gradient at :math:`{\omega=0}`) 
    into account when fitting 
++++
  ``CONTINUE 2*``   Allow parameters :math:`{b_i}` with positive imaginary parts 
    (should be negative) to contribute. (Default with Pade method.) 
++++

The second method is ``contour integration``, in which the frequency integration is performed explicitly, however not along the real frequency axis but on a deformed integration contour that avoids the real frequency axis as well as possible. This integration contour starts from :math:`{\infty}`, describes an infinite quarter circle to :math:`{i\infty}`, then runs along the imaginary frequency axis to :math:`{i\infty}`, and finishes, again with an infinite quarter circle, to :math:`{\infty}`. The two quarter circles do not contribute to the integral (because the integrand behaves as :math:`{\propto \omega^{2}}`). Furthermore, depending on the frequency argument :math:`{\omega}` of the selfenergy, one has to add a few residues coming from the poles of the Green function in the interval [:math:`{0,\omega\epsilon_\mathrm{F}}`] if :math:`{\omega>\epsilon_\mathrm{F}}` and [:math:`{\epsilon_\mathrm{F}\omega,0}`] otherwise, which requires the correlation screened interaction :math:`{W^\mathrm{c}(\omega)}` to be evaluated on this interval of the real axis. As a consequence, the calculations are more demanding in terms of computational complexity and cost (time and memory). Also, contour integration requires additional input parameters and is therefore somewhat more difficult to apply. However, the results are more accurate. In particular, they are not affected by the "illdefinedness" of the analytic continuation.
+The keyword ``CONTINUE`` (section ``SENERGY``) chooses the analytic continuation method.
+It can have optional parameters. If given without parameters (or if not specified at all),
+PadÃ© approximants are used. An integer number, such as ``CONTINUE 2``,
+lets Spex perform a fit to an *n*\ pole function (here, *n*\ =2) with the Newton method.
+The latter approach is somewhat obsolete by now and recommended only for test calculations.
+The argument can have additional flags ``+``, ``c``, and ``*``. Any combination is possible. They are explained in the examples.
+
+.. listtable:: Examples
+ :widths: 22 100
+
+ *  ``CONTINUE``
+  Use PadÃ© approximants (Default).
+ *  ``CONTINUE 2``
+  Fit to the twopole function :math:`{a_1/(\omegab_1)+a_2/(\omegab_2)}`.
+ *  ``CONTINUE 2+``
+  Include a constant in the fit function :math:`{a_1/(\omegab_1)+a_2/(\omegab_2)+c}`.
+ *  ``CONTINUE 2c``
+  Take constraints (continuity of value and gradient at :math:`{\omega=0}`) into account when fitting.
+ *  ``CONTINUE 2*``
+  Allow parameters :math:`{b_i}` with positive imaginary parts (should be negative) to contribute. (Default with PadÃ© method.)
CONTOUR (SENERGY)

The corresponding keyword is called ``CONTOUR`` and belongs to the section ``SENERGY``. Obviously, ``CONTOUR`` and ``CONTINUE`` must not be given at the same time. The keyword ``CONTOUR`` expects two arguments. The first defines the frequencies :math:`{\omega}`, for which the selfenergy :math:`{\Sigma^\mathrm{xc}(\omega)}` is to be evaluated. At least, two frequencies are needed to approximate the selfenergy as a linear function in :math:`{\omega}` and, thus, to calculate the ''linearized'' solution of the quasiparticle equation (see above). For this, a single value suffices (example ``0.01``), with which the selfenergy for a state :math:`{\mathbf{k}n}` is evaluated at two frequencies (:math:`{\epsilon_{\mathbf{k}n}0.01}` and :math:`{\epsilon_{\mathbf{k}n}+0.01}`). The more accurate ''direct'' solution is only available if we specify a range of frequencies for the selfenergy instead of a single number. This is possible by an argument such as :math:`{0.1:0.15,0.01}`. Here, the range of values is relative to :math:`{\epsilon_{\mathbf{k}n}}`. Note that the range is flipped (to :math:`{0.15:0.1,0.01}` in the example) for occupied states :math:`{\mathbf{k}n}` to reflect the fact that occupied and unoccupied states tend to shift in opposite directions by the renormalization. One can also specify an absolute frequency mesh by [:math:`{...}`], i.e., relative to the Fermi energy. This is mandatory for ``FULL`` calculations. It is sometimes a bit inconvenient to determine suitable values for the upper and lower bound of [:math:`{...}`].
Therefore, Spex allows the usage of wildcards for one of the boundaries or both (see below). The second argument to ``CONTOUR`` gives the increment of the (equidistant) realfrequency mesh for :math:`{W(\omega')}`. The lower and upper boundaries of this mesh are predetermined already by the first argument. As a third method, Spex also allows to omit the second argument altogether. Then, it uses a ''hybrid'' method where the screened interaction is analytically continued from the imaginary (where it has to be known for the integral along this axis, see above) to the real axis, thereby obviating the need of calculating and storing ''W'' on a realfrequency mesh. The disadvantage is that the badly controlled Pade extrapolation introduces an element of randomness (also see keyword ``SMOOTH`` below). Our experience so far is that the ''hybrid'' method is inbetween the two other methods in terms of both computational cost and numerical accuracy.

++++
 Examples  ``CONTOUR 0.01 0.005``   Use contour integration to obtain the two 
    values :math:`{\Sigma^\mathrm{xc}(\epsilon\pm 0.01)}`, giving :math:`{\Sigma^\mathrm{xc}}` as a linear 
    function. 
++++
  ``CONTOUR {0.1:0.15,0.01} 0.005``   Calculate :math:`{\Sigma^\mathrm{xc}(\omega)}` on a frequency mesh relative to 
    the KS energy :math:`{\epsilon}`. 
++++
  ``CONTOUR [{*:*,0.01}] 0.005``   Use an absolute frequency mesh 
    (relative to the Fermi energy) instead. 
    Wildcards are used for the upper and lower bounds. 
++++
  ``CONTOUR [{*:*,0.01}]``  Use ``hybrid`` method. 
++++

``FREQINT`` (``SENERGY``)
(*) Independently of whether :math:`{\Sigma^\mathrm{xc}(i\omega_n)}` (``CONTINUE``) or :math:`{\Sigma^\mathrm{xc}(\omega_n)}` (``CONTOUR``) is evaluated, an important step in the calculation of the selfenergy is to perform the frequency convolution :math:`{\int_{\infty}^\infty G(z+i\omega')W(i\omega')d\omega'}` with :math:`{z=i\omega}` or :math:`{z=\omega}`. For this frequency integration, we interpolate ``W`` and then perform the convolution with the Green function analytically. The keyword ``FREQINT`` determines how the interpolation should be done. It can take two values, ``SPLINE`` and ``PADE`` for spline [after the transformation :math:`{\omega\rightarrow \omega/(1+\omega)}`] and Pade interpolation. The default is ``SPLINE``. In the case of ``GT`` calculations, there is a similar frequency integration with the ``T`` matrix replacing ``W``. There, the default is ``PADE``.

++++
 Examples  ``FREQINT SPLINE``   Use spline interpolation for ``W`` (or ``T``) in the frequency 
    convolution :math:`{G\cdot W}` (:math:`{G\cdot T}`). 
++++
  ``FREQINT PADE``  Use Pade interpolation. 
++++

+
+The second method is *contour integration*, in which the frequency integration is performed explicitly,
+however not along the real frequency axis but on a deformed integration contour that avoids the real frequency axis as best as possible.
+This integration contour starts from :math:`{\infty}`, describes an infinite quarter circle to :math:`i\infty`,
+then runs along the imaginary frequency axis to :math:`{i\infty}`, and finishes, again with an infinite quarter circle, to :math:`\infty`.
+The two quarter circles do not contribute to the integral (because the integrand behaves as :math:`{\propto \omega^{2}}`).
+Furthermore, depending on the frequency argument :math:`{\omega}` of the selfenergy,
+one has to add a few residues coming from the poles of the Green function in the interval :math:`[0,\epsilon_\mathrm{F}\omega]`
+if :math:`{\omega<\epsilon_\mathrm{F}}` and :math:`[\epsilon_\mathrm{F}\omega,0]` otherwise,
+which requires the correlation screened interaction :math:`{W^\mathrm{c}(\omega)}` to be evaluated on this interval of the real axis.
+As a consequence, the calculations are more demanding in terms of computational complexity and cost (time and memory).
+Also, contour integration requires additional input parameters and is therefore somewhat more difficult to use.
+However, the results are more accurate. In particular, they are not affected by the "illdefinedness" of the analytic continuation.
+
+The corresponding keyword is called ``CONTOUR`` and belongs to the section ``SENERGY``.
+Obviously, ``CONTOUR`` and ``CONTINUE`` must not be given at the same time. The keyword ``CONTOUR`` expects two arguments.
+The first defines the frequencies :math:`{\omega}`, for which the selfenergy :math:`{\Sigma^\mathrm{xc}(\omega)}`
+is to be evaluated. At least, two frequencies are needed to approximate the selfenergy as a linear function in :math:`{\omega}` and,
+thus, to calculate the ''linearized'' solution of the quasiparticle equation [Eq. :eq:`qppert`]). For this, a single value suffices
+(example ``0.01``), with which the selfenergy for a state :math:`{\mathbf{k}n}` is evaluated at two frequencies
+(:math:`{\epsilon_{\mathbf{k}n}0.01}` htr and :math:`{\epsilon_{\mathbf{k}n}+0.01}` htr).
+The more accurate ''direct'' solution is only available if we specify a range of frequencies for the selfenergy instead of a single number,
+making an interpolation of the selfenergy beyond a linear function possible.
+This is possible by an argument such as ``{0.1:0.15,0.01}``. Here, the range of values is relative to :math:`\epsilon_{\mathbf{k}n}>\epsilon_\mathrm{F}`.
+The range is reversed (to {0.15:0.1,0.01} in the example) for occupied states (:math:`\epsilon_{\mathbf{k}n}<\epsilon_\mathrm{F}`)
+to reflect the fact that occupied and unoccupied states tend to shift in opposite directions by the renormalization.
+(Thus, it makes sense to choose a frequency range shifted to larger energies instead of a symmetrical one.)
+One can also specify an absolute frequency mesh by ``[{...}]``, where the energies are given relative to the Fermi energy.
+This is mandatory for calculations with ``FULL`` (:numref:`beyondpert`).
+It is sometimes a bit inconvenient to determine suitable values for the upper and lower bound of ``[{...}]``.
+Therefore, Spex allows the usage of wildcards for one of the boundaries or both (see examples).
+The second argument to ``CONTOUR`` gives the increment of the (equidistant) realfrequency mesh for :math:`{W(\omega')}`.
+The lower and upper boundaries of this mesh are determined automatically from the first argument. As a third method,
+Spex also allows the second argument to be omitted altogether. Then, it uses a *hybrid* method where the screened interaction is
+analytically continued from the imaginary axis (where it has to be known for the evaluation of the integral along this axis, see above) to the real axis,
+thereby obviating the need of calculating and storing *W* on a realfrequency mesh.
+The disadvantage is that the badly controlled PadÃ© extrapolation introduces an element of randomness
+(also see keyword ``SMOOTH`` below). Our experience so far is that the *hybrid* method is inbetween the two other methods
+in terms of both computational cost and numerical accuracy.
+
+.. listtable:: Examples
+ :widths: 70 100
+
+ *  ``CONTOUR 0.01 0.005``
+  Use contour integration to obtain the two values :math:`{\Sigma^\mathrm{xc}(\epsilon\pm 0.01\,\mathrm{htr})}` with the KS energy :math:`\epsilon`, giving :math:`{\Sigma^\mathrm{xc}}` as a linear function.
+ *  ``CONTOUR {0.1:0.15,0.01} 0.005``
+  Calculate :math:`{\Sigma^\mathrm{xc}(\omega)}` on a frequency mesh relative to the KS energy :math:`{\epsilon}`.
+ *  ``CONTOUR [{*:*,0.01}] 0.005``
+  Use an absolute frequency mesh (relative to the Fermi energy) instead. Wildcards are used for the upper and lower bounds.
+ *  ``CONTOUR [{*:*,0.01}]``
+  Use *hybrid* method.
+
+FREQINT (SENERGY)
+
+(*) Independently of whether :math:`{\Sigma^\mathrm{xc}(i\omega_n)}` (``CONTINUE``) or :math:`{\Sigma^\mathrm{xc}(\omega_n)}` (``CONTOUR``) is evaluated,
+an important step in the calculation of the selfenergy is to perform the frequency convolution
+:math:`{\int_{\infty}^\infty G(z+i\omega')W(i\omega')d\omega'}` with :math:`{z=i\omega_n}` or :math:`{z=\omega_n}`.
+For this frequency integration, we interpolate *W* and then perform the convolution with the Green function analytically.
+The keyword ``FREQINT`` determines how the interpolation should be done. It can take the two arguments ``SPLINE`` (default) and ``PADE``
+for spline [after the transformation :math:`{\omega'\rightarrow \omega'/(1+\omega')}`] and PadÃ© interpolation, respectively.
+In the case of *GT* calculations, there is a similar frequency integration with the *T* matrix replacing *W*.
+There, the default is ``PADE``.
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``FREQINT SPLINE``
+  Use spline interpolation for *W* (or *T*) in the frequency convolution :math:`{G\cdot W}` (:math:`{G\cdot T}`).
+ *  ``FREQINT PADE``
+  Use PadÃ© interpolation.
+
+.. _SMOOTH:
SMOOTH (SENERGY)

(*) We have already mentioned that the Pade approximation can lead to spurious features in the extrapolated or interpolated quantity, e.g., the selfenergy or the screened interaction (also see ``FREQINT``), if one of the Pade poles happen to lie close to the real (or imaginary) frequency axis. To avoid such unphysical results, we can make use of an averaging procedure where, for each set of values, :math:`{\{\Sigma^\mathrm{xc}(i\omega_n)\}}` or :math:`{\{W(i\omega_n)\}}`, we evaluate a large number of Pade approximants with randomly shifted frequency points according to :math:`{\omega_n\rightarrow\omega_n(1+10^{7}r_n)}` with random numbers :math:`{r\in[0.5,0.5]}` and average over them. This simple but effective way to smoothen the functions is activated with the keyword ``SMOOTH`` in section ``SENERGY``. It can have up to two arguments, giving the number of Pade approximants to average over, the first for the selfenergy, the second for the screened interaction or the ``T`` matrix in the case of ``GT`` calculations (``FREQINT PADE``). Smoothening is usually unnecessary for ``GW`` calculations but is helpful in the case of the ``GT`` approach. By default, no smoothening is applied. Note that this feature leads to a small degree of randomness in the results.

++++
 Examples  ``SMOOTH 500``  Average over 500 Pade approximants to smoothen the selfenergy. 
++++
  ``SMOOTH 500 250``   In addition, average over 250 approximants to 
    smoothen ``W`` (or the ``T`` matrix). (Requires ``FREQINT PADE``.) 
++++

+(*) We have already mentioned that the PadÃ© approximation can lead to spurious features in the extrapolated or interpolated quantity,
+i.e., the selfenergy or the screened interaction (also see ``FREQINT``), if one of the PadÃ© poles happen to lie close to the real
+(or imaginary) frequency axis. To avoid such unphysical results, we can make use of an averaging procedure where,
+for each set of values, :math:`{\{\Sigma^\mathrm{xc}(i\omega_n)\}}` or :math:`{\{W(i\omega_n)\}}`,
+we evaluate a large number of PadÃ© approximants with randomly shifted frequency points according to
+:math:`{\omega_n\rightarrow\omega_n(1+10^{7}r_n)}` with random numbers :math:`{r_n\in[0.5,0.5]}` and average over them.
+This simple but effective way to smoothen the functions is activated with the keyword ``SMOOTH`` in section ``SENERGY``.
+It can have up to two arguments, giving the number of PadÃ© approximants over which to average, the first for the selfenergy,
+the second for the screened interaction or the *T* matrix in the case of *GT* calculations (``FREQINT PADE``).
+Smoothening is usually unnecessary for *GW* calculations but is helpful in the case of the *GT* approach.
+By default, no smoothening is applied. Note that this feature leads to a small degree of randomness in the results.
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``SMOOTH 500``
+  Average over 500 PadÃ© approximants to smoothen the selfenergy.
+ *  ``SMOOTH 500 250``
+  In addition, average over 250 approximants to smoothen *W* (or the *T* matrix). (Requires ``FREQINT PADE``.)
+
NOZERO (SENERGY)

Both the exchange (HF) and the correlation selfenergy contain terms of the form :math:`{\langle...\rangle\langle...\rangle I(\mathbf{k})}` with :math:`{I=v}` and :math:`{I=W}`, respectively. The quantities :math:`{\langle...\rangle}` are the projections discussed in the section :ref:`mbp` and behave as :math:`{\propto k}` in the longwavelength limit :math:`{k\rightarrow 0}`. The prefactor is obtained from :math:`{k.p}` perturbation theory. Multiplying with :math:`{I\propto k^{2}}` gives rise to zeroorder terms involving that prefactor. These zeroorder terms (which only appear at the :math:`{\Gamma}`; point) can improve the kpoint convergence. However, in the case of systems with small band gaps, the zeroorder terms can become unnaturally large (because the integrand is strongly varying in the neighborhood of k=0 in these systems), impeding a favourable convergence with respect to the kpoint sampling. Therefore, the keyword ``NOZERO`` in section ``SENERGY`` allows the zeroorder terms to be neglected. They are included by default.
+Both the exchange (HF) and the correlation selfenergy contain terms of the form
+:math:`\langle...\rangle\langle...\rangle V(\mathbf{k})`
+with :math:`{V=v}` and :math:`{V=W}`, respectively. The quantities :math:`{\langle...\rangle}` are the "projections" discussed in :numref:`mpb`.
+One of the involved projections acquires the form
+:math:`\langle e^{i\mathbf{k}r} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle`, which
+behaves as :math:`{\propto k}` in the longwavelength limit :math:`{k\rightarrow 0}` and
+:math:`n\neq n'`. The prefactor is obtained from :math:`k\cdot p`
+perturbation theory. Multiplying two of these projections with :math:`V\propto k^{2}` (or one with :math:`V\propto k^{1}`,
+which is valid for the "wings" of *W*) gives rise to zeroorder terms.
+These zeroorder terms (which only appear at the :math:`{\Gamma}` point) can improve the kpoint convergence.
+However, in the case of systems with small band gaps, the zeroorder terms can become unnaturally large
+(because the integrand is strongly varying in the neighborhood of k=0 in such systems),
+impeding a favourable convergence with respect to the kpoint sampling.
+Therefore, the keyword ``NOZERO`` (section ``SENERGY``) allows the zeroorder terms to be neglected. They are included by default.
ALIGNVXC (SENERGY)

(*) The equation :eq:`qpeq` depends explicitly on the meanfield starting point through the eigensolutions :math:`{\{\phi_{\mathbf{k}n}\}}`. This dependence is more obvious in the equation :eq:`qppert` quasiparticle equation, which contains the exchangecorrelation potential explicitly. Thus, a constant shift :math:`{v^\mathrm{xc}\rightarrow v^\mathrm{xc}+\Delta}` will not only shift the quasiparticle energies as a whole (as it would the KS energies) but also individually so that energy differences (e.g., the quasiparticle band gap) can depend on this constant shift. (Again, this is different from the KS gap, which would not be affected.) The underlying reason is that quasiparticle energies represent totalenergy differences (between the ``N`` and the ``N``:math:`{\pm}` 1 particle system) and are thus defined as absolute energies, which depend individually on the "energy zero" set by :math:`{v^\mathrm{xc}}`. One can utilize this dependence to simulate the selfconsistency condition that the input and output ionization potentials (valenceband maximum) be identical; in other words, the ionization potential should not change by the quasiparticle renormalization. The keyword ``ALIGNVXC`` in section ``SENERGY`` enforces this condition. To this end, the valenceband maximum should be included in the job definition. (Spex uses the highest occupied one among the states defined after ``JOB`` for the correction.) One can also specify the shift explicitly as a realvalued argument after ``ALIGNVXC``.

++++
 Examples  ``ALIGNVXC``   Align exchangecorrelation potential in such a way that the 
    ionization potential remains unchanged by the quasiparticle 
    correction. 
++++
  ``ALIGNVXC 0.2eV``   Apply a constant positive shift of 0.2 eV to the 
    exchangecorrelation potential. 
++++

+(*) Equation :eq:`qpeq` depends explicitly on the meanfield starting point through the eigensolutions :math:`{\{\phi_{\mathbf{k}n}\}}`.
+This dependence is more obvious from Eq. :eq:`qppert`, which contains the exchangecorrelation potential explicitly.
+Thus, a constant shift :math:`{v^\mathrm{xc}\rightarrow v^\mathrm{xc}+\Delta}` will not only shift the quasiparticle
+energies as a whole (as it would the KS energies) but also individually so that energy differences (e.g., the quasiparticle band gap)
+can depend on this constant shift. (Again, this is different from the KS gap, which would not be affected.)
+The underlying reason is that quasiparticle energies represent totalenergy differences (between the :math:`N` and the :math:`N{\pm}1` particle system)
+and are thus defined as absolute energies, which depend individually on the "energy zero" set by :math:`{v^\mathrm{xc}}`.
+One can utilize this dependence to simulate the selfconsistency condition that the input and output ionization potentials
+(valenceband maximum) be identical; in other words, the ionization potential should not change by the quasiparticle renormalization.
+The keyword ``ALIGNVXC`` in section ``SENERGY`` enforces this condition. To this end, the valenceband maximum should be included in the job definition.
+(Spex uses the highest occupied one among the states defined after ``JOB`` for the correction.)
+One can also specify the shift explicitly as a realvalued argument after ``ALIGNVXC``. This is useful if a different criterion is to be used,
+for example, requiring the Fermi energy (metallic case) or core state energies (core spectroscopy) to remain unchanged.
+In this case, the shift has to be determined manually by repeatedly trying different shift values until the constancy of the desired quantity
+is achieved (because the shift affects the results, so there is a feedback effect). In this trialanderror approach, it is very helpful
+to use ``RESTART 2`` (:numref:`RESTART_GW`), in which case the selfenergy is read from harddisc and the calculation takes very little time.
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``ALIGNVXC``
+  Align exchangecorrelation potential in such a way that the ionization potential remains unchanged by the quasiparticle correction.
+ *  ``ALIGNVXC 0.2eV``
+  Apply a constant positive shift of 0.2 eV to the exchangecorrelation potential.
+
VXC (SENERGY)

(*) The matrix elements of the exchangecorrelation potential are needed in the solution of the equation :eq:`qppert`. By default, these matrix elements are taken from the input data provided by the DFT code (e.g., in the file "vxcfull" written by Fleur). Alternatively, Spex can calculate the matrix using the potentials from the input data. (This is required for PBE0.) For that, the keyword ``VXC`` in section ``SENERGY`` can take the arguments ``READ`` and ``CALC``.
+(*) The matrix elements of the exchangecorrelation potential are needed in the solution of Eq. :eq:`qppert`.
+By default, these matrix elements are taken from the input data provided by the DFT code (e.g., in the file "vxcfull" written by Fleur). Alternatively, Spex can calculate the matrix using the potentials from the input data. (This is required for PBE0.) For that, the keyword ``VXC`` in section ``SENERGY`` can take the arguments ``READ`` (default) and ``CALC``.
+
+.. listtable:: Examples
+ :widths: 16 100
+
+ *  ``VXC READ``
+  Read :math:`{v^\mathrm{xc}}` expectation values from harddisc.
+ *  ``VXC CALC``
+  Calculate :math:`{v^\mathrm{xc}}` expectation values using the xc potential defined in stars and lattice harmonics.
++++
 Examples  ``VXC READ``  Read :math:`{v^\mathrm{xc}}` expectation values from harddisc. 
++++
  ``VXC CALC``   Calculate :math:`{v^\mathrm{xc}}` expectation values 
    using the xc potential defined in stars and lattice harmonics.
++++
+.. _RESTART_GW:
RESTART

Spex can reuse data from a previous ``GW`` run that has finished successfully, crashed, or has been stopped by the user. A ``GW`` calculation consists mainly of a loop over the irreducible k points. For each k point, Spex (a) calculates the matrix ``W``(``'k``') and (b) updates the selfenergy matrix (or expectation values) :math:`{\Sigma^{xc}}` with the contribution of the current k point (and its symmetryequivalent k points). After completion of step (b), the current selfenergy is always written to the (binary) files "spex.sigx" and "spex.sigc". If the ``RESTART`` option is specified (independently of its argument), Spex also writes the matrix ``W(k)`` (in addition to some other data) to the (HDF5) file "spex.cor" unless it is already present in that file. If it is present, the corresponding data is read instead of being calculated. In this way, the keyword ``RESTART`` enables reusage of the calculated ``W(k)`` from a previous run. The matrix ``W(k)`` does not depend on the k points and band indices defined after ``JOB``. So, these parameters can be changed before a run with ``RESTART``, in which the ``W`` data is then reused. For example, band structures can be calculated efficiently in this way (see below). Especially for long calculations, it is recommended to use the ``RESTART`` option. Spex can also restart a calculation using selfenergy data contained in "spex.sigx" and "spex.sigc". To this end, an argument is added: ``RESTART 2``. Spex then starts with the k point, at which the calculation was interrupted before. In contrast to "spex.cor", the files "spex.sigx" and "spex.sigc" do depend on the job definition, which must therefore ``not`` be changed before a run with ``RESTART 2``. However, there are few parameters (see below) that may be changed before a rerun with ``RESTART 2``. These concern details of solving the quasiparticle equation, which follows after completion of the selfenergy calculation. The following logical table gives an overview.

++++
  ``spex.cor``  ``spex.sigx/c`` 
+===============+==============+=================+
     write 
++++
 ``RESTART``  readwrite  write 
++++
 ``RESTART 2``  readwrite  readwrite 
++++

The different rules for "spex.cor" and "spex.sigx/c" are motivated by the facts that (a) the file "spex.cor" is much bigger than "spex.sigx/c" (so, writing of "spex.cor" to harddisc should not be the default), (b) the files "spex.sigx/c" include the updated selfenergy (requiring more computation than for ``W``, thus representing "more valuable" data).

There are other (binary) files that are written during a program run and that may be reused later. The ones relevant for ``GW`` (same rules as for "spex.cor") are
+Spex can reuse data from a previous *GW* run that has finished successfully, crashed, or has been stopped by the user.
+A *GW* calculation consists mainly of a loop over the irreducible k points. For each k point, Spex (a) calculates the matrix :math:`W(\mathbf{k})` and (b) updates the selfenergy matrix (or expectation values) :math:`{\Sigma^\mathrm{xc}}` with the contribution of the current k point (and its symmetryequivalent k points). After completion of step (b), the current selfenergy is always written to the (binary) files "spex.sigx" and "spex.sigc". If the ``RESTART`` option is specified (independently of its argument), Spex also writes the matrix :math:`W(\mathbf{k})` (in addition to some other data) to the (HDF5) file "spex.cor" unless it is already present in that file. If it is present, the corresponding data is read instead of being calculated. In this way, the keyword ``RESTART`` enables reusage of the calculated :math:`W(\mathbf{k})` from a previous run. The matrix :math:`W(\mathbf{k})` does not depend on the k points and band indices defined after ``JOB``. So, these parameters can be changed before a run with ``RESTART``, in which the *W* data is then reused.
+For example, band structures can be calculated efficiently in this way (:numref:`GWbandstr`).
+Especially for long calculations, it is recommended to use the ``RESTART`` option.
+Spex can also restart a calculation using selfenergy data contained in "spex.sigx" and "spex.sigc".
+To this end, the argument 2 is added: ``RESTART 2``. Spex then starts with the k point, at which the calculation was interrupted before.
+In contrast to "spex.cor", the files "spex.sigx" and "spex.sigc" do depend on the job definition,
+which must therefore **not** be changed before a run with ``RESTART 2``.
+However, there are few parameters (see below) that may be changed before a rerun with ``RESTART 2``.
+These concern details about how the quasiparticle equation is solved after completion of the selfenergy calculation. The following logical table gives an overview.
+
+.. listtable::
+ :widths: 30 30 30
+ :headerrows: 1
+
+ * 
+  ``spex.cor``
+  ``spex.sigx/c``
+ *  
+  
+  write
+ *  ``RESTART``
+  readwrite
+  write
+ *  ``RESTART 2``
+  readwrite
+  readwrite
+
+The different rules for "spex.cor" and "spex.sigx/c" are motivated by the facts that (a) the file "spex.cor" is much bigger than "spex.sigx/c" (so, writing of "spex.cor" to harddisc should not be the default), (b) the files "spex.sigx/c" include the updated selfenergy (requiring more computation than for *W*, thus representing "more valuable" data).
+
+There are other (binary) files that are written during a program run and that may be reused later.
+The ones relevant for *GW* (same rules as for "spex.cor") are
* ``spex.mb``: contains MT part of product basis,
* ``spex.ccou``: contains ``contracted`` screened interaction (i.e., summed over ``'k``'); this contraction is needed for the selfenergy core contribution (``CORES``), ``COSX``, and ``IBC``,
* ``spex.core``: contains core contribution to ``W`` (keyword ``CORES``).

The ``RESTART`` option can be exploited to customize the calculations. For example, it is possible to perform a selfconsistent QS ``GW`` while keeping ``W`` fixed to the LDA level throughout the iterations. (The screened interaction is then always read from "spex.cor" and the MT product basis from "spex.mb".) As another example, when ``CORES`` is used, the MT product basis additionally contains products of core states with LAPW basis functions. To avoid that, one can run a calculation without ``CORES`` first (which can be stopped just after the product basis is written to "spex.mb"), followed by a calculation with ``RESTART``. The product basis is then read from "spex.mb" (excluding the corebasis products) instead of being constructed anew.

As already mentioned above, when preliminary data is read from one of the restart files, changing input parameters might conflict with the data in the files. Spex tries to detect such conflicts. In case of conflict, Spex either stops the calculation ("spex.cor") or recalculates the data and overwrites the files ("spex.sigx/c"). It is, however, not guaranteed that Spex can detect all possible conflicts. So, care has to be taken if parameters are changed in calculations with ``RESTART``.

The following is a (incomplete) list of parameters that may be changed for a calculation with ``RESTART 2``, in which the files "spex.sigx/c" are read: ``ALIGNVXC``, ``VXC``, ``SPECTRAL``, ``SMOOTH``, ``CONTINUE``, ``WRITE``, all of section ``WANNIER`` (e.g., for Wannier interpolation).
In the case of ``RESTART 1`` (or ``RESTART``) (i.e, when "spex.cor" is to be read), one may change most parameters except: ``MESH``, 2nd argument of ``CONTOUR``, ``JOB`` types (most cases). Many changes of parameters will be overridden with values from the restart files (e.g., parameters in the section ``MBASIS`` by data from the file "spex.mb"). Other changed parameters might affect only the selfenergy but not the ``W``, for example, ``NBAND``. (The latter allows to check the band convergence of ``W`` and :math:`{\Sigma^{xc}}`, separately.) Again, this can be used efficiently to customize a calculation.

++++
 Examples  ``RESTART``   Spex tries to read the screened interaction ``W`` 
    (and other data) from the file "spex.cor". 
    If the required data does not exist, the calculation proceeds 
    normally and appends ``W`` to the file (reusable in a later restart). 
++++
  ``RESTART 1``   Same as ``RESTART``. 
++++
  ``RESTART 2``   Implies ``RESTART``. Spex reads the selfenergy from the 
    files "spex.sigx" and "spex.sigc" if they exist. If not, the 
    calculation proceeds normally. 
++++


GW band structures

In this section, three ways of calculating ``GW`` band structures are discussed.

* Using the :ref:`kpath` label (single Spex run)
* With :ref:`kptadd` additional :math:`q` points (multiple Spex run)
* Wannier interpolation (single ; but long ; Spex run).

The calculation of a ``GW`` band structure is not as straightforward as in the case of KSDFT, where the local nature of the effective potential (e.g., LDA or GGA) allows the KS Hamiltonian to be diagonalized independently for each q point along a qpoint path. The ``GW`` selfenergy, on the other hand, is nonlocal. As a consequence, its evaluation involves a convolution in reciprocal space :math:`{\Sigma^\mathrm{xc}(\mathbf{q})=i\sum_\mathbf{q}G(\mathbf{q+k})W(\mathbf{k})}` (in simplified notation), requiring summations over the full Brillouin zone for any arbitrarily chosen :math:`q` point. So, in addition to :math:`q`, one would need to include all shifted points :math:`q+k`, too, where :math:`k` stands for the k points defined in ``BZ``.

If we restrict ourselves to the set of k points defined in ``BZ``, the band structure along a certain highsymmetry line (e.g., :math:`{\GammaX}`) can consist only of the few k points that happen to lie on this line. Sometimes this is already sufficient. For example, if the ``GW`` renormalization affects the band dispersion only little, then a few points along the highsymmetry line are enough to shift/modify the KS bands accordingly and produce, in this way, a quasiparticle band structure. Such a calculation can be performed in a single Spex run.

.. _kpath:
+* ``spex.ccou``: contains "contracted" screened interaction (i.e., summed over **k**); this contraction is needed for the selfenergy core contribution (``CORES``), for ``COSX``, and for ``IBC``,
+* ``spex.core``: contains core contribution to *W* (keyword ``CORES``, :numref:`CORES_GW`).
+
+The ``RESTART`` option can be exploited to customize the calculations. For example, it is possible to perform a selfconsistent QSGW calculation
+while keeping *W* fixed to the LDA level throughout the iterations. (The screened interaction is then always read from "spex.cor" and the MT product basis from "spex.mb".) As another example, when ``CORES`` is used, the MT product basis additionally contains products of core states with LAPW basis functions. To avoid that, one can run a calculation without ``CORES`` first (which can be stopped just after the product basis is written to "spex.mb"), followed by a calculation with ``RESTART``. The product basis is then read from "spex.mb" (excluding the corebasis products) instead of being constructed anew.
+
+As already mentioned above, changing input parameters in "spex.inp" and running a calculation with ``RESTART`` might lead to a conflict
+between the changed parameters and the data read from the restart files.
+Spex tries to detect such conflicts. In case of conflict, Spex either stops the calculation ("spex.cor") or recalculates the data and overwrites the files ("spex.sigx/c"). It is, however, not guaranteed that Spex can detect all possible conflicts. So, care has to be taken if parameters are changed in calculations with ``RESTART``.
+
+The following is a (incomplete) list of parameters that may be changed for a calculation with ``RESTART 2``, in which the files "spex.sigx/c" are read:
+``ALIGNVXC``, ``VXC``, ``SPECTRAL``, ``SMOOTH``, ``CONTINUE``, ``WRITE``, all of section ``WANNIER`` (e.g., for Wannier interpolation).
+In the case of ``RESTART 1`` (or just ``RESTART``) (i.e, when "spex.cor" is to be read), one may change most parameters except: ``MESH``, 2nd argument of ``CONTOUR``, ``JOB`` types (most cases). Many changes of parameters will be overridden with values from the restart files (e.g., parameters in the section ``MBASIS`` by data from the file "spex.mb"). Other changed parameters might affect only the selfenergy but not the
+screened interaction *W*, for example, ``NBAND``. (The latter allows the band convergence of *W* and :math:`{\Sigma^\mathrm{xc}}` to be checked separately.) Again, this can be used efficiently to customize a calculation.
+
+.. listtable:: Examples
+ :widths: 18 100
+
+ *  ``RESTART``
+  Spex tries to read the screened interaction *W* (and other data) from the file "spex.cor". If the required data does not exist, the calculation proceeds normally and appends *W* to the file (reusable in a later restart).
+ *  ``RESTART 1``
+  Same as ``RESTART``.
+ *  ``RESTART 2``
+  Implies ``RESTART``. Spex reads the selfenergy from the files "spex.sigx" and "spex.sigc" if they exist. If not, the calculation proceeds normally.
+
+.. _GWbandstr:
+
+*GW* band structures
+
+In this section, three ways of calculating *GW* band structures are discussed.
+
+* Using the ``KPTPATH`` label (single Spex run, :numref:`KPTPATH`),
+* With the special ``+`` label in ``KPT`` for additional q points (multiple Spex runs, :numref:`KPTadd`),
+* Wannier interpolation (single  but long  Spex run, :numref:`INTERPOL_GW`).
+
+The calculation of a *GW* band structure is not as straightforward as in the case of KSDFT, where the local nature of the effective potential (e.g., LDA or GGA) allows the KS Hamiltonian to be diagonalized independently for each q point along a qpoint path.
+The *GW* selfenergy, on the other hand, is nonlocal. As a consequence,
+its evaluation involves a convolution in reciprocal space
+:math:`{\Sigma^\mathrm{xc}(\mathbf{q})=i\sum_\mathbf{q}G(\mathbf{q+k})W(\mathbf{k})}` (in simplified notation),
+requiring summations over the full Brillouin zone for any arbitrarily chosen q point. So, in addition to :math:`\mathbf{q}`,
+one would need to include all shifted points :math:`\mathbf{q}+\mathbf{k}`, as well, where :math:`\mathbf{k}` stands for the k points defined in ``BZ``.
+
+If we restrict ourselves to the set of k points defined in ``BZ``, the band structure along a certain highsymmetry line (e.g., :math:`{\GammaX}`)
+can consist only of the few k points that happen to lie on this line. Sometimes this is already sufficient. For example, if the *GW* renormalization
+affects the band dispersion only little, then a few points along the highsymmetry line are enough to shift/modify the KS bands accordingly and produce,
+in this way, a quasiparticle band structure. Such a calculation can be performed in a single Spex run.
+
+.. _KPTPATH:
KPTPATH

As an example, we want to plot a band structure for the path from L over :math:`{\Gamma}`; to :math:`X` for our Si example. In principle, we would have to identify all kpoint indices along this path and set the job definition accordingly, but Spex can do this for you. The kpoint path can be defined with a line ``KPTPATH (L,1,X)``. In the :ref:`job` job definition, one can then use the special label ``PATH`` to address the corresponding list of k points, e.g., ``JOB GW PATH:(110)``. Spex will automatically choose the k points that lie on the path defined by ``KPTPATH`` and calculate the ``GW`` quasiparticle energies. As usual, the results are written to the Spex :ref:`spex.out` file. The data can be extracted from the output file with the spex.extr utility: ``spex.extr g b spex.out > bandstr`` (assuming that the output has been piped into the file "spex.out"). The file "bandstr" then contains a plottable list of xy values (momentum/energy) for each quasiparticle band. However, in most cases the band structure will not be smooth enough because of the small number of k points. One solution is, as already mentioned, to modify a KS band structure. Another possibility is to interpolate the energies with splines; the spex.extr utility has this capability. Finally, one could simply use much denser kpoint sets, but this would entail very expensive calculations.

++++
 Examples  ``KPTPATH (L,1,X)``   Define kpoint path from L over the 
    :math:`\Gamma`; point (k index is always 1) to X. 
    The labels L and X must be defined with [[#KPT``KPT``]] 
++++
  ``JOB GW PATH:(110)``   Run ``GW`` calculation for the first ten bands at all 
    k points defined by ``KPTPATH``. 
++++

.. _kptadd:
+As an example, we want to plot a band structure for the path from L over :math:`\Gamma` to :math:`X` for our Si example. In principle, we would have to identify all kpoint indices along this path and set the job definition accordingly, but Spex can do this for you. The kpoint path can be defined with a line ``KPTPATH (L,1,X)``. In the job definition, one can then use the special label ``PATH`` to address the corresponding list of k points, e.g., ``JOB GW PATH:(110)``. Spex will automatically choose the k points that lie on the path defined by ``KPTPATH`` and calculate the *GW* quasiparticle energies. As usual, the results are written to the output file. The data can be extracted from the output file with the spex.extr utility: ``spex.extr g b spex.out > bandstr`` (assuming that the output has been piped into the file "spex.out"). The file "bandstr" then contains a plottable list of xy values (momentum/energy) for each quasiparticle band. However, in most cases the band structure will not be smooth enough because of the small number of k points. One solution is, as already mentioned, to modify a KS band structure. Another possibility is to interpolate the energies with splines; the spex.extr utility has this capability. Finally, one could simply use much denser kpoint sets, but this would entail very expensive calculations.
+
+.. listtable:: Examples
+ :widths: 36 100
+
+ *  ``KPTPATH (L,1,X)``
+  Define kpoint path from L over the :math:`\Gamma` point (k index is always 1) to X. The labels L and X must be defined with ``KPT``.
+ *  ``JOB GW PATH:(110)``
+  Run *GW* calculation for the first ten bands at all k points defined by ``KPTPATH``.
+
+.. _KPTadd:
KPT +=...

This section explains an extension to the keyword :ref:`kpt` (already described before), which will later enable the calculation of smooth band structures. The extension allows adding arbitrary q points, one at a time, allowing one to investigate the quasiparticle spectrum at momenta that are not element of the original kpoint set (defined by ``BZ``). For example, in Si the conductionband minimum is not at a highsymmetry k point but at around 3/4 along the line :math:`{\Gamma  X}`. This particular q point can be defined using the special k label "+" by ``+=[0,0,0.75]`` (alternatively, ``+=3/8*(1,1,0)`` in internal coordinates) after the keyword :ref:`kpt`. As explained above, each additional point :math:`q` must be accompanied with all shifted points :math:`q+k`. So, before running Spex, we must let the DFT code generate the wavefunctions and energies at the shifted kpoint set ``'q+k``'. This is done in the usual way by creating the kpoint file :ref:`wrtkpt` or ``w`` flag) and running the DFT code. (In the case of Fleur, the shifted k points are appended to the list in the file "kpts".) The additional q point can be addressed in the job definition with the special ``+`` label, e.g., ``JOB GW 1:(15) +:(15)``, which will yield the bands 15 for the Γ point and the additional :math:`q`. The energy difference between the fourth state at Γ and the fifth state at :math:`q` yields the fundamental ``GW`` gap.


++++
 Examples  ``KPT +=1/7*[1,0,0]``   Define special k point at 1/7 along the line :math:`{\GammaX}` 
    (here, the X point in x direction) of an fcc structure. 
++++
  ``KPT +=1/14*(0,1,1)``   The same in internal coordinates. 
++++
  ``JOB GW +:(110)``   Run ``GW`` calculation for the states 110 at the added q point. 
++++

The possibility of adding arbitrary q points enables the calculation of smooth band structures. To this end, a ``GW`` run (together with the generation of the eigenstates with the DFT code) has to be performed for each q point in a list of q points that make up a path in the Brillouin zone, e.g., from :math:`{\Gamma}`; to :math:`X`. This list is defined with the keyword ``KPTPATH`` as explained above. Fortunately, the screened interaction has to be calculated only once and can be reused for the other points. For this, the keyword ``RESTART`` should be given in the input file. Of course, we also need a job definition, such as ``JOB GW +:(110)``. Here, only the additional q point should be specified. (It would be difficult to extract the data for the band structure plot otherwise.) Now, running Spex with ``KPTPATH`` defined in the input file and :ref:`wrtkpt` (or with the option ``w``) creates two files, containing the usual list of k points (Fleur: "kpts") and the list of q points (Fleur: "qpts"). For each of the q points, we would have to run, in this order, Spex, Fleur, and Spex again. To simplify this task, there is a shell script called "spex.band". This shell script performs all the necessary steps automatically (it uses the same environment variables as `selfc.sh` and produces, for each q point, one output file named ``spex_???.out`` (where ``???`` is a threedigit counting index, or more digits if ``???>999``).

From these files the bandstructure data is extracted with ``spex.extr g b spex_???.out > bandstr``. The data is written to the file ``bandstr``. The band structure can then be plotted with xmgrace or gnuplot. If you extract the bandstructure data, instead, with ``spex.extr g b c spex_??.out > bandstr``, you get an additional column, which contains the imaginary parts of the quasiparticle energies. These imaginary parts are proportional to the line width and inversely proportional to the excitation lifetime. The line widths can be plotted together with the band structure in gnuplot with ``plot "bandstr" using 1:2:(10*abs(3)) with points ps var``. You will see that the bands get wider the further they are away from the Fermi energy.

As a third alternative, we can use Wannier interpolation to interpolate between the few ``'k``' points of the original k mesh. The construction of Wannier orbitals is explained in a :ref:`wannier`. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to ``spex.inp``:
+This section explains an extension to the keyword ``KPT`` (:numref:`KPT`), which will later enable the calculation of smooth band structures.
+The extension enables adding arbitrary q points, one at a time, allowing one to investigate the quasiparticle spectrum at momenta that
+are not element of the original kpoint set (defined by ``BZ``).
+For example, in Si the conductionband minimum is not at a highsymmetry k point but at around 3/4 along the line :math:`{\Gamma  \mathrm{X}}`.
+This particular q point can be defined using the special k label "+" by ``+=[0,0,0.75]`` (alternatively, ``+=3/8*(1,1,0)`` in internal coordinates)
+after the keyword ``KPT``. As explained above, each additional point :math:`\mathbf{q}` must in principle be accompanied with all shifted points :math:`\mathbf{q}+\mathbf{k}`.
+So, before running Spex, we must let the DFT code generate the wavefunctions and energies at the shifted kpoint set :math:`\{\mathbf{q}+\mathbf{k}\}`.
+This is done in the usual way by creating the kpoint file (:numref:`WRTKPT`) and running the DFT code.
+(In the case of Fleur, the shifted k points are appended to the list in the file "kpts".)
+The additional q point can be addressed in the job definition with the special "+" label, e.g., ``JOB GW 1:(15) +:(15)``,
+which will yield the bands 15 for the :math:`\Gamma` point and the additional q point.
+The energy difference between the fourth state at :math:`\Gamma` and the fifth state at :math:`\mathbf{q}` yields the fundamental *GW* gap.
+
+.. listtable:: Examples
+ :widths: 36 100
+
+ *  ``KPT +=1/7*[1,0,0]``
+  Define special k point at 1/7 along the line :math:`{\Gamma  \mathrm{X}}` (here, the X point in *x* direction) of an fcc structure.
+ *  ``KPT +=1/14*(0,1,1)``
+  The same in internal coordinates.
+ *  ``JOB GW +:(110)``
+  Run *GW* calculation for the states 110 at the added q point.
+
+The possibility of adding arbitrary q points enables the calculation of smooth band structures.
+To this end, a *GW* run (together with the generation of the eigenstates with the DFT code)
+has to be performed for each q point in a list of q points that make up a path in the Brillouin zone,
+e.g., from :math:`{\Gamma}` to :math:`X`. This list is defined with the keyword ``KPTPATH`` as explained above.
+Fortunately, the screened interaction has to be calculated only once and can be reused for the other points.
+For this, the keyword ``RESTART`` should be given in the input file. Of course, we also need a job definition,
+such as ``JOB GW +:(110)``. Here, only the additional q point should be specified. It would be difficult to
+extract the data for the band structure plot otherwise.) Now, running Spex with ``KPTPATH`` defined in the input
+file and ``WRTKPT`` (or with the command line option ``w``) creates two files, containing the usual list of
+k points (Fleur: "kpts") and the list of q points (Fleur: "qpts"). For each of the q points, we would have to run,
+in this order, Spex, Fleur, and Spex again. To simplify this task, there is a shell script called "spex.band".
+This shell script performs all the necessary steps automatically. It uses the same environment variables as "spex.selfc" (:numref:`beyondpert`)
+and produces, for each q point, one output file named "spex_NNN.out", where NNN is a threedigit counting index, or more digits if the index exceeds 999.
+
+From these files the bandstructure data is extracted with ``spex.extr g b spex_???.out > bandstr``.
+The data is written to the file ``bandstr``. The band structure can then be plotted with "xmgrace" or "gnuplot".
+If you extract the bandstructure data, instead, with ``spex.extr g b c spex_???.out > bandstr``,
+you get an additional column, which contains the imaginary parts of the quasiparticle energies.
+These imaginary parts are proportional to the line width and inversely proportional to the excitation lifetime.
+The line widths can be plotted together with the band structure in "gnuplot" with ``plot "bandstr" using 1:2:(10*abs(3)) with points ps var``.
+You will see that the bands get wider the further they are away from the Fermi energy.
+
+If "spex.inp" contains a ``SPECTRAL`` definition (:numref:`SPECTRAL`), the script "spex.band" also renames the file "spectral",
+written in each Spex run, to "spectralNNN" with the same counting index NNN as above. The spectral data can
+be compiled with the "spex.extr" utility (e.g., ``spex.extr s b spectral??? > spec``)
+into a single data file containing a list of xyz value triples: the Bloch momenta as the first column, the frequency as
+the second, and the value of the spectral function as the third. With "gnuplot", for example,
+a momentumfrequency plot of the spectral function with color coding is then prepared with the command
+``plot "spec" using 1:2:3 palette with lines``. (The respective color range is set with ``set cbrange [...]``,
+see the Gnuplot manual.)
+
+.. _INTERPOL_GW:
+
+INTERPOL
+
+As a third alternative, we can use Wannier interpolation to interpolate between the few k points of the original k mesh.
+The construction of Wannier orbitals is explained in :numref:`wannier`. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to ``spex.inp``:
.. codeblock:: bash
@@ 329,16 +583,57 @@ As a third alternative, we can use Wannier interpolation to interpolate between
MAXIMIZE
INTERPOL
END

Please also remove the entry ``+=[0,0,0.75]`` from the ``KPT`` line. The first line lets SPEX calculate quasiparticle energies in the whole irreducible Brillouin zone (IBZ). (The Wannier interpolation can be understood as a backandforth Fourier transformation with a realspace truncation of matrix elements inbetween. The Fourier transformation from momentum to real space involves the band energies in the whole BZ.) The other lines define a set of maximally localized Wannier functions (MLWFs) of bonding :math:`sp^3` hybrid orbitals. These four hybrid orbitals are generated from the four lowest valence bands (first two arguments after ``ORBITALS``). Running Spex will construct the MLWFs and use them to interpolate the band structure. The band structure data is written to the files "bands0" for the KS energies and "bands1" for the ``GW`` quasiparticle energies and can be plotted with xmgrace or gnuplot. The file "bands1" has one data column more than "bands0". This additional column contains the interpolated imaginary parts of the quasiparticle energies. You can plot the band structure together with the variable line thickness as described above.

Wannier interpolation can also be used in conjunction with the keyword :ref:`spectral`. The interpolated spectral function is written to the file "spectralw" or, if the system is spindependent, to the files "spectralw1" and "spectralw2" for spin up and spin down, respectively.

Beyond perturbative GW and QSGW

One can also go beyond the perturbative solution of the equation :eq:`qpeq` and solve it by iterative diagonalization in the basis of the singleparticle states. This requires the full selfenergy matrix including offdiagonal elements to be calculated. The respective line in the input file would read
``JOB GW FULL 1:(110) X:(110) L:(110)`` giving the selfenergy as 10x10 matrices at the Γ, X, and L points. (In principle, it is also possible to have a band definition like ``(15,9,10)``, which would yield a 7x7 matrix. Degenerate states, not included in the definition, are automatically included.) In ``FULL`` calculations, the Spex code uses irreducible representations (irreps) to speed up the calculation. As a result, the bands are grouped in terms of these irreps (called "blocks") in the output. For example, in the example below (for the X point) the block is composed of the bands 3, 4, 9, and 10. If a block contains only a single band or a single subgroup of degenerate states, the full solution is omitted because it would be identical to the perturbative (diagonal) solution. The full solutions are tabulated in the output::
+Please also remove the entry ``+=[0,0,0.75]`` from the ``KPT`` line. The first line lets SPEX calculate quasiparticle energies
+in the whole irreducible Brillouin zone (IBZ). (The Wannier interpolation can be understood as a backandforth Fourier
+transformation with a realspace truncation of matrix elements inbetween.
+The Fourier transformation from momentum to real space involves the band energies in the whole BZ.)
+The section ``WANNIER`` defines a set of maximally localized Wannier functions (MLWFs) of bonding sp\ :sup:`3` hybrid orbitals. These four hybrid orbitals are generated from the four lowest valence bands (first two arguments after ``ORBITALS``). Running Spex will construct the MLWFs and use them to interpolate the band structure. The band structure data is written to the files "bands0" for the KS energies and "bands1" for the *GW* quasiparticle energies and can be plotted with "xmgrace" or "gnuplot". The file "bands1" has one data column more than "bands0". This additional column contains the interpolated imaginary parts of the quasiparticle energies.
+You can plot the band structure with the line thickness scaled by the imaginary part in the same way as described before.
+
+Wannier interpolation can also be used in conjunction with the keyword ``SPECTRAL`` (:numref:`spectral`).
+The interpolated spectral function is written to the file "spectralw" or, if the system is spindependent,
+to the files "spectralw1" and "spectralw2" for the spin up and spin down channels, respectively. These files
+can be plotted in the same way as described in :numref:`KPTadd`.
+
+.. _CORES_GW:
+
+CORES
+
+By default, the *n* summation in the Green function
+
+.. math:: G_0(\mathbf{r},\mathbf{r}';\omega)=\sum_\mathbf{k}\sum_n \frac{\phi_{\mathbf{k}n}(\mathbf{r})\phi^*_{\mathbf{k}n}(\mathbf{r}')}{\omega\epsilon_{\mathbf{k}n}+ i\eta\,\mathrm{sgn}(\epsilon_{\mathbf{k}n}\epsilon_\mathrm{F})}
+ :label: lehmann
+
+for the selfenergy [Eq. :eq:`selfene`] includes only the valence and conduction states represented in the LAPW basis but excludes the core states.
+In theory, of course, the core states should be included. Their contribution is, however, numerically small. Spex allows selective core states
+to be included in the evaluation of the selfenergy with the keyword ``CORES``. This keyword also affects the calculation of the polarization function, see
+:numref:`CORES_P`. As arguments, Spex expects a comma separated list of core states between round brackets, e.g., ``(2s,2p,3d)``, for each atom type.
+*Empty* definitions ``()`` are allowed. If ``CORESOC`` is specified, one can also address the SOCsplitted states individually, e.g., ``(2p1/2)``.
+
+.. listtable:: Examples
+ :widths: 32 100
+
+ *  ``CORES () (1s)``
+  Include the 1s state of the second atom type.
+ *  ``CORES (2s,2p1/2)``
+  Include the 2s and the 2p\ :sup:`1/2` states but exclude the 2p\ :sup:`3/2` states.
+
+.. _beyondpert:
+
+Beyond perturbative *GW* and QSGW
+
+One can also go beyond the perturbative solution of Eq. :eq:`qpeq` and solve it by iterative diagonalization in the basis of the singleparticle states.
+This requires the full selfenergy matrix including offdiagonal elements to be calculated.
+The respective line in the input file would read ``JOB GW FULL 1:(110) X:(110) L:(110)``
+giving the selfenergy as :math:`10\times 10` matrices at the :math:`\Gamma`, X, and L points.
+(In principle, it is also possible to have a band definition like ``(15,9,10)``, which would yield a :math:`7\times 7` matrix.
+Degenerate states, not included in the definition, are automatically included.)
+In ``FULL`` calculations, the Spex code uses irreducible representations (irreps) to speed up the calculation.
+As a result, the bands are grouped in terms of these irreps (called "blocks") in the output.
+For example, in the example below (for the X point) the block is composed of the bands 3, 4, 9, and 10.
+If a block contains only a single band or a single subgroup of degenerate states,
+the full solution is omitted because it would be identical to the perturbative (diagonal) solution. The full solutions are tabulated in the output::
Bd KS olap HF olap GW
3 2.92306 0.99971 0.12683 0.99985 2.50496 0.06006
@@ 346,29 +641,34 @@ One can also go beyond the perturbative solution of the equation :eq:`qpeq` and
9 16.77136 0.99732 23.25660 0.99962 16.97051 0.30294
10 16.77136 0.99732 23.25660 0.99962 16.97051 0.30294
``HF`` are the HartreeFock values (from diagonalization). They are ``'not``' ordered according to energy but related to the KS states with a condition of maximal overlap (also making sure that there is a onetoone correspondence). The overlap is given on the left. The full ``GW`` energies are related to the KS states as well. This is due to the fact that the quasiparticle equation is nonlinear in energy and must, therefore, be solved iteratively for each quasiparticle state. The iterations are started with the KS energy :math:`{\epsilon_{\mathbf{k}n}}` of that state, giving :math:`{H^\mathrm{KS}+\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})v^\mathrm{xc}}` as the (complex) quasiparticle Hamiltonian. The diagonalization then yields a spectrum of states, of which the one with the largest overlap is picked for the next iteration. In the first few iterations, the offdiagonal elements are switched on adiabatically to alleviate a smooth convergence. The overlap of the final quasiparticle wave function with the original KS state is given, too. Although, this approach usually yields welldefined quasiparticle solutions, it cannot be ruled out that the iterative procedure, when started at two different energies, might result in the same quasiparticle energy. This can happen especially for highlying states, which do not have a welldefined quasiparticle character anymore.
+The values with the label ``HF`` are the HartreeFock values (from diagonalization).
+They are **not** necessarily ordered according to energy but related to the KS states with a condition of maximal overlap
+(and at the same time making sure that a onetoone correspondence can be established).
+The corresponding overlap is given on the left of the HF column.
+The full *GW* energies (column labeled ``GW``; the last column lists the imaginary parts) are related to the KS states as well.
+This is due to the fact that the quasiparticle equation is nonlinear in energy and must, therefore, be solved iteratively for each quasiparticle state. The iterations are started with the KS energy :math:`{\epsilon_{\mathbf{k}n}}` of that state, giving :math:`{H^\mathrm{KS}+\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})v^\mathrm{xc}}` as the (complex) quasiparticle Hamiltonian. The diagonalization then yields a spectrum of states, of which the one with the largest overlap is picked for the next iteration. In the first few iterations, the offdiagonal elements are switched on adiabatically to alleviate a smooth convergence. The overlap of the final quasiparticle wavefunction with the original KS state is given, too. Although, this approach usually yields welldefined quasiparticle solutions, it cannot be ruled out that the iterative procedure, when started at two different energies, might result in the same quasiparticle energy. This can happen especially for highlying states, which do not have a welldefined quasiparticle character anymore.
If the job definition contains ``FULL`` and IBZ, the full ``GW`` selfenergy matrix is evaluated for the whole IBZ, which enables selfconsistent calculations in the framework of the ``quasiparticle selfconsistent GW`` (QSGW) approach. In this approach, one creates a meanfield system from the ``GW`` selfenergy whose singleparticle energies are as close as possible to the quasiparticle energies. This meanfield system is subsequently solved to selfconsistency in a DFT code. The resulting solution can then serve as a starting point for a new oneshot ``GW`` calculation, which constitutes the second iteration and so on until the quasiparticle energies are converged. The construction of the meanfield system is, to some extent, arbitrary. We use the following definition, which is slightly modified from the original work [PRL 93, 126406]:
+If the job definition contains ``FULL`` and ``IBZ``, the full *GW* selfenergy matrix is evaluated for the whole IBZ, which enables selfconsistent calculations in the framework of the quasiparticle selfconsistent *GW* (QSGW) approach. In this approach, one creates a meanfield system from the *GW* selfenergy whose singleparticle energies are as close as possible to the quasiparticle energies. This meanfield system is subsequently solved to selfconsistency in a DFT code. The resulting solution can then serve as a starting point for a new oneshot *GW* calculation, which constitutes the second iteration and so on until the quasiparticle energies are converged. The construction of the meanfield system is, to some extent, arbitrary. We use the following definition, which is slightly modified from the original work [Phys. Rev. Lett. 93, 126406]:
.. math::
\displaystyle A_{\mathbf{k}nn}=Z_{\mathbf{k}n}^{1} \langle \phi_{\mathbf{k}n}  \Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})  \phi_{\mathbf{k}n} \rangle
for diagonal elements and
.. math::
+.. math::
\displaystyle A_{\mathbf{k}nn'}=\langle \phi_{\mathbf{k}n}  \Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})+\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n'})  \phi_{\mathbf{k}n'} \rangle
for offdiagonal 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:
+for offdiagonal 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:
* ``rm fleur.qsgw``  remove any previous version of the *hermitianized* matrix.
+* ``rm fleur.qsgw``  remove the file "fleur.qsgw" containing the *hermitianized* matrix used internally by Fleur from previous QSGW iterations.
* ``rm broyd*``  remove Broyden information about previous iterations because this information is inconsistent with the new Hamiltonian (the SCF calculation does not converge otherwise).
* Set ``gw=3`` in the Fleur input file.
* Run Fleur.
 Fleur translates the file "spex.qsgw" (in basis of eigenstates) to a file "fleur.qsgw" (in LAPW basis). Then, a SCF run is performed where instead of the KS Hamiltonian :math:`{H^\mathrm{KS}}` the Hamiltonian :math:`{H^\mathrm{KS}+(\Sigma^{\mathrm{xc,H}}v^\mathrm{xc})}` is employed. After the run finishes, the ``gw=2`` step has to be performed:
+ Fleur translates the file "spex.qsgw" (in basis of eigenstates) to a file "fleur.qsgw" (in LAPW basis). Then, a SCF run is performed where, instead of the KS Hamiltonian :math:`{H^\mathrm{KS}}`, the Hamiltonian :math:`{H^\mathrm{KS}+(\Sigma^{\mathrm{xc,H}}v^\mathrm{xc})}` is employed. After the run finishes, the ``gw=2`` step has to be performed.
* Set ``gw=2`` in the Fleur input file.
* Run Fleur.
Apart from the usual output data for a ``GW`` calculation, Fleur writes, in addition, the file "qsgw", which contains the matrix elements of :math:`{\Sigma^{\mathrm{xc,H}}v^\mathrm{xc}}` in the basis of the new eigenstates. The file "qsgw" is later read by Spex to subtract correctly the double counting. Of course, doing these steps manually is tedious and also errorprone. Therefore, we provide the shell script ``spex.selfc``, which performs the steps needed for a selfconsistent calculation automatically. For example, ``spex.selfc 5`` will run five iterations. The shell script assumes that Spex and Fleur are run with ``spex.x`` and ``fleur.x``, respectively. Environment variables can be used to change this:
+Apart from the usual output data for a *GW* calculation, Fleur writes, in addition, the file "qsgw", which contains the matrix elements of :math:`{\Sigma^{\mathrm{xc,H}}v^\mathrm{xc}}` in the basis of the new eigenstates. The file "qsgw" is later read by Spex to subtract correctly the double counting. Of course, doing these steps manually is tedious and also errorprone. Therefore, we provide the shell script ``spex.selfc``, which performs the steps needed for a selfconsistent calculation automatically. For example, ``spex.selfc 5`` will run five iterations. The shell script assumes that Spex and Fleur are run with ``spex.x`` and ``fleur.x``, respectively. Environment variables can be used to change this:
* ``SPEX_EXEC``  Spex executable,
* ``FLEUR_EXEC``  Fleur executable,
@@ 376,103 +676,353 @@ Apart from the usual output data for a ``GW`` calculation, Fleur writes, in addi
* ``FLEUR_QUEUE``  Same for Fleur,
* ``QUEUE``  Default queue command to be used for Spex and Fleur.
The job types ``HF``, ``PBE0``, ``SX``, and ``COSX`` also allow ``FULL`` in the definition, i.e., a full matrix is evaluated, enabling a selfconsistent calculation in the same way as QSGW (see above). The file names used are the same: ``spex.qsgw``, ``fleur.qsgw``, ``qsgw``, even though it might actually be a HF calculation. The job ``GT`` (also ``GWT``) does not allow for selfconsistency at the moment.

+The job types ``HF``, ``PBE0``, ``SX``, and ``COSX`` also allow ``FULL`` in the definition, i.e., a full matrix is evaluated, enabling a selfconsistent calculation in the same way as QSGW (see above). The file names used are the same: "spex.qsgw", "fleur.qsgw", "qsgw", even though it might actually be a HF calculation. The job ``GT`` (also ``GWT``) does not allow for selfconsistency at the moment.
GW with MPI


RESTART...


MPIBLK


MPISYM

.. _polar:
Polarization function
=====================
The polarization function gives the linear change in the electronic density of a noninteracting system with respect to changes in the effective potential. It is, thus, a fundamental quantity in the calculation of screening properties of a manyelectron 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)=1P(\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
+The polarization function gives the linear change in the electronic density of a noninteracting system with respect to changes in the effective potential. It is, thus, a fundamental quantity in the calculation of screening properties of a manyelectron 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)=1P(\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
.. math::
 \scriptstyle P_{\mu\nu}(\mathbf{k},\omega)=2\sum_{\mathbf{q}}^{\mathrm{BZ}}\sum_{n}^{\mathrm{occ}}\sum_{n'}^{\mathrm{unocc}}\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle\langle \phi_{\mathbf{k+q}n'}  \phi_{\mathbf{q}n} M_{\mathbf{k}\nu} \rangle \cdot\left(\frac{1}{\omega+\epsilon_{\mathbf{q}n}\epsilon_{\mathbf{q}+\mathbf{k}n'}+i\eta}\frac{1}{\omega\epsilon_{\mathbf{q}n}+\epsilon_{\mathbf{q}+\mathbf{k}n'}i\eta}\right) =\int_{\infty}^\infty \frac{S_{\mu\nu}(\mathbf{k},\omega')}{\omega\omega'+i\eta\mathrm{sgn}(\omega')}d\omega'\,.
+ P_{\mu\nu}(\mathbf{k},\omega)&=2\sum_{\mathbf{q}}^{\mathrm{BZ}}\sum_{n}^{\mathrm{occ}}\sum_{n'}^{\mathrm{unocc}}\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle\langle \phi_{\mathbf{k+q}n'}  \phi_{\mathbf{q}n} M_{\mathbf{k}\nu} \rangle \\
+ & \cdot \left(\frac{1}{\omega+\epsilon_{\mathbf{q}n}\epsilon_{\mathbf{q}+\mathbf{k}n'}+i\eta}\frac{1}{\omega\epsilon_{\mathbf{q}n}+\epsilon_{\mathbf{q}+\mathbf{k}n'}i\eta}\right) \\
+ & = \int_{\infty}^\infty \frac{S_{\mu\nu}(\mathbf{k},\omega')}{\omega\omega'+i\eta\mathrm{sgn}(\omega')}d\omega'\,.
:label: eqP
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.)
+A system without spin polarization has been assumed, hence the factor 2. In case of spin polarization, the formula has a
+spin summation instead of the factor 2 and corresponding spin quantum numbers. In the case of spinorbit coupling (and no spin polarization),
+two spin summation labels are added, one in front of each "projection" :math:`\langle...\rangle`.
+
+We have implicitly defined the spectral function *S* in Eq. :eq:`eqP`, 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.)
The expression for equation :eq:`eqP` involves an infinite band summation over :math:`{n'}`. In practice, this summation is truncated by the number of bands (keyword `nbabd`, which thus becomes an important convergence parameter. The selfenergy contains an infinite band summation as well.
+The expression for Eq. :eq:`eqP` involves an infinite emptyband summation over :math:`{n'}`. In practice, this summation is truncated by the number of bands (keyword ``NBAND``), which thus becomes an important convergence parameter. The selfenergy contains an infinite band summation as well.
In principle, the k summation is infinite, too: In an infinite crystal, there are infinitely many k points, and the k summation formally turns into a k integration. In practice, of course, there can only by a discreet sampling of the Brillouin zone. The tetrahedron method enables a geometrical interpolation between the explicit k points and, thus, an approximate k integration. Therefore, the tetrahedron method is the default and recommended method. However, there are two other implemented methods: Simple summation over the kpoint mesh (see ``HILBERT NONE``) and the Gaussian method (see ``GAUSS``).
+In principle, the k summation is infinite, too: In an infinite crystal, there are infinitely many k points, and the k summation formally turns into a k integration. In practice, of course, there can only be a discrete sampling of the Brillouin zone. The tetrahedron method enables a geometrical interpolation between the explicit k points and, thus, an approximate k integration. Therefore, the tetrahedron method is the default and recommended method. However, there are two other implemented methods: Simple summation over the kpoint mesh (see ``HILBERT NONE``) and the Gaussian method (see ``GAUSS``).
+
+HILBERT (SUSCEP)
+
+The most important keyword in the calculation of the polarization function is ``HILBERT``, which defines a (Hilbert) mesh for the frequency integration in Eq. :eq:`eqP`. (The old keyword ``FSPEC`` still works but is deprecated.) The Hilbert mesh, used in conjunction with the tetrahedron and the Gaussian method, extends from the lowest (roughly the band gap) to the highest virtual electron transition energy. (Note that the spectral function is symmetric :math:`{S(\omega)=S(\omega)}`, which makes a mapping to :math:`{\int_0^\infty...}` possible.) The lower and upper bounds are, thus, predetermined by the input data. As the spectral function usually shows more variation at low energies (also, these energies are more important from a perturbation theory point of view), we employ an exponential mesh, which is dense at low and coarse at high energies. Two parameters are necessary to define the Hilbert mesh :math:`{\{w_i\}}`, for example, the first step size :math:`{\omega_2\omega_1}` and the stretching factor :math:`{a=(\omega_{i+1}\omega_i)/(\omega_i\omega_{i1})}`. These two parameters are accepted as realvalued arguments (example: ``HILBERT 0.01 1.04``). However, with this definition it is not straightforward to increase the density of the mesh without changing the exponential form of the mesh. Therefore, it is recommended to use a different definition, namely using an integer as first argument (without a period!) and a real number as second (example: ``HILBERT 30 40.0``; the second argument is always interpreted as realvalued, so ``HILBERT 30 40`` is equivalent). The first argument gives the number of mesh points between 0 and 5 htr and the second argument is the "accumulated" stretching factor at 5 htr, i.e., by what factor the step size at 5 htr has increased from :math:`{\omega_2\omega_1}`. (The reason for defining an arbitrary but fixed energy of 5 htr is that modifying the energy range, e.g., by changing ``NBAND``, then leaves the form of the Hilbert mesh unchanged.) In this way, increasing the first argument makes the mesh denser but does not affect the exponential form of the mesh, which, in turn, is governed by the second argument. So, there are two possible definitions (distinguished by whether the first argument has a period or not). For convenience, Spex always writes the equivalent parameters of the other definition to the output. By default, Spex uses a Hilbert mesh defined by ``HILBERT 50 30.0`` for bandgap materials and ``HILBERT 150 50.0`` for metals, unless a spectral frequency mesh is defined (e.g., ``DIELEC {0:1,0.01}``), in which case Spex adjusts the Hilbert mesh to this frequency mesh. An idea of how well a given Hilbert mesh resolves the spectral function can be got by specifying ``WRITE INFO`` in the input file. Then, a file "spex.sf.nnn" (where nnn is a threedigit counting index) is written for each k point, and it contains the spectral function for the head element of *P* on the Hilbert mesh as plottable data. One can also define ``HILBERT NONE`` as a special option, which implements the k summation as simple sums over the k points without interpolation or broadening. ``HILBERT NONE`` is only available if the polarization function is to be calculated for purely imaginary frequencies (including zero).
+
+.. listtable:: Examples
+ :widths: 34 100
+
+ *  ``HILBERT NONE``
+  Simple summation over k points.
+ *  ``HILBERT 50 30.0``
+  Use Hilbert mesh with fifty frequencies (between 0 and 5 htr) and an accumulated stretching factor of 30 at 5 htr.
+ *  ``HILBERT 0.01 1.05``
+  Use Hilbert mesh with a first step size of :math:`{\omega_2\omega_1=}` 0.01 htr and a stretching factor of :math:`{a=1.05}`. (First argument is realvalued.)
+
+MULTDIFF (SUSCEP)
+
+(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of Eq. :eq:`eqP` 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 :math:`k\cdot 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 :math:`{\Gamma}` point (:math:`k=0`), Spex multiplies the linear term with this energy difference, resulting in smooth and nondivergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the frequency integration of Eq. :eq:`eqP`, 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 at :math:`k=0`. The behavior can be changed with the keyword ``MULTDIFF`` in the section ``SUSCEP``.
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``MULTDIFF OFF``
+  Never separate (divergent) energy differences.
+ *  ``MULTDIFF ON``
+  Always separate energy differences (i.e., for all k points).
+ *  ``MULTDIFF INT``
+  Use default behavior (separate for :math:`k=0`, do not for :math:`k\neq 0`) but stick energy differences into integration weights.
+ *  ``MULTDIFF INTON``
+  Always separate energy differences by sticking them into integration weights.
+
+.. _PLASMA:
+
+PLASMA (SUSCEP)
+
+(*) In the case of metals, the polarization function contains an additional term, the socalled Drude term, which gives rise to metallic screening (diverging static dielectric constant, shortrange static *W*). The Drude term stems from virtual intraband transitions across the Fermi surface. It can be formulated with the square of the plasma frequency, which in turn is evaluated by an integration over the Fermi surface. The Drude term can
+be treated analytically and, as long as the Fermi surface is sufficiently big, it normally does not pose a numerical problem. However, a very small Fermi surface eventually leads to a very sharp "Drude peak" in the *GW* selfenergy, impeding a straightforward numerical solution of the nonlinear quasiparticle equation. One could also say that, while the Drude term is actually treated correctly, it gets too much weight because of the coarseness of the kpoint mesh. It is, therefore, sometimes helpful to modify the plasma frequency. This is possible with the keyword ``PLASMA``. Its argument replaces the plasma frequency calculated by the code. Setting ``PLASMA 0`` switches the Drude term off altogether. Of course, ``PLASMA 0`` thus also disables metallic screening, which might be unphysical. Therefore, there is a special option, ``PLASMA METAL``, which scales the head element of :math:`{W(\mathbf{k},\omega)}` in the limit :math:`{k\rightarrow 0}` to enforce metallic screening. This latter option can be helpful, for example, in the case of semimetallic systems with a tiny Fermi surface.
+
+.. listtable:: Examples
+ :widths: 24 100
+
+ *  ``PLASMA 1.5eV``
+  Set plasma frequency manually to 1.5 eV.
+ *  ``PLASMA 0``
+  Set plasma frequency to zero. Disables Drude term.
+ *  ``PLASMA METAL``
+  Disable Drude term but enforce metallic screening by scaling the head element of *W*.
+
+.. _CORES_P:
+
+CORES
+
+By default, the *n* summation of Eq. :eq:`eqP` extends only over the valence states (represented in the LAPW basis) and excludes the core states.
+The core states can be included in this summation in the same way as for the *GW* selfenergy (see :numref:`CORES_GW`). If ``CORESOC`` is specified,
+the SOCsplitted core states are included, otherwise the SOCaveraged core states (see :numref:`CORESOC`). In the former case, the Weiss field
+of a spinpolarized system can further lift the degeneracies. This is taken into account as well. See :numref:`CORES_GW` for details and examples.
TETRAF (SUSCEP)

(*) In rare cases, especially for very large kpoint sets, the determination of the tetrahedron weights (Timing 'wghts') can become computationally expensive. This is because Spex calculates tetrahedron weights for all k points by default, which makes it possible to average over equivalent k points, thereby avoiding a (slight) symmetry breaking connected with the spatial form of the tetrahedra. (The tetrahedra are not unique!) One can disable this averaging by setting ``TETRAF`` in the section ``SUSCEP``. If set, tetrahedron weights are only calculated for the irreducible part of the BZ. This accelerates the calculation of the weights but introduces slight deviations due to symmetry breaking.
+(*) In rare cases, especially for very large kpoint sets, the determination of the tetrahedron weights (Timing ``wghts``) can become computationally expensive. This is because Spex calculates tetrahedron weights for all k points by default, which makes it possible to average over equivalent k points, thereby avoiding a (slight) symmetry breaking connected with the spatial form of the tetrahedra. (The tetrahedra are not unique!) One can disable this averaging by setting ``TETRAF`` in the section ``SUSCEP``. If set, tetrahedron weights are only calculated for the irreducible part of the BZ. This accelerates the calculation of the weights but introduces slight deviations due to symmetry breaking.
WGHTTHR (WFPROD)

(*) In both tetrahedron and Gaussian method, integration weights are calculated to perform the k summation. To be more precise, an integration weight is calculated for each virtual transition, i.e., for each combination of band index pair (occnocc) and k point. One can reduce the number of terms by introducing a threshold value below which the respective weight is neglected. The corresponding keyword is ``WGHTTHR`` in the section ``SUSCEP``. The default value is :math:`{10^{10}}`. It is usually not necessary to change this value in practice.
+(*) In both tetrahedron and Gaussian methods, integration weights are calculated to perform the k summation. To be more precise, an integration weight is calculated for each virtual transition, i.e., for each combination of band index pair (occnocc) and k point. One can reduce the number of terms by introducing a threshold value below which the respective weight is neglected. The corresponding keyword is ``WGHTTHR`` in the section ``SUSCEP``. The default value is :math:`{10^{10}}`. It is usually not necessary to change this value in practice.
GAUSS

(*) The Gaussian method is included mostly for testing. It does not only affect the polarization function but also all other k summations (such as for the HartreeFock exchange potential and the determination of the Fermi energy). Therefore, ``GAUSS`` is a global keyword. It effectively replaces each singleparticle eigenvalue by a Gaussian function with a certain width so that the density of states becomes a smooth function. In the calculation of the polarization function, one additionally needs a finite width for the Fermi edge. The keyword ``GAUSS`` needs the two width parameters as arguments given as real positive values. The advantage of the Gaussian method is its relative mathematical simplicity compared to the tetrahedron method and, in particular, that the Gaussian method cannot give rise to symmetry breaking.
+(*) The Gaussian method is included mostly for testing. It does not only affect the polarization function but also all other k summations (such as for the HartreeFock exchange potential and the determination of the Fermi energy). Therefore, ``GAUSS`` is a global keyword. It effectively replaces each singleparticle eigenvalue by a Gaussian function with a certain width so that the density of states becomes a smooth function. In the calculation of the polarization function, one additionally needs a finite width for the Fermi edge. The keyword ``GAUSS`` needs the two width parameters as arguments given as real positive values. The advantage of the Gaussian method is its relative mathematical simplicity compared to the tetrahedron method and, in particular, that the Gaussian method cannot give rise to symmetry breaking.
++++
 Example  ``GAUSS 0.001 0.01``   Use Gaussian kintegration method with 
    widths 0.001 and 0.01 htr. (Mostly used for testing.) 
++++

HILBERT (SUSCEP)

The most important keyword in the calculation of the polarization function is ``HILBERT``, which defines a (Hilbert) mesh for the frequency integration in equation :eq:`eqP`. (The old keyword ``FSPEC`` still works but is deprecated.) The Hilbert mesh, used in conjunction with the tetrahedron and the Gaussian method, extends from the lowest (roughly the band gap) to the highest virtual electron transition energy. (Note that the spectral function is symmetric :math:`{S(\omega)=S(\omega)}`, which makes a mapping to :math:`{\int_0^\infty...}` possible.) The lower and upper bounds are, thus, predetermined by the input data. As the spectral function usually shows more variation at low energies (also, these energies are more important from a perturbation theory point of view), we employ an exponential mesh, which is dense at low and coarse at high energies. Two parameters are necessary to define the Hilbert mesh :math:`{\{w_i\}}`, for example, the first step size :math:`{\omega_2\omega_1}` and the stretching factor :math:`{a=(\omega_{i+1}\omega_i)/(\omega_i\omega_{i1})}`. These two parameters are accepted as realvalued arguments (example: ``HILBERT 0.01 1.04``). However, with this definition it is not straightforward to increase the density of the mesh without changing the exponential form of the mesh. Therefore, it is recommended to use a different definition, namely using an integer as first argument (without a period!) and a real number as second (example: ``HILBERT 30 40.0``; the second argument is always interpreted as realvalued, so ``HILBERT 30 40`` is equivalent). The first argument gives the number of mesh points between 0 and 5 htr and the second argument is the "accumulated" stretching factor at 5 htr, i.e., by what factor the step size at 5 htr has increased from :math:`{\omega_2\omega_1}`. (The reason for defining an arbitrary but fixed energy of 5 htr is that modifying the energy range, e.g., by changing ``NBAND``, then leaves the form of the Hilbert mesh unchanged.) In this way, increasing the first argument makes the mesh denser but does not affect the exponential form of the mesh, which, in turn, is governed by the second argument. So, there are two possible definitions (distinguished by whether the first argument has a period or not). For convenience, Spex always writes the equivalent parameters of the other definition to the output. By default, Spex uses a Hilbert mesh defined by ``HILBERT 50 30.0`` for bandgap materials and ``HILBERT 150 50.0`` for metals, unless a spectral frequency mesh is defined (e.g., ``DIELEC {0:1,0.01}``), in which case Spex adjusts the Hilbert mesh to this frequency mesh. An idea of how well a given Hilbert mesh resolves the spectral function can be got by specifying ``WRITE INFO`` in the input file. Then, a file "spex.sf.nnn" (where nnn is a threedigit counting index) is written for each k point, and it contains the spectral function for the head element of ``P`` on the Hilbert mesh as plottable data. One can also define ``HILBERT NONE`` as a special option, which implements the k summation as simple sums over the k points without interpolation or broadening. ``HILBERT NONE`` is only available if the polarization function is to be calculated for purely imaginary frequencies (including zero).

++++
 Examples  ``HILBERT NONE``  Simple summation over k points. 
++++
  ``HILBERT 50 30.0``   Use Hilbert mesh with fifty frequencies (between 0 and 5 htr) 
    and an accumulated stretching factor of 30 at 5 htr. 
++++
  ``HILBERT 0.01 1.05``   Use Hilbert mesh with a first step size of :math:`{\omega_2\omega_1=}` 0.01 htr 
    and a stretching factor of :math:`{a=1.05}`. 
    (First argument is realvalued.) 
++++
+.. listtable:: Example
+ :widths: 32 100
+
+ *  ``GAUSS 0.001 0.01``
+  Use Gaussian kintegration method with widths 0.001 and 0.01 htr. (Mostly used for testing.)
MULTDIFF (SUSCEP)
+DISORDER (SUSCEP)

(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of equation :eq:`eqP` 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 :math:`{\Gamma}`; point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and nondivergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the equation :eq:`eqP` 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``.

++++
 Examples  ``MULTDIFF OFF``  Never separate (divergent) energy difference. 
++++
  ``MULTDIFF ON``  Always separate energy difference (i.e., for all k points). 
++++
  ``MULTDIFF INT``   Use default behavior (separate for k=0, do not for k≠0) 
    but stick energy difference into integration weights.a 
++++
  ``MULTDIFF INTON``   Always separate energy difference by sticking them into 
    integration weights. 
++++
+(*) Mathematically, the parameter :math:`{\eta}` in :eq:`eqP` is a positive infinitesimal ensuring the correct time order of the response quantity. Spex effectively treats it as an infinitesimally small positive parameter. Sometimes, however, it can be useful to use a finite value for :math:`{\eta}`, e.g., to simulate disorder in a material. This is possible with the keyword ``DISORDER`` in the section ``SUSCEP``. Note that this keyword has rarely been used so far.
PLASMA (SUSCEP)

(*) In the case of metals, the polarization function contains an additional term, the socalled Drude term, which gives rise to metallic screening (diverging static dielectric constant, shortrange static ``W``). The Drude term stems from virtual intraband transitions across the Fermi surface. It can be formulated with the square of the plasma frequency, which in turn is evaluated by an integration over the Fermi surface. The Drude term can
be treated analytically and, as long as the Fermi surface is sufficiently big, it normally does not pose a numerical problem. However, a very small Fermi surface eventually leads to a very sharp "Drude peak" in the ``GW`` selfenergy, impeding a straightforward numerical solution of the nonlinear quasiparticle equation. One could also say that, while the Drude term is actually treated correctly, it gets too much weight because of the coarseness of the kpoint mesh. It is, therefore, sometimes helpful to modify the plasma frequency. This is possible with the keyword ``PLASMA``. Its argument replaces the plasma frequency calculated by the code. Setting ``PLASMA 0`` switches the Drude term off altogether. Of course, ``PLASMA 0`` would also disable metallic screening, which might be unphysical. Therefore, there is a special option, ``PLASMA METAL``, which scales the head element of :math:`{W(\mathbf{k},\omega)}` in the limit :math:`{k\rightarrow 0}` to enforce metallic screening. This latter option can be helpful, for example, in the case of semimetallic systems with a tiny Fermi surface.

++++
 Examples  ``PLASMA 1.5eV``  Set plasma frequency manually to 1.5 eV. 
++++
  ``PLASMA 0``  Set plasma frequency to zero. Disables Drude term. 
++++
  ``PLASMA METAL``   Disable Drude term but enforce metallic 
    screening by scaling the head element of ``W``. 
++++
+.. listtable:: Example
+ :widths: 26 100
+
+ *  ``DISORDER 1000``
+  Use a finite value :math:`{\eta=1/(2\cdot 1000)\,\mathrm{htr}}`.
DISORDER (SUSCEP)

(*) Mathematically, the parameter :math:`{\eta}` in :eq:`eqP` is a positive infinitesimal ensuring the correct time order of the response quantity. Spex effectively treats it as an infinitesimally small positive parameter. Sometimes, however, it can be useful to use a finite value for :math:`{\eta}`, e.g., to simulate disorder in a material. This is possible with the keyword ``DISORDER`` in the section ``SUSCEP``. Note that the keyword has rarely been used so far.
+.. _spectra:
+
+Spectra
+===============
+The dielectric function :math:`\epsilon_{\mathbf{GG}'}(\mathbf{k},\omega)` is related to several spectroscopic
+experimental techniques where
+the system is perturbed (excited) by some incoming beam of particles but does not lose or gain particles
+(in contrast to photoemission, for example, where electrons are lost from the sample).
+
+In optical spectroscopy measurements, the absorption cross section is proportional to the
+imaginary part of the longwavelength macroscopic dielectric function :math:`\varepsilon_{M}(\omega)`,
+which can be obtained from the microscopic dielectric function by
+
+.. math:: \varepsilon_{M}(\omega)=\lim_{\mathbf{k}\rightarrow 0}\frac{1}{\varepsilon^{1}_\mathbf{00}(\mathbf{k},\omega)}\,.
+
+Here, :math:`\varepsilon^{1}` refers to the matrix inverse.
+
+Another example is electronenergy loss spectroscopy (EELS), where the scattering cross section can be shown to be proportional
+to the head element (:math:`\mathbf{G}=\mathbf{G}'=0`) of the inverse dielectric function
+
+.. math:: \varepsilon^{1}_\mathbf{00}(\mathbf{k},\omega)\,.
+
+Spex can be used to calculate the dielectric function. The respective job is called ``DIELEC``. Arguments are expected to define
+the k point index or label and a frequency range. For example, ``JOB DIELEC X:{0:1,0.001}`` would yield the dielectric function
+evaluated at the X point between 0 and 1 htr in steps of 0.001 htr. Two files are written, "dielec" and "dielecR", containing the
+head elements of the (bare) dielectric function :math:`\varepsilon_\mathbf{00}(\mathbf{k},\omega)` and of the renormalized
+dielectric function :math:`1/\varepsilon^{1}_\mathbf{00}(\mathbf{k},\omega)`. The imaginary part of the latter can directly be
+identified with an optical absorption spectrum (excluding excitonic effects).
+For the EELS spectrum, one would have to plot the imaginary part of the respective
+reciprocal value (for example, in "gnuplot": ``plot "dielecR" using 1:(imag(1/($2+{0,1}*$3)))``).
+
+In the same way, other quantities can be calculated: the polarization function with ``SUSCEP`` (written to the file "suscep"), the renormalized polarization
+function (:math:`R=P+PvR`) with ``SUSCEPR`` (written to the file "suscepR"; "suscep" is also written), the screened interaction with
+``SCREEN`` (written to the file "screen"). In the former case (``SUSCEP``), one can also specify spin indices. For example,
+``JOB SUSCEP X:ud{0:1,0.001}`` restricts the the polarization function to virtual up\ :math:`\rightarrow`\ down transitions.
+Other valid spin labels are ``uu``, ``dd``, ``du``, and ``+`` (spin summed, default).
+
+By default, only the head element is written to the output files. More elements can be written using, e.g.,
+``JOB SUSCEP 1:{0.5eV:2eV,0.001eV},OUT=4``, which would yield :math:`4\times 4` matrices instead of only the head element.
+One can also write the full matrix to a binary file with ``...,OUT=BIN``.
+
+The special ``+`` kpoint label can be used as well (e.g., ``JOB SUSCEP +:{0..50eV,0.01eV}``). This enables using the script
+"spex.band" to prepare a whole list of spectra for a series of Bloch vectors on a highsymmetry line. As explained in
+:numref:`KPTadd` for the file "spectral", the script appends a threedigit counting index to the name of the output file
+("suscep" in the present example, which is thus renamed to "suscep001", "suscep002", ...). The data can then be compiled
+into a single file using the spex.extr utility. (See :numref:`KPTadd` for details).
+
+Total xc energies
+==================
+Spex can calculate total exchange and/or correlation energies of the manyelectron system. The total exchange energy is defined
+by the HF expression
+
+.. math:: E^\mathrm{x}=\frac{1}{2} \sum_{\mathbf{k}n} \sum_{\mathbf{k}'n'} \iint \frac{ \varphi^*_{\mathbf{k}n}(\mathbf{r}) \varphi_{\mathbf{k}'n'}(\mathbf{r}) \varphi^*_{\mathbf{k}'n'}(\mathbf{r}') \varphi_{\mathbf{k}n}(\mathbf{r'}) } { \mathbf{r}\mathbf{r}' } d^3r\,d^3r'\,.
+ :label: hfene
+
+The corresponding job is called ``HFENE``. Internally, it works very much like a HF calculation. So, all parameters pertaining to ``HF``
+are available to ``HFENE``, as well. However, the final result is a single number, :math:`E^\mathrm{x}`.
+
+The total correlation energy can be calculated using the adiabaticconnection fluctuationdissipation theorem (ACFDT) with the
+randomphase approximation (RPA) for the (partially renormalized) response function :math:`\chi_\lambda` (:math:`\lambda=[0,1]`).
+This response function is identical to the polarization function :math:`P` (:numref:`polar`) for :math:`\lambda=0` and
+to the "renormalized polarization function" :math:`R` mentioned in :numref:`spectra` for :math:`\lambda=1`.
+
+The ACFDT method gives the following form of the correlation energy functional
+
+.. math:: E_{\mathrm{c}}[n]=\int_{0}^{1}d\lambda\iint d^{3}rd^{3}r'\frac{1}{\vec{r}\vec{r}'}\times\left[\int_{0}^{\infty}\frac{du}{2\pi}\chi_{\lambda}(\vec{r},\vec{r}',iu)\chi_{0}(\vec{r},\vec{r}',iu)\right]\,
+ :label: ACFDT
+
+where :math:`\chi_0` corresponds to the (bare) polarization function [Eq. :eq:`eqP`] and :math:`\chi_\lambda` is a
+partially renormalized response function (i.e., with a scaled electronelectron interaction) fulfilling
+
+.. math:: \chi_{\lambda}(\vec{r},\vec{r}',\omega)=\chi_{\mathrm{0}}(\vec{r},\vec{r}',\omega)+\iint d^{3}r''d^{3}r'''\chi_{\mathrm{0}}(\vec{r},\vec{r}'',\omega)\\\hspace{0.8cm}\times\left[\frac{\lambda}{\vec{r}''\vec{r}'''}+f_{\mathrm{xc,\lambda}}(\vec{r}'',\vec{r}''',\omega)\right]\chi_{\lambda}(\vec{r}''',\vec{r}',\omega)\,,
+
+here written in the general form including the (scaled) exchangecorrelation kernel :math:`f_{\mathrm{xc},\lambda}`.
+For RPA, :math:`f_{\mathrm{xc},\lambda}=0`.
+
+The evaluation of Eq. :eq:`ACFDT` requires the polarization function to be calculated for all
+(irreducible) k points on an imaginary frequency mesh. The procedure is thus similar to a oneshot *GW* calculation
+and needs a similar computation time. Many of the parameters discussed in connection with the *GW* method
+(e.g., details related to the calculation of the polarization function, :numref:`polar`)
+can be used for ACFDTRPA in the same way. In particular, the same rules using the ``RESTART`` option (:numref:`RESTART_GW`) apply.
+It is, for example, possible to reuse the file "spex.cor" obtained from a previous *GW* calculation for an ACFDTRPA calculation.
+A run with ``RPAENE`` comprises the calculation of the total exchange energy using Eq. :eq:`hfene`.
+
+.. _hubbardu:
+
+Interaction parameters (Hubbard *U*)
+=====================================
+In a manyelectron system, the motion of all electrons are correlated with each other
+through the Coulomb interaction. The Coulomb interaction is thus responsible for the
+fact that we have do deal with a highly complex manybody problem. DFT
+takes electronic correlation into account, but most of the available functionals
+are suitable only for weak to moderate correlation. One way to go beyond standard
+DFT is by manybody perturbation theory, for example, using the *GW*
+approximation. However, this approach has its limits for very strong electronic
+correlations as in Mott insulators, for example.
+
+As an alternative, one can resort to methods with a special treatment of local
+correlations. The first one is the socalled LDA+\ *U*
+scheme, in which the localdensity approximation (LDA) of DFT is
+augmented by an onsite Coulomb repulsion term and an exchange term with the Hubbard
+*U* and Hund exchange *J* parameters, respectively.
+
+A more elaborate computional scheme, which combines manybody model Hamiltonian methods with DFT, is the socalled LDA+DMFT method. In this scheme,
+the interacting manybody system is mapped onto the subspace of localized states,
+for example, formed by d orbitals. The subspace contains much less electrons
+and can thus be treated with highlevel quantum manybody
+approaches such as quantum Monte Carlo, numerical
+renormalization group, diagrammatic techniques, exact diagonalization, et cetera.
+
+The other electrons, not included in the subspace, cannot simply be ignored.
+For example, their screening processes lead to a renormalization of the
+effective electronelectron interaction that acts in the subspace, again
+giving rise to the parameters *U* and *J*. Thus, the Coulomb interaction parameters
+play a crucial role in the study of the correlation effects in solids.
+
+There are two ways to obtain the Hubbard *U* parameter from first principles,
+constrained localdensity approximation (cLDA) and constrained randomphase
+approximation (cRPA). The latter is implemented in the Spex code and relies on the
+fact that, since the electrons outside the localized subspace can be assumed delocalized,
+the randomphase approximation should be appropriate for describing their screening
+contribution. The cRPA offers an efficient way to calculate the effective Coulomb interaction
+parameters in solids. Moreover, cRPA allows individual Coulomb matrix
+elements to be determined, e.g., onsite, offsite, intraorbital, interorbital, and exchange, as well
+as their frequency dependence.
+
+JOB SCREENW ...
+
+In general, the bare and (fully or partially) screened interaction is calculated in reciprocal
+space and represented in the mixed product basis (:numref:`mpb`). The job definition ``SCREENW`` tells Spex
+to project this interaction potential onto the set of local Wannier functions defined in the section ``WANNIER``
+(:numref:`wannier`). The resulting interaction matrix elements
+
+.. math::
+ V_{n_1n_2,n_3n_4}(\omega) = = \iint w_{n_1}^*(\mathbf{r}) w_{n_3}(\mathbf{r}) V(\mathbf{r},\mathbf{r}';\omega) w_{n_2}^*(\mathbf{r}') w_{n_4}^*(\mathbf{r}') d^3r d^3r'
+ :label: screenwdef
+
+are written to the file "spex.cou" for the bare interaction (:math:`V=v`) and the screened interaction
+(:math:`V=U` if ``HUBBARD`` is specified, otherwise :math:`V=W`, see below).
+In addition, if :math:`\omega=0` is included in the frequency mesh,
+effective interaction parameters (HubbardHund and/or Kanamori parameterization) are
+written to the output file for both the bare (*v*) and the screened interaction (*U* or *W*).
+(Note again that, without the keyword ``HUBBARD``, the effective parameters are obtained from the fully screened interaction
+and would thus be unsuitable for a Hubbard model Hamiltonian.)
+The HubbardHund parameters are given for each complete
+Wannier :math:`l` shell (e.g., specifying ``(d)`` in ``ORBITALS`` defines the complete d shell, see :numref:`wannier`),
+here written for the partially screened interaction *U*,
+
+.. math:: U_l = \frac{1}{(2l+1)^2} \sum_{m=l}^l \sum_{m'=l}^l U_{mm',mm'}(0)
+.. math:: J_l = U_l  \frac{1}{2l(2l+1)} \sum_{m=l}^l \sum_{m'=l}^l \left[ U_{mm',mm'}(0)  U_{mm',m'm}(0) \right]
+
+The Kanamori parameters are given for groups of t\ :sub:`2g` (:math:`N=3`) or e\ :sub:`g` orbitals (:math:`N=2`) as
+
+.. math:: U_l = \frac{1}{N} \sum_{m=l}^l U_{mm,mm}(0) \\
+.. math:: U'_l = \frac{1}{N(N1)} \sum_{m=l}^l \sum_{m'=l,m'\neq m}^l U_{mm',mm'}(0) \\
+.. math:: J_l = \frac{1}{N(N1)} \sum_{m=l}^l \sum_{m'=l,m'\neq m}^l U_{mm',m'm}(0)
+
+In the case of a spinpolarized system, the spinup and spindown Wannier functions differ slightly, and the respective spinindices
+have to be specified in the job definition using one of the labels ``uu``, ``dd``, ``du``, or ``+`` (spin summed). There is no default.
+
+.. listtable:: Examples
+ :widths: 60 100
+
+ *  ``JOB SCREENW {0:40eV,0.1eV}``
+  Calculate :math:`W(\mathbf{r},\mathbf{r}';\omega)` and project it onto the Wannier functions; {...} defines the frequency mesh.
+ *  ``JOB SCREENW {0}``
+  The same only for the static case :math:`\omega=0`. (Here, ``{0}`` is short for ``{0:0,1}``.)
+ *  ``JOB SCREENW uu{0}``
+  The same for the spinpolarized case: the upup combination of Wannier functions is selected.
+
+
+HUBBARD (SUSCEP)
+
+With the keyword ``HUBBARD`` (section ``SUSCEP``), the screening that takes place in the correlated subspace is eliminated from the screened interaction, thus
+yielding the partially screened interaction and, with ``JOB SCREENW ...``, the effective interaction parameters
+for a Hubbard model Hamiltonian. The correlated subspace is spanned by the Wannier functions. So, the Wannier definition affects the
+calculation in two ways: (1) It defines the local Wannier basis on which the (bare and) partially screened interaction is
+projected, and (2) it spans the subspace whose screening is eliminated. Sometimes it makes sense to separate these two
+effects, for example, if the subspace is spanned by Wannier functions located at several (possibly symmetryequivalent) atoms,
+but the effective parameters are only needed for one of the atoms. Such a calculation is performed in two steps
+with a file name as argument (example: ``HUBBARD hub``). For the first run (i.e, the file does not exist yet),
+one specifies the full set of Wannier functions. When Spex is run, it just writes information about the
+screening channels to the specified file ("hub") and stops. Before the second run, one can modify the Wannier definition. In the
+present example, one would reduce the Wannier set (retaining the orbitals at only one atom). The second Spex run then calculates
+the partially screened interaction and projects it onto the reduced set of Wannier functions.
+
+.. listtable:: Examples
+ :widths: 28 100
+
+ *  ``HUBBARD``
+  Calculate Hubbard *U* interaction matrix including effective parameters.
+ *  ``HUBBARD hub.dat``
+  Write information about subspace screening to "hub.dat", if the file does not exist, and calculate the Hubbard *U* interaction matrix using data from "hub.dat" otherwise.
+
+.. _RESTART_U:
+
+RESTART
+
+Since ``SCREENW`` calculations can be quite costly, the keyword ``RESTART`` enables a restart of a calculation that has stopped before
+finishing. The usage is similar to *GW* calculations (:numref:`RESTART_GW`). The file "spex.cor" contains the (fully or partially) screened
+interaction for each (completed) k point, and the file "spex.wcou" holds the last update of the projected interaction matrix.
+The logical table is similar to :numref:`RESTART_GW`.
+
+.. listtable::
+ :widths: 30 30 30
+ :headerrows: 1
+
+ * 
+  ``spex.cor``
+  ``spex.wcou``
+ *  
+  
+  write
+ *  ``RESTART``
+  readwrite
+  write
+ *  ``RESTART 2``
+  readwrite
+  readwrite
+
+DISENTGL (WANNIER)
+
+(*) There are different ways to eliminate the screening channels in the case of an entangled band structure, i.e., if the
+correlated subspace is not made up of isolated bands. In this case, it is not possible to unambiguously define the screening
+that takes place only in the subspace. The default approach of Spex,
+described in [Phys. Rev. B 83, 121101 (2011)], scales each virtual transition :math:`\phi\rightarrow\phi'` by the probability that the
+transition takes place in the subspace. (To be more precise, it is the probability that the electronhole pair
+occupying the states :math:`\phi` and :math:`\phi'` happens to be in the subspace.) There is another method,
+called *disentanglement* method [Phys. Rev. 80, 155134 (2009)],
+in which the reference meanfield system is modified in such a way that each state :math:`\phi` is either fully contained in
+the subspace or completely outside (corresponding to probabilities 1 or 0). This method can be used by specifying ``DISENTGL``
+in the section ``WANNIER``.
+
+.. note:: Using this keyword actually modifies the meanfield system! (The keyword has not been tested thoroughly yet.)
++++
 Example  ``DISORDER 1000``  Use a finite value :math:`{\eta=1/(2\cdot 1000)}`. 
++++
diff git a/source/version b/source/version
new file mode 100644
index 0000000000000000000000000000000000000000..71d6fae8cd433e45de4668ab696995cfae0195c4
 /dev/null
+++ b/source/version
@@ 0,0 +1 @@
+05.00pre29