Part 1: Introduction¶

Content¶

 0 A short introduction to GW
 1 Band gap of Si within the $G_0W_0$ approximation
 2 $GW_0$ band structure using Wannier90 for Si

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:

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

1.1 Task¶

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¶

1.2.1 Input files¶

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!

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?

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

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

1.3.2 Run the calculation¶

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!

Does this calculation include local-field effects?

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!

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

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

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.

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

2.1 Task¶

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.

2.1.1 Input files¶

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!

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

2.2 Task¶

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.

2.2.1 Input files¶

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!

Provided the $TUTORIALS/GW/e01_Si-G0W0/vaspout.h5 file exists, you can add the DFT-band structure to your $GW_0$-band-structure plot.

First execute the following python code!

Next, generate a new band-structure 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¶

1. Why do we use wannierization to compute the band structure for a GW calculation?
2. How many bands should be included for a GW calculation?

Good job! You have finished this tutorial about the GW approximation!¶

Go to Top $\uparrow$