tutorials.rst 87 KB
 anoop chandran committed Dec 16, 2018 1 2 3 4 5 .. _tutorials: =================== Tutorials ===================  anoop chandran committed Dec 16, 2018 6 7 In this section, we provide detailed tutorials for the main features of the Spex code.  Christoph Friedrich committed Jan 11, 2019 8 9 10 *GW* calculations ================== The Spex code was started as a pure *GW* code. The *GW* approach remains the main computational method implemented in the code.  anoop chandran committed Dec 16, 2018 11   Christoph Friedrich committed Jan 11, 2019 12 13 14 15 One-shot *GW* approach ----------------------- A simple input file for a one-shot *GW* calculation for bulk silicon was already presented in :numref:getting_started. The following, equivalent input file is even more minimalistic.  anoop chandran committed Dec 16, 2018 16 17 18 19 20 21  .. code-block:: bash BZ 4 4 4 JOB GW 1:(1,2,5) 7:(1,3,5) 3:(1-3,5)  Christoph Friedrich committed Jan 11, 2019 22 23 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  Christoph Friedrich committed Jan 30, 2019 24 for the bands (1,2,5) at the first, (1,3,5) at the seventh, and (1,2,3,5) at the third k point. The band indices are defined  Christoph Friedrich committed Jan 11, 2019 25 26 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.  anoop chandran committed Dec 16, 2018 27   Christoph Friedrich committed Jan 11, 2019 28 29 30 In the case of a spin-polarized 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).  anoop chandran committed Dec 16, 2018 31   Christoph Friedrich committed Jan 30, 2019 32 33 34 35 36 The same definition of bands is used for the job types HF, PBE0, SX, COSX, GT, GWT, and KS. For example, JOB COSX 1:(1,2,5) would run a COHSEX calculation. Many of the keywords described in this section work identically for the other job types, others can only be used for *GW* calculations. If a keyword is inapplicable to the defined job, it is simply ignored by Spex (e.g., SPECTRAL for a PBE0 calculation).  anoop chandran committed Dec 16, 2018 37 38 .. _kpt:  Christoph Friedrich committed Jan 11, 2019 39 KPT  anoop chandran committed Dec 16, 2018 40 ----  Christoph Friedrich committed Jan 11, 2019 41 42 43 44 45 46 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 :math:{\Gamma} (index 1), X (index 7), and L point (index 3). Clearly, for a different k-point set, e.g., :math:8\times 8\times 8, the k-point indices would change (except for the index 1, which always corresponds to the :math:{\Gamma} point). Therefore, there is the possibility of defining k-point labels by  anoop chandran committed Dec 16, 2018 47 48  KPT X=1/2*(0,1,1) L=1/2*(0,0,1)  Christoph Friedrich committed Jan 11, 2019 49 or  anoop chandran committed Dec 16, 2018 50 51 KPT X=[1,0,0] L=1/2*[1,1,-1]  Christoph Friedrich committed Jan 11, 2019 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 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 :math:a. In the case of the fcc lattice of silicon, the lattice basis vectors are :math:{a_1=a(0,1,1)/2} , :math:a_2=a(1,0,1)/2 , and :math:a_3=a(1,1,0)/2 . The reciprocal basis vectors are thus :math:{b_1=2\pi/a(-1,1,1)} et cetera. For the X point, e.g., one then gets :math:{(a_2+a_3)/2=2\pi/a(1,0,0)} . In order for the cartesian coordinates (square brackets) to be interpreted correctly, Spex obviously needs to know the lattice constant :math: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 described in :numref:getting_started. 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 Wannier interpolation of *GW* quasiparticle energies or a self-consistent *GW* calculation is to be performed, which require the self-energy to be calculated for the whole irreducible Brillouin zone.  anoop chandran committed Dec 16, 2018 70 71 72 73 74 75 76 77 78 79 80 81 82 83  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,  Christoph Friedrich committed Jan 11, 2019 84 85 86 87 88 89 90 * vxc expectation value of the exchange-correlation potential :math:v^\mathrm{xc}, e.g., -12.15897 eV, * sigmax expectation value of the exchange self-energy :math:\Sigma^x, e.g., -19.20415 eV, * sigmac expectation value of the correlation self-energy :math:\Sigma^c, e.g., (6.60649+1.05993\ *i*\ ) eV (note that the self-energy is complex), * Z renormalization factor :math:Z, e.g., (0.65146-0.10556\ *i*\ ) 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.73681\ *i*\ ) eV (linearized solution) and (-6.36806+0.72704\ *i*\ ) eV (direct solution).  anoop chandran committed Dec 16, 2018 91   Christoph Friedrich committed Jan 11, 2019 92 The two *GW* values for the quasiparticle energy follow two common methods to approximately solve the quasiparticle equation  anoop chandran committed Dec 16, 2018 93   Anoop Chandran committed Dec 18, 2018 94 .. math::  Christoph Friedrich committed Jan 11, 2019 95  \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})  Anoop Chandran committed Dec 18, 2018 96  :label: qpeq  anoop chandran committed Dec 16, 2018 97   Christoph Friedrich committed Jan 11, 2019 98 99 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  anoop chandran committed Dec 16, 2018 100   Anoop Chandran committed Dec 18, 2018 101 .. math::  Christoph Friedrich committed Jan 11, 2019 102  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  Anoop Chandran committed Dec 18, 2018 103  :label: qppert  anoop chandran committed Dec 16, 2018 104   Christoph Friedrich committed Jan 11, 2019 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 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). (The keyword FULL is explained in :numref:beyondpert.) 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. .. _SPECTRAL: SPECTRAL (SENERGY)  anoop chandran committed Dec 16, 2018 120 121 ------------------  Christoph Friedrich committed Jan 11, 2019 122 123 (*) 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  anoop chandran committed Dec 16, 2018 124   Anoop Chandran committed Dec 18, 2018 125 .. math::  Christoph Friedrich committed Jan 11, 2019 126  G(\omega)=G_0(\omega)+G_0(\omega)[\Sigma^{\mathrm{xc}}(\omega)-v^{\mathrm{xc}}]G(\omega)  anoop chandran committed Dec 16, 2018 127   Christoph Friedrich committed Jan 11, 2019 128 129 130 which links the interacting Green function :math:G to the non-interacting KS one :math:{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  anoop chandran committed Dec 16, 2018 131   Christoph Friedrich committed Jan 11, 2019 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 .. 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 keyword SPECTRAL in the section SENERGY. 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 (third 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. The SPECTRAL keyword can be used in conjunction with band-structure calculations, see :numref:GWbandstr for details. .. list-table:: Examples :widths: 95 100 * - SPECTRAL - Write spectral function :math:{\text{Im}G(\omega)} to the file "spectral"; frequency mesh automatically defined. * - SPECTRAL {-10eV:1eV,0.01eV} - Write spectral function on the specified frequency mesh. * - SPECTRAL {-10eV:1eV,0.01eV} 0.002eV - Bound imaginary part from below by 0.002 eV, preventing peak widths to become too small to be plotted. *GW* basics ----------------------- 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. To lay the basis for the subsequent sections, we briefly recapitulate the basics of the *GW* method. The *GW* self-energy .. math:: \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'\,. :label: selfene can be understood as a scattering potential that contains the exact exchange potential and correlation effects through the inclusion of :math: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:: 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:: 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 mixed product basis (:numref:mpb), 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 :numref: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 Eq. :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, :math: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 in the badly controlled extrapolation of the Padé 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. This method requires the screened interaction :math:W to be calculated in a small range on the real frequency axis, another contribution comes from an integration along the imaginary frequency axis. As a third alternative, there is a *hybrid* method that is nearly as accurate as the "explicit" contour integration method, but it is faster and avoids the calculation of :math:W for real frequencies altogether.  anoop chandran committed Dec 16, 2018 216 217 218 219  MESH (SENERGY) --------------  Christoph Friedrich committed Jan 11, 2019 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 All methods to calculate the self-energy 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 dense 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). .. list-table:: Examples :widths: 28 100 * - MESH 12 15.0 - Use a mesh containing twelve frequencies for the range :math:[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".  anoop chandran committed Dec 16, 2018 242 243 CONTINUE (SENERGY) -------------------  Christoph Friedrich committed Jan 11, 2019 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 The keyword CONTINUE (section SENERGY) chooses the analytic continuation method. It can have optional parameters. If given without parameters (or if not specified at all), Padé 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, and *. Any combination is possible. They are explained in the examples. .. list-table:: Examples :widths: 22 100 * - CONTINUE - Use Padé 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 Padé method.)  anoop chandran committed Dec 16, 2018 264 265  CONTOUR (SENERGY)  Christoph Friedrich committed Jan 11, 2019 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 ----------------- 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 best 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,\epsilon_\mathrm{F}-\omega] 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 use. 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 [Eq. :eq:qppert]). 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} htr and :math:{\epsilon_{\mathbf{k}n}+0.01} htr). The more accurate ''direct'' solution is only available if we specify a range of frequencies for the self-energy instead of a single number, making an interpolation of the self-energy beyond a linear function possible.  Christoph Friedrich committed Mar 24, 2019 289 290 This is possible by an argument such as +-{-0.1:0.15,0.01}. Here, the range of values is relative to :math:\epsilon_{\mathbf{k}n}>\epsilon_\mathrm{F}. The range is reversed (to -0.15:0.1 in the example) for occupied states (:math:\epsilon_{\mathbf{k}n}<\epsilon_\mathrm{F})  Christoph Friedrich committed Jan 11, 2019 291 292 to reflect the fact that occupied and unoccupied states tend to shift in opposite directions by the renormalization. (Thus, it makes sense to choose a frequency range shifted to larger energies instead of a symmetrical one.)  Christoph Friedrich committed Mar 24, 2019 293 One can also specify an absolute frequency mesh by {...}, where the energies are given relative to the Fermi energy.  Christoph Friedrich committed Jan 11, 2019 294 This is mandatory for calculations with FULL (:numref:beyondpert).  Christoph Friedrich committed Mar 24, 2019 295 It is sometimes a bit inconvenient to determine suitable values for the upper and lower bound of {...}.  Christoph Friedrich committed Jan 11, 2019 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 Therefore, Spex allows the usage of wildcards for one of the boundaries or both (see examples). 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 determined automatically from the first argument. As a third method, Spex also allows the second argument to be omitted altogether. Then, it uses a *hybrid* method where the screened interaction is analytically continued from the imaginary axis (where it has to be known for the evaluation of 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 Padé 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. .. list-table:: Examples :widths: 70 100 * - CONTOUR 0.01 0.005 - Use contour integration to obtain the two values :math:{\Sigma^\mathrm{xc}(\epsilon\pm 0.01\,\mathrm{htr})} with the KS energy :math:\epsilon, giving :math:{\Sigma^\mathrm{xc}} as a linear function.  Christoph Friedrich committed Mar 24, 2019 311  * - CONTOUR +-{-0.1:0.15,0.01} 0.005  Christoph Friedrich committed Jan 11, 2019 312  - Calculate :math:{\Sigma^\mathrm{xc}(\omega)} on a frequency mesh relative to the KS energy :math:{\epsilon}.  Christoph Friedrich committed Mar 24, 2019 313  * - CONTOUR {*:*,0.01} 0.005  Christoph Friedrich committed Jan 11, 2019 314  - Use an absolute frequency mesh (relative to the Fermi energy) instead. Wildcards are used for the upper and lower bounds.  Christoph Friedrich committed Mar 24, 2019 315  * - CONTOUR {*:*,0.01}  Christoph Friedrich committed Jan 11, 2019 316 317  - Use *hybrid* method.  Christoph Friedrich committed Mar 24, 2019 318 319 .. note:: Until version 05.00pre31, the syntax was {...} for +-{...} and [{...}] for {...}.  Christoph Friedrich committed Jan 11, 2019 320 321 322 323 324 325 326 327 328 FREQINT (SENERGY) ------------------ (*) 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_n} or :math:{z=\omega_n}. 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 the two arguments SPLINE (default) and PADE for spline [after the transformation :math:{\omega'\rightarrow \omega'/(1+\omega')}] and Padé interpolation, respectively. In the case of *GT* calculations, there is a similar frequency integration with the *T* matrix replacing *W*.  Christoph Friedrich committed Mar 07, 2019 329 330 331 332 333 334 335 An experimental option is NONLIN, which invokes a new tetrahedron k integration method that uses weight functions instead of weight factors. This enables smooth integrations over strongly varying (highly nonlinear) functions, which is particularly useful for the magnetic *T* matrix with its very sharp spin-wave peaks. Therefore, the default for *GT* calculations is NONLIN. FREQINT NONLIN, in particular, affects the calculation of the residues part in a CONTOUR calculation but also the integration along the imaginary frequency axis (both for CONTINUE and CONTOUR), where it uses Padé interpolation. We note that FREQINT NONLIN and, to a lesser extent, also FREQINT PADE tend to break energy degeneracies slightly.  Christoph Friedrich committed Jan 11, 2019 336 337 338 339 340  .. list-table:: Examples :widths: 28 100 * - FREQINT SPLINE  Christoph Friedrich committed Mar 07, 2019 341  - Use spline interpolation for *W* (or *T*) in the frequency convolution :math:{G\cdot W} (:math:{G\cdot T}) (default for *GW*).  Christoph Friedrich committed Jan 11, 2019 342 343  * - FREQINT PADE - Use Padé interpolation.  Christoph Friedrich committed Mar 07, 2019 344 345 346  * - FREQINT NONLIN - Use new tetrahedron k integration method (default for *GT*).  Christoph Friedrich committed Jan 11, 2019 347 .. _SMOOTH:  anoop chandran committed Dec 16, 2018 348 349 350  SMOOTH (SENERGY) -----------------  Christoph Friedrich committed Jan 11, 2019 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 (*) We have already mentioned that the Padé approximation can lead to spurious features in the extrapolated or interpolated quantity, i.e., the self-energy or the screened interaction (also see FREQINT), if one of the Padé 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 Padé approximants with randomly shifted frequency points according to :math:{\omega_n\rightarrow\omega_n(1+10^{-7}r_n)} with random numbers :math:{r_n\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 Padé approximants over which to average, 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. .. list-table:: Examples :widths: 28 100 * - SMOOTH 500 - Average over 500 Padé 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.)  Christoph Friedrich committed Mar 24, 2019 371 ZERO (SENERGY)  anoop chandran committed Dec 16, 2018 372 ----------------  Christoph Friedrich committed Mar 24, 2019 373 (*) Both the exchange (HF) and the correlation self-energy contain terms of the form  Christoph Friedrich committed Jan 11, 2019 374 375 376 377 378 379 380 381 382 383 384 385 :math:\langle...\rangle\langle...\rangle V(\mathbf{k}) with :math:{V=v} and :math:{V=W}, respectively. The quantities :math:{\langle...\rangle} are the "projections" discussed in :numref:mpb. One of the involved projections acquires the form :math:\langle e^{i\mathbf{k}r} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle, which behaves as :math:{\propto k} in the long-wavelength limit :math:{k\rightarrow 0} and :math:n\neq n'. The prefactor is obtained from :math:k\cdot p perturbation theory. Multiplying two of these projections with :math:V\propto k^{-2} (or one with :math:V\propto k^{-1}, which is valid for the "wings" of *W*) gives rise to zero-order terms. 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 such systems), impeding a favourable convergence with respect to the k-point sampling.  Christoph Friedrich committed Mar 24, 2019 386 387 388 389 The keyword ZERO (section SENERGY) includes the zero-order terms. They are neglected by default. .. note:: Keyword ZERO replaces NOZERO of older versions (until 05.00pre31), which switched off zero-order corrections. So, starting from version 05.00pre32, the default has changed to excluding zero-order corrections.  anoop chandran committed Dec 16, 2018 390 391 392  ALIGNVXC (SENERGY) ------------------  Christoph Friedrich committed Jan 11, 2019 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 (*) Equation :eq:qpeq depends explicitly on the mean-field starting point through the eigensolutions :math:{\{\phi_{\mathbf{k}n}\}}. This dependence is more obvious from Eq. :eq:qppert, 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 :math:N and the :math:N{\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. This is useful if a different criterion is to be used, for example, requiring the Fermi energy (metallic case) or core state energies (core spectroscopy) to remain unchanged. In this case, the shift has to be determined manually by repeatedly trying different shift values until the constancy of the desired quantity is achieved (because the shift affects the results, so there is a feedback effect). In this trial-and-error approach, it is very helpful to use RESTART 2 (:numref:RESTART_GW), in which case the self-energy is read from harddisc and the calculation takes very little time. .. list-table:: Examples :widths: 28 100 * - 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.  anoop chandran committed Dec 16, 2018 418 419 420  VXC (SENERGY) --------------  Christoph Friedrich committed Jan 11, 2019 421 422 423 424 425 426 427 428 429 430 (*) The matrix elements of the exchange-correlation potential are needed in the solution of Eq. :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 (default) and CALC. .. list-table:: Examples :widths: 16 100 * - 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.  anoop chandran committed Dec 16, 2018 431   Christoph Friedrich committed Jan 11, 2019 432 .. _RESTART_GW:  anoop chandran committed Dec 16, 2018 433 434 435  RESTART ---------  Christoph Friedrich committed Jan 11, 2019 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 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 :math:W(\mathbf{k}) and (b) updates the self-energy matrix (or expectation values) :math:{\Sigma^\mathrm{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 :math:W(\mathbf{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 :math:W(\mathbf{k}) from a previous run. The matrix :math:W(\mathbf{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 (:numref:GWbandstr). 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, the argument 2 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 about how the quasiparticle equation is solved after completion of the self-energy calculation. The following logical table gives an overview. .. list-table:: :widths: 30 30 30 :header-rows: 1 * - - 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  anoop chandran committed Dec 16, 2018 468 469  * spex.mb: contains MT part of product basis,  Christoph Friedrich committed Jan 11, 2019 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 * spex.ccou: contains "contracted" screened interaction (i.e., summed over **k**); this contraction is needed for the self-energy core contribution (CORES), for COSX, and for IBC, * spex.core: contains core contribution to *W* (keyword CORES, :numref:CORES_GW). The RESTART option can be exploited to customize the calculations. For example, it is possible to perform a self-consistent QSGW calculation 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, changing input parameters in "spex.inp" and running a calculation with RESTART might lead to a conflict between the changed parameters and the data read from the restart 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 just 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 screened interaction *W*, for example, NBAND. (The latter allows the band convergence of *W* and :math:{\Sigma^\mathrm{xc}} to be checked separately.) Again, this can be used efficiently to customize a calculation. .. list-table:: Examples :widths: 18 100 * - 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. .. _GWbandstr: *GW* band structures --------------------- In this section, three ways of calculating *GW* band structures are discussed. * Using the KPTPATH label (single Spex run, :numref:KPTPATH), * With the special + label in KPT for additional q points (multiple Spex runs, :numref:KPTadd), * Wannier interpolation (single -- but long -- Spex run, :numref:INTERPOL_GW). 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 q point. So, in addition to :math:\mathbf{q}, one would need to include all shifted points :math:\mathbf{q}+\mathbf{k}, as well, where :math:\mathbf{k} stands for the k points defined in BZ. If we restrict ourselves to the set of k points defined in BZ, the band structure along a certain high-symmetry line (e.g., :math:{\Gamma-X}) can consist only of the few k points that happen to lie on this line. Sometimes this is already sufficient. For example, if the *GW* renormalization affects the band dispersion only little, then a few points along the high-symmetry line are enough to shift/modify the KS bands accordingly and produce, in this way, a quasiparticle band structure. Such a calculation can be performed in a single Spex run. .. _KPTPATH:  anoop chandran committed Dec 16, 2018 518 519 520  KPTPATH --------  Christoph Friedrich committed Jan 30, 2019 521 522 523 524 525 526 527 528  .. .. figure:: figures/band1.png (Figures are disabled because Sphinx does not place them correctly.) :scale: 55% :figwidth: 50% :align: right :alt: band structure Si (1) Band structure for bulk silicon from DFT (solid line) and *GW* (symbols) calculated with JOB GW PATH:(1-12).  529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 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 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 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. The keyword KPTPATH can also be used for the other two methods to calculate a band structure (see :numref:KPTadd and :numref:INTERPOL_GW). In these methods, one is not restricted to the k points of the k set and can define k-point paths with arbitrary step sizes. To this end, an integer argument can be specified at the end of the k-path definition. For example, KPTPATH (L,1,X),50 would give a k path with a step size of (roughly) :math:2\pi/(50\sqrt[3]{\Omega}). The k-points given (L, :math:\Gamma, and X in the example) must lie on the path. The actual step size is adjusted accordingly. Here, :math:\Omega is the unit-cell volume, so :math:8\pi^3/\Omega is the volume of the Brillouin zone and :math:2\pi/\sqrt[3]{\Omega} the average reciprocal lattice constant. Larger integer arguments give denser k meshes. The default is 20 for the method explained in :numref:KPTadd and 100 for the one in :numref:INTERPOL_GW.  Christoph Friedrich committed Jan 11, 2019 549 550 551 552 553 554 555 556  .. list-table:: Examples :widths: 36 100 * - KPTPATH (L,1,X) - Define k-point path from L over the :math:\Gamma point (k index is always 1) to X. The labels L and X must be defined with KPT. * - JOB GW PATH:(1-10) - Run *GW* calculation for the first ten bands at all k points defined by KPTPATH.  557 558  * - KPTPATH (1,N),100 - Define k-point path from :math:\Gamma to N and use step size of :math:2\pi/(100\sqrt[3]{\Omega}).  Christoph Friedrich committed Jan 11, 2019 559 560  .. _KPTadd:  anoop chandran committed Dec 16, 2018 561 562 563  KPT +=... ---------  Christoph Friedrich committed Jan 11, 2019 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 This section explains an extension to the keyword KPT (:numref:KPT), which will later enable the calculation of smooth band structures. The extension enables 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 - \mathrm{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. As explained above, each additional point :math:\mathbf{q} must in principle be accompanied with all shifted points :math:\mathbf{q}+\mathbf{k}. So, before running Spex, we must let the DFT code generate the wavefunctions and energies at the shifted k-point set :math:\{\mathbf{q}+\mathbf{k}\}. This is done in the usual way by creating the k-point file (:numref:WRTKPT) 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 :math:\Gamma point and the additional q point. The energy difference between the fourth state at :math:\Gamma and the fifth state at :math:\mathbf{q} yields the fundamental *GW* gap. .. list-table:: Examples :widths: 36 100 * - KPT +=1/7*[1,0,0] - Define special k point at 1/7 along the line :math:{\Gamma - \mathrm{X}} (here, the X point in *x* direction) of an fcc structure. * - KPT +=1/14*(0,1,1) - The same in internal coordinates. * - JOB GW +:(1-10) - Run *GW* calculation for the states 1-10 at the added q point.  Christoph Friedrich committed Jan 30, 2019 587 588 589 590 591 592 593 .. .. figure:: figures/band2.png :scale: 55% :figwidth: 50% :align: right :alt: band structure Si (2) Band structure for bulk silicon from DFT (magenta) and *GW* (green) calculated with multiple runs of Spex using the "spex.band" shell script.  Christoph Friedrich committed Jan 11, 2019 594 595 596 597 598 599 600 601 602 603 604 605 606 607 The possibility of adding arbitrary q points enables the calculation of smooth band structures. To this end, a *GW* run (together with the generation of the eigenstates with the DFT code) has to be performed for each q point in a list of q points that make up a path in the Brillouin zone, e.g., from :math:{\Gamma} to :math:X. This list is defined with the keyword KPTPATH as explained above. Fortunately, the screened interaction has to be calculated only once and can be reused for the other points. For this, the keyword RESTART should be given in the input file. Of course, we also need a job definition, such as JOB GW +:(1-10). Here, only the additional q point should be specified. It would be difficult to extract the data for the band structure plot otherwise.) Now, running Spex with KPTPATH defined in the input file and WRTKPT (or with the command line 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 "spex.selfc" (:numref:beyondpert) and produces, for each q point, one output file named "spex_NNN.out", where NNN is a three-digit counting index, or more digits if the index exceeds 999.  Christoph Friedrich committed Jan 30, 2019 608 609 610 611 612 613 614 .. .. figure:: figures/band3.png :scale: 55% :figwidth: 50% :align: right :alt: band structure Si (3) Band structure for bulk silicon from DFT (magenta) and *GW* (green) with lifetime broadening.  Christoph Friedrich committed Jan 11, 2019 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 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. If "spex.inp" contains a SPECTRAL definition (:numref:SPECTRAL), the script "spex.band" also renames the file "spectral", written in each Spex run, to "spectralNNN" with the same counting index NNN as above. The spectral data can be compiled with the "spex.extr" utility (e.g., spex.extr s -b spectral??? > spec) into a single data file containing a list of xyz value triples: the Bloch momenta as the first column, the frequency as the second, and the value of the spectral function as the third. With "gnuplot", for example, a momentum-frequency plot of the spectral function with color coding is then prepared with the command plot "spec" using 1:2:3 palette with lines. (The respective color range is set with set cbrange [...], see the Gnuplot manual.) .. _INTERPOL_GW: INTERPOL --------- As a third alternative, we can use Wannier interpolation to interpolate between the few k points of the original k mesh.  637 638 The construction of Wannier orbitals is explained in :numref:wannier. Here, we use a minimalistic definition for demonstration purposes. We have to modify and add some lines to spex.inp:  anoop chandran committed Dec 16, 2018 639 640 641 642 643 644 645 646 647  .. code-block:: bash JOB GW IBZ:(1-4) SECTION WANNIER ORBITALS 1 4 (sp3) MAXIMIZE INTERPOL END  anoop chandran committed Dec 16, 2018 648   Christoph Friedrich committed Jan 30, 2019 649 650 651 652 653 654 655 .. .. figure:: figures/band4.png :scale: 60% :figwidth: 55% :align: right :alt: band structure Si (4) Band structure for bulk silicon from DFT (magenta) and *GW* (green) calculated with Wannier interpolation (only filled bands).  Christoph Friedrich committed Jan 11, 2019 656 657 658 659 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.)  660 661 662 663 664 665 The section WANNIER defines a set of maximally localized Wannier functions (MLWFs) of bonding sp\ :sup: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.  Christoph Friedrich committed Jan 11, 2019 666 667 You can plot the band structure with the line thickness scaled by the imaginary part in the same way as described before.  668 669 670 671 The band structure is calculated along the k-point path defined by KPTPATH (:numref:KPTPATH). By default, Spex uses a step size of (roughly) :math:2\pi/(100\sqrt[3]{\Omega}). Alternatively, one can specify a file after INTERPOL, e.g., INTERPOL qpts. The k-point path is then taken from that file. (The file structure is the same as of "kpts" or "qpts".)  Christoph Friedrich committed Jan 11, 2019 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 Wannier interpolation can also be used in conjunction with the keyword SPECTRAL (:numref: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 the spin up and spin down channels, respectively. These files can be plotted in the same way as described in :numref:KPTadd. .. _CORES_GW: CORES ------ By default, the *n* summation in the Green function .. math:: G_0(\mathbf{r},\mathbf{r}';\omega)=\sum_\mathbf{k}\sum_n \frac{\phi_{\mathbf{k}n}(\mathbf{r})\phi^*_{\mathbf{k}n}(\mathbf{r}')}{\omega-\epsilon_{\mathbf{k}n}+ i\eta\,\mathrm{sgn}(\epsilon_{\mathbf{k}n}-\epsilon_\mathrm{F})} :label: lehmann for the self-energy [Eq. :eq:selfene] includes only the valence and conduction states represented in the LAPW basis but excludes the core states. In theory, of course, the core states should be included. Their contribution is, however, numerically small. Spex allows selective core states to be included in the evaluation of the self-energy with the keyword CORES. This keyword also affects the calculation of the polarization function, see :numref:CORES_P. As arguments, Spex expects a comma separated list of core states between round brackets, e.g., (2s,2p,3d), for each atom type. *Empty* definitions () are allowed. If CORESOC is specified, one can also address the SOC-splitted states individually, e.g., (2p1/2). .. list-table:: Examples :widths: 32 100 * - CORES () (1s) - Include the 1s state of the second atom type. * - CORES (2s,2p1/2) - Include the 2s and the 2p\ :sup:1/2 states but exclude the 2p\ :sup:3/2 states. .. _beyondpert: Beyond perturbative *GW* and QSGW ---------------------------------- One can also go beyond the perturbative solution of Eq. :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 :math:10\times 10 matrices at the :math:\Gamma, X, and L points. (In principle, it is also possible to have a band definition like (1-5,9,10), which would yield a :math:7\times 7 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::  anoop chandran committed Dec 16, 2018 715 716 717 718 719 720 721  Bd KS olap HF olap GW 3 2.92306 0.99971 0.12683 0.99985 2.50496 -0.06006 4 2.92306 0.99971 0.12683 0.99985 2.50496 -0.06006 9 16.77136 0.99732 23.25660 0.99962 16.97051 -0.30294 10 16.77136 0.99732 23.25660 0.99962 16.97051 -0.30294  Christoph Friedrich committed Jan 11, 2019 722 723 724 725 726 727 The values with the label HF are the Hartree-Fock values (from diagonalization). They are **not** necessarily ordered according to energy but related to the KS states with a condition of maximal overlap (and at the same time making sure that a one-to-one correspondence can be established). The corresponding overlap is given on the left of the HF column. The full *GW* energies (column labeled GW; the last column lists the imaginary parts) 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 wavefunction 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.  anoop chandran committed Dec 16, 2018 728   Christoph Friedrich committed Jan 11, 2019 729 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 [Phys. Rev. Lett. 93, 126406]:  anoop chandran committed Dec 16, 2018 730   Anoop Chandran committed Dec 18, 2018 731 732 .. 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  anoop chandran committed Dec 16, 2018 733 734 735  for diagonal elements and  Christoph Friedrich committed Jan 11, 2019 736 .. math::  Anoop Chandran committed Dec 18, 2018 737  \displaystyle A_{\mathbf{k}nn'}=\langle \phi_{\mathbf{k}n} | \Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n})+\Sigma^\mathrm{xc}(\epsilon_{\mathbf{k}n'}) | \phi_{\mathbf{k}n'} \rangle  anoop chandran committed Dec 16, 2018 738   Christoph Friedrich committed Jan 11, 2019 739 for off-diagonal elements. The *hermitianized* QSGW operator is then obtained from :math:{\Sigma^\mathrm{xc,H}=(A+A^\dagger)/2}. The difference to the original definition is the inclusion of the renormalization factor to better reproduce the *GW* quasiparticle energies. The *hermitianized* matrix, or rather the difference :math:{\Sigma^{\mathrm{xc,H}}-v^\mathrm{xc}}, is written to the file "spex.qsgw", which is later read by the DFT code. In Fleur, the following steps are required:  anoop chandran committed Dec 16, 2018 740   Christoph Friedrich committed Jan 11, 2019 741 * rm fleur.qsgw - remove the file "fleur.qsgw" containing the *hermitianized* matrix used internally by Fleur from previous QSGW iterations.  anoop chandran committed Dec 16, 2018 742 743 744 * rm broyd* - remove Broyden information about previous iterations because this information is inconsistent with the new Hamiltonian (the SCF calculation does not converge otherwise). * Set gw=3 in the Fleur input file. * Run Fleur.  Christoph Friedrich committed Jan 11, 2019 745  Fleur translates the file "spex.qsgw" (in basis of eigenstates) to a file "fleur.qsgw" (in LAPW basis). Then, a SCF run is performed where, instead of the KS Hamiltonian :math:{H^\mathrm{KS}}, the Hamiltonian :math:{H^\mathrm{KS}+(\Sigma^{\mathrm{xc,H}}-v^\mathrm{xc})} is employed. After the run finishes, the gw=2 step has to be performed.  anoop chandran committed Dec 16, 2018 746 747 748 * Set gw=2 in the Fleur input file. * Run Fleur.  Christoph Friedrich committed Jan 11, 2019 749 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:  anoop chandran committed Dec 16, 2018 750 751 752 753 754 755 756  * SPEX_EXEC - Spex executable, * FLEUR_EXEC - Fleur executable, * SPEX_QUEUE - Queue command for parallel execution of Spex (e.g., mpiexec -np 20), * FLEUR_QUEUE - Same for Fleur, * QUEUE - Default queue command to be used for Spex and Fleur.  Christoph Friedrich committed Jan 11, 2019 757 The job types HF, PBE0, SX, and COSX also allow FULL in the definition, i.e., a full matrix is evaluated, enabling a self-consistent calculation in the same way as QSGW (see above). The file names used are the same: "spex.qsgw", "fleur.qsgw", "qsgw", even though it might actually be a HF calculation. The job GT (also GWT) does not allow for self-consistency at the moment.  anoop chandran committed Dec 16, 2018 758 759   anoop chandran committed Dec 16, 2018 760 761 .. _polar:  anoop chandran committed Dec 16, 2018 762 763 Polarization function =====================  Christoph Friedrich committed Jan 11, 2019 764 The polarization function gives the linear change in the electronic density of a non-interacting system with respect to changes in the effective potential. It is, thus, a fundamental quantity in the calculation of screening properties of a many-electron systems. For example, the dielectric function, instrumental in the calculation of spectroscopic quantities (e.g. EELS) and the screened interaction needed in *GW*, is related to the polarization matrix through :math:{\varepsilon(\mathbf{k},\omega)=1-P(\mathbf{k},\omega)v(\mathbf{k})}, here given in matrix notation. The corresponding explicit formula for matrix elements of *P* in the mixed product basis is  anoop chandran committed Dec 16, 2018 765   Anoop Chandran committed Dec 18, 2018 766 .. math::  Christoph Friedrich committed Jan 11, 2019 767 768 769  P_{\mu\nu}(\mathbf{k},\omega)&=2\sum_{\mathbf{q}}^{\mathrm{BZ}}\sum_{n}^{\mathrm{occ}}\sum_{n'}^{\mathrm{unocc}}\langle M_{\mathbf{k}\mu} \phi_{\mathbf{q}n} | \phi_{\mathbf{k+q}n'} \rangle\langle \phi_{\mathbf{k+q}n'} | \phi_{\mathbf{q}n} M_{\mathbf{k}\nu} \rangle \\ & \cdot \left(\frac{1}{\omega+\epsilon_{\mathbf{q}n}-\epsilon_{\mathbf{q}+\mathbf{k}n'}+i\eta}-\frac{1}{\omega-\epsilon_{\mathbf{q}n}+\epsilon_{\mathbf{q}+\mathbf{k}n'}-i\eta}\right) \\ & = \int_{-\infty}^\infty \frac{S_{\mu\nu}(\mathbf{k},\omega')}{\omega-\omega'+i\eta\mathrm{sgn}(\omega')}d\omega'\,.  Anoop Chandran committed Dec 18, 2018 770  :label: eqP  anoop chandran committed Dec 16, 2018 771   Christoph Friedrich committed Jan 11, 2019 772 773 774 775 776 A system without spin polarization has been assumed, hence the factor 2. In case of spin polarization, the formula has a spin summation instead of the factor 2 and corresponding spin quantum numbers. In the case of spin-orbit coupling (and no spin polarization), two spin summation labels are added, one in front of each "projection" :math:\langle...\rangle. We have implicitly defined the spectral function *S* in Eq. :eq:eqP, 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.)  anoop chandran committed Dec 16, 2018 777   Christoph Friedrich committed Jan 11, 2019 778 The expression for Eq. :eq:eqP involves an infinite empty-band summation over :math:{n'}. In practice, this summation is truncated by the number of bands (keyword NBAND), which thus becomes an important convergence parameter. The self-energy contains an infinite band summation as well.  anoop chandran committed Dec 16, 2018 779   Christoph Friedrich committed Jan 11, 2019 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 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 be a discrete 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). 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. :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). .. list-table:: Examples :widths: 34 100 * - HILBERT NONE - Simple summation over k points. * - HILBERT 50 30.0 - Use Hilbert mesh with fifty frequencies (between 0 and 5 htr) and an accumulated stretching factor of 30 at 5 htr. * - HILBERT 0.01 1.05 - Use Hilbert mesh with a first step size of :math:{\omega_2-\omega_1=} 0.01 htr and a stretching factor of :math:{a=1.05}. (First argument is real-valued.) MULTDIFF (SUSCEP) ------------------- (*) In the limit :math:{k\rightarrow 0}, the projections in the numerator of Eq. :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 :math:k\cdot 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.  Christoph Friedrich committed Mar 07, 2019 804 805 806 807 808 809 810 Therefore, when treating the :math:{\Gamma} point (:math: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 frequency integration of Eq. :eq:eqP, thereby avoiding the numerical difficulties. By default, Spex does that only at :math:k=0. The behavior can be changed with the keyword MULTDIFF in the section SUSCEP. 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.  Christoph Friedrich committed Jan 11, 2019 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848  .. list-table:: Examples :widths: 28 100 * - MULTDIFF OFF - Never separate (divergent) energy differences. * - MULTDIFF ON - Always separate energy differences (i.e., for all k points). * - MULTDIFF INT - Use default behavior (separate for :math:k=0, do not for :math:k\neq 0) but stick energy differences into integration weights. * - MULTDIFF INTON - Always separate energy differences by sticking them into integration weights. .. _PLASMA: PLASMA (SUSCEP) --------------- (*) In the case of metals, the polarization function contains an additional term, the so-called Drude term, which gives rise to metallic screening (diverging static dielectric constant, short-range static *W*). The Drude term stems from virtual intraband transitions across the Fermi surface. It can be formulated with the square of the plasma frequency, which in turn is evaluated by an integration over the Fermi surface. The Drude term can be treated analytically and, as long as the Fermi surface is sufficiently big, it normally does not pose a numerical problem. However, a very small Fermi surface eventually leads to a very sharp "Drude peak" in the *GW* self-energy, impeding a straight-forward numerical solution of the non-linear quasiparticle equation. One could also say that, while the Drude term is actually treated correctly, it gets too much weight because of the coarseness of the k-point mesh. It is, therefore, sometimes helpful to modify the plasma frequency. This is possible with the keyword PLASMA. Its argument replaces the plasma frequency calculated by the code. Setting PLASMA 0 switches the Drude term off altogether. Of course, PLASMA 0 thus also disables metallic screening, which might be unphysical. Therefore, there is a special option, PLASMA METAL, which scales the head element of :math:{W(\mathbf{k},\omega)} in the limit :math:{k\rightarrow 0} to enforce metallic screening. This latter option can be helpful, for example, in the case of semimetallic systems with a tiny Fermi surface. .. list-table:: Examples :widths: 24 100 * - PLASMA 1.5eV - Set plasma frequency manually to 1.5 eV. * - PLASMA 0 - Set plasma frequency to zero. Disables Drude term. * - PLASMA METAL - Disable Drude term but enforce metallic screening by scaling the head element of *W*. .. _CORES_P: CORES ------ By default, the *n* summation of Eq. :eq:eqP extends only over the valence states (represented in the LAPW basis) and excludes the core states. The core states can be included in this summation in the same way as for the *GW* self-energy (see :numref:CORES_GW). If CORESOC is specified, the SOC-splitted core states are included, otherwise the SOC-averaged core states (see :numref:CORESOC). In the former case, the Weiss field of a spin-polarized system can further lift the degeneracies. This is taken into account as well. See :numref:CORES_GW for details and examples.  anoop chandran committed Dec 16, 2018 849 850 851  TETRAF (SUSCEP) -----------------  Christoph Friedrich committed Jan 11, 2019 852 (*) In rare cases, especially for very large k-point sets, the determination of the tetrahedron weights (Timing wghts) can become computationally expensive. This is because Spex calculates tetrahedron weights for all k points by default, which makes it possible to average over equivalent k points, thereby avoiding a (slight) symmetry breaking connected with the spatial form of the tetrahedra. (The tetrahedra are not unique!) One can disable this averaging by setting TETRAF in the section SUSCEP. If set, tetrahedron weights are only calculated for the irreducible part of the BZ. This accelerates the calculation of the weights but introduces slight deviations due to symmetry breaking.  anoop chandran committed Dec 16, 2018 853 854 855  WGHTTHR (WFPROD) -----------------  Christoph Friedrich committed Jan 11, 2019 856 (*) In both tetrahedron and Gaussian methods, integration weights are calculated to perform the k summation. To be more precise, an integration weight is calculated for each virtual transition, i.e., for each combination of band index pair (occ-nocc) and k point. One can reduce the number of terms by introducing a threshold value below which the respective weight is neglected. The corresponding keyword is WGHTTHR in the section SUSCEP. The default value is :math:{10^{-10}}. It is usually not necessary to change this value in practice.  anoop chandran committed Dec 16, 2018 857 858 859  GAUSS -------  Christoph Friedrich committed Jan 11, 2019 860 (*) The Gaussian method is included mostly for testing. It does not only affect the polarization function but also all other k summations (such as for the Hartree-Fock exchange potential and the determination of the Fermi energy). Therefore, GAUSS is a global keyword. It effectively replaces each single-particle eigenvalue by a Gaussian function with a certain width so that the density of states becomes a smooth function. In the calculation of the polarization function, one additionally needs a finite width for the Fermi edge. The keyword GAUSS needs the two width parameters as arguments given as real positive values. The advantage of the Gaussian method is its relative mathematical simplicity compared to the tetrahedron method and, in particular, that the Gaussian method cannot give rise to symmetry breaking.  anoop chandran committed Dec 16, 2018 861   Christoph Friedrich committed Jan 11, 2019 862 863 864 865 866 .. list-table:: Example :widths: 32 100 * - GAUSS 0.001 0.01 - Use Gaussian k-integration method with widths 0.001 and 0.01 htr. (Mostly used for testing.)  anoop chandran committed Dec 16, 2018 867   Christoph Friedrich committed Jan 11, 2019 868 DISORDER (SUSCEP)  anoop chandran committed Dec 16, 2018 869 -------------------  Christoph Friedrich committed Jan 11, 2019 870 (*) 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 this keyword has rarely been used so far.  anoop chandran committed Dec 16, 2018 871   Christoph Friedrich committed Jan 11, 2019 872 873 874 875 876 .. list-table:: Example :widths: 26 100 * - DISORDER 1000 - Use a finite value :math:{\eta=1/(2\cdot 1000)\,\mathrm{htr}}.  anoop chandran committed Dec 16, 2018 877   Christoph Friedrich committed Jan 11, 2019 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 .. _spectra: Spectra =============== The dielectric function :math:\epsilon_{\mathbf{GG}'}(\mathbf{k},\omega) is related to several spectroscopic experimental techniques where the system is perturbed (excited) by some incoming beam of particles but does not lose or gain particles (in contrast to photoemission, for example, where electrons are lost from the sample). In optical spectroscopy measurements, the absorption cross section is proportional to the imaginary part of the long-wavelength macroscopic dielectric function :math:\varepsilon_{M}(\omega), which can be obtained from the microscopic dielectric function by .. math:: \varepsilon_{M}(\omega)=\lim_{\mathbf{k}\rightarrow 0}\frac{1}{\varepsilon^{-1}_\mathbf{00}(\mathbf{k},\omega)}\,. Here, :math:\varepsilon^{-1} refers to the matrix inverse. Another example is electron-energy loss spectroscopy (EELS), where the scattering cross section can be shown to be proportional to the head element (:math:\mathbf{G}=\mathbf{G}'=0) of the inverse dielectric function .. math:: \varepsilon^{-1}_\mathbf{00}(\mathbf{k},\omega)\,. Spex can be used to calculate the dielectric function. The respective job is called DIELEC. Arguments are expected to define the k point index or label and a frequency range. For example, JOB DIELEC X:{0:1,0.001} would yield the dielectric function evaluated at the X point between 0 and 1 htr in steps of 0.001 htr. Two files are written, "dielec" and "dielecR", containing the head elements of the (bare) dielectric function :math:\varepsilon_\mathbf{00}(\mathbf{k},\omega) and of the renormalized dielectric function :math:1/\varepsilon^{-1}_\mathbf{00}(\mathbf{k},\omega). The imaginary part of the latter can directly be identified with an optical absorption spectrum (excluding excitonic effects). For the EELS spectrum, one would have to plot the imaginary part of the respective reciprocal value (for example, in "gnuplot": plot "dielecR" using 1:(imag(1/($2+{0,1}*$3)))). In the same way, other quantities can be calculated: the polarization function with SUSCEP (written to the file "suscep"), the renormalized polarization function (:math:R=P+PvR) with SUSCEPR (written to the file "suscepR"; "suscep" is also written), the screened interaction with  911 912 SCREEN (written to the file "screen"). In the first case (SUSCEP), one can also specify spin indices. For example, JOB SUSCEP X:ud{0:1,0.001} restricts the polarization function to virtual up\ :math:\rightarrow\ down transitions.  Christoph Friedrich committed Jan 11, 2019 913 914 915 916 917 918 919 920 921 922 Other valid spin labels are uu, dd, du, and + (spin summed, default). By default, only the head element is written to the output files. More elements can be written using, e.g., JOB SUSCEP 1:{0.5eV:2eV,0.001eV},OUT=4, which would yield :math:4\times 4 matrices instead of only the head element. One can also write the full matrix to a binary file with ...,OUT=BIN. The special + k-point label can be used as well (e.g., JOB SUSCEP +:{0..50eV,0.01eV}). This enables using the script "spex.band" to prepare a whole list of spectra for a series of Bloch vectors on a high-symmetry line. As explained in :numref:KPTadd for the file "spectral", the script appends a three-digit counting index to the name of the output file ("suscep" in the present example, which is thus renamed to "suscep001", "suscep002", ...). The data can then be compiled  923 into a single file using the "spex.extr" utility. (See :numref:KPTadd for details).  Christoph Friedrich committed Jan 11, 2019 924 925 926  Total xc energies ==================  927 Spex can calculate total exchange and correlation energies of the many-electron system. The total exchange energy is defined  Christoph Friedrich committed Jan 11, 2019 928 929 by the HF expression  930 .. math:: E^\mathrm{x}=\frac{1}{2} \sum_{\mathbf{k}n} \sum_{\mathbf{k}'n'} \iint \frac{ \phi^*_{\mathbf{k}n}(\mathbf{r}) \phi_{\mathbf{k}'n'}(\mathbf{r}) \phi^*_{\mathbf{k}'n'}(\mathbf{r}') \phi_{\mathbf{k}n}(\mathbf{r'}) } { \mathbf{r}-\mathbf{r}' } d^3r\,d^3r'\,.  Christoph Friedrich committed Jan 11, 2019 931 932 933 934 935 936 937 938 939 940 941 942  :label: hfene The corresponding job is called HFENE. Internally, it works very much like a HF calculation. So, all parameters pertaining to HF are available to HFENE, as well. However, the final result is a single number, :math:E^\mathrm{x}. The total correlation energy can be calculated using the adiabatic-connection fluctuation-dissipation theorem (ACFDT) with the random-phase approximation (RPA) for the (partially renormalized) response function :math:\chi_\lambda (:math:\lambda=[0,1]). This response function is identical to the polarization function :math:P (:numref:polar) for :math:\lambda=0 and to the "renormalized polarization function" :math:R mentioned in :numref:spectra for :math:\lambda=1. The ACFDT method gives the following form of the correlation energy functional  Christoph Friedrich committed Jan 30, 2019 943 .. math:: E_{\mathrm{c}}[n]=-\int_{0}^{1}d\lambda\iint d^{3}rd^{3}r'\frac{1}{|\mathbf{r}-\mathbf{r}'|}\times\left[\int_{0}^{\infty}\frac{du}{2\pi}\chi_{\lambda}(\mathbf{r},\mathbf{r}',iu)-\chi_{0}(\mathbf{r},\mathbf{r}',iu)\right]\,  Christoph Friedrich committed Jan 11, 2019 944 945 946 947 948  :label: ACFDT where :math:\chi_0 corresponds to the (bare) polarization function [Eq. :eq:eqP] and :math:\chi_\lambda is a partially renormalized response function (i.e., with a scaled electron-electron interaction) fulfilling  Christoph Friedrich committed Jan 30, 2019 949 .. math:: \chi_{\lambda}(\mathbf{r},\mathbf{r}',\omega)=\chi_{\mathrm{0}}(\mathbf{r},\mathbf{r}',\omega)+\iint d^{3}r''d^{3}r'''\chi_{\mathrm{0}}(\mathbf{r},\mathbf{r}'',\omega)\\\hspace{0.8cm}\times\left[\frac{\lambda}{|\mathbf{r}''-\mathbf{r}'''|}+f_{\mathrm{xc,\lambda}}(\mathbf{r}'',\mathbf{r}''',\omega)\right]\chi_{\lambda}(\mathbf{r}''',\mathbf{r}',\omega)\,,  Christoph Friedrich committed Jan 11, 2019 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000 1001 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011  here written in the general form including the (scaled) exchange-correlation kernel :math:f_{\mathrm{xc},\lambda}. For RPA, :math:f_{\mathrm{xc},\lambda}=0. The evaluation of Eq. :eq:ACFDT requires the polarization function to be calculated for all (irreducible) k points on an imaginary frequency mesh. The procedure is thus similar to a one-shot *GW* calculation and needs a similar computation time. Many of the parameters discussed in connection with the *GW* method (e.g., details related to the calculation of the polarization function, :numref:polar) can be used for ACFDT-RPA in the same way. In particular, the same rules using the RESTART option (:numref:RESTART_GW) apply. It is, for example, possible to reuse the file "spex.cor" obtained from a previous *GW* calculation for an ACFDT-RPA calculation. A run with RPAENE comprises the calculation of the total exchange energy using Eq. :eq:hfene. .. _hubbard-u: Interaction parameters (Hubbard *U*) ===================================== In a many-electron system, the motion of all electrons are correlated with each other through the Coulomb interaction. The Coulomb interaction is thus responsible for the fact that we have do deal with a highly complex many-body problem. DFT takes electronic correlation into account, but most of the available functionals are suitable only for weak to moderate correlation. One way to go beyond standard DFT is by many-body perturbation theory, for example, using the *GW* approximation. However, this approach has its limits for very strong electronic correlations as in Mott insulators, for example. As an alternative, one can resort to methods with a special treatment of local correlations. The first one is the so-called LDA+\ *U* scheme, in which the local-density approximation (LDA) of DFT is augmented by an on-site Coulomb repulsion term and an exchange term with the Hubbard *U* and Hund exchange *J* parameters, respectively. A more elaborate computional scheme, which combines many-body model Hamiltonian methods with DFT, is the so-called LDA+DMFT method. In this scheme, the interacting many-body system is mapped onto the subspace of localized states, for example, formed by d orbitals. The subspace contains much less electrons and can thus be treated with high-level quantum many-body approaches such as quantum Monte Carlo, numerical renormalization group, diagrammatic techniques, exact diagonalization, et cetera. The other electrons, not included in the subspace, cannot simply be ignored. For example, their screening processes lead to a renormalization of the effective electron-electron interaction that acts in the subspace, again giving rise to the parameters *U* and *J*. Thus, the Coulomb interaction parameters play a crucial role in the study of the correlation effects in solids. There are two ways to obtain the Hubbard *U* parameter from first principles, constrained local-density approximation (cLDA) and constrained random-phase approximation (cRPA). The latter is implemented in the Spex code and relies on the fact that, since the electrons outside the localized subspace can be assumed delocalized, the random-phase approximation should be appropriate for describing their screening contribution. The cRPA offers an efficient way to calculate the effective Coulomb interaction parameters in solids. Moreover, cRPA allows individual Coulomb matrix elements to be determined, e.g., on-site, off-site, intra-orbital, inter-orbital, and exchange, as well as their frequency dependence. JOB SCREENW ... ---------------- In general, the bare and (fully or partially) screened interaction is calculated in reciprocal space and represented in the mixed product basis (:numref:mpb). The job definition SCREENW tells Spex to project this interaction potential onto the set of local Wannier functions defined in the section WANNIER (:numref:wannier). The resulting interaction matrix elements .. math::  1012  V_{n_1n_2,n_3n_4}(\omega) = \langle n_1 n_2|V(\omega)|n_3 n_4\rangle = \iint w_{n_1}^*(\mathbf{r}) w_{n_3}(\mathbf{r}) V(\mathbf{r},\mathbf{r}';\omega) w_{n_2}^*(\mathbf{r}') w_{n_4}^*(\mathbf{r}') d^3r d^3r'  Christoph Friedrich committed Jan 11, 2019 1013 1014 1015 1016 1017 1018 1019 1020 1021 1022 1023 1024 1025 1026 1027 1028 1029 1030 1031 1032 1033 1034 1035 1036 1037 1038 1039 1040 1041 1042 1043 1044 1045 1046 1047 1048 1049 1050 1051 1052 1053 1054 1055 1056 1057 1058 1059 1060 1061 1062 1063 1064 1065 1066 1067 1068 1069 1070 1071 1072 1073 1074 1075 1076 1077 1078 1079 1080 1081 1082 1083 1084 1085 1086 1087 1088 1089 1090 1091 1092 1093 1094 1095 1096 1097 1098 1099 1100 1101 1102 1103 1104 1105 1106 1107 1108 1109 1110  :label: screenwdef are written to the file "spex.cou" for the bare interaction (:math:V=v) and the screened interaction (:math:V=U if HUBBARD is specified, otherwise :math:V=W, see below). In addition, if :math:\omega=0 is included in the frequency mesh, effective interaction parameters (Hubbard-Hund and/or Kanamori parameterization) are written to the output file for both the bare (*v*) and the screened interaction (*U* or *W*). (Note again that, without the keyword HUBBARD, the effective parameters are obtained from the fully screened interaction and would thus be unsuitable for a Hubbard model Hamiltonian.) The Hubbard-Hund parameters are given for each complete Wannier :math:l shell (e.g., specifying (d) in ORBITALS defines the complete d shell, see :numref:wannier), here written for the partially screened interaction *U*, .. math:: U_l = \frac{1}{(2l+1)^2} \sum_{m=-l}^l \sum_{m'=-l}^l U_{mm',mm'}(0) .. math:: J_l = U_l - \frac{1}{2l(2l+1)} \sum_{m=-l}^l \sum_{m'=-l}^l \left[ U_{mm',mm'}(0) - U_{mm',m'm}(0) \right] The Kanamori parameters are given for groups of t\ :sub:2g (:math:N=3) or e\ :sub:g orbitals (:math:N=2) as .. math:: U_l = \frac{1}{N} \sum_{m=-l}^l U_{mm,mm}(0) \\ .. math:: U'_l = \frac{1}{N(N-1)} \sum_{m=-l}^l \sum_{m'=-l,m'\neq m}^l U_{mm',mm'}(0) \\ .. math:: J_l = \frac{1}{N(N-1)} \sum_{m=-l}^l \sum_{m'=-l,m'\neq m}^l U_{mm',m'm}(0) In the case of a spin-polarized system, the spin-up and spin-down Wannier functions differ slightly, and the respective spin-indices have to be specified in the job definition using one of the labels uu, dd, du, or + (spin summed). There is no default. .. list-table:: Examples :widths: 60 100 * - JOB SCREENW {0:40eV,0.1eV} - Calculate :math:W(\mathbf{r},\mathbf{r}';\omega) and project it onto the Wannier functions; {...} defines the frequency mesh. * - JOB SCREENW {0} - The same only for the static case :math:\omega=0. (Here, {0} is short for {0:0,1}.) * - JOB SCREENW uu{0} - The same for the spin-polarized case: the up-up combination of Wannier functions is selected. HUBBARD (SUSCEP) ----------------- With the keyword HUBBARD (section SUSCEP), the screening that takes place in the correlated subspace is eliminated from the screened interaction, thus yielding the partially screened interaction and, with JOB SCREENW ..., the effective interaction parameters for a Hubbard model Hamiltonian. The correlated subspace is spanned by the Wannier functions. So, the Wannier definition affects the calculation in two ways: (1) It defines the local Wannier basis on which the (bare and) partially screened interaction is projected, and (2) it spans the subspace whose screening is eliminated. Sometimes it makes sense to separate these two effects, for example, if the subspace is spanned by Wannier functions located at several (possibly symmetry-equivalent) atoms, but the effective parameters are only needed for one of the atoms. Such a calculation is performed in two steps with a file name as argument (example: HUBBARD hub). For the first run (i.e, the file does not exist yet), one specifies the full set of Wannier functions. When Spex is run, it just writes information about the screening channels to the specified file ("hub") and stops. Before the second run, one can modify the Wannier definition. In the present example, one would reduce the Wannier set (retaining the orbitals at only one atom). The second Spex run then calculates the partially screened interaction and projects it onto the reduced set of Wannier functions. .. list-table:: Examples :widths: 28 100 * - HUBBARD - Calculate Hubbard *U* interaction matrix including effective parameters. * - HUBBARD hub.dat - Write information about subspace screening to "hub.dat", if the file does not exist, and calculate the Hubbard *U* interaction matrix using data from "hub.dat" otherwise. .. _RESTART_U: RESTART -------- Since SCREENW calculations can be quite costly, the keyword RESTART enables a restart of a calculation that has stopped before finishing. The usage is similar to *GW* calculations (:numref:RESTART_GW). The file "spex.cor" contains the (fully or partially) screened interaction for each (completed) k point, and the file "spex.wcou" holds the last update of the projected interaction matrix. The logical table is similar to :numref:RESTART_GW. .. list-table:: :widths: 30 30 30 :header-rows: 1 * - - spex.cor - spex.wcou * - -- - -- - write * - RESTART - read-write - write * - RESTART 2 - read-write - read-write DISENTGL (WANNIER) ---------------------- (*) There are different ways to eliminate the screening channels in the case of an entangled band structure, i.e., if the correlated subspace is not made up of isolated bands. In this case, it is not possible to unambiguously define the screening that takes place only in the subspace. The default approach of Spex, described in [Phys. Rev. B 83, 121101 (2011)], scales each virtual transition :math:\phi\rightarrow\phi' by the probability that the transition takes place in the subspace. (To be more precise, it is the probability that the electron-hole pair occupying the states :math:\phi and :math:\phi' happens to be in the subspace.) There is another method, called *disentanglement* method [Phys. Rev. 80, 155134 (2009)], in which the reference mean-field system is modified in such a way that each state :math:\phi is either fully contained in the subspace or completely outside (corresponding to probabilities 1 or 0). This method can be used by specifying DISENTGL in the section WANNIER.  1111 .. warning:: Using this keyword actually modifies the mean-field system! (The keyword has not been tested thoroughly yet.)  anoop chandran committed Dec 16, 2018 1112