tutorials.rst 39.5 KB
Newer Older
1 2 3 4 5
.. _tutorials:

6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270
In this section, we provide detailed tutorials for the main features of the Spex code.

GW calculations
The Spex code was started as a pure GW code. The GW approach remains the main computational method implemented in the code.

One-shot approach
A simple input file for a one-shot GW calculation for bulk silicon was already presented in :ref:`getting_started`. The following, equivalent input file is even more minimalistic.

.. code-block:: bash

  BZ 4 4 4
  JOB GW 1:(1,2,5) 7:(1,3,5) 3:(1-3,5)

It only contains the two essential keywords that any Spex input file must contain: ``BZ`` and ``JOB``, specifying, respectively, the k mesh for the Brillouin-zone sampling and the requested type of calculation. Here, we calculate quasiparticle corrections for the bands (1,2,5) of the first, (1,3,5) of the seventh, and (1,2,3,5) of the third k point. The band indices are defined according to the ordering in the mean-field input data. In the example, we have left out some band indices because of energy degeneracy. For example, bands 3 and 4 (of the first k point) are degenerate with band 2.

In the case of magnetic systems, the quasiparticle energies for both spins are calculated by default. Alternatively, one can choose the spin index by using u (spin up) or d (spin down), e.g., ``JOB GW 1:u(1,2,5) 1:d(1,2,5)`` is identical to ``JOB GW 1:(1,2,5)``.

The k-points follow an internal ordering: first the irreducible parent k points, then the remaining equivalent k points. A list of the irreducible k-vectors is written to the output file. In the present example, there are eight irreducible k points and 64 k points in total. All 64 k-point indices can be used in the job definition. We have chosen the ones from the irreducible set. They correspond to the Γ (index 1), X (index 7), and L point (index 3). Clearly, for a different k-point set, e.g., 8x8x8, the k-point indices would change (except for the index 1, which always corresponds to the Γ point). Therefore, there is the possibility of defining k-point labels with 

``KPT X=1/2*(0,1,1) L=1/2*(0,0,1)``
``KPT X=[1,0,0] L=1/2*[1,1,-1]``

The labels must be single (upper- or lower-case) characters. The two definitions are equivalent. The first gives the k vectors in internal coordinates, i.e., in the basis of the reciprocal basis vector, the second in cartesian coordinates in units of :math:`{2\pi/a}`  with the lattice constant a . In the case of the fcc lattice of silicon, the lattice basis vectors are :math:`{a_1=a(011)/2}` , :math:`a_2=a(101)/2` , and :math:`a_3=a(110)/2` . The reciprocal basis vectors are thus :math:`{b_1=2\pi/a(−111)}`  et cetera. For the X point, e.g., one then gets :math:`{(a_2+a_3)/2=2\pi/a(100)}` . In order for the cartesian coordinates (square brackets) to be interpreted correctly, Spex must obviously know the lattice constant a. (In Fleur, the lattice constant is taken to be the global scaling factor for the lattice vectors.) The so-defined k points must be elements of the k mesh defined by ``BZ``. We will later discuss how points outside the k mesh can be considered using the special label +. With the k labels, the above job definition can be written as ``JOB GW 1:(1,2,5) X:(1,3,5) L:(1-3,5)`` as in the input file above There are two more special k-point labels: ``IBZ`` and ``BZ`` (e.g., ``JOB GW IBZ:(1,2,5)``) stand for all k points in the irreducible and the full k set, respectively. (The label ``BZ`` is included for completeness but is not needed in practice.) The ``IBZ`` label is helpful when a self-consistent GW calculation is to be performed, which requires the self-energy to be calculated for the whole irreducible Brillouin zone.

The quasiparticle energies are written to standard output in tabular form::

    Bd     vxc    sigmax    sigmac         Z        KS        HF          GW lin/dir
    1 -12.15897 -19.20415   6.60649   0.65146  -6.19011 -13.23529  -6.36401   0.73681
                            1.05993  -0.10556                      -6.36806   0.72704
    2 -13.30408 -14.45814   0.69809   0.76876   5.77669   4.62263   5.42615   0.00070
                            0.00039  -0.00088                       5.42695   0.00081
    5 -11.49631  -7.27160  -3.86726   0.76385   8.34013  12.56485   8.61313  -0.00731
                           -0.00708  -0.00533                       8.61239  -0.00749

