# Part 1: Introduction¶

### 0 A short introduction to GW$\uparrow$¶

The GW approximation is an advanced way to account for electronic correlation in electronic structure calculations. In the GW approximation the electrons are not considered to be independent particles, instead an electron that moves through a material interacts with the other electrons and polarizes its surroundings. The electrons are said to be dressed by this polarization cloud they induce in their immediate environment. These dressed electrons are called quasiparticles.

These quasiparticles (eigenstates and eigenenergies) are obtained by solving the following one-electron equation:

$$(T+V_{ext}+V_h)\psi_{n{\bf k}}({\bf r})+\int d{\bf r}\Sigma({\bf r},{\bf r}', E_{n{\bf k}})\psi_{n{\bf k}}({\bf r}') = E_{n{\bf k}}\psi_{n{\bf k}}({\bf r})$$

Compared to the usual Kohn-Sham (KS) equation of density-functional theory (DFT), one sees that in this so-called quasiparticle equation effects of "exchange and correlation" are introduced through a non-local, orbital- and energy-dependent function, $\Sigma({\bf r},{\bf r}',E_{n{\bf k}})$, called the self-energy.

In GW the self-energy $\Sigma$, is approximated as $\Sigma=GW$, the product of a Green's function $G$ and the screened Coulomb interaction $W$, hence the name. The dielectric screening of the Coulomb interaction is commonly calculated within the random-phase approximation (RPA): in this approximation the electronic interactions between an electron that travels through a medium and its environment are strictly limited to polarization events.

Without going into detail, we note that in principle both $G$ as well as $W$ depend on the quasiparticle eigenstates and eigenenergies ($\psi_{n{\bf k}}$ and $E_{n{\bf k}}$). This means that solving the quasiparticle equation is in fact a self-consistency problem. In praxis, however, one hardly ever iterates the solution of the quasiparticle equation to full self-consistency. Instead most GW calculations fall in either one of the following categories:

• $G_0W_0$: single shot GW

In single shot GW calculations the Green's function and the screened Coulomb interaction are calculated from KS eigenstates and eigenenergies, and the quasiparticle equation is solved a single time (a "single shot"), i.e., without any self-consistency.

• $GW_0$: partially self-consistent GW (in the Green's function only)

In the $GW_0$ approximation the Green's function is iterated only with respect to the eigenenergies, i.e., after the initial single-shot GW calculation the Green's function is recomputed using the quasiparticle eigenenergies but with the initial KS eigenstates. This updated Green's function is then used to set up a new quasiparticle equation, which is then solved for a new set of quasiparticle energies. This procedure is repeated several times.

For more details on the theoretical background of GW, check out the VASP Wiki.

### 1 Band gap of Si within the $G_0W_0$ approximation$\uparrow$¶

By the end of this tutorial, you will be able to:

• run a single-shot GW calculation (G$_0$W$_0$)
• obtain the quasiparticle energies and renormalization factor

Compute the quasiparticle energies, renormalization factor and band gap of Si with a single-shot GW, i.e., G$_0$W$_0$, calculation.

A single-shot GW calculation (G$_0$W$_0$) is performed in three steps:

1. The first step is a DFT-ground-state calculation.
2. In the second step one restarts from the self-consistent solution (WAVECAR) of (1) and computes additional unoccupied KS eigenstates. This is done by a single exact diagonalization of the Hamilton matrix. These additional unoccupied KS states are needed to compute the dielectric screening properties. Furthermore, in this second step we compute the change of the cell periodic part of the KS orbitals with respect to the Bloch vector $\bf{k}$ (this information is stored on the WAVEDER file).
3. The third and final step is the actual $G_0W_0$ calculation. This step needs the WAVECAR and WAVEDER from (2) as a starting point.

#### 1.2 Determine the DFT ground state¶

The input files to run the DFT grounstate calculation are prepared at $TUTORIALS/GW/e01_Si-G0W0. Check them out! POSCAR: Si in the cubic-diamond structure. cd Si 5.430 0.5 0.5 0.0 0.0 0.5 0.5 0.5 0.0 0.5 Si 2 Cartesian 0.00 0.00 0.00 Si 0.25 0.25 0.25 Si  INCAR: A DFT run using the default generalized-gradient-approximation functional of Perdew, Burke, and Ernzerhof (GGA-PBE). Use the tetrahedron method with Blöchl correction by setting ISMEAR=-5 to obtain a smooth density of states (DOS). SYSTEM = cd Si DFT ground state ISMEAR = -5 ! tetrahedron method w Blöchl corrections EDIFF = 1E-8 ! convergence criterion ENCUT = 400 ! energy cutoff  KPOINTS: An equally spaced (6$\times$6$\times$6) Γ-centered grid of$\bf{k}$points. k mesh 0 Gamma centered 6 6 6 0 0 0  KPOINTS_OPT: Specifies$\bf{k}$points along several high-symmetry lines in the first Brillouin zone. After self-consistency is reached, a non-selfconsistent calculation is performed to obtain the band structure (dispersion of the eigenenergies) at these$\bf{k}$points. k points for band structure 10 ! intersections line-mode reciprocal 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.5000000000 0.0000000000 0.5000000000 X 0.5000000000 0.0000000000 0.5000000000 X 0.6250000000 0.2500000000 0.6250000000 U 0.3750000000 0.3750000000 0.7500000000 K 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.0000000000 0.0000000000 0.0000000000 GAMMA 0.5000000000 0.5000000000 0.5000000000 L 0.5000000000 0.5000000000 0.5000000000 L 0.5000000000 0.2500000000 0.7500000000 W 0.5000000000 0.2500000000 0.7500000000 W 0.5000000000 0.0000000000 0.5000000000 X  POTCAR: For GW caluclations, one should use the POTCAR files with _GW appended to their name. These potentials were constructed to reproduce the atomic scattering properties over a larger energy range than the conventional potentials. This is needed when one needs to compute a large number of unoccupied states. You can check the TITEL field in your POTCAR file.  PAW_PBE Si_GW 04May2012 4.00000000000000 parameters from PSCTR are: VRHFIN =Si: s2p2 LEXCH = PE EATOM = 103.0669 eV, 7.5752 Ry TITEL = PAW_PBE Si_GW 04May2012 LULTRA = F use ultrasoft PP ? ... ...  ##### 1.2.2 Run the calculation¶ Open a terminal, navigate to this example's directory and run VASP by entering the following: cd$TUTORIALS/GW/e01_*
vasp_std


This calculation computes the DFT ground state. What is the size of the band gap?
Hint: Execute the following cell to obtain the difference between the highest-occupied-molecular-orbital (HOMO) and lowest-unoccupied-molecular-orbital (LUMO) eigenenergies and to plot the DOS using py4vasp!

In [2]:
import py4vasp
import numpy as np

my_calc = py4vasp.Calculation.from_path("./e01_Si-G0W0")

def gap(calc):
homo = np.amax(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] > 0.5])
lumo = np.amin(calc.band.to_dict()['bands'][calc.band.to_dict()['occupations'] < 0.5])
return lumo - homo

