(hf-on-hbn-tutorial-target)= # Hartree-Fock on hBN In this tutorial you will learn how to calculate the quasiparticle corrections in the Hartree-Fock (HF) approximation and how to choose the related input parameters to obtain meaningful converged results. We will use bulk hBN as an example system. ```{admonition} Requirements - `SAVE` folder of bulk hBN: download and extract `hBN.tar.gz` from the [Tutorial files](#tutorial-files-target) page - `yambo` executable - `gnuplot` for plotting the results ``` In general we want to obtain quasiparticle correction to energy levels using many-body perturbation theory (MBPT). The general non-linear quasiparticle equation reads: ```{math} E_{n\boldsymbol{k}}^{\mathrm{QP}} = \varepsilon_{n\boldsymbol{k}} + \langle \psi_{n\boldsymbol{k}} | \Sigma(E_{n\boldsymbol{k}}^{\mathrm{QP}}) - V_{xc} | \psi_{n\boldsymbol{k}} \rangle \, , ``` where {math}`E_{n\boldsymbol{k}}^{\mathrm{QP}}` are the quasiparticle-corrected energies, {math}`\varepsilon_{n\boldsymbol{k}}` and {math}`\psi_{n\boldsymbol{k}}` are the DFT energies and wavefunctions, {math}`V_{xc}` the DFT exchange-correlation potential and {math}`\Sigma` the self-energy. In practice, what Yambo does is to compute the self-energy to obtain the corrections to the one-particle (Kohn-Sham DFT) energies. In the _GW_ approach the self-energy can be separated into two components: a static term called the exchange self-energy ({math}`\Sigma^x`), and a dynamical term (energy dependent) called the correlation self-energy ({math}`\Sigma^c`), that is ```{math} \Sigma_{n\boldsymbol{k}}(\omega) = \Sigma_{n\boldsymbol{k}}^x + \Sigma_{n\boldsymbol{k}}^c(\omega) \, . ``` In this tutorial we only consider the exchange self-energy and compute the corresponding quasiparticle energies, which basically consists in a HF calculation. In reciprocal space the exchange part contribution to the self-energy for a generic state ({math}`n\boldsymbol{k}`) reads ```{math} :label: sigma_x_nk \Sigma_{n\boldsymbol{k}}^x = \langle n\boldsymbol{k} | \Sigma^x | n\boldsymbol{k} \rangle = - \sum_m \int_{\mathrm{BZ}} \frac{d\boldsymbol{q}}{(2\pi)^3} \sum_{\boldsymbol{G}} v(\boldsymbol{q}+\boldsymbol{G}) |\rho_{nm}(\boldsymbol{k},\boldsymbol{q},\boldsymbol{G})|^2 f_{m\boldsymbol{k-q}} \, , ``` where $v$ is the Coulomb potential, $ \rho_{mn}(\boldsymbol{k},\boldsymbol{q},\boldsymbol{G}) =\langle n\boldsymbol{k} | e^{i(\boldsymbol{q}+\boldsymbol{G})\cdot\boldsymbol{r}} | m\boldsymbol{k-q}\rangle $ are the dipole matrix elements and $f$ are the occupation factors. ````{note} In this way we are adding the HF contribution in a perturbative way to previously calculated DFT energies, i.e. ```{math} :label: E_nk_HF E_{n\boldsymbol{k}}^{\mathrm{HF}} = \varepsilon_{n\boldsymbol{k}}^{\mathrm{DFT}} + (\Sigma_{n\boldsymbol{k}}^x - V_{n\boldsymbol{k}}^{xc}) \, , ``` so they will differ from a standard self-consistent HF calculation. ```` Before starting go to the `hBN/YAMBO` directory and initialize the `SAVE` folder: ```{code-block} none $ cd hBN/YAMBO $ ls SAVE $ yambo ``` ## Input file and variables Let's start by building up an input file for a HF calculation. From `yambo -h` we can see that the option to be used is `-hf` (or `-x`): ```{code-block} none $ yambo -x -F hf.in ``` The input file `hf.in` contains the following variables: ```{code-block} none HF_and_locXC # [R] Hartree-Fock EXXRLvcs= 3187 RL # [XX] Exchange RL components VXCRLvcs= 3187 RL # [XC] XCpotential RL components %QPkrange # [GW] QP generalized Kpoint/Band indices 1|14|1|100| % ``` ```{note} :class: dropdown `HF_and_locXC` is a **runlevel** variable (marked with `[R]`). These variables tell Yambo which parts of the code should be executed. Each runlevel enables its own set of input variables. In particular `HF_and_locXC` is needed to compute {math}`\Sigma^x`. ``` - The variable `EXXRLvcs` controls the number of G vectors to be summed in the expression of the exchange self-energy (see Eq.{eq}`sigma_x_nk`). - The variable `VXCRLvcs` sets the number of G vectors used in the evaluation of the density for the $V_{xc}$ matrix element. ```{important} A large number of G vectors is needed in order to have a precise cancellation with the ground state calculation, in particular when GGA potentials are used. A good practice is to leave `VXCRLvcs` at the maximum number (default)[^lessVXCRL]. Also notice that {math}`\Sigma_{n\boldsymbol{k}}^x` and {math}`V_{n\boldsymbol{k}}^{xc}` are to be subtracted, it is important to keep the same value of `EXXRLvcs` and `VXCRLvcs` in order to compute them with the same accuracy. In this tutorial we will ignore these recommendations just to show the convergence behaviour of a HF calculation. ``` [^lessVXCRL]: In case you need to reduce `VXCRLvcs` (e.g. for memory reason), the accuracy can be checked by comparing the `E_xc` values printed in the Yambo report file with those from the DFT calculation (usually printed in the output file). - The `QPkrange` name list it is a generalized index needed to select the first and last kpoints as well as the bands for which we want to calculate the QP correction. In general, we do not want to calculate the corrections for the entire band structure, but instead we are interested in the gap of the system or in the energies around the Fermi level. Let's start by editing the `hf.in` file by selecting the last occupied and first unoccupied bands. Looking into the `r_setup` file we find that these are the bands number 8 and 9, so we set the `QPkrange` variable to ```{code-block} none %QPkrange 1|14|8|9| % ``` ````{tip} :class: dropdown If you want to compute the quasiparticle correction for only one eigenvalue, for example the one with band index 22 and k point number 3, set ```{code-block} none %QPkrange 3|3|22|22| % ``` If you want to compute two specific direct gaps between bands 8 and 9 located at k points 7 and 10, set ```{code-block} none %QPkrange 7|7|8|9| 10|10|8|9| % ``` ```` In the expression of the exchange self-energy (Eq.{eq}`sigma_x_nk`), we have a summation over bands, an integral over the Brillouin Zone and a sum over the G vectors. The presence of the occupation factors $f$ imply that only occupied states enter in the expression, so we do not need to worry about the bands as Yambo will include all the occupied bands by default and thus the sum over $m$ is automatically restricted to valence bands only. We can however study the convergence of the G vectors sum. In the following we will perform different calculations varying the kinetic energy cutoff governing the number of G vectors entering the sum. ## Calculations Let's start by setting ```{code-block} none EXXRLvcs= 10 Ry ``` in `hf.in`. Now run the calculation with ```{code-block} none $ yambo -F hf.in -J 3D ``` ````{dropdown} HF_and_locXC log ```{code-block} none __ __ ____ ___ ___ ____ ___ | | |/ | | | \ / \ | | | o | _ _ | o ) | | ~ | | \_/ | | O | |___, | _ | | | O | | | | | | | | | | |____/|__|__|___|___|_____|\___/ <---> [01] MPI/OPENMP structure, Files & I/O Directories <---> MPI Cores-Threads : 1(CPU)-48(threads) <---> [02] CORE Variables Setup <---> [02.01] Unit cells <---> [02.02] Symmetries <---> [02.03] Reciprocal space <---> [02.04] K-grid lattice <---> Grid dimensions : 6 6 2 <---> [02.05] Energies & Occupations <---> [03] Transferred momenta grid and indexing <---> [04] Local Exchange-Correlation + Non-Local Fock <---> [FFT-HF/Rho] Mesh size: 18 18 44 <---> EXS |########################################| [100%] --(E) --(X) <---> [xc] Functional : Slater exchange(X)+Perdew & Zunger(C) <---> [xc] LIBXC used to calculate xc functional <---> [04.01] Hartree-Fock occupations report <---> [05] Timing Overview <---> [06] Memory Overview <---> [07] Game Over & Game summary ``` ```` Let's inspect the output file `o-3D.hf` and the report file `r-3D_HF_and_locXC`. In the output we have: ``````{tab-set} `````{tab-item} v5 ````{code-block} none $ grep -A7 '# K-point' o-3D.hf ```` ````{dropdown} grep output :open: ```{literalinclude} snippets/o-3D.hf_v5.txt :language: none ``` ```` ````` `````{tab-item} v4.5 (or older) ````{code-block} none $ grep -A7 '# K-point' o-3D.hf ```` ````{dropdown} grep output ```{literalinclude} snippets/o-3D.hf_v4.txt :language: none ``` ```` ````` `````` Looking at the definition of the HF energy (Eq.{eq}`E_nk_HF`), the third column (`Eo`) is the DFT eigenvalue, `Ehf` is the HF energy, while columns 5 and 6 report the two different contributions to the unperturbed DFT eigenvalue: `Vxc` (i.e. {math}`V_{n\boldsymbol{k}}^{xc}` to be subtracted) and `Vnlxc` (i.e. {math}`\Sigma_{n\boldsymbol{k}}^x` to be added). From these data we can calculate the HF gap. In the report file we find that the direct gap is located at k-point number 7 and its value is 4.29 eV at DFT level and 11.95 eV in the HF approximation. ````{tip} :class: dropdown These values are printed in their respective sections of the report file: `Energies & Occupations` for DFT and `Hartree-Fock occupations report` for HF. ```{code-block} none $ grep "Direct Gap" r-3D_HF_and_locXC [X] Direct Gap : 4.289853 [eV] [X] Direct Gap localized at k : 7 [Hartree-Fock] Direct Gap : 11.94549 [eV] [Hartree-Fock] Direct Gap localized at k : 7 ``` ```` ## Convergence of the direct gap In order to converge the direct gap we will focus the k point where it is located. Inspecting either the setup report file or any other report file we can see that the direct minimum gap is found at the k-point number 7. ```{note} In this example of bulk hBN the 7-th k point corresponds to the M point: you can see this by searching the DFT "Direct Gap" and looking at the eigenvalues reported for each k point. Remember that the zero of the energy is fixed at the top of the valence band, which in this case is located at the k point number 14 (H point). ``` Thus we restrict the convergence calculations to the occupied and empty bands at point M by modifying `QPkrange` in the input file as follows: ```{code-block} none %QPkrange 7|7|8|9| % ``` ```{tip} :class: dropdown Even if in this example the calculations are not at all expensive, reducing the calculations to few bands when performing convergence tests allows to save memory and time. ``` We now run a series of calculations changing the value of `EXXRLvcs`, let's say to 10, 20, 30 and 40 Ry. Create four different input files which only differ for the `EXXRLvcs` value and name them `hf_10Ry.in`, `hf_20Ry.in`, etc. ```{code-block} none $ cp hf.in hf_10Ry.in $ cp hf.in hf_20Ry.in ... ``` After you have set `EXXRLvcs` to the corresponding value in each `hf_*.in` file, run the calculations one by one: ```{code-block} none $ yambo -F hf_10Ry.in -J 3D $ yambo -F hf_20Ry.in -J 3D ... ``` Once you have run all the calculations collect the results found in the output files, for example putting in a file named `hf.dat` the following data: cutoff energy (10, ..., 40 Ry), number of G vector in the sum (found in the report), HF energy for the occupied state (`Ehf` for band 8), HF energy for the unoccupied state (`Ehf` for band 9). You can use these two grep commands to extract `Ehf` for band 8 ```{code-block} none $ grep "^ *7 *8" o-3D.hf_* | awk '{print $5}' -2.870097 -3.308719 -3.402065 -3.433922 ``` and band 9 ```{code-block} none $ grep "^ *7 *9" o-3D.hf_* | awk '{print $5}' 9.075388 9.018848 9.012869 9.011314 ``` and the following one to extract the number of G vectors ```{code-block} none $ grep ' Exchange RL vectors' r-3D_HF_and_locXC_0* | awk '{print $6}' 151 371 651 1009 ``` Your `hf.dat` file should look like this: ```{code-block} none # Energy cutoff [Ry] N. G vectors Ehf_8 [eV] Ehf_9 [eV] 10 151 -2.870097 9.075388 20 371 -3.308719 9.018848 30 651 -3.402065 9.012869 40 1009 -3.433922 9.011314 ``` You can plot the results like this ```{code-block} none $ gnuplot gnuplot> plot "hf.dat" u 1:($4-$3) w lp t "HF gap vs Energy cutoff" gnuplot> plot "hf.dat" u 2:($4-$3) w lp t "HF gap vs Number of G vector" ``` ````{dropdown} or use a slightly more fancy gnuplot script ```{code-block} none set key c set grid y set grid x set xtics nomirror set x2tics set xlabel "Cutoff energy [Ry]" set x2label 'Number of RL vectors in {/Symbol S}^x' set ylabel "HF direct gap [eV]" plot "hf.dat" u 1:($4-$3):x2tic(2) w p pt 7 ps 1 t "" rep "hf.dat" u 1:($4-$3) w lp lw 3 pt 7 ps 2 lc 'red' t "HF gap" ``` (put these lines in a `hf.gnu` file and then type `load 'hf.gnu'` in gnuplot). ```` ```{figure} ./images/hBN-HF_gap-conv.png :width: 85% Convergence of the Hartree-Fock direct gap of bulk hBN with respect to the cutoff on the sum over G vectors (see Eq.{eq}`sigma_x_nk`). ``` You can see that at 40 Ry the energy gap is very close to convergence. If you plot the occupied and unoccupied bands separately you can see that for this system the unoccupied band converges much faster than the occupied one. ### k points convergence (_optional_) The last important convergence test deals with the accuracy of the integral over the Brillouin zone in the expression of the Exchange Self-Energy (Eq.{eq}`sigma_x_nk`), which translates into k point convergence. To test the convergence with respecto to the k point sampling, you should: 1. Perform a new non-scf calculation with a bigger k point grid; 2. Convert wave functions and electronic structure to Yambo databases in a different directory; 3. Initialize the new Yambo databases; 4. Redo the steps explained in this tutorial for each k point grid.