The columns are

* ``Bd``  band index,
* ``vxc`` expectation value of the exchange-correlation potential :math:`v^{xc}`, e.g., -12.15897 eV, 
* ``sigmax``  expectation value of the exchange self-energy :math:`\sum^x`, e.g., -19.20415 eV, 
* ``sigmac``  expectation value of the correlation self-energy :math:`\sum^c`, e.g., (6.60649+1.05993i) eV (note that the self-energy is complex), 
* ``Z`` renormalization factor Z, e.g., (0.65146-0.10556i) eV, 
* ``KS``  mean-field energy eigenvalue (e.g., from the previous KS-DFT calculation), e.g., -6.19011 eV, 
* ``HF``  Hartree-Fock value (one-shot), e.g., -13.23529 eV, 
* ``GW``  GW quasiparticle values, e.g., (-6.36401+0.73681i ) eV (linearized solution) and (-6.36806+0.72704i ) eV (direct solution). 

The two GW values for the quasiparticle energy follow two common methods to approximately solve the quasiparticle equation

:math:`{\displaystyle\left\{-\frac{\nabla}{2}+v^\mathrm{ext}(\mathbf{r})+v^\mathrm{H}(\mathbf{r})\right\}\psi_{\mathbf{k}n}(\mathbf{r})+\int \Sigma^\mathrm{xc}(\mathbf{r},\mathbf{r}';E_{\mathbf{k}n})\psi_{\mathbf{k}n}(\mathbf{r}')d^3 r=E_{\mathbf{k}n}\psi_{\mathbf{k}n}(\mathbf{r})}`

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

:math:`{\displaystyle E_{\mathbf{k}n}=\epsilon_{\mathbf{k}n}+\langle\phi_{\mathbf{k}n}|\Sigma^\mathrm{xc}(E_{\mathbf{k}n})-v^\mathrm{xc}|\phi_{\mathbf{k}n}\rangle\approx\epsilon_{\mathbf{k}n}+Z_{\mathbf{k}n}\langle\phi_{\mathbf{k}n}|\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})-v^\mathrm{xc}|\phi_{\mathbf{k}n}\rangle}`

with the single-particle wavefunction :math:`{\phi_{\mathbf{k}n}}` and the frequency-independent potential :math:`{v^{\mathrm{xc}}}`, which in the case of a KS solution would correspond to the local exchange-correlation potential; the nonlocal Hartree-Fock exchange potential and the ''hermitianized'' self-energy of QSGW (see below) are other examples. :math:`{Z_{\mathbf{k}n}=[1-\partial\Sigma^{\mathrm{xc}}/\partial\omega(\epsilon_{\mathbf{k}n})]^{-1}}` is called the renormalization factor. The two expressions on the right-hand side correspond to the "linearized" and "direct" (iterative) solutions given in the output. The direct solution takes into account the non-linearity of the quasiparticle equation and is thus considered the more accurate result. However, there is usually only little difference between the two values.
Up to this point, the job syntax for Hartree Fock (``JOB HF``), PBE0 (``JOB PBE0``), screened exchange (``JOB SX``), COHSEX (``JOB COSX``), and ''GT'' (``JOB GT`` and ``JOB GWT``) calculations is identical to the one of ``GW`` calculations, e.g., ``JOB HF FULL X:(1-10)``. Except the latter (``GT``), all of these methods are mean-field approaches, so one only gets one single-particle energy (instead of a ''linearized'' and a ''direct'' solution) for each band.


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

::math:`{\displaystyle G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)-v^{\mathrm{xc}}]G(\omega)}`

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

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

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


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

