Hartree Fock
Prerequisites
- You must first complete the "How to use Yambo" modules (Initialization, input file generation, and command line options)
You will need:
- The
SAVEdatabases for bulk hBN - The
yamboexecutable gnuplotfor plotting spectra
Background
The exchange part contribution of the self-energy for a generic state (nk) in reciprocal space reads:
It is important to note that in this way we are adding the HF contribution in a perturbative way to previously calculated DFT energies:
and hence they will differ from a standard self-consistent HF calculation.
Input file generation
Let's start by building up an extended input file for a HF calculation. From yambo -H you should understand that the correct option is yambo -x:
$ cd YAMBO_TUTORIALS/hBN/YAMBO $ yambo (Initialization) $ yambo -x -V all -F hf.in
Looking inside the input file you will find:
EXXRLvcs = 1491 RL # [XX] Exchange RL components %QPkrange # [GW] QP generalized Kpoint/Band indices 1| 14| 1|100| %
The variable EXXRLvcs governs the number of G vector to be summed in the expression of the exchange self-energy reported above. The QPkrange name list it is a generalized index needed to select the first and last kpoints and bands we want to calculate the QP correction. In general, we are interested in the gap of the system or in the band structure across the Fermi Energy. Let's start by editing the hf.in file by selecting the last occupied and first unoccupied bands These are the bands 8 and 9 as reported in the r_setup file. So we change the QPkrange variable in:
%QPkrange # [GW] QP generalized Kpoint/Band indices 1| 14| 8| 9| %
and run calculations to converge the expression above. In the expression of the Exchange Self Energy, we have a summation over bands, an integral over the Brillouin Zone and a Sum over the G vectors. Looking at the occupation factor we realize that occupied states only enters in the expression, so we do not need to worry about the bands as Yambo will include all the occupied bands by default. In order to check convergence of the G vectors in the summation in the Σx expression, we will perform different calculation varying the kinetic energy cutoff governing the number of G vectors entering in the sum. Let's start by setting in hf.in:
EXXRLvcs= 10 Ry
Calculation and analysis of output
Let's now run yambo:
$ yambo -F hf.in -J 3D
The result can be found int the output o-3D.hf file and in the r-3D_HF_and_locXC file. Looking inside the output file we have:
# K-point Band Eo Ehf DFT HF # 1.00000 8.00000 -1.29642 -4.78877 -18.03791 -21.53026 1.000000 9.000000 4.832399 9.755138 -9.592175 -4.669436 2.00000 8.00000 -1.33551 -4.80180 -18.08289 -21.54919 2.00000 9.00000 7.56742 13.44057 -11.55955 -5.68641 ........
Looking at the definition of the HF energy, the third column is the DFT energy (KS eigenvalue), Ehf is the HF energy, and column 4 and 5 report the different contribution to be subtracted Vx (DFT) and added Σx (HF) to the unperturbed DFT eigenvalue. From these data we can calculate the HF gap: you can recognize that the direct gap is found at k-point number 7, being 4.29 eV at DFT level and 11.95 eV in HF approximation. This information can be also easily found in the report file r-D3_HF_and_locXC searching for Direct Gaps:
Before the HF computation:
States summary : Full Metallic Empty
0001-0008 0009-0100
Indirect Gaps [ev]: 3.877976 7.279063
Direct Gaps [ev]: 4.28985 11.35417
... after the HF computation:
[05.01] HF occupations report =============================
States summary : Full Metallic Empty
0001-0008 0009-0100
Indirect Gaps [ev]: 11.31653 16.10704
Direct Gaps [ev]: 11.94662 21.62986
or simply by typing:
$ grep "Direct Gaps" r-D3_HF_and_locXC
and the gap value before and after the HF correction is shown:
Direct Gaps [ev]: 4.28985 11.35417 Direct Gaps [ev]: 11.94662 21.62986
In the report file, the first value indicates the minimum gap while the second one is the maximum energy gap between the same bands.
In the following, in order to converge the direct gap, we will focus the K point where the direct band is found. Inspecting the r_setup file, or any other report file it can be recognized that the direct minimum gap is found at the K-point number 7 (M point): this can be seen by searching the DFT "Direct Gap" and looking at the eigenvalues reported for each k point. Note that the zero energy is fixed at the top of the valence band, K point number 14 (H point). So we can restrict the convergence calculations to the occupied and empty bands at point M, by modifying the input file as:
%QPkrange # [GW] QP generalized Kpoint/Band indices 7| 7| 8| 9| %
Even if for this case the calculations are not at all expensive, reducing the calculations to few bands when performing convergence test allows to save memory and time. Now we can run different calculations changing the value of EXXRLvcs, let's say 10,20,30 and 40 Ry. A possible way is to build four different input files differing for the EXXRLvcs values naming them hf_10Ry.in, hf_20Ry.in, etc..
$ yambo -F hf_10Ry.in -J 3D
next,
$ yambo -F hf_20Ry.in -J 3D .....
Once the calculations are terminated, collect the results found in the output, for instance, putting in a file named hf.dat the following data: Cutoff Energy, Number of G vector in the Sum, HF energy for the occupied state (Ehf for the state 8), HF energy for the unoccupied state (Ehf for the state 9). You can use the grep command to extract these data from the output e.g by typing.
grep 8.000 o-3D.hf_*
and
grep 9.000 o-3D.hf_*
and looking at the column corresponding to Ehf value (column 4). For the number of G vectors corresponding to the cutoff energy cutoff in input search Exchange RL vectors in the report files. Note that it is not possible anymore to extract directly the HF gap from the report file as we did above, as in this case, we did not include all the BZ in the calculation. Now it possible to plot the energy gap in function of the energy cutoff or number of Gvectors:
$gnuplot plot "hf.dat" u 1:($4-$3) w lp t "HF gap vs Energy cutoff" plot "hf.dat" u 2:($4-$3) w lp t "HF gap vs Number of G vector"
or you can plot the gap in function of both quantities using the hf.gnu gnuplot script:
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 Exchange' p "hf.dat" u 1:($4-$3) w lp lw 2 lt -1 t "HF gap", using 1:($4-$3):x2tic(2) lt 7 t
and next
$gnuplot lo "hf.gnu"
and you can see that at 40Ry the energy gap is very close to convergence. Plotting the occupied and unoccupied bands separately you can recognize that for this system the unoccupied band does converge much faster.
The last important convergence test deals with the accuracy of the integral over the Brillouin zone in the expression of the Exchange Self-Energy, which translates into K point convergence. In order to test the K point sampling, you should:
- perform a new non-scf calculation with a bigger k point grid,
- convert wave functions and electronic structure to Yambo databases in a different directory as explained in the DFT and p2y module,
- Initialize the Yambo databases,
- Redo the steps explained in this section.
For now, however, you can skip this convergence test and continue the tutorial at How to obtain the quasiparticle band structure of a bulk material