diff --git a/source/tutorials.rst b/source/tutorials.rst
index 9f39a0e2be378dca61dad070a0b0446aae911ce9..97c6f3c3cac64be74ba33e6ab0293a5d97a72ef5 100644
--- a/source/tutorials.rst
+++ b/source/tutorials.rst
@@ -119,7 +119,7 @@ This integral equation turns into a matrix equation
if the quantities :math:`{W}`, :math:`{\varepsilon}`, and :math:`{v}` ; are represented in the :ref:`mbp`, which thus has to be converged properly. The dielectric function, in turn, describes the change of the internal potential through screening processes and is related to the polarization matrix by :math:`{\varepsilon(\mathbf{k},\omega)=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 :ref:`polar`.
-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.)
+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 equation :eq:`selfene` 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.
@@ -206,7 +206,7 @@ Both the exchange (HF) and the correlation self-energy contain terms of the form
ALIGNVXC (SENERGY)
------------------
-(*) 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``.
+(*) The equation :eq:`qpeq` depends explicitly on the mean-field starting point through the eigensolutions :math:`{\{\phi_{\mathbf{k}n}\}}`. This dependence is more obvious in the equation :eq:`qppert` 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 |
@@ -220,7 +220,7 @@ ALIGNVXC (SENERGY)
VXC (SENERGY)
--------------
-(*) 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``.
+(*) The matrix elements of the exchange-correlation potential are needed in the solution of the equation :eq:`qppert`. 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. |
@@ -277,8 +277,8 @@ GW band structures
In this section, three ways of calculating ``GW`` band structures are discussed.
* Using the :ref:`kpath` label (single Spex run)
-* With [[#KPTADD|additional ``'q``']] points (multiple Spex run)
-* [[#WanInt|Wannier interpolation]] (single --; but long --; Spex run).
+* With :ref:`kptadd` additional :math:`q` points (multiple 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``.
@@ -288,7 +288,7 @@ If we restrict ourselves to the set of k points defined in ``BZ``, the band stru
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.
+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 :ref:`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 :ref:`spex.out` 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 |
@@ -303,7 +303,7 @@ As an example, we want to plot a band structure for the path from L over :math:`
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 [[#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.
+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.
+----------+-------------------------+---------------------------------------------------------------------+
@@ -315,11 +315,11 @@ This section explains an extension to the keyword :ref:`kpt` (already described
| | ``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 :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|"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``).
+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.
-[[#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
@@ -332,12 +332,12 @@ From these files the band-structure data is extracted with ``spex.extr g -b spex
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.
+Wannier interpolation can also be used in conjunction with the keyword :ref:`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
+One can also go beyond the perturbative solution of the equation :eq:`qpeq` 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
@@ -348,7 +348,7 @@ One can also go beyond the perturbative solution of the [[#Eq:qpeq|quasiparticle
``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]:
+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]:
.. 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
@@ -368,7 +368,7 @@ for off-diagonal elements. The *hermitianized* QSGW operator is then obtained fr
* 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:
+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,
* ``FLEUR_EXEC`` - Fleur executable,
@@ -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.)
-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.
+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``).
@@ -426,7 +426,7 @@ GAUSS
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).
+The most important keyword in the calculation of the polarization function is ``HILBERT``, which defines a (Hilbert) mesh for the frequency integration in equation :eq:`eqP`. (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. |
@@ -441,7 +441,7 @@ The most important keyword in the calculation of the polarization function is ``
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 :math:`k` becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}-\epsilon_{\mathbf{q}n})^{-1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the Γ point (k=0), Spex multiplies the linear term with this energy difference, resulting in smooth and 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``.
+(*) In the limit :math:`{k\rightarrow 0}`, the projections in the numerator of equation :eq:`eqP` approach linearly to zero. However, when calculating the dielectric function, one has to multiply with :math:`{\sqrt{4\pi}/k}` (square root of Coulomb matrix) in this limit. So, the first order of :math:`{\langle e^{i\mathbf{kr}} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle}` (corresponding to :math:`{\mu=1}`) in :math:`k` becomes relevant. Using k·p perturbation theory, one can show that the linear term is proportional to :math:`{(\epsilon_{\mathbf{q}n'}-\epsilon_{\mathbf{q}n})^{-1}}`. This can lead to numerical problems if the two energies are very close to each other. Therefore, when treating the :math:`{\Gamma}`; 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 equation :eq:`eqP` 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. |
@@ -471,7 +471,7 @@ be treated analytically and, as long as the Fermi surface is sufficiently big, i
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.
+(*) Mathematically, the parameter :math:`{\eta}` in :eq:`eqP` 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)}`. |