print("PBE band gap for Si:", gap(my_calc))

my_calc.dos.plot()


PBE band gap for Si: 0.7078575386433759


Now, plot the band structure along the $\bf{k}$ points specified in the KPOINTS_OPT file using py4vasp! What is the size of the band gap here?

In [3]:
import py4vasp

my_calc = py4vasp.Calculation.from_path("./e01_Si-G0W0")

my_calc.band.plot("kpoints_opt")


Do the band gaps you deduce from the DOS and from the band structure plot agree? If not, why?

The band gap you have deduced from the difference between the HOMO and LUMO eigenenergies (or the DOS) are based upon the eigenenergies at the $\bf{k}$ points of the regular mesh specified in the KPOINTS file, whereas in the band-structure plot one obtains the dispersion along the lines specified in the KPOINTS_OPT file. The latter includes more $\bf{k}$ points along the path between $\Gamma - X$. In Si, the band gap is indirect: the HOMO is found at $\Gamma$, and the LUMO is located slight away from the $X$ point along $\Gamma - X$. This point is not included in the regular grid of $\bf{k}$ points specified in the KPOINTS file!

#### 1.3 Compute additional unoccupied Kohn-Sham orbitals¶

In step two of this exercise, we compute the unoccupied KS orbitals and the frequency-dependent dielectric function in the independent particle picture. It is necessary to compute roughly 3 to 4 times as many KS orbitals as VASP would compute by default. In the previous DFT calculation, check the value of NBANDS written to the OUTCAR file! How many of these bands are occupied? Also, compare this to the value of NBANDS in the unoccupied-states/INCAR below!

