diff --git a/source/basis_sets.rst b/source/basis_sets.rst
index d9860ff608ac331a8be2561eb5f4b896ba2f73ed..2c76a9f5db5984054d3f38786511e50fa599b91c 100644
--- a/source/basis_sets.rst
+++ b/source/basis_sets.rst
@@ -4,4 +4,148 @@
Basis sets
===================
-TODO
+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 non-overlapping 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 high-lying 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.
+
+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 spin-polarized calculations. The files contain sets of data for all l quantum numbers, which can be plotted separately with gnuplot (e.g., ``plot "spex.mt.1" i 3`` for :math:`{l=3}`).
+
++----------+------------------+--------------------------------------+
+| Examples | ``MTACCUR -1 2`` | | Calculate MT representation error |
+| | | | between -1 and 2 Hartree. |
++----------+------------------+--------------------------------------+
+| | ``MTACCUR`` | Automatic energy range. |
++----------+------------------+--------------------------------------+
+
+
+Mixed product basis [[#MPB]]
+============================
+
+The methods implemented in Spex distinguish themselves from conventional DFT with local or semilocal functionals by the fact that individual two-particle 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.
+
+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 (MBASIS) [[#LCUT]]
+-----------------------
+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=|l-l'|,...,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. |
++---------+--------------+------------------------------------------------------------+
+
+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}` (so-called :math:`{u}`) and no function of :math:`{p=1}` (so-called :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 wave-function products that must be represented consist mostly of occupied-unoccupied 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.
+
+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). |
++---------+-----------------+-------------------------------------------------------------------------------+
+
+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 basis-set size further by introducing a Coulomb cutoff :math:`{v_\mathrm{min}}`. The reduced basis set is then used for the polarization function, the dielectric function, and the screened interaction. The parameter :math:`{v_\mathrm{min}}` can be specified after the keyword ``OPTIMIZE MB`` in three ways: first, as a "pseudo" reciprocal cutoff radius :math:`{\sqrt{4\pi/v_\mathrm{min}}}` (which derives from the plane-wave Coulomb eigenvalues :math:`{v_\mathbf{G}=4\pi/G^2}), second, directly as the parameter :math:`{v_\mathrm{min}}` by using a negative real number, and, finally, as the number of basis functions that should be retained when given as an integer. The so-defined basis functions are mathematically close to plane waves. For testing purposes, one can also enforce the usage of plane waves (or rather projections onto plane waves) with the keyword ``OPTIMIZE PW``, in which case the Coulomb matrix is known analytically. No optimization of the basis is applied, if ``OPTIMIZE`` is omitted.
+
++----------+-----------------------+-------------------------------------------------------------------------------------------------------------+
+| 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 ``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``.
+
+Finally, we show a section ``MBASIS`` for a system with three atom types as an example
+
+.. code-block:: bash
+
+ SECTION MBASIS
+ GCUT 3.0
+ LCUT 5 4 4
+ SELECT 3,2;4;2 2;3 2;3
+ TOL 0.0001
+ OPTIMIZE MB 4.5
+ END
+
+NOAPW (MBASIS)
+--------------
+(*) By default, Spex constructs augmented plane waves from the mixed product basis in a similar way as in the LAPW approach. In this way, the unphysical "step" at the MT sphere boundaries is avoided. It also leads to a further reduction of the basis set without compromising accuracy. This APW construction can be avoided with the keyword ``NOAPW``.
+
+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.
+
+CHKPROD (MBASIS)
+-----------------
+(*) Checks the accuracy of the MT product basis.
+
+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.)
+
+FFT (WFPROD)
+--------------
+When the interaction potential is represented in the mixed product basis, the coupling to the single-particle 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 real-valued 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}}`. |
++----------+---------------+-------------------------------------------------------------------------------+
+
+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.)
+
+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 [[#LCUT|``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 wave-function coefficients are neglected. By default, the two options are unused.
+
+MINCPW (MBASIS)
+----------------
+(**) The keyword ``MINCPW`` modifies the plane-wave wave-function coefficients such that their absolute values become as small as possible. This reduces the error in the evaluation of the wave-function 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 orbitals [[#Wannier]]
+=============================
+
+...
+
diff --git a/source/keywords.rst b/source/common_keywords.rst
similarity index 99%
rename from source/keywords.rst
rename to source/common_keywords.rst
index 20d08d330a77914c28be8e2f6c40fd3d8f0a728c..d920a2389ad25e0bc344bf8407d26382434572f5 100644
--- a/source/keywords.rst
+++ b/source/common_keywords.rst
@@ -1,4 +1,4 @@
-.. _keywords:
+.. _common_keywords:
****************
Common Keywords
diff --git a/source/index.rst b/source/index.rst
index e41ac41ad9e26b60da682aa2396490a7f382fc26..2659b9c6ccf5bb5fbe19de264953c3b278b2c37f 100644
--- a/source/index.rst
+++ b/source/index.rst
@@ -19,12 +19,14 @@ User's guide
++++++++++++
.. toctree::
- :maxdepth: 3
+ :maxdepth: 2
installation
getting_started
general_information
- keywords
+ common_keywords
tutorials
basis_sets
+ keyword_reference
+ spex_faq
diff --git a/source/keyword_reference.rst b/source/keyword_reference.rst
new file mode 100644
index 0000000000000000000000000000000000000000..b8b76a4e385e1f89102ce7c7ebd691965d3a6e5e
--- /dev/null
+++ b/source/keyword_reference.rst
@@ -0,0 +1,174 @@
+.. _keyword_reference:
+
+========================
+Keyword Reference
+========================
+This document gives a reference of available keywords in the input file ``spex.inp``. The syntax of the input file is as follows.
+
+* 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 and may be given only once.
+* 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.
+* Some keywords are only defined inside sections. A section starts with ``SECTION sectionname`` and ends with ``END``.
+
+In the following the keywords are listed and explained alphabetically and for each section separately. If there is no default value, the keyword must be given. **Experimental, badly tested, and obsolete keywords as well as keywords for testing are shown in round brackets.**
+
+.. warning:: The list can be partly out of date or obsolete.
+
+Global Keywords
+-----------------
+
+(``ALIGNBD``): Aligns the degenerate states on the unshifted k-point set (see ``KPT +=...``) to the states on the shifted k points. Corresponding transformation matrices are written to the file ``spex.abd`` or read from it if it exists. Normally used to calculate response quantities of metals for small k vectors, e.g., in order to test the expansion around the point '''k=0'''. (Default false)
+
+``BANDINFO``: Prints information about bands: weight of s,p,d,f,g character of each state, (P)DOS plot in separate file (``spex.dos.NNN``), square integral of wave-function products. (Default false)
+
+(``BANDOMIT (bands) ABC``): Omits bands defined in ``bands`` in exchange (``A=1``), susceptibility (``B=1``), and/or correlation self-energy (``C=1``). (Otherwise ``A=0`` etc.) (Default none)
+
+``BZ k l m``: Defines a k-point set ``k``x``l``x``m``.
+
+``CHKMISM``: Checks the mismatch of the DFT wave functions at the MT spheres. (Default false)
+
+``DIPOLE``: Calculates dipole moments for all direct (optical) transitions. (Default false)
+
+``CHKOLAP``: Checks the overlap of the DFT wave functions. (Default false)
+
+(``CORES cores``): Includes core states into the calculation of the susceptibility and self-energy. The core states are defined separately for each atom type as a list in brackets, e.g., ``(4f,5s,5p)``. Empty lists ``()`` are allowed. (Default none)
+
+``CORESOC``: Treats core states as { j,m_j } Dirac spinors. Otherwise, the Dirac states are averaged in each { l,m } channel. (Default false)
+
+(``GAUSS x y``): Use Gauss integration with widths ``x`` and ``y`` instead of tetrahedron integration. Only for debugging. (Default false)
+
+(``IBC n``): Uses the incomplete-basis-set correction (IBC). Experimental! (Default false)
+
+(``ITERATE opt e``): Constructs a (``opt=NR``: non-relativistic, ``opt=SR``: scalar-relativistic, ``opt=FR``: fully relativistic, i.e., including SOC) LAPW(+LO) Hamiltonian from the converged DFT potential and calculates wave functions and energies from it, which are then used in the calculation. This spares you the last DFT run. The parameter ``e`` is the lower energy bound for the electron states. The calculation is fairly slow. (Default none)
+
+``JOB job1 job2 ...``: Defines what Spex should do. Consult the [[Spex_Tutorial|tutorial]] for details.
+
+``KPT label1 label2 ...``: Defines labels (single characters) for k points. Consult the [[Spex_Tutorial|tutorial]] for details. (Default none)
+
+``KPTPATH (A,B,C,...)``: Comma-separated list of k points that defines a k-point path. Spex automatically determines all k points (according to the ``BZ`` definition) that lie on the path. The path can be used (e.g., in the ``JOB`` definition) with ``PATH``. If Wannier functions are defined, Spex performs a Wannier interpolation along this path. (default none)
+
+``MEM m``: Restricts the memory usage to ``m`` MB (per node). (Does not work exactly. The code actually needs more memory.) (Default 200)
+
+``MTACCUR emin emax``: Writes plots for the MT representation error in the energy range {``emin``,``emax``} to the files ``spex.mt.NNN``. (Default none)
+
+``NBAND n``: Restricts the number of bands to ``n`` bands or, if ``n`` is a real number, to all bands with energies below that value. (:comment All bands that are not treated explicitly can be taken into account in an approximate way by placing them at a fixed energy above all other energies. For example with ``NBAND 100+5`` this energy would lie 5 ha above the maximum energy (see PRB 78, 085125). :) (Default all bands)
+
+(``PLUSSOC``): Adds spin-orbit coupling to the input data in second variation. (Default false)
+
+``RESTART``: Enables to continue an older calculation (e.g., if the latter has crashed). Intermediate results are written to files (``spex.mb`` - Mixed product basis, ``spex.cou(s)`` - Coulomb matrix, ``spex.cor`` - Dielectric function, ``spex.spc`` - (e.g.) polarization function, ``spex.uwan`` - Wannier U matrices, ``spex.kolap`` - Overlap matrices for Wannier construction) or read from the files if they exist. With ``RESTART SIG``, Spex reads the self-energy from the files ``spex.sigx`` and ``spex.sigc``, which are always written in ``GW/HF FULL`` calculations. (Default false)
+
+``STOREIBZ``: Wave functions are only stored in the IBZ rather than in the full BZ. This decreases the memory demand considerably but slows down the calculation somewhat because wave functions must be regenerated when needed. (Default false)
+
+``WRITE q1 q2 ...``: Tells Spex to write the head element of the specified quantities ``q1/2=`` ``SUSCEP``, ``DIELEC`` or ``SCREEN`` to a file. Optionally a job number can be specified in brackets, as in ``WRITE SUSCEP(2)``. If one of these quantities are explicitly defined in the job definition, they are automatically written to files. With ``q1=INFO``, Spex writes preliminary results to several output files. Consult the [[Spex_Tutorial|tutorial]] for details. (Default none)
+
+``WRTKPT``: Tells Spex to write out the k-point set and exit. (Default false)
+
+Section ``MBASIS``
+-------------------
+
+``ADDBAS``: Augments mixed basis by the radial {u} functions (according to first part of SELECT entry). This improves the basis convergence for COHSEX calculations. (Default false)
+
+``CHKPROD``: Checks the accuracy of the MT product basis. (Default false)
+
+``GCUT G``: G cutoff for the IPWs of the mixed basis.
+
+``LCUT l``: l cutoffs (one for each atom type) for the mixed basis.
+
+``NOAPW``: Switches off construction of "product APWs", which are continuous at the MT sphere boundary. (Default false)
+
+``OPTIMIZE``: Optimize mixed basis for correlation self-energy. Consult the [[Spex_Tutorial|tutorial]] for details. (Default none)
+
+``SELECT l1,l3,lo1;l2,l4,lo2``: Selection of the MT radial functions from which the mixed basis is constructed. Consult the [[Spex_Tutorial|tutorial]] for details.
+
+``TOL t``: Tolerance value for the removal of linear dependencies from the MT part of the mixed basis. (Default 0.0001)
+
+(``WFADJUST``): 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.) (Default false)
+
+Section ``COULOMB``
+--------------------
+
+``CHKCOUL``: Performs checks on the Coulomb matrix: (1) Project on plane waves and see whether we get :math:`{v_{\mathbf{GG'}}=4\pi \delta_{\mathbf{GG'}}/ |\mathbf{k+G}|^2 }`, (2) compare largest eigenvalue of Coulomb matrix to :math:`{4\pi / k^2}`. (Default false)
+
+(``CUTZERO``): Removes linear dependencies in the overlap matrix of the mixed basis. Otherwise, in the case of a large ``GCUT`` the diagonalization of the Coulomb matrix might fail due to numerical rounding errors in LAPACK. (Default false)
+
+``LEXP``: l cutoff for the expansion of plane waves inside MT spheres. (Default 14)
+
+``NOSTORE``: The default is to precalculate and keep the Coulomb matrices for all k points in memory. With this keyword the matrices are calculated whenever they are needed. This reduces the memory demand but somewhat slows down the calculation. Unless ``APPROXPW`` is chosen, this also applies to the overlap-correction matrix. (Only affects HF and QP calculations.) (Default false)
+
+(``STEPRAD K``): Uses the Fourier transform of the step function as an alternative to the Rayleigh expansion for representing the IPWs and prints the resulting difference. The parameter ``K`` is the corresponding reciprocal cutoff radius. Only for testing. (Default none)
+
+(``TSTCOUL``): Similar to ``CHKCOUL`` but regenerates MT part of the mixed basis from spherical Bessel functions and stops after the test. (Default false)
+
+Section ``WFPROD``
+-------------------
+
+(``APPROXPW``): Since version 02.03 the plane-wave part of the wave-function products are calculated exactly leading to much better 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. (Default false)
+
+``LCUT l``:Optional l cutoffs (one for each atom type) for the wave functions in the calculation of the products. (Default none)
+
+``MTTHR d``: Threshold value below which wave-function coefficients are neglected. (1E-10 is a safe choice.) (Default 0)
+
+(``MINCPW n``): Modifies the IPW wave-function coefficients such that their absolute values become as small as possible. This reduces the error in the evaluation of the wave-function 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``. (Default 0)
+
+``FFT gcut``: Employs Fast Fourier transformation for the calculation of the plane-wave part of the wave-function products. Scales much better with system size than the direct (default) evaluation. Requires a reciprocal cutoff ``gcut``. The exact limit is :math:`gcut=GCUT{+2k_{\mathrm{max}}}`, where :math:`{k_{\mathrm{max}}}` is the PW cutoff of the DFT calculation, but much lower cutoffs are enough; try with ``gcut=8`` for example. If no ``gcut`` is given, Spex uses the exact limit. (Default none)
+
+
+Section ``SUSCEP``
+------------------
+
+``DISORDER`` :math:`{\tau}`: Replaces the infinitesimal :math:`{i\eta}` by the finite frequency :math:`{i/2\tau}` in the response function. (Default none)
+
+(``FSPEC N x``): Deprecated keyword. Use ``HILBERT`` instead. (Default none)
+
+``HILBERT d x``: Defines the exponential frequency mesh for the Hilbert transformation in one of two ways: (1) ``d`` is the difference between the first (set at the band gap value or zero) mesh point :math:`{\omega_1}` and the second mesh point :math:`{\omega_2}`, and :math:`x` is the factor :math:`{x=(\omega_{i+1}-\omega_i)/(\omega_i-\omega_{i-1})}`; (2) if ``d`` is an integer, it is interpreted as the number of mesh points in the range :math:`{{\omega_1}..5 htr}`, and ``x`` is interpreted as the stretching factor "acquired at 5 htr", i.e., :math:`{x^{1/N}}` is the corresponding factor defined in (1). (Default 0.01 1)
+:comment Alternatively, the first argument ``N`` can be a real number in which case it is interpreted as the first nonzero mesh point. Latest version: A negative second argument means the stretching factor that is reached after {N} points, i.e., :math:`{x^N}`. (Default 30 1.05):)
+
+``HUBBARD file``: Calculates Hubbard parameters. If the optional filename ``file`` is given, the Hubbard parameters are read from (written to) ``file`` if it exists (does not exist). (Default false)
+
+``MULTDIFF ON`` or ``OFF``: If set "ON" (and ``MB/PW=`` is used), the BZ integrand in the calculation of the head and wing elements of the polarization function (not for the magnetic susceptibility) is multiplied by KS energy differences. (The integration weights are correspondingly adjusted.) For a small (or zero) '''k''' vector this leads to better numerics (because the integrand is exactly inversely proportional to the KS energy differences in the case '''k=0'''). (Default ``ON `` for '''k=0''', ``OFF`` else)
+
+``PLASMA p``: Sets the plasma frequency at ``p``. ``p=METAL``: use metallic screening. (Default none)
+
+(``TETRAF``): If set, tetrahedron weights are only calculated in the (extended) irreducible part of the BZ. This accelerates the calculation of the weights considerably but introduces slight deviations due to symmetry breaking. (Default false)
+
+``WGHTTHR w``: Sets a threshold value for the (tetrahedron) weights. All transition with weights below ``w`` are neglected. (Default :math:`{10^{-10}}`)
+
+Section ``SENERGY``
+----------------------
+
+(``ALIGNVXC``): Aligns chemical potential to quasiparticle state of highest occupied state. (Default false)
+
+``CONTINUE (n)``: If the optional (and now somewhat obsolete) parameter ``n`` is not given, a pade continuation is used. Otherwise, ``n`` specifies the number of poles used in the fit function for the analytic continuation of the self-energy. The fitting procedure can be influenced further by adding single characters to ``n``; character "``c``": employ additional constraints (i.e., integrability of the self-energy and continuity of its derivative at {\omega=0}), "``+``": add a constant to the fit function as a parameter, "``*``": allow poles above the real frequency axis. (Default Pade)
+
+``CONTOUR D d``: Specifies contour integration for the evaluation of the self-energy. The first argument ``D``, if it is a number, defines the interval for the numerical derivative :math:`{d/d\epsilon \Sigma(\epsilon) = [\Sigma(\epsilon+D)-\Sigma(\epsilon-D)]/(2D)}` where :math:`{\epsilon}` is the corresponding KS energy. It can also be defined as a range, e.g. ``{-0.5..0.5,0.01}`` or ``[{-0.5..0.5,0.01}]``, which defines an equidistant mesh (relative to the Kohn-Sham energy or in absolute energies, the latter is the only viable definition for a ``GW FULL`` calculation) on which the self-energy is evaluated and which is used to solve the nonlinear quasiparticle equation (neglecting off-diagonal elements). The second argument ``d`` defines the step size of an equidistant mesh for the screened interaction. If ``d`` is omitted, a Pade continuation is used for the screened interaction (Keywords ``CONTINUE`` and ``CONTOUR`` are mutually exclusive.) (Default none)
+
+``FREQINT SPLINE/PADE``: Algorithm for frequency integration in :math:`{\int G(\omega+\omega')W(\omega')d\omega'}`: the screened interaction is interpolated with a splines or with a Pade approximation. (Default SPLINE)
+
+``LFRAC l``: l cutoff for the treatment of the anisotropy at the point :math:`{k=0}`. For isotropic systems ``l=1`` is enough. (Default 1)
+
+``MESH n x``: Defines the imaginary frequency mesh where ``n`` is the number of points and ``x`` the frequency of the last point. Instead, one can also specify a file that contains the frequency mesh as a list with ``MESH filename``. (Default 10 10.0)
+
+``NOZERO``: Neglect zero order terms in the self-energy coming from combinations of divergent terms math:`{1/k}` and math:`{1/k^2}` in the (bare or screened) interaction with linear and quadratic terms of the other quantities. This is necessary to converge the k-point sampling for materials with small band gaps. (Default false)
+
+``RPAENE``: Calculate RPA correlation energy "on the fly". (Default false)
+
+``VXC READ/CALC``: Specifies whether the vxc matrix elements are calculated or read from the file ``vxcfull``. (Default READ)
+
+.. note:: ``IMAGTIME``: Since version 02.00 the correlation self-energy is evaluated by a direct frequency convolution of {G} and {W} rather than a multiplication in imaginary time. With the keyword ``IMAGTIME`` the latter (older) algorithm can still be used. (In some cases the old algorithm is considerably faster.) (Default false) :)
+
+Section ``WANNIER``
+--------------------
+
+``FROZEN e``: Use frozen energy window. Lower bound is defined automatically, upper bound can be defined with the optional parameter ``e``. If ``e`` is not given, Spex guesses a reasonable value. (Default false)
+
+``MAXIMIZE``: Use maximally localized Wannier functions instead of the first-guess ones. (Default false)
+
+``ORBITALS (..) (..) ..``: Defines the orbitals for the first guess. Round brackets ``()`` for atoms, square brackets ``[]`` for atom types. In the brackets, one can use ;-separated lists of s,p,d,f,px,py,pz,dxy,dyz,dxz,dz2,dx2y2,t2g,eg,fz3,fxz2,fyz2,fzxy,fxyz,fxxy,fyxy,sp,sp2,sp3,sp3d,sp3d2. Euler angles can be specified with, e.g., ``sp,120,0,0``. (Default none)
+
+``RSITE latvec``: Calculate off-site Hubbard interaction parameters (in addition to the on-site ones) with ``latvec`` being a lattice vector, e.g., ``(0,0,1)`` or several lattice vectors, e.g., ``(0,0,1) (0,0,3)`` or ``(0,0,1)-(0,0,3)`` or ``(0,1..2,1..3)``. The latter two define a list of lattice vectors: (0,0,1), (0,0,2), (0,0,3) and (0,1,1), (0,2,1), (0,1,2), (0,2,2), (0,1,3), (0,2,3), respectively. (Default none)
+
+(``UREAD``): Wannier U matrix elements are read from the Fleur output file. (Default false)
+
+(``WSCALE s``): Scaling factor ``s`` for the screened Coulomb interaction for spin-wave calculations. ``s`` can be a range ...... to be completed! (Default none)
diff --git a/source/tutorials.rst b/source/tutorials.rst
index 8605a6e0fe09729626cc8bf101e8015b6157478e..6835dc233cfd64f2d0886cc83349f0d753a847c4 100644
--- a/source/tutorials.rst
+++ b/source/tutorials.rst
@@ -188,7 +188,7 @@ SMOOTH (SENERGY)
NOZERO (SENERGY)
----------------
-Both the exchange (HF) and the correlation self-energy 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 [[#MPB|"Mixed Product Basis"]] and behave as :math:`{\propto k}` in the long-wavelength limit :math:`{k\rightarrow 0}`. The prefactor is obtained from :math:`{k\middot;p}` perturbation theory. Multiplying with :math:`{I\propto k^{-2}}` gives rise to zero-order terms involving that prefactor. These zero-order terms (which only appear at the :math:`{\Gamma}`; point) can improve the k-point convergence. However, in the case of systems with small band gaps, the zero-order 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 k-point sampling. Therefore, the keyword ``NOZERO`` in section ``SENERGY`` allows the zero-order terms to be neglected. They are included by default.
+Both the exchange (HF) and the correlation self-energy 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 [[#MPB|"Mixed Product Basis"]] and behave as :math:`{\propto k}` in the long-wavelength 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 zero-order terms involving that prefactor. These zero-order terms (which only appear at the :math:`{\Gamma}`; point) can improve the k-point convergence. However, in the case of systems with small band gaps, the zero-order 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 k-point sampling. Therefore, the keyword ``NOZERO`` in section ``SENERGY`` allows the zero-order terms to be neglected. They are included by default.
ALIGNVXC (SENERGY)
------------------
@@ -215,7 +215,7 @@ VXC (SENERGY)
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 self-energy matrix (or expectation values) :math:`{\Sigma^{xc}}` with the contribution of the current k point (and its symmetry-equivalent k points). After completion of step (b), the current self-energy 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 self-energy 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 self-energy calculation. The following logical table gives an overview.
+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 self-energy matrix (or expectation values) :math:`{\Sigma^{xc}}` with the contribution of the current k point (and its symmetry-equivalent k points). After completion of step (b), the current self-energy 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 self-energy 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 self-energy calculation. The following logical table gives an overview.
+--+---------------+-----------------------+
| | "spex.cor" | "spex.sigx/c" |
@@ -258,14 +258,194 @@ In the case of ``RESTART 1`` (or ``RESTART``) (i.e, when "spex.cor" is to be rea
GW band structures
-------------------
+In this section, three ways of calculating ``GW`` band structures are discussed.
+
+* Using the [[#KPTPATH|``PATH``]] label (single Spex run)
+* With [[#KPTADD|additional ``'q``']] points (multiple Spex run)
+* [[#WanInt|Wannier interpolation]] (single --; but long --; Spex run).
+
+The calculation of a ``GW`` band structure is not as straightforward as in the case of KS-DFT, 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 q-point path. The ``GW`` self-energy, on the other hand, is non-local. 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 high-symmetry line (e.g., Γ–X) 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 high-symmetry 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
+--------
+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 k-point indices along this path and set the job definition accordingly, but Spex can do this for you. The k-point path can be defined with a line ``KPTPATH (L,1,X)``. In the [[#JOB|job definition]], one can then use the special label ``PATH`` to address the corresponding list of k points, e.g., ``JOB GW PATH:(1-10)``. 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 [[#spex.out|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 k-point sets, but this would entail very expensive calculations.
+
++-----------+------------------------+------------------------------------------------------------+
+| Examples | ``KPTPATH (L,1,X)`` | | Define k-point 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:(1-10)`` | | Run ``GW`` calculation for the first ten bands at all |
+| | | | k points defined by ``KPTPATH``. |
++-----------+------------------------+------------------------------------------------------------+
+
+
+KPT +=...
+---------
+This section explains an extension to the keyword [[#KPT|``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 k-point set (defined by ``BZ``). For example, in Si the conduction-band minimum is not at a high-symmetry 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 [[#KPT|``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 k-point set ``'q+k``'. This is done in the usual way by creating the k-point file ([[#WRTKPT|``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:(1-5) +:(1-5)``, which will yield the bands 1-5 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:`{\Gamma--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 +:(1-10)`` | | Run ``GW`` calculation for the states 1-10 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 +:(1-10)``. 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|``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|"spex.selfc"]]) and produces, for each q point, one output file named ``spex_???.out`` (where ``???`` is a three-digit counting index, or more digits if ``???>999``).
+
+From these files the band-structure 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 band-structure 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.
+
+[[#WanInt]]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 [[#Wannier|separate section]]. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to ``spex.inp``:
+
+.. code-block:: bash
+
+ JOB GW IBZ:(1-4)
+ SECTION WANNIER
+ ORBITALS 1 4 (sp3)
+ 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 back-and-forth Fourier transformation with a real-space truncation of matrix elements in-between. 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 [[#SPECTRAL|``SPECTRAL``]]. The interpolated spectral function is written to the file "spectralw" or, if the system is spin-dependent, 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 [[#Eq:qpeq|quasiparticle equation]] and solve it by iterative diagonalization in the basis of the single-particle states. This requires the full self-energy matrix including off-diagonal elements to be calculated. The respective line in the input file would read
+``JOB GW FULL 1:(1-10) X:(1-10) L:(1-10)`` giving the self-energy as 10x10 matrices at the Γ, X, and L points. (In principle, it is also possible to have a band definition like ``(1-5,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::
+
+ Bd KS olap HF olap GW
+ 3 2.92306 0.99971 0.12683 0.99985 2.50496 -0.06006
+ 4 2.92306 0.99971 0.12683 0.99985 2.50496 -0.06006
+ 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 Hartree-Fock 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 one-to-one 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 non-linear 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 off-diagonal 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 well-defined 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 high-lying states, which do not have a well-defined quasiparticle character anymore.
+
+[[#QSGW]]If the job definition contains ``FULL`` and [[#IBZlabel|``IBZ``]], the full ``GW`` self-energy matrix is evaluated for the whole IBZ, which enables self-consistent calculations in the framework of the ``quasiparticle self-consistent GW`` (QSGW) approach. In this approach, one creates a mean-field system from the ``GW`` self-energy whose single-particle energies are as close as possible to the quasiparticle energies. This mean-field system is subsequently solved to self-consistency in a DFT code. The resulting solution can then serve as a starting point for a new one-shot ``GW`` calculation, which constitutes the second iteration and so on until the quasiparticle energies are converged. The construction of the mean-field system is, to some extent, arbitrary. We use the following definition, which is slightly modified from the original work [PRL 93, 126406]:
+
+: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}`
+
+for off-diagonal elements. The ``hermitianized`` QSGW operator is then obtained from :math:`{\Sigma^\mathrm{xc,H}=(A+A^\dagger)/2}`. The difference to the original definition is the inclusion of the renormalization factor to better reproduce the ``GW`` quasiparticle energies. The ``hermitianized`` matrix, or rather the difference :math:`{\Sigma^{\mathrm{xc,H}}-v^\mathrm{xc}}`, is written to the file "spex.qsgw", which is later read by the DFT code. In Fleur, the following steps are required:
+
+* ``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.
+ 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 error-prone. [[#selfc.sh]]Therefore, we provide the shell script ``spex.selfc``, which performs the steps needed for a self-consistent 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,
+* ``SPEX_QUEUE`` - Queue command for parallel execution of Spex (e.g., ``mpiexec -np 20``),
+* ``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 self-consistent 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 self-consistency at the moment.
+
GW with MPI
-------------
+RESTART...
+----------
+
+MPIBLK
+---------
+
+MPISYM
+--------
Polarization function
=====================
+The polarization function gives the linear change in the electronic density of a non-interacting system with respect to changes in the effective potential. It is, thus, a fundamental quantity in the calculation of screening properties of a many-electron systems. For example, the dielectric function, instrumental in the calculation of spectroscopic quantities (e.g. EELS) and the screened interaction needed in ``GW``, is related to the polarization matrix through :math:`{\varepsilon(\mathbf{k},\omega)=1-P(\mathbf{k},\omega)v(\mathbf{k})}`, here given in matrix notation. The corresponding explicit formula for matrix elements of ``P`` in the mixed product basis is
+
+:math:`{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]]
+
+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.)
+
+The expression for [[#Eq:P|P]] involves an infinite band summation over :math:`{n'}`. In practice, this summation is truncated by the number of bands (keyword [[#NBAND|``NBAND``]]), which thus becomes an important convergence parameter. The self-energy 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 k-point mesh (see ``HILBERT NONE``) and the Gaussian method (see ``GAUSS``).
+
+TETRAF (SUSCEP)
+-----------------
+(*) In rare cases, especially for very large k-point 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 (occ-nocc) 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 Hartree-Fock exchange potential and the determination of the Fermi energy). Therefore, ``GAUSS`` is a global keyword. It effectively replaces each single-particle 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 k-integration 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 [[#Eq:P|``P``]]. (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_{i-1})}`. These two parameters are accepted as real-valued 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 real-valued, 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 band-gap 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 three-digit 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 real-valued.) |
++----------+------------------------+----------------------------------------------------------------------------------------+
+
+MULTDIFF (SUSCEP)
+-------------------
+(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of [[#Eq:P|``P``]] approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in ``'k``' becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}-\epsilon_{\mathbf{q}n})^{-1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and non-divergent values, and takes the energy difference into account by replacing :math:`{S(\omega)\rightarrow S(\omega)/\omega}` in the [[#Eq:P|frequency integration]], thereby avoiding the numerical difficulties. (As an alternative, the energy differences can also be incorporated into the integration weights, which is arguably even more stable numerically, see option ``INT`` below.) By default, Spex does that only for k=0. The behavior can be changed with the keyword ``MULTDIFF`` in the section ``SUSCEP``.
+
++-------------------+--------------------+------------------------------------------------------------------------------+
+| 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. |
++-------------------+--------------------+------------------------------------------------------------------------------+
+
+PLASMA (SUSCEP)
+---------------
+(*) In the case of metals, the polarization function contains an additional term, the so-called Drude term, which gives rise to metallic screening (diverging static dielectric constant, short-range 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`` self-energy, impeding a straight-forward numerical solution of the non-linear 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 k-point 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``. |
++----------+------------------+----------------------------------------------------+
+
+DISORDER (SUSCEP)
+-------------------
+(*) Mathematically, the parameter :math:`{\eta}` in [[#Eq:P|''P'']] 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.
++---------+-------------------+----------------------------------------------------+
+| Example | ``DISORDER 1000`` | Use a finite value :math:`{\eta=1/(2\cdot 1000)}`. |
++---------+-------------------+----------------------------------------------------+