Part 1: Introduction¶
0 A short introduction to GW ¶
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 oneelectron 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 KohnSham (KS) equation of densityfunctional theory (DFT), one sees that in this socalled quasiparticle equation effects of "exchange and correlation" are introduced through a nonlocal, orbital and energydependent function, $\Sigma({\bf r},{\bf r}',E_{n{\bf k}})$, called the selfenergy.
In GW the selfenergy $\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 randomphase 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 selfconsistency problem. In praxis, however, one hardly ever iterates the solution of the quasiparticle equation to full selfconsistency. Instead most GW calculations fall in either one of the following categories:

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 selfconsistency.
$GW_0$: partially selfconsistent 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 singleshot 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.
By the end of this tutorial, you will be able to:
 run a singleshot GW calculation (G$_0$W$_0$)
 obtain the quasiparticle energies and renormalization factor
1.1 Task¶
Compute the quasiparticle energies, renormalization factor and band gap of Si with a singleshot GW, i.e., G$_0$W$_0$, calculation.
A singleshot GW calculation (G$_0$W$_0$) is performed in three steps:
 The first step is a DFTgroundstate calculation.
 In the second step one restarts from the selfconsistent 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).
 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¶
1.2.1 Input files¶
The input files to run the DFT grounstate calculation are prepared at $TUTORIALS/GW/e01_SiG0W0
. Check them out!
POSCAR:
Si in the cubicdiamond 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 generalizedgradientapproximation functional of Perdew, Burke, and Ernzerhof (GGAPBE). 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 = 1E8 ! 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 highsymmetry lines in the first Brillouin zone. After selfconsistency is reached, a nonselfconsistent calculation is performed to obtain the band structure (dispersion of the eigenenergies) at these $\bf{k}$ points.
k points for band structure
10 ! intersections
linemode
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 highestoccupiedmolecularorbital (HOMO) and lowestunoccupiedmolecularorbital (LUMO) eigenenergies and to plot the DOS using py4vasp!
import py4vasp
import numpy as np
my_calc = py4vasp.Calculation.from_path("./e01_SiG0W0")
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()
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?
import py4vasp
my_calc = py4vasp.Calculation.from_path("./e01_SiG0W0")
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?
Click to see the answer!
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 bandstructure 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 KohnSham orbitals¶
In step two of this exercise, we compute the unoccupied KS orbitals and the frequencydependent 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 unoccupiedstates/
INCAR below!
1.3.1 Input files¶
unoccupiedstates/
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
1.3.2 Run the calculation¶
cd $TUTORIALS/GW/e01_*/unoccupiedstates
cp ../WAVECAR .
mpirun np 2 vasp_std
As a byproduct, you are now able to plot the electronic dielectric function using py4vasp!
import py4vasp
optics_calc = py4vasp.Calculation.from_path("./e01_SiG0W0/unoccupiedstates")
optics_calc.dielectric_function.plot()
Does this calculation include localfield effects?
Click to see the answer!
No, this is the dielectric function in the independentparticle 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¶
singleshot/
INCAR:
To perform a singleshot 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_*/singleshot
cp ../unoccupiedstates/WAVECAR ../unoccupiedstates/WAVEDER .
mpirun np 2 vasp_std
You can find the quasiparticle energies QPenergies
and the renormalization factor Z
for each ${\bf k}$point and band in the OUTCAR file. Search for the corresponding lines!
Click to see the answer!
...
QP shifts <psi_nk G(iteration)W_0 psi_nk>: iteration 1
for scGW calculations column KSenergies equals QPenergies in previous step
and V_xc(KS)= KSenergies  (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)
kpoint 1 : 0.0000 0.0000 0.0000
band No. KSenergies QPenergies 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
kpoint 2 : 0.1667 0.0000 0.0000
band No. KSenergies QPenergies 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 selfenergy $\langle \psi _{{n{{\bf {k}}}}}^{{(0)}}\Sigma (\omega =E_{{n{{\bf {k}}}}}^{{(0)}})\psi _{{n{{\bf {k}}}}}^{{(0)}}\rangle$ , the exchangecorrelation 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?
Click to see the answer!
...
QP shifts <psi_nk G(iteration)W_0 psi_nk>: iteration 1
for scGW calculations column KSenergies equals QPenergies in previous step
and V_xc(KS)= KSenergies  (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)
kpoint 1 : 0.0000 0.0000 0.0000
band No. KSenergies QPenergies 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
...
...
kpoint 13 : 0.5000 0.5000 0.0000
band No. KSenergies QPenergies 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 (kpoint 1
), and the LUMO at $X$ (kpoint 13
): $E_{\rm LUMO}E_{\rm HOMO} = 6.2826  5.1235 = 1.1591$ eV.
Compare the dielectric function obtained within the independentparticle approximation (IPA) to the dielectric function in the randomphase 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
import numpy as np
from py4vasp import plot
filename = "./e01_SiG0W0/singleshot/dielectric_function_IPA.dat"
# uncomment the following line to use the RPA dielectric function
#filename = "./e01_SiG0W0/singleshot/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.
from py4vasp import Calculation
calc = Calculation.from_path("./e01_SiG0W0/singleshot")
calc.dielectric_function.plot("Re(IPA, RPA)")
1.5 Questions¶
 What quantity is approximated by $GW$?
 Where does VASP write the quasiparticle energies and renormalization factors?
 At what level of approximation does VASP compute the screened Coulomb interaction in GW per default?