##### 1.3.1 Input files¶

unoccupied-states/INCAR:
Increase NBANDS with respect to the default and perform a single exact diagonalization of the Hamiltonian (ALGO=Exact and NELM=1). The LOPTICS tag is set so that the derivative of the Bloch states with respect to the Bloch wavevector $\bf{k}$ is calculated and written to the WAVEDER file. The WAVECAR and WAVEDER files of this calculation will be the starting point for the actual GW calculation.

SYSTEM  = Si DFT unoccupied states

ISMEAR  = 0      ! Gaussian smearing
SIGMA   = 0.05   ! smearing in eV, small sigma is required to avoid partial occupancies
ENCUT   = 400    ! energy cutoff

NBANDS  = 64     ! maybe even more bands

ALGO    = Exact  ! use exact diagonalization of the Hamiltonian
NELM    = 1      ! since we are already converged stop after one step

LOPTICS = T      ! write WAVEDER


cd $TUTORIALS/GW/e01_*/unoccupied-states cp ../WAVECAR . mpirun -np 2 vasp_std  As a byproduct, you are now able to plot the electronic dielectric function using py4vasp! In [5]: import py4vasp optics_calc = py4vasp.Calculation.from_path("./e01_Si-G0W0/unoccupied-states") optics_calc.dielectric_function.plot()  Does this calculation include local-field effects? Click to see the answer! No, this is the dielectric function in the independent-particle approximation. See the documentation of the LOPTICS tag. #### 1.4 The G$_0$W$_0$calculation¶ In the third and final step, we perform the actual G$_0$W$_0$calculation! ##### 1.4.1 Input files¶ single-shot/INCAR: To perform a single-shot GW calculation, specify ALGO=EVGW0 ("eigenvalue GW") and request a single electronic update (NELMGW=1).Do not forget to set NBANDS to the same value that was used to generate the WAVECAR and WAVEDER files in the previous step. SYSTEM = Si G0W0 ISMEAR = 0 ! Gaussian smearing SIGMA = 0.05 ! smearing in eV, small sigma is required to avoid partial occupancies ENCUT = 400 ! energy cutoff NBANDS = 64 ! maybe even more bands ! G0W0 ALGO = EVGW0 ! use "GW0" for VASP.5.X NELMGW = 1 ! use "NELM" prior to VASP.6.3 NOMEGA = 50 ! default  ##### 1.4.2 Run the calculation¶ cd$TUTORIALS/GW/e01_*/single-shot
cp ../unoccupied-states/WAVECAR ../unoccupied-states/WAVEDER .
mpirun -np 2 vasp_std


You can find the quasiparticle energies QP-energies and the renormalization factor Z for each ${\bf k}$-point and band in the OUTCAR file. Search for the corresponding lines!


...

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

k-point   1 :       0.0000    0.0000    0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.4909      -6.8379     -10.9930     -10.4578     -17.5178       0.6483       2.0000       1.1943
2       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
3       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
4       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
5       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
6       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
7       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
8       8.8415       9.2390     -10.5240     -11.0570      -6.0614       0.7458       0.0000      -0.0773

k-point   2 :       0.1667   -0.0000   -0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.1297      -6.4866     -11.0371     -10.4913     -17.3967       0.6539       2.0000       1.1417
2       3.0917       2.7016     -11.4416     -10.9098     -13.1355       0.7334       2.0000       0.0940
3       5.0228       4.6678     -11.7041     -11.2334     -12.6637       0.7542       2.0000       0.0629
4       5.0228       4.6678     -11.7041     -11.2334     -12.6637       0.7542       2.0000       0.0629
5       7.8246       8.0877      -9.8624     -10.2108      -5.8624       0.7552       0.0000      -0.0595
6       8.6837       8.9664      -9.8785     -10.2580      -5.6701       0.7450       0.0000      -0.0721
7       8.6837       8.9664      -9.8785     -10.2580      -5.6701       0.7450       0.0000      -0.0721
8      10.9292      11.2599     -10.4796     -10.9314      -5.5612       0.7320       0.0000      -0.1380
...
...



