@@ -278,7 +278,7 @@ In this section, three ways of calculating ``GW`` band structures are discussed.

...

@@ -278,7 +278,7 @@ In this section, three ways of calculating ``GW`` band structures are discussed.

* Using the :ref:`kpath` label (single Spex run)

* Using the :ref:`kpath` label (single Spex run)

* With :ref:`kptadd` additional :math:`q` points (multiple Spex run)

* With :ref:`kptadd` additional :math:`q` points (multiple Spex run)

* :ref:`_wanint` Wannier interpolation (single --; but long --; Spex run).

* Wannier interpolation (single --; but long --; Spex run).

The calculation of a ``GW`` band structure is not as straightforward as in the case of 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``.

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``.

...

@@ -303,7 +303,7 @@ As an example, we want to plot a band structure for the path from L over :math:`

...

@@ -303,7 +303,7 @@ As an example, we want to plot a band structure for the path from L over :math:`

KPT +=...

KPT +=...

---------

---------

This section explains an extension to the keyword :ref:`kpt` (already described before), which will later enable the calculation of smooth band structures. The extension allows adding arbitrary q points, one at a time, allowing one to investigate the quasiparticle spectrum at momenta that are not element of the original 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 :ref:`_kpt`. As explained above, each additional point :math:`q` must be accompanied with all shifted points :math:`q+k`. So, before running Spex, we must let the DFT code generate the wavefunctions and energies at the shifted k-point set ``'q+k``'. This is done in the usual way by creating the k-point file :ref:`wrtkpt` or ``-w`` flag) and running the DFT code. (In the case of Fleur, the shifted k points are appended to the list in the file "kpts".) The additional q point can be addressed in the job definition with the special ``+`` label, e.g., ``JOB GW 1:(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.

This section explains an extension to the keyword :ref:`kpt` (already described before), which will later enable the calculation of smooth band structures. The extension allows adding arbitrary q points, one at a time, allowing one to investigate the quasiparticle spectrum at momenta that are not element of the original 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 :ref:`kpt`. As explained above, each additional point :math:`q` must be accompanied with all shifted points :math:`q+k`. So, before running Spex, we must let the DFT code generate the wavefunctions and energies at the shifted k-point set ``'q+k``'. This is done in the usual way by creating the k-point file :ref:`wrtkpt` or ``-w`` flag) and running the DFT code. (In the case of Fleur, the shifted k points are appended to the list in the file "kpts".) The additional q point can be addressed in the job definition with the special ``+`` label, e.g., ``JOB GW 1:(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.

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 :ref:`wrtkpt` (or with the option ``-w``) creates two files, containing the usual list of k points (Fleur: "kpts") and the list of q points (Fleur: "qpts"). For each of the q points, we would have to run, in this order, Spex, Fleur, and Spex again. To simplify this task, there is a shell script called "spex.band". This shell script performs all the necessary steps automatically (it uses the same environment variables as :ref:`selfc.sh` and produces, for each q point, one output file named ``spex_???.out`` (where ``???`` is a three-digit counting index, or more digits if ``???>999``).

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 :ref:`wrtkpt` (or with the option ``-w``) creates two files, containing the usual list of k points (Fleur: "kpts") and the list of q points (Fleur: "qpts"). For each of the q points, we would have to run, in this order, Spex, Fleur, and Spex again. To simplify this task, there is a shell script called "spex.band". This shell script performs all the necessary steps automatically (it uses the same environment variables as `selfc.sh` and produces, for each q point, one output file named ``spex_???.out`` (where ``???`` is a 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.

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.

:ref:`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 :ref:`wannier`. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to ``spex.inp``:

As a third alternative, we can use Wannier interpolation to interpolate between the few ``'k``' points of the original k mesh. The construction of Wannier orbitals is explained in a :ref:`wannier`. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to ``spex.inp``:

.. code-block:: bash

.. code-block:: bash

...

@@ -348,7 +348,7 @@ One can also go beyond the perturbative solution of the equation :eq:`qpeq` and

...

@@ -348,7 +348,7 @@ One can also go beyond the perturbative solution of the equation :eq:`qpeq` and

``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.

``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.

:ref:`qsgw` If the job definition contains ``FULL`` and :ref:`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]:

If the job definition contains ``FULL`` and 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]:

@@ -368,7 +368,7 @@ for off-diagonal elements. The *hermitianized* QSGW operator is then obtained fr

...

@@ -368,7 +368,7 @@ for off-diagonal elements. The *hermitianized* QSGW operator is then obtained fr

* Set ``gw=2`` in the Fleur input file.

* Set ``gw=2`` in the Fleur input file.

* Run Fleur.

* 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. :ref:`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:

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. 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,

* ``SPEX_EXEC`` - Spex executable,

* ``FLEUR_EXEC`` - Fleur executable,

* ``FLEUR_EXEC`` - Fleur executable,

...

@@ -403,7 +403,7 @@ The polarization function gives the linear change in the electronic density of a

...

@@ -403,7 +403,7 @@ The polarization function gives the linear change in the electronic density of a

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.)

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 equation :eq:`eqP` involves an infinite band summation over :math:`{n'}`. In practice, this summation is truncated by the number of bands (keyword :ref:`nbabd`, which thus becomes an important convergence parameter. The self-energy contains an infinite band summation as well.

The expression for equation :eq:`eqP` involves an infinite band summation over :math:`{n'}`. In practice, this summation is truncated by the number of bands (keyword `nbabd`, which thus becomes an important convergence parameter. The 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``).

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``).