By the end of this tutorial, you will be able to:
 perform a partially selfconsistent $GW_0$ calculation
 aproximate the band structure using Wannierization
2.1 Task¶
Compute the band structure of cubicdiamond Si within the partially selfconsistent $GW_0$ approximation.
In this example, we perform a partially selfconsistent $GW_0$ calculation (partially selfconsistent 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 singleshot $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.
2.1.1 Input files¶
The input files can be found in $TUTORIALS/GW/e02_SiGW0band
. 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 selfconsistent 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_*/unoccupiedstates/WAVECAR ../e01_*/unoccupiedstates/WAVEDER .
mpirun np 4 vasp_std
You can find the quasiparticle energies QPenergies
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!
Click to see the answer!
...
QP shifts <psi_nk G(iteration)W_0 psi_nk>: iteration 1
for scGW calculations column KSenergies equals QPenergies in previous step
and V_xc(KS)= KSenergies  (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)
kpoint 1 : 0.0000 0.0000 0.0000
band No. KSenergies QPenergies 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
...
...
kpoint 1 : 0.0000 0.0000 0.0000
band No. old QPenery QPenergies 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
...
...
kpoint 1 : 0.0000 0.0000 0.0000
band No. old QPenery QPenergies 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?
Click to see the answer!
...
QP shifts <psi_nk G(iteration)W_0 psi_nk>: iteration 1
for scGW calculations column KSenergies equals QPenergies in previous step
and V_xc(KS)= KSenergies  (<T + V_ion + V_H > + <T+V_H+V_ion>^1 + <V_x>^1)
kpoint 1 : 0.0000 0.0000 0.0000
band No. old QPenery QPenergies 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
...
...
kpoint 13 : 0.5000 0.5000 0.0000
band No. old QPenery QPenergies 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.
2.2 Task¶
Approximate the $GW_0$ of cubicdiamond Si using Wannier90.
In a bandstructure plot we vizualize dependence of the oneelectron 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 reciprocalspace 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 GWband 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 Fourierinterpolation between the $GW_0$ quasiparticle energies we have computed on the $\Gamma$centered MonkhorstPack mesh specified in the KPOINTS file. The quasiparticle energies have been stored in the WAVECAR file of the previous task.
2.2.1 Input files¶
The input files can be found in $TUTORIALS/GW/e02_SiGW0band/bandstructure
. 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 postprocessing phase, skipping any actual work (hence "NONE").
If LWANNIER90=T, then VASP will write input files for Wannier90 in this postprocessing 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 1764
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_*/bandstructure
cp ../WAVECAR .
vasp_std
This produces a Wannier90related 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$bandstructure plot. Execute the following commands in your terminal!
echo """
# Bandstructure 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!
Provided the $TUTORIALS/GW/e01_SiG0W0/vaspout.h5
file exists, you can add the DFTband structure to your $GW_0$bandstructure plot.
First execute the following python code!
import py4vasp
e02_calc = py4vasp.Calculation.from_path("./e01_SiG0W0")
d = e02_calc.band.to_dict(selection="kpoints_opt")
# rescaling, the last value of the first column in wannier90_band.dat
x = 0.51924171E+01
with open("./e02_SiGW0band/bandstructure/e01_kpoints_opt_band.dat", "w") as f:
for i in range(8):
for k in range(60):
f.write(str(d['kpoint_distances'][k]/d['kpoint_distances'][1]*x) + " " + str(d['bands'][k, i]+d['fermi_energy']) + "\n")
f.write("\n")
print("DFT Efermi: ", d['fermi_energy'])
e02_calc.band.plot(selection="kpoints_opt")
Next, generate a new bandstructure plot!
echo """
set out 'band_DFT_vs_GW.png'
plot 'wannier90_band.dat', 'e01_kpoints_opt_band.dat' w l, 5.628
""" >> band.gp
gnuplot band.gp
Open the band_DFT_vs_GW.png file to see the result!
2.3 Questions¶
 Why do we use wannierization to compute the band structure for a GW calculation?
 How many bands should be included for a GW calculation?