| Examples | ``SPECTRAL {-10eV:1eV,0.01eV}``       | | Write spectral function :math:`{\text{Im}G(\omega)}` on the                                          |
|          |                                       | | specified frequency mesh to the file "spectral".                                                     |
|          | ``SPECTRAL {-10eV:1eV,0.01eV} 0.002`` | | Bound imaginary part from below by 0.02, preventing                                                  |
|          |                                       | | peak widths to become too small to be plotted.                                                       |

Alhough ``GW`` calculations can be readily performed with the default settings, the user should be familiar to some degree with the details of the computation and how he/she can influence each step of the program run. Also note that the default settings might work well for one physical system but be unsuitable for another. 

The ``GW`` self-energy 

: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'\,.}`

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

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

This integral equation turns into a matrix equation 

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

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

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

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

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

| Examples | ``MESH 12 15.0``   | Use a mesh containing twelve frequencies for [0,i15] htr                                         |
|          | ``MESH 12+2 15.0`` | | Add points at low frequencies, two (one) between :math:`{\omega_1}` and :math:`{\omega_2}`     |
|          |                    | | (:math:`{\omega_2}`  and :math:`{\omega_3}` )                                                  |
|          | ``MESH 12 -15.0``  | Use alternative mesh definition.                                                                 |
|          | ``MESH mesh.dat``  | Read frequency mesh from file "mesh.dat"                                                         |

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

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

The second method is ``contour integration``, in which the frequency integration is performed explicitly, however not along the real frequency axis but on a deformed integration contour that avoids the real frequency axis as well as possible. This integration contour starts from :math:`{-\infty}`, describes an infinite quarter circle to :math:`{-i\infty}`, then runs along the imaginary frequency axis to :math:`{i\infty}`, and finishes, again with an infinite quarter circle, to :math:`{\infty}`. The two quarter circles do not contribute to the integral (because the integrand behaves as :math:`{\propto \omega^{-2}}`). Furthermore, depending on the frequency argument :math:`{\omega}` of the self-energy, one has to add a few residues coming from the poles of the Green function in the interval [:math:`{0,\omega-\epsilon_\mathrm{F}}`] if :math:`{\omega>\epsilon_\mathrm{F}}` and [:math:`{\epsilon_\mathrm{F}-\omega,0}`] otherwise, which requires the correlation screened interaction :math:`{W^\mathrm{c}(\omega)}` to be evaluated on this interval of the real axis. As a consequence, the calculations are more demanding in terms of computational complexity and cost (time and memory). Also, contour integration requires additional input parameters and is therefore somewhat more difficult to apply. However, the results are more accurate. In particular, they are not affected by the "ill-definedness" of the analytic continuation. 

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

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

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

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

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

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

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.

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

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

(*) The matrix elements of the exchange-correlation potential are needed in the solution of the [[#eq:qppert|quasiparticle equation]]. By default, these matrix elements are taken from the input data provided by the DFT code (e.g., in the file "vxcfull" written by Fleur). Alternatively, Spex can calculate the matrix using the potentials from the input data. (This is required for PBE0.) For that, the keyword ``VXC`` in section ``SENERGY`` can take the arguments ``READ`` and ``CALC``.

| Examples | ``VXC READ`` | Read :math:`{v^\mathrm{xc}}` expectation values from harddisc.  |
|          | ``VXC CALC`` | | Calculate :math:`{v^\mathrm{xc}}` expectation values          |
|          |              | | using the xc potential defined in stars and lattice harmonics.|

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"         |
|  |     --        | write                 |
|  | ``RESTART``   | read-write write      |
|  | ``RESTART 2`` | read-write read-write |

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

There are other (binary) files that are written during a program run and that may be reused later. The ones relevant for ``GW`` (same rules as for "spex.cor") are

* ``spex.mb``: contains MT part of product basis,
* ``spex.ccou``: contains ``contracted`` screened interaction (i.e., summed over ``'k``'); this contraction is needed for the self-energy core contribution (``CORES``), ``COSX``, and ``IBC``,
* ``spex.core``: contains core contribution to ``W`` (keyword ``CORES``).

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

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

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

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

GW band structures

Beyond perturbative GW and QSGW

GW with MPI

Polarization function