The first column is the band index of the KS orbital. The third column lists the quasiparticle energies $E_{{n{{\bf {k}}}}}$.

Column two, four, five and seven are the KS energies $E_{{n{{\bf {k}}}}}^{{(0)}}$, the diagonal matrix elements of the self-energy $\langle \psi _{{n{{\bf {k}}}}}^{{(0)}}|\Sigma (\omega =E_{{n{{\bf {k}}}}}^{{(0)}})|\psi _{{n{{\bf {k}}}}}^{{(0)}}\rangle$ , the exchange-correlation potential energy, and the renormalization factor $Z_{{n{{\bf {k}}}}}$, respectively.

How large is the quasiparticle band gap in the G$_0$W$_0$ approximation?


...

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

k-point   1 :       0.0000    0.0000    0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.4909      -6.8379     -10.9930     -10.4578     -17.5178       0.6483       2.0000       1.1943
2       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
3       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
4       5.4736       5.1235     -11.8740     -11.4114     -12.7301       0.7568       2.0000       0.0584
5       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
6       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
7       8.0356       8.3142      -9.7373     -10.1085      -5.7301       0.7504       0.0000      -0.0603
8       8.8415       9.2390     -10.5240     -11.0570      -6.0614       0.7458       0.0000      -0.0773

...
...
k-point  13 :       0.5000    0.5000    0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -2.3453      -2.7160     -11.4052     -10.8687     -16.0347       0.6909       2.0000       0.5746
2      -2.3453      -2.7160     -11.4052     -10.8687     -16.0347       0.6909       2.0000       0.5746
3       2.6195       2.1679     -11.2874     -10.6657     -13.2787       0.7264       2.0000       0.0980
4       2.6195       2.1679     -11.2874     -10.6657     -13.2787       0.7264       2.0000       0.0980
5       6.1815       6.2826      -8.8671      -8.9982      -5.2301       0.7712       0.0000      -0.0495
6       6.1815       6.2826      -8.8671      -8.9982      -5.2301       0.7712       0.0000      -0.0495
7      15.5850      15.8769     -10.3675     -10.7774      -3.6978       0.7122       0.0000      -0.5730
8      15.5850      15.8769     -10.3675     -10.7774      -3.6978       0.7122       0.0000      -0.5730
...
...



The HOMO is found at the $\Gamma$ point (k-point 1), and the LUMO at $X$ (k-point 13): $E_{\rm LUMO}-E_{\rm HOMO} = 6.2826 - 5.1235 = 1.1591$ eV.

Compare the dielectric function obtained within the independent-particle approximation (IPA) to the dielectric function in the random-phase approximation (RPA). The latter accounts for the fact that an electron that travels through a medium interacts (polarizes) its surroundings. Per default VASP computes the screened Coulomb interaction in GW calculations within the RPA.

To plot the aforementioned dielectric functions, first execute the get_dielectric_function.sh bash script and then the python code below.

bash get_dielectric_function.sh

In [8]:
import numpy as np
from py4vasp import plot

filename = "./e01_Si-G0W0/single-shot/dielectric_function_IPA.dat"
# uncomment the following line to use the RPA dielectric function
#filename = "./e01_Si-G0W0/single-shot/dielectric_function_RPA.dat"

# read the dielectric function from file
w, eps_imag, eps_real = np.loadtxt(filename, unpack=True)

# plot the dielectric function as a function of frequency
# DFPT
plot(
(w, eps_real, "Dielectric function (real)"),
(w, eps_imag, "Dielectric function (imaginary)"),
xlabel = "Frequency (eV)",
)


How does the dielectric function change from the IPA to the RPA? Look for instance to the change in the real part of the dielectric function at $\omega=0$ eV. This can also be done conveniently with py4vasp.

In [9]:
from py4vasp import Calculation

calc = Calculation.from_path("./e01_Si-G0W0/single-shot")
calc.dielectric_function.plot("Re(IPA, RPA)")


#### 1.5 Questions¶

1. What quantity is approximated by $GW$?
2. Where does VASP write the quasiparticle energies and renormalization factors?
3. At what level of approximation does VASP compute the screened Coulomb interaction in GW per default?

### 2 The $GW_0$ band structure of Si$\uparrow$¶

