diff git a/.gitignore b/.gitignore
index 207410a41a7aaf45f8b5f63ef2bd0ad55e9a138c..01c3d69278d24be8005b2a4c2c42b10744717739 100644
 a/.gitignore
+++ b/.gitignore
@@ 1,4 +1,2 @@
build
_static
_templates
.gitlabci.yml
diff git a/source/_static/css/custom.css b/source/_static/css/custom.css
new file mode 100644
index 0000000000000000000000000000000000000000..8ccd9b889c1098c9974b36653c15faeeaa35b95c
 /dev/null
+++ b/source/_static/css/custom.css
@@ 0,0 +1,3 @@
+span.eqno {
+ float: right;
+}
diff git a/source/basis_sets.rst b/source/basis_sets.rst
index cee2a6fe32e292c1b95c1785ddb793c9d6540911..072ef08c8758bd8c7402bf41bc31f6196951a632 100644
 a/source/basis_sets.rst
+++ b/source/basis_sets.rst
@@ 33,10 +33,10 @@ 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. 
+++++
+++++
.. _lcut:
@@ 54,14 +54,17 @@ In the LAPW basis, the matching conditions at the MT boundaries require large :m
++++
 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}}`). 
+    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,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. 
+    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. 
++++
@@ 73,21 +76,25 @@ 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). 
+ Example  ``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.
+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``. 
+ 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 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 \mathrm{Bohr}^{1}`` (for testing only, can be quite slow). 
+  ``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``.
@@ 124,7 +131,7 @@ 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^'.
+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^'. 
diff git a/source/common_keywords.rst b/source/common_keywords.rst
index 8fd7af4d08be42ad028ad732dfe0d3bf1692a9be..77655f82512444ebdbbe5bb0da09d088ae260e75 100644
 a/source/common_keywords.rst
+++ b/source/common_keywords.rst
@@ 41,13 +41,14 @@ 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 
++++
+++++
+ 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

@@ 68,6 +69,7 @@ This is an important keyword. It enables (a) continuing a calculation that has,
++++
 Examples  ``RESTART``  Read/write restart file 
+++++
  ``RESTART 2``  Try to reuse data from standard output files 
++++
@@ 104,7 +106,7 @@ CHKMISM
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}`
+(*) 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 
@@ 116,19 +118,19 @@ BANDINFO

(*) 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`` 
++++
+++++
+ Example  ``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) 
++++
+++++
+ Example  ``DELTAEX 0.2eV``   Increase the exchange splitting by 0.2eV 
+    (spinup/down energies are decreased/increased by 0.1eV) 
+++++
PLUSSOC

@@ 138,17 +140,18 @@ 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 
++++
+++++
+ 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
====================
diff git a/source/conf.py b/source/conf.py
index 39caa8430bcdb73ba1d94fa77d7ebba5a634d120..d990991a1e6c7ceeb9d08de8078793b24831a5c5 100644
 a/source/conf.py
+++ b/source/conf.py
@@ 15,6 +15,8 @@
# import os
# import sys
# sys.path.insert(0, os.path.abspath('.'))
+def setup(app):
+ app.add_stylesheet('css/custom.css')
#  Project information 
diff git a/source/index.rst b/source/index.rst
index 4b2718c088da0749953c081e34ac17e238749125..eb2793c1c1978069aa48bd243e2b15fdfcd247b8 100644
 a/source/index.rst
+++ b/source/index.rst
@@ 13,7 +13,7 @@ It needs input from a converged DFT calculation, which can be generated by Fleur
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
diff git a/source/tutorials.rst b/source/tutorials.rst
index f367e9625e40438a0f629fd10ce460fd8a3b525d..9f39a0e2be378dca61dad070a0b0446aae911ce9 100644
 a/source/tutorials.rst
+++ b/source/tutorials.rst
@@ 57,17 +57,17 @@ The columns are
The two GW values for the quasiparticle energy follow two common methods to approximately solve the quasiparticle equation
.. _qpeq:

: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})}`
+.. 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})
+ :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
.. _qppert:

: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}`
+.. 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
+ :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.
+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:
@@ 77,11 +77,13 @@ SPECTRAL (SENERGY)
(*) It should be pointed out that the quasiparticle energies given in the output rely on the quasiparticle approximation. The more fundamental equation, as it were, is the Dyson equation
::math:`{\displaystyle G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)v^{\mathrm{xc}}]G(\omega)}`
+.. 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 {$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
+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\}\,,}`
+.. 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
@@ 91,25 +93,25 @@ in the section ``SENERGY`` of the input file. This option does not require the `
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. 
++++
+++++
+ 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. 
+++++
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.
The ``GW`` selfenergy
.. _selfene:

: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'\,.}`
+.. 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''\,.}``
+.. 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
@@ 117,13 +119,13 @@ This integral equation turns into a matrix equation
if the quantities :math:`{W}`, :math:`{\varepsilon}`, and :math:`{v}` ; are represented in the :ref:`mbp`, which thus has to be converged properly. The dielectric function, in turn, describes the change of the internal potential through screening processes and is related to the polarization matrix by :math:`{\varepsilon(\mathbf{k},\omega)=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 [[#Eq:Selfeneselfenergy]] with :math:{W(\omega)} replaced by :math:`{W^\mathrm{c}(\omega)=W(\omega)v}`. The frequency integration is carried out analytically for the exchange part (by summing over the residues). The correlation part is more complex to evaluate because of the frequency dependence of the interaction. (Fast electrons experience a different potential than slow electrons.)
+The 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 [[#Eq:Selfeneselfenergy]] 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, ``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).
+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 
@@ 140,19 +142,19 @@ 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.) 
++++
+++++
+ 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.
@@ 161,18 +163,20 @@ 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. 
++++
+++++
+ 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``.
@@ 206,9 +210,11 @@ ALIGNVXC (SENERGY)
++++
 Examples  ``ALIGNVXC``   Align exchangecorrelation potential in such a way that the 
    ionization potential remains unchanged by the quasiparticle correction. 
+    ionization potential remains unchanged by the quasiparticle 
+    correction. 
++++
  ``ALIGNVXC 0.2eV``  Apply a constant positive shift of 0.2 eV to the exchangecorrelation potential.
+  ``ALIGNVXC 0.2eV``   Apply a constant positive shift of 0.2 eV to the 
+    exchangecorrelation potential. 
++++
@@ 227,15 +233,15 @@ 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 
++++
+++++
+  ``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).
@@ 344,15 +350,17 @@ One can also go beyond the perturbative solution of the [[#Eq:qpeqquasiparticle
[[#QSGW]]If the job definition contains ``FULL`` and [[#IBZlabel``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]:
: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}`
+.. 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:`{\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}`
+.. 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 any previous version of the *hermitianized* matrix.
* ``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.
@@ 389,7 +397,9 @@ 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
: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'\,.}` [[#Eq:P]]
+.. 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'\,.
+ :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.)
@@ 425,12 +435,13 @@ The most important keyword in the calculation of the polarization function is ``
    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.) 
+    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:P``P``]] approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in ``'k``' becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}\epsilon_{\mathbf{q}n})^{1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and nondivergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the [[#Eq:Pfrequency integration]], thereby avoiding the numerical difficulties. (As an alternative, the energy differences can also be incorporated into the integration weights, which is arguably even more stable numerically, see option ``INT`` below.) By default, Spex does that only for k=0. The behavior can be changed with the keyword ``MULTDIFF`` in the section ``SUSCEP``.
+(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of [[#Eq:P``P``]] approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n}  \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in :math:`k` becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}\epsilon_{\mathbf{q}n})^{1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and nondivergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the [[#Eq:Pfrequency 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. 
@@ 440,7 +451,8 @@ MULTDIFF (SUSCEP)
  ``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. 
+  ``MULTDIFF INTON``   Always separate energy difference by sticking them into 
+    integration weights. 
++++
PLASMA (SUSCEP)