diff --git a/source/specifics.rst b/source/specifics.rst
new file mode 100644
index 0000000000000000000000000000000000000000..d52bf83269e27e6005197dfb0bf1f7bc65e44d00
--- /dev/null
+++ b/source/specifics.rst
@@ -0,0 +1,146 @@
+.. _specifics:
+
+===================
+Specifics
+===================
+
+Spex and Fleur
+---------------
+In the Fleur [[input]] file there is a flag for writing out the necessary data for a Spex calculation. It is called ``gw`` and is appended to line 23 of the Fleur ``inp`` file.
+[`23|vchk=F,cdinf=F,pot8=F,gw=n,numbands=N`]
+``n`` can take four values:
+
+* ``gw=0``: No output files are generated (default; the same as omitting ``gw`` and ``numbands`` altogether).
+* ``gw=1``: Some basic parameters are written to several files. Otherwise Fleur runs as usual.
+* ``gw=2``: All output files necessary for Spex are created and Fleur stops after one iteration (without updating the potential). In this case you should also set ``pot8=T`` and the maximal number of bands ``N``. The latter overrides the default ``neigd`` value in [[fl7para]].
+* ``gw=3``: Self-consistent cycle for QSGW. Same as ``gw=1`` but adds {$ \Sigma^{\mathrm{QSGW}}_{\mathrm{xc}}(\mathbf{r},\mathbf{r}')-v_{\mathrm{xc}}(\mathbf{r}) $} to the Hamiltonian in each iteration.
+
+The output files are
+
+* ``gwa``: Basic parameters including the atomic numbers and positions in the unit cell, lattice parameters and basis vectors, FLAPW l cutoff, and local-orbital parameters,
+* ``LATTC``: FLAPW G cutoff,
+* ``radfun``: radial basis functions,
+* ``ecore``: core-electron functions,
+* ``eig``: k points, wave-function coeffients (interstitial), energies,
+* ``abcoeff``: wave-function coefficients (muffin tins),
+* ``vxc``: expectation values of the xc potential (diagonal elements),
+* ``vxcfull``: matrix of the xc potential,
+* ``qsgw``: matrix of the QSGW self-energy (only for QSGW calculations).
+
+Furthermore Spex needs the Fleur files [[fl7para]] and ``sym.out``.
+
+Some comments about the DFT calculation
+---------------------------------------
+
+Accurate description of conduction states
+"""""""""""""""""""""""""""""""""""""""""
+As in a ''GW'' calculation the conduction states enter both the Green function and the screened interaction, they should be accurately described already in the KS system which is our starting point for the perturbation treatment. The LAPW basis, however, only guarantees well described occupied states. In order to obtain well-converged ''GW'' results, one must make the basis set more flexible, especially in the MT regions, by introducing local orbitals either at high energy parameters or as higher-order energy derivatives [see C. Friedrich et al., Phys. Rev. B 74, 045104 (2006)].
+
+Role of semicore states
+""""""""""""""""""""""""
+High-lying semicore states can appear as badly described "ghost bands" in the valence band region leading to a wrong DFT ground state. A possible solution is to use local orbitals and extend the energy window, such that the semicore states are treated as valence states. On the other hand if the ghost bands are above the Fermi energy, they usually pose no problem in a DFT calculation. In a ''GW'' calculation, however, the conduction states are also relevant, and one must make sure that there are no ghost bands at all in the energy spectrum. Therefore, it might be necessary to treat more (and deeper) semicore states as valence states than are needed in the DFT calculation.
+As an indication for such a case the overlap of core states with basis and wave functions calculated in the routine "checkinput" should be checked. As an example, here is the output for Strontium Titanate with a Titanium 3s ghost band::
+
+ Overlap
+ Atom type 1
+ u(s) udot(s)
+ 1s -0.001431 0.000644
+ 2s -0.000830 0.000414
+ 3s 0.287789 -0.878671
+ u(p) udot(p) ulo(p) ...
+ 2p -0.000665 0.000112 -0.000494
+ ...
+ Maximum overlap at (band/kpoint)
+ Atom type 1
+ 1s 0.000611 (031/ 001)
+ 2s 0.000343 (031/ 001)
+ 3s 0.664265 (047/ 001)
+ 2p 0.000340 (003/ 036) 0.000481 (003/ 042) 0.000340 (003/ 036)
+ ...
+
+(:OBSOLETE: !!! Overcompleteness of the interstitial plane waves
+The set of interstitial plane waves (IPW) (i.e., plane waves where the MT spheres are cut out) is overcomplete in the sense that a function (defined only in the interstitial region) has a nonunique representation in this set: one can always add the Fourier coefficients of a function, which is restricted to the MT spheres and zero outside. This leads to the fact that the IPW coefficients of the wave functions might not go quickly to zero for large G vectors, especially for a large '''G''' cutoff in the DFT calculation and large MT radii. In the ''GW'' calculation this gives rise to bad convergence with respect to the parameter ``GCUT`` ('''G''' cutoff in Spex) as well as jumps in the band structure. In this case Spex can exploit the nonuniqueness of the coefficients to make them approach zero for large '''G''' vectors as quickly as possible. Use the keyword ``MINCPW`` for this.:)
+
+Setting up ``spex.inp``
+-----------------------
+
+Mixed Basis (section ``MBASIS``)
+""""""""""""""""""""""""""""""""""""
+The bare Coulomb scattering {$ v(\mathbf r,\mathbf r') $} couples two electrons whose incoming and outgoing waves give rise to wave-function products (occupied/unoccupied). The mixed basis is optimized for representing wave-function products and derived from the FLAPW basis set. It is a combined set of interstitial plane waves (IPWs) and MT functions and is defined in section ``MBASIS`` of the input file ``spex.inp`` with the keywords
+
+* ``GCUT``: '''G''' cutoff for the interstitial plane waves,
+* ``LCUT``: L cutoffs (one for every atom type) for the MT functions,
+* ``SELECT``: determines the pairs of radial functions from whose products the MT functions are constructed. For every atom type two l cutoffs ``l1`` and ``l2`` must be specified for the functions {$u_l(r)$}. Optionally l cutoffs ``l3`` and ``l4`` can also be defined for {$\dot u_l(r)$}, which are neglected by default. The local orbitals, on the other hand, are always used in the construction of the products by default, but can also be chosen individually by series of zeros (neglect) and ones (select), e.g., ``10100``, ``lo1`` and ``lo2`` for occupied and unoccupied states. (Instead of ``0000111``, one can also use ``4031``.) Examples are --- ``2;3`` --- ``2,2;3,2`` --- ``2,2,10100;3,2,1110`` --- ``2,,10100;3,2,1110``
+ The so-defined set of product MT functions is usually highly linearly dependent. Therefore, in a subsequent optimization step one diagonalizes the overlap matrix and retains only those eigenfunctions with eigenvalues exceeding a tolerance value which is defined by the keyword ``TOL``.
+* ``OPTIMIZE``: Truncation of the polarization and related matrices to ``n``x``n`` by ``MB=n`` where ``n`` is an integer number between 1 and the size of the full mixed basis. Alternatively, one can use ``PW=n`` where plane waves are used instead of the mixed basis (or better: projections of plane waves onto the mixed basis) in which case the Coulomb matrix is analytically known. In both cases (``MB`` and ``PW``), you can also specify the minimum Coulomb eigenvalue instead of ``n`` as a negative real number, e.g. ``MB=-1.0`` . Alternatively, you can specify a "pseudo" reciprocal cutoff radius as a positive real number, e.g. ``MB=4.5``, which is more suitable for convergence tests. This cutoff value {$K$} is related to the minimum Coulomb eigenvalue {$c$} by {$K= \sqrt{4\pi/c}$}. The options ``MB`` and ``PW`` can speed up the calculation considerably.
+
+If any of the keywords (or the whole section) is not defined, Spex guesses reasonable values according to the FLAPW basis parameters.
+
+Coulomb Matrix (section ``COULOMB``)
+""""""""""""""""""""""""""""""""""""
+In many-body perturbation theory the perturbing force is the Coulomb interaction. Thus its representation in the mixed basis is one of the basic ingredients in the implementation. The IPWs are plane waves {$\exp[i(\mathbf{k+G})\mathbf{r}]$} from which the MT spheres are cut away. In a MT sphere a plane wave is given by the Rayleigh expansion over l quantum numbers. In a practical implementation we must introduce an l cutoff here (keyword ``LEXP``). Multipole interactions between the periodic MT functions and their periodic images are calculated in an Ewald summation involving a real and a Fourier-space sum. The cutoff radii of these are determined during the program run. There is a scaling parameter (``SCALE``) which must be adjusted if one of the sums take much longer to evaluate than the other one. In most calculations, the default values work reliably, though.
+
+Labels for k points (keyword ``KPT``)
+"""""""""""""""""""""""""""""""""""""
+The keyword ``KPT`` allows you to define labels (single characters) for k points, e.g. ``KPT X=[0,0,1]``. Here, ``X`` is the label and the square brackets denote that cartesian coordinates are used instead of internal coordinates (round brackets). Note: In order for the cartesian coordinates (square brackets) to be interpreted correctly, Spex must know the lattice constant. (In Fleur, this is the global scaling factor for the lattice vectors.) You can also add a prefactor by ``L=0.5*(0,0,1)`` or ``L=1/2*(0,0,1)``. The so-defined k points must be elements of the k mesh defined by ``BZ``. A single k point outside that mesh can be specified by the plus sign ``+`` as a special label. Each time you change this single k point you have to regenerate the k-point set and rerun the DFT code.
+
+Job definition (keyword ``JOB``)
+""""""""""""""""""""""""""""""""
+The most important entry in ``spex.inp`` is the job definition. It defines what Spex should do. Each job begins with one of the keywords ``HF, QP, SUSCEP, SUSCEPR, DIELEC`` or ``SCREEN`` followed by additional parameters:
+
+* ``HF``: Hartree-Fock correction, e.g., ``HF A:(1-3,5)`` where A is a k-point label or index and ``1-3,5`` are band indices. A spin index (up or down) can be specified by ``A:u(...)`` and ``A:d(...)``. More than one k-point definition is possible in one job, e.g., ``HF A:(1-3,5) 12:(1,2,6)``. Normally, only diagonal elements are calculated, i.e., {$ \langle \phi_{n\mathbf{k}} | \Sigma_{\mathrm{x}} | \phi_{n'\mathbf{k}} \rangle $} with {$ n=n' $}. Full calculations of the self-energy matrix (i.e., {$ n \neq n' $}) are possible by adding ``FULL`` (e.g., ``HF FULL A:(1-10)``). This allows a subsequent self-consistent HF calculation.
+* ``GW``: Quasiparticle correction with the ''GW'' approximation. Spin, k points and bands are defined as in ``HF``. A full self-energy calculation is also possible with ``GW FULL``.
+* ``SUSCEP``: Calculation of the susceptibility, e.g., ``SUSCEP A:{0..2,0.01}`` where A is a k-point label or index and a frequency mesh from 0 to 2 Ha in 0.02 Ha steps is specified. A more general frequency mesh can also be specified as a list of real or complex frequencies in a file. Then simply write the file name into the brackets, e.g., ``{foo}``. A spin index (``uu, dd, ud, du`` or ``+``) can be specified with ``A:ud{...}``. The spin index ``+`` denotes addition of ``uu`` and ``dd`` and is the default. Furthermore you can add ``MB`` or ``PW`` (see above) and ``OUT=n`` (separated with commas) to the job definition. With ``OUT=n`` an ``n``x``n`` matrix is written instead of only the head element. You can also write the full matrix to a binary file with ``OUT=BIN``. An expression like ``A:u{0..2,0.01},MB=100,OUT=4`` defines one job which can be followed by additional expressions each defining one job.
+* ``SUSCEPR``: Calculation of the renormalized susceptibility. Spin, k points and frequencies are defined as in ``SUSCEP``. Additionally you can specify a TDDFT kernel as in ``SUSCEPR A:{0..2,0.01},KERN=ALDA,MB=100``. As spin labels only ``ud``, ``du`` and ``+`` are allowed where the latter is the default. The former spin indices require a TDDFT kernel.
+* ``DIELEC``: Calculation of the dielectric function. K points and frequencies are defined as in ``SUSCEP``. A spin index cannot be specified.
+* ``SCREEN``: Calculation of the screened interaction. Syntax as in ``DIELEC``.
+
+In the case of an "empty" job definition (``JOB``), Spex will do basic checkings and then stop. (This is useful, e.g., for a simple Wannier interpolation or P/DOS calculation.)
+
+Output files
+------------
+The main results of the Spex calculation is written to standard output ([[SpexOutputFile|output example]]). You should pipe the output to a file (e.g. ``spex > spex.out``). The spectra are written to separate files:
+
+* ``suscep``: Susceptibility head element.
+* ``suscep.XXX.YYY``: Susceptibility ``XXX/YYY`` matrix element (key: ``OUT=NNN`` and ``XXX,YYY`` {$\le$} ``NNN``).
+* ``suscep.bin``: Binary output of the susceptibility matrix (key: ``OUT=BIN``).
+* ``suscepr``: Renormalized susceptibility head element.
+* ``suscepr.XXX.YYY``: Corresponding matrix element.
+* ``suscepr.bin``: Binary output.
+* ``dielec``: Head element of the dielectric matrix (dielectric function).
+* ``dielec.XXX.YYY``: Corresponding matrix element.
+* ``dielecLF``: Dielectric function with local-field effects.
+* ``screen``: Screened interaction head element.
+* ``screen.XXX.YYY``: Corresponding matrix element.
+* ``screen.bin``: Binary output.
+
+If more than one job is defined, the job number is appended to the file names (before the dot), e.g. ``suscep001.bin``.
+
+There are a number of additional files which can be written during the program run (``WRITE INFO``) and which tell you more about intermediate results. Some of them are numbered, such that files are never overwritten. From time to time you should remove them from disk, because the highest file count is ``NNN=999``, and they clog up directory listings.
+
+* ``spex.sf.NNN``: If projections ``MB=...`` or ``PW=...`` are used, these files contain the spectral function on the mesh given by ``FSPEC`` for any '''k''' point projected on {$exp(ikr)$}. Column 1: Frequency, column 2: value. The files are useful to check whether the mesh is dense enough and whether the spectral function approaches zero quickly. If it does not, ``GCUT`` might be underconverged or you must use ``MINCPW``.
+* ``spex.sew.NNN``: Correlation self-energy on the imaginary frequency axis (only for ``CONTINUE``). Column 1: imaginary frequency, column 2: fitted value (real part), column 3: fitted value (imaginary part), column 4: numerical value (real part), column 5: numerical value (imaginary part).
+* ``spex.sec.NNN``: Correlation self-energy on the real frequency axis. Column 1: real frequency, column 2: (fitted) value (real part), column 3: (fitted) value (imaginary part).
+* ``spex.head``: Head element of screened interaction on imaginary axis.
+
+Other output files written during the run are
+
+* ``spex.dos.NNN``: (Partial) density of states (keyword ``BANDINFO``).
+* ``bands``: Wannier-interpolated band structure (written after the Wannier construction and if ``KPTPATH`` is defined).
+* ``hamiltonian_r``: Wannier Hamiltonian at all R points defined by ``BZ`` (if ``bands`` is written).
+* ``hamiltonian_k``: Wannier-interpolated Hamiltonian at k points defined in the file ``klist`` (if ``bands`` is written).
+
+For certain keywords in ``spex.inp`` there are special output files that can be used for subsequent calculations.
+
+* ``spex.mb``: Mixed product basis (keyword ``RESTART``).
+* ``spex.core``: Core susceptibility (keyword ``CORES``, unmaintained).
+* ``spex.cou``: Coulomb matrix (keyword ``RESTART``).
+* ``spex.cous``: Special Coulomb matrix for calculation of spectra (keywords ``RESTART`` and ``NOSTORE``).
+* ``spex.cor``: Renormalized dielectric matrix for correlation self-energy (keyword ``RESTART``).
+* ``spex.spc``: Susceptibility/dielectric/etc. matrix for calculation of spectra (keyword ``RESTART``).
+* ``spex.sigx``: Exchange self-energy (HF) matrix elements (keyword ``JOB ... FULL``, can be used with ``RESTART SIG``).
+* ``spex.sigc``: Correlation self-energy (GW) matrix elements (keyword ``JOB ... FULL``, can be used with ``RESTART SIG``).
+* ``spex.qsgw``: QSGW (or HF) self-energy matrix (keyword ``JOB ... FULL``) for a subsequent self-consistent-field run (``gw=3`` in the FLEUR inp file).
+* ``spex.qp``: Quasiparticle energies and amplitudes (keyword ``JOB ... FULL``).
+
+.. note:: ``spex.set.NNN``: Correlation self-energy for any quasiparticle state on the imaginary time axis (only for ``CONTINUE`` and ``IMAGTIME``). Column 1: imaginary time, column 2: interpolated value, column 3: numerical value. :)
diff --git a/source/spex_faq.rst b/source/spex_faq.rst
new file mode 100644
index 0000000000000000000000000000000000000000..5e4a607c9d56cb03c62bef1f076284463d4d7262
--- /dev/null
+++ b/source/spex_faq.rst
@@ -0,0 +1,91 @@
+.. _spex_faq:
+
+===============
+Spex FAQ
+===============
+
+Configuration
+---------------
+
+**Configuration stops with "MPI-3 compiler requested, but couldn't use MPI-3."**
+
+.. highlights:: Your MPI compiler does not support the MPI-3 standard, which is required for the parallelized version.
+
+**Configuration stops with "Could not link HDF5 library."**
+
+.. highlights:: Either (1) you do not have the HDF5 library, or (2) you have HDF5 but it is not installed in a standard directory (then, use the environment variables ``LDFLAGS`` and ``FCFLAGS``), or (3) you have HDF5 but it has been compiled with a different compiler.
+
+You can get further information from the file ``config.log``.
+
+Installation
+-------------
+
+**I want to use my own DFT program instead of Fleur.**
+
+.. highlights:: You will have to rewrite or modify the routines in ``read_write.f`` according to your output files and recompile the source. If you use a different mesh for radial integration than Fleur/Spex, you must either interpolate the values to the new mesh (recommended) or write new integration routines in ``numerics.f`` (``intgr, intgrf_init, intgrf, primitive, primitivef``).
+
+**My compiler fails to compile the source.**
+
+.. highlights:: We have successfully compiled and run the code with the Intel Fortran compiler (since version 12.1.3), Gfortran (since version 5), and NAG Fortran (since version 6.2 6214). As of May 2018, PGI compilers have a bug (TPR 23745) that prevents compilation of the Spex code. For the parallelized version, we have used IntelMPI and MPICH. In case of configuration or compilation problems, please inform us.
+
+Usage
+-------
+
+**Spex needs too much memory.**
+
+.. highlights:: Try the options ``STOREIBZ`` and/or ``NOSTORE``. You can also reduce the maximally allowed memory storage with the keyword ``MEM``.
+
+**The program stops with a segmentation fault.**
+
+.. highlights:: If this is caused by a too high memory demand you should try the option ``STOREIBZ`` (and ``NOSTORE``) or reduce the maximally allowed storage with the keyword ``MEM``. With some compilers this can also happen, if the stack size is limited. In this case set the stack size to "unlimited" (e.g., ``limit stacksize unlimited`` in csh). However, this is not possible with Mac computers. In this case add "Wl,-stack_size,0x10000000,-stack_addr,0xc0000000" to LDFLAGS or "-heap-arrays" to FFLAGS in Makefile.
+
+**Degenerate quasiparticle states have slightly different energies or I get slightly different results at equivalent k points**
+
+.. highlights:: There are several possible reasons for this:
+
+ * In your DFT output groups of degenerate states are cut at the top. Then you should specify ``NBAND`` explicitly.
+ * In the calculation of the susceptibility groups of degenerate states are split into different packets. This can only be avoided with a larger ``MEM`` if possible.
+ * The tetrahedron method breaks the symmetry. You can use ``GAUSS`` to test this.
+
+**My ``GW`` calculations do not converge wrt the k-point set.**
+
+.. highlights:: In systems with small band gaps (e.g. GaAs) the default treatment of the point k=0 might lead to bad k-point convergence. In this case try the option ``NOZERO`` which guarantees a smooth (but slower) convergence.
+
+**My calculations do not converge wrt the number of bands.**
+
+.. highlights:: The parameters are not completely independent. If you increase the number of bands, you might need a higher G cutoff (``GCUT`` in section ``MBASIS``). You might also have to increase the number of frequencies for sampling the spectral function. Instead, you can also specify the value of the first nonzero frequency point as the first argument to ``FSPEC``. Then the number of frequencies automatically increases as you increase the number of bands.
+
+**The real-space (reciprocal-space) summation in the routine structureconstant takes very long.**
+
+.. highlights:: Increase (decrease) ``SCALE`` in section ``COULOMB``.
+
+**There are jumps in the band structure.**
+
+.. highlights:: Band structures should not be calculated without the keyword ``NOZERO``. In the case of a metal a jagged band structure (especially in Hartree-Fock) might also be caused by the tetrahedron method. Then you must improve the k-point sampling. Alternatively, you may try the Gauss integration method (``GAUSS``).
+
+Error Messages
+--------------
+.. highlights:: This is not a complete list of errors. In many error messages a possible solution is already indicated. These are not listed. The program might also stop at places where this is not supposed to happen. Then the message ends with ``...(bug?)``. In this case you should contact the programmers.
+
+**error while loading shared libraries: lib...: cannot open shared object file: No such file or directory**
+
+.. highlights:: The library "lib..." is not found in standard directories. In this case, the location of the library has to be defined in the environment variable ``LD_LIBRARY_PATH``.
+
+**correlation: Solution of quasiparticle equation did not converge after 1000 iterations.**
+
+.. highlights:: You probably have poles very close to the real frequency axis (normally with a small weight). Usually a slight change of parameters solves the problem.
+
+**diagonalize_coulomb: Negative/Residual eigenvalue.**
+
+.. highlights:: Your Coulomb matrix seems to be underconverged. Try to increase the parameter ``LEXP`` or adjust the scaling with ``SCALE``. If this does not solve the problem, the reason might be that the LAPACK routine fails to solve the general eigenvalue problem because of a nearly singular overlap matrix. Then you might try the option ``CUTZERO`` which removes linear dependencies from the overlap matrix and reattempts the diagonalization.
+
+**intgrf: Negative exponent x in extrapolation a+c*r**x.**
+
+.. highlights:: The first (nonzero) point of your radial MT mesh might be too far away from zero. Modify the mesh accordingly. You can also try to increase ``TOL`` in section ``MBASIS``.
+
+**modulo1: argument not an element of k-point set.**
+
+.. highlights:: This message is caused by an error in the k-point set. Either the k-point sets defined by ``BZ`` and in the DFT output file are different or the k points defined by ``KPT`` are not elements of the set. The latter happens, if you use ``KPT ?=[?,?,?]`` and the lattice constant is not properly defined since ``[?,?,?]`` is interpreted with respect to the lattice constant.
+
+To be continued ...
+
diff --git a/source/theory.rst b/source/theory.rst
new file mode 100644
index 0000000000000000000000000000000000000000..fba6d5aa5369b075bedf73803c63e4c24e76bca2
--- /dev/null
+++ b/source/theory.rst
@@ -0,0 +1,81 @@
+.. _theory:
+
+========
+Theory
+========
+
+``GW`` calculations [[#GW]]
+---------------------------
+
+The ''GW'' approach has become a routine method to calculate quasiparticle spectra from first principles. It is based on many-body perturbation theory, in which the Dyson equation (here in simplified notation)
+
+:math:`{\displaystyle G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)-v^{\mathrm{xc}}]G(\omega)}`
+
+links the interacting Green function :math:`{G}` to the non-interacting one :math:`{G_0}`. (A non-magnetic system is considered in this section for simplicity.) The latter is obtained from a mean-field solution with frequency-independent potential :math:`{v^{\mathrm{xc}}}`, which in the case of a KS system would correspond to the local exchange-correlation potential (the nonlocal Hartree-Fock exchange potential and the ''hermitianized'' self-energy of QSGW are other examples). The self-energy :math:`{\Sigma^{\mathrm{xc}}}` connects the non-interacting mean-field system to the real interacting system and thus accounts for the many-body electronic exchange and correlation effects. It is not known exactly and must be approximated in practice. The most popular approximation to the electronic self-energy is the GW approximation,
+
+: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'\,}`
+
+which has been shown to yield accurate band structures for a wide range of materials. Here, :math:`{G}` is the single-particle Green function (usually :math:`{G_0}` is employed) and :math:`{W}` is the screened Coulomb interaction, which is the effective interaction potential in an interacting system, incorporating the dynamical screening processes. The screened Coulomb interaction is the solution of a Dyson-type equation
+
+:math:`{\displaystyle W(\omega)=v+vP(\omega)W(\omega)\,}`
+
+with the random-phase approximation for the polarization function :math:`{P}`
+
+:math:`{\displaystyle P(\omega)=-\frac{i}{\pi}\int_{-\infty}^\infty G(\omega-\omega')G(\omega')d\omega'\,.}` [[#eq:P]]
+
+Usually, the Green function :math:`{G_0}` of a mean-field system is taken for :math:`{G}`, in which case it takes the form
+
+:math:`{\displaystyle G_0(\mathbf{r},\mathbf{r}';\omega)=\sum_\mathbf{k}\sum_n \frac{\phi_{\mathbf{k}n}(\mathbf{r})\phi^*_{\mathbf{k}n}(\mathbf{r}')}{\omega-\epsilon_{\mathbf{k}n}+ i\eta\mathrm{sgn}(-\epsilon_{\mathbf{k}n}\mathrm{F}}\,,}` [[#eq:G0]]
+
+where :math:`{\phi_{n\mathbf{k}}(\mathbf{r})}` and :math:`{\epsilon_{n\mathbf{k}}}` are the eigenfunctions and eigenvalues, respectively, of the mean-field system. Here, a natural choice is to employ the solution of the KS equations
+
+:math:`{\displaystyle\left\{-\frac{\nabla}{2}+v^\mathrm{ext}(\mathbf{r})+v^\mathrm{H}(\mathbf{r})+v^\mathrm{xc}(\mathbf{r})\right\}\phi_{\mathbf{k}n}(\mathbf{r})=\epsilon_{\mathbf{k}n}\phi_{\mathbf{k}n}(\mathbf{r})\,,}`
+
+often with the local-density approximation (LDA) for the exchange-correlation potential. In this respect, one says that LDA is the ''starting point'' for the ''GW'' calculation. Other starting points are generalized gradient approximation (GGA), LDA+''U'', hybrid functionals, Hartree-Fock, QSGW, et cetera. The explicit form of :math:`{G_0}` enables us to perform the frequency integration in [[#eq:P|:math:`{P}`]], which yields
+
+:math:`{\displaystyle P\left(\mathbf{r},\mathbf{r}',\omega\right)=2\sum_{\mathbf{q},\mathbf{k}}^{\mathrm{BZ}}\sum_{n}^{\mathrm{occ}}\sum_{n'}^{\mathrm{unocc}}\phi_{\mathbf{q}n}\left(\mathbf{r}\right)\phi_{\mathbf{q}n}^{*}\left(\mathbf{r}'\right)\phi_{\mathbf{q}+\mathbf{k}n'}\left(\mathbf{r}'\right)\phi_{\mathbf{q}+\mathbf{k}n'}^{*}\left(\mathbf{r}\right)\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(\mathbf{r},\mathbf{r}';\omega')}{\omega-\omega'+i\eta\mathrm{sgn}(\omega')}d\omega'\,,}` [[#eq:Pexp]]
+
+where the last equality implicitly defines the spectral function :math:`{S(\omega)}`.
+
+The Dyson equation can be rewritten in the form of an effective equation of motion of a single particle, the quasiparticle equation
+
+:math:`{\displaystyle\left\{-\frac{\nabla}{2}+v^\mathrm{ext}(\mathbf{r})+v^\mathrm{H}(\mathbf{r})\right\}\psi_{\mathbf{k}n}(\mathbf{r})+\int \Sigma^\mathrm{xc}(\mathbf{r},\mathbf{r}';E_{\mathbf{k}n})\psi_{\mathbf{k}n}(\mathbf{r}')d^3 r=E_{\mathbf{k}n}\psi_{\mathbf{k}n}(\mathbf{r})}` [[#eq:qpeq]]
+
+with the (non-orthonormal) quasiparticle wavefunctions and (complex-valued) quasiparticle energies. It can be shown that the interacting Green function takes the same form as [[#eq:G0|:math:`{G_0}`]] (Lehmann representation) if the eigenfunctions and eigenvalues are replaced by the latter. The fact that this equation of motion is formally similar to the KS equation motivates to use perturbation theory of first order and write the quasiparticle energies as
+
+:math:`{\displaystyle E_{n\mathbf{k}}=\epsilon_{n\mathbf{k}}+\langle\phi_{n\mathbf{k}}|\Sigma^\mathrm{xc}(E_{n\mathbf{k}})-v^\mathrm{xc}|\phi_{n\mathbf{k}}\rangle\approx\epsilon_{n\mathbf{k}}+Z_{n\mathbf{k}}\langle\phi_{n\mathbf{k}}|\Sigma^\mathrm{xc}(\epsilon_{n\mathbf{k}})-v^\mathrm{xc}|\phi_{n\mathbf{k}}\rangle}` [[#eq:qppert]]
+
+with the renormalization constant :math:`{Z_{n\mathbf{k}}=[1-\partial\Sigma^{\mathrm{xc}}/\partial\omega(\epsilon_{n\mathbf{k}})]^{-1}}`. Both expressions can be taken to evaluate the quasiparticle energy, the first one is a non-linear equation in :math:`{E_{n\mathbf{k}}}` and requires an iterative solution. We thus have three ways to calculate the quasiparticle energies: (a) by solving the full [[#eq:qpeq|quasiparticle equation]], (b) with the [[#eq:qppert|non-linear equation]], and (c) using the [[#eq:qppert|linearized solution]]. The mathematical complexity, the computational cost, and the accuracy decreases in this order. For example, solution (a) requires the full self-energy matrix (i.e., including the off-diagonal elements) to be evaluated. Once the full matrix :math:`{\Sigma^\mathrm{xc}_{\mathbf{k},nn'}(\omega)}` is known, one can proceed to perform a self-consistent calculation within the QSGW approach. In this approach, the self-energy matrix is ''hermitianized'' and made frequency independent. This self-energy can then be used to replace :math:`{v^\mathrm{xc}}` in the KS equation, defining a new mean-field system. A self-consistent solution of this mean-field system (using a DFT code) is then employed as a new starting point and so on until self-consistency in the quasiparticle energies is achieved.
+
+
+
+Interaction parameters (Hubbard ``U``)
+--------------------------------------
+
+
+Wannier functions
+--------------------
+
+
+RPA exchange-correlation energy
+-------------------------------
+
+
+
+Spectra (EELS, spin excitation spectra etc.)
+--------------------------------------------
+
+The spectral function of response quantities (e.g., [[#eq:Pexp|:math:`{S(\omega)}`]]) are often related to spectra measured in experiments. The imaginary part of the head element (:math:`{\varepsilon^{-1}_{\mathbf{k},\mathbf{GG}'}, \mathbf{G}=\mathbf{G}'=0}`) of the inverse of the dielectric matrix
+
+:math:`{\displaystyle \varepsilon_\mathbf{k}(\omega)=I - P_\mathbf{k}(\omega) v_\mathbf{k}}`
+
+with [[#eq:Pexp|:math:`{P(\omega)}`]], the identity matrix :math:`{I}`, and the bare Coulomb interaction :math:`{v}` is measured in valence electron-energy loss spectroscopy (EELS). Note that :math:`{W_\mathbf{k}(\omega)=\varepsilon_\mathbf{k}^{-1}(\omega)v_\mathbf{k}}`. Optical absorption is related to the macroscopic dielectric function :math:`{\varepsilon_\mathrm{M}(\omega)=1/\varepsilon^{-1}_{\mathbf{k=0},\mathbf{00}'}(\omega)}`. We have adopted a pure plane-wave representation. The dielectric function can also be used to calculate the density response function :math:`{R_\mathbf{k}(\omega)=\varepsilon_\mathbf{k}^{-1}(\omega)P_\mathbf{k}(\omega)}`, where :math:`{R(\mathbf{r},\mathbf{r}';\omega)=\delta n(\mathbf{r},\omega)/\delta v_\mathrm{ext}(\mathbf{r}',\omega)}` gives the dynamical linear response of the electron number density with respect to changes in the external potential. In this respect, :math:`{P(\mathbf{r},\mathbf{r}';\omega)}` gives the linear response of the density with respect to changes in the effective potential, e.g., the effective KS potential.
+
+In photoelectron spectroscopy, one measures the imaginary part of the Green function, which, using the Dyson equation, can be written as
+
+:math:`{\displaystyle A_\mathbf{k}(\omega)=\mathrm{Im} \sum_n \frac{1}{\omega-\epsilon_{\mathbf{k}n}-\Sigma_\mathbf{k}^{\mathrm{xc}}(\omega)}.}`
+
+
+``GT`` self-energy
+------------------
+