By the end of this tutorial, you will be able to:

• perform a partially self-consistent $GW_0$ calculation
• aproximate the band structure using Wannierization

Compute the band structure of cubic-diamond Si within the partially self-consistent $GW_0$ approximation.

In this example, we perform a partially self-consistent $GW_0$ calculation (partially self-consistent in the Green's function). In the $GW_0$ approximation, the Green's function is updated only with respect to the eigenenergies, i.e., after the initial single-shot $GW$ calculation the Green's function is recomputed using the quasiparticle eigenenergies, but with the initial KS orbitals. This updated Green's function is then used to set up a new quasiparticle equation, which is solved for a new set of quasiparticle energies. This procedure is repeated several times, either until the quasiparticles energies are converged or a specified number of iterations has been reached.

The input files can be found in $TUTORIALS/GW/e02_Si-GW0-band. Check them out! As you will notice the KPOINTS, POSCAR, and POTCAR files are identical to the ones in the previous example. The INCAR file for the$GW_0$calculation is slightly different from the$G_0W_0$one: INCAR: To perform a partially self-consistent GW$_0$calculation, specify ALGO=EVGW0 ("eigenvalue GW") and request a number of electronic updates (NELMGW>1). Do not forget to set NBANDS to the same value that was used to generate the WAVECAR and WAVEDER files from which the GW0 calculation will restart. In addition we have increased the number of quasiparticle states that will be computed (NBANDSGW=16). This was done in order to facilitate the construction of a bandstructucture plot in the next task. SYSTEM = cd Si GW0 ISMEAR = 0 SIGMA = 0.05 ENCUT = 400 NBANDS = 64 NBANDSGW = 16 # GW0 ALGO = EVGW0 NELMGW = 3  ##### 2.1.2 Run the calculation¶ In addition to the input files mentioned above, you will need the WAVECAR and WAVEDER files from task 1.3! In case you do not have these readily available you need to (re)run task 1.2 and 1.3. Open a terminal, navigate to this example's directory and run VASP by entering the following: cd$TUTORIALS/GW/e02_*
cp ../e01_*/unoccupied-states/WAVECAR ../e01_*/unoccupied-states/WAVEDER .
mpirun -np 4 vasp_std


You can find the quasiparticle energies QP-energies and the renormalization factor Z for each $\bf{k}$-point and band, written into the OUTCAR file for each iteration of the eigenenergies in $G$. Search for the corresponding lines!


...

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

k-point   1 :       0.0000    0.0000    0.0000
band No.  KS-energies  QP-energies   sigma(KS)   V_xc(KS)     V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.4909      -6.7562     -10.8777     -10.4578     -17.5178       0.6317       2.0000       1.0614
2       5.4736       5.1331     -11.8594     -11.4114     -12.7301       0.7600       2.0000       0.0086
3       5.4736       5.1331     -11.8594     -11.4114     -12.7301       0.7600       2.0000       0.0086
4       5.4736       5.1331     -11.8594     -11.4114     -12.7301       0.7600       2.0000       0.0086
5       8.0356       8.3077      -9.7483     -10.1085      -5.7301       0.7555       0.0000      -0.0428
6       8.0356       8.3077      -9.7483     -10.1085      -5.7301       0.7555       0.0000      -0.0428
7       8.0356       8.3077      -9.7483     -10.1085      -5.7301       0.7555       0.0000      -0.0428
8       8.8415       9.2299     -10.5366     -11.0570      -6.0614       0.7464       0.0000      -0.0517

...
...

k-point   1 :       0.0000    0.0000    0.0000
band No. old QP-enery  QP-energies   sigma(KS)   T+V_ion+V_H  V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.7562      -6.8635     -10.8915       3.9669     -17.5178       0.6370       2.0000       1.0485
2       5.1331       5.0860     -11.8138      16.8851     -12.7301       0.7615       2.0000       0.0095
3       5.1331       5.0860     -11.8138      16.8851     -12.7301       0.7615       2.0000       0.0095
4       5.1331       5.0860     -11.8138      16.8851     -12.7301       0.7615       2.0000       0.0095
5       8.3077       8.3382      -9.7960      18.1441      -5.7301       0.7546       0.0000      -0.0442
6       8.3077       8.3382      -9.7960      18.1441      -5.7301       0.7546       0.0000      -0.0442
7       8.3077       8.3382      -9.7960      18.1441      -5.7301       0.7546       0.0000      -0.0442
8       9.2299       9.2733     -10.6104      19.8985      -6.0614       0.7460       0.0000      -0.0530

...
...

k-point   1 :       0.0000    0.0000    0.0000
band No. old QP-enery  QP-energies   sigma(KS)   T+V_ion+V_H  V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.8635      -6.8952     -10.8803       3.9669     -17.5178       0.6354       2.0000       1.0531
2       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
3       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
4       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
5       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
6       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
7       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
8       9.2733       9.2787     -10.6179      19.8985      -6.0614       0.7454       0.0000      -0.0535
...
...



The first column is the band index of the KS orbital. The third column lists the quasiparticle energies $E_{{n{{\bf {k}}}}}$.

Note that in the first iteration, column two lists the KS eigenergies, whereas in the subsequent iterations this column shows the quasiparticle energies of the previous iteration.

How large is the quasiparticle band gap in the $GW_0$ approximation after three iterations?


...

QP shifts <psi_nk| G(iteration)W_0 |psi_nk>: iteration 1
for sc-GW calculations column KS-energies equals QP-energies in previous step
and V_xc(KS)=  KS-energies - (<T + V_ion + V_H > + <T+V_H+V_ion>^1  + <V_x>^1)

k-point   1 :       0.0000    0.0000    0.0000
band No. old QP-enery  QP-energies   sigma(KS)   T+V_ion+V_H  V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -6.8635      -6.8952     -10.8803       3.9669     -17.5178       0.6354       2.0000       1.0531
2       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
3       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
4       5.0860       5.0779     -11.8096      16.8851     -12.7301       0.7623       2.0000       0.0098
5       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
6       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
7       8.3382       8.3417      -9.8012      18.1441      -5.7301       0.7545       0.0000      -0.0443
8       9.2733       9.2787     -10.6179      19.8985      -6.0614       0.7454       0.0000      -0.0535

...
...
k-point  13 :       0.5000    0.5000    0.0000
band No. old QP-enery  QP-energies   sigma(KS)   T+V_ion+V_H  V^pw_x(r,r')   Z            occupation Imag(sigma)

1      -2.7418      -2.7648     -11.2986       8.5234     -16.0347       0.6893       2.0000       0.4684
2      -2.7418      -2.7648     -11.2986       8.5234     -16.0347       0.6893       2.0000       0.4684
3       2.1118       2.0970     -11.1937      13.2852     -13.2787       0.7278       2.0000       0.0344
4       2.1118       2.0970     -11.1937      13.2852     -13.2787       0.7278       2.0000       0.0344
5       6.2991       6.2999      -8.8795      15.1796      -5.2301       0.7775       0.0000      -0.0355
6       6.2991       6.2999      -8.8795      15.1796      -5.2301       0.7775       0.0000      -0.0355
7      15.8806      15.8902     -10.4683      26.3624      -3.6978       0.7094       0.0000      -0.5165
8      15.8806      15.8902     -10.4683      26.3624      -3.6978       0.7094       0.0000      -0.5165
...
...



The HOMO is found at the $\Gamma$-point ($\bf{k}$-point 1), and the LUMO at $X$ ($\bf{k}$-point 13). The bandgap:

• After the first iteration ($G_0W_0$): $E_{\rm LUMO}-E_{\rm HOMO} = 6.2826 - 5.1235 = 1.1591$ eV.
• After three iterations ($GW_0$): $E_{\rm LUMO}-E_{\rm HOMO} = 6.2999 - 5.0779 = 1.2220$ eV.

Approximate the $GW_0$ of cubic-diamond Si using Wannier90.

In a band-structure plot we vizualize dependence of the one-electron eigenenergies $\epsilon_{n{\bf k}}$ on the Bloch vector $\bf{k}$ (aka, the dispersion relation of the material). This is commonly done along particular lines of high symmetry in the first Brillouin zone of the reciprocal-space lattice of the material. Unfortunately, within the $GW$ approximation, VASP can only calculate quasiparticle energies at $\bf{k}$ points that are part of a regular mesh, and not for an arbitrary set of $\bf{k}$ points along a line in reciprocal space. The GW-band structure may, however, be approximated using Wannier90.

It is important to realize this really only is an approximate band structure. Essentially we obtain the quasiparticle energies $E_{n{\bf k}}$ at arbitrary $\bf{k}$ from a Fourier-interpolation between the $GW_0$ quasiparticle energies we have computed on the $\Gamma$-centered Monkhorst-Pack mesh specified in the KPOINTS file. The quasiparticle energies have been stored in the WAVECAR file of the previous task.

The input files can be found in $TUTORIALS/GW/e02_Si-GW0-band/band-structure. Check them out! As you will notice, the KPOINTS, POSCAR, and POTCAR files are identical to the ones in the previous example. Only the INCAR file is different: INCAR: The combination of ALGO=NONE and NELM=1 tells VASP to read its input files and immediately go to the post-processing phase, skipping any actual work (hence "NONE"). If LWANNIER90=T, then VASP will write input files for Wannier90 in this post-processing phase. NUM_WANN will inform Wannier90 how many Wannier functions are to be constructed and the contents of the WANNIER90_WIN tag will be copied verbatim into the wannier90.win input file of Wannier90. We will restart this calculation from the WAVECAR file that was produced by the$GW_0$calculation. Do not forget to set NBANDS to the same value that was used in that calculation. SYSTEM = cd Si GW0 ENCUT = 400 NBANDS = 64 NELM = 1 ALGO = NONE LWAVE = F LCHARG = F # write input files for Wannier90 LWANNIER90 = T NUM_WANN = 8 WANNIER90_WIN = " exclude_bands 17-64 Begin Projections Si:sp3 End Projections # Disentanglement dis_win_min = -7 dis_win_max = 16 dis_num_iter = 100 guiding_centres = true "  ##### 2.2.2 Run the calculation¶ cd$TUTORIALS/GW/e02_*/band-structure
cp ../WAVECAR .
vasp_std


This produces a Wannier90-related input files (wannier90.*). Next, you can call Wannier90 to contruct maximally localized Wannier functions:

wannier90.x wannier90.win


Check whether Wannier90 succeeded in creating the aforementioned Wannier functions by opening the wannier90.wout file. Near the end of of the file, you should find something very similar to:

 Final State
WF centre and spread    1  (  0.208905,  0.208905,  0.208905 )     2.96906743
WF centre and spread    2  (  0.208905, -0.208905, -0.208905 )     2.96906716
WF centre and spread    3  ( -0.208905,  0.208905, -0.208905 )     2.96906733
WF centre and spread    4  ( -0.208905, -0.208905,  0.208905 )     2.96906752
WF centre and spread    5  (  1.813343,  1.813343,  1.813343 )     2.52490771
WF centre and spread    6  (  1.813343,  0.901657,  0.901656 )     2.52490775
WF centre and spread    7  (  0.901657,  1.813343,  0.901656 )     2.52490773
WF centre and spread    8  (  0.901657,  0.901657,  1.813343 )     2.52490772
Sum of centres and spreads (  5.430000,  5.430000,  5.430000 )    21.97590034

Spreads (Ang^2)       Omega I      =    16.261745330
================       Omega D      =     0.187447192
Omega OD     =     5.526707818
Final Spread (Ang^2)       Omega Total  =    21.975900340


Finally, you can generate the $GW_0$-band-structure plot. Execute the following commands in your terminal!

echo """
# Band-structure plot
restart         =  plot
bands_plot      =  true
begin kpoint_path
G 0.000 0.000 0.000 X 0.500 0.000 0.500
X 0.500 0.000 0.500 U 0.625 0.250 0.625
K 0.375 0.375 0.750 G 0.000 0.000 0.000
G 0.000 0.000 0.000 L 0.500 0.500 0.500
L 0.500 0.500 0.500 W 0.500 0.250 0.750
W 0.500 0.250 0.750 X 0.500 0.000 0.500
end kpoint_path
bands_num_points 40
bands_plot_format gnuplot
""" >> wannier90.win

wannier90.x wannier90.win

echo """
set term png
set out 'band.png'
set ytics 5
set mytics 5
set grid ytics mytics
""" > band.gp
cat wannier90_band.gnu >> band.gp

gnuplot band.gp


and open the band.png file from the file browser!