Electronic Bandstructure from GW

by Dmitrii Nabok for exciting carbon

Purpose: In this tutorial you will learn how to perform a basic G$_{0}$W$_{0}$ calculation. As an example, the electronic bandstructure of bulk Si is calculated. Notice that self-consistent GW is not yet implemented in exciting.


0. Define relevant environment variables

Read the following paragraphs before starting with the rest of this tutorial!

Before starting, be sure that relevant environment variables are already defined as specified in How to set environment variables for tutorials scripts. Here, you can find a list of the scripts which are relevant for this tutorial, with a short description.

  • PLOT-gwbands.py: Python visualization script for plotting and comparing energy bands.
  • PLOT-gwdos.py: Python visualization script for plotting and comparing density of states.

From now on the symbol $ will indicate the shell prompt.

Requirements: Bash shell. Python numpy, lxml, matplotlib.pyplot, and sys libraries.

Important note: All input parameters are in atomic units!


1. Theoretical background: G0W0 approximation.

The quasiparticle energies in the G$_{0}$W$_{0}$ approximation are given as a solution of the linearized quasi-particle (QP) equation [Hybertsen1986] as

(1)
\begin{align} \epsilon_{n\mathbf{k}}^{QP} = \epsilon_{n\mathbf{k}} + Z_{n\mathbf{k}} \left[ {\rm Re}~\Sigma_{n\mathbf{k}}(\epsilon_{n\mathbf{k}}) - V_{n\mathbf{k}}^{xc} \right]\,, \end{align}

where $\epsilon_{n\mathbf{k}}$ are the Kohn-Sham (KS) eigenvalues, $\Sigma_{n\mathbf{k}}$ and $V_{n\mathbf{k}}^{xc}$ are, respectively, the diagonal matrix elements of the self-energy and the exchange-correlation (xc) potential that is employed in the single-particle KS Hamiltonian. The QP renormalization factor $Z_{n\mathbf{k}}$ accounts for the energy-dependence of the self-energy. The self-energy $\Sigma(\mathbf{r},\mathbf{r}';\omega)$ is given by [Hedin1965]

(2)
\begin{align} \Sigma(\mathbf{r},\mathbf{r}';\omega) = \frac{i}{2\pi} \int G_0(\mathbf{r},\mathbf{r}';\omega+\omega')~W_0(\mathbf{r},\mathbf{r}';\omega')~e^{i\omega'\eta}~d\omega', \end{align}

where $G_0$ is the noninteracting single-particle Green function obtained from the Kohn-Sham states and $W_0$ is the dynamically screened Coulomb potential

(3)
\begin{align} W_0(\mathbf{r},\mathbf{r}';\omega) = \int \epsilon^{-1}(\mathbf{r},\mathbf{r}'';\omega)~v_{\rm C}(\mathbf{r}'',\mathbf{r}') d\mathbf{r}''. \end{align}

Here, $v_{\rm C}(\mathbf{r},\mathbf{r}') = 1/|\mathbf{r}-\mathbf{r}'|$ is the Coulomb potential and $\epsilon(\mathbf{r},\mathbf{r}';\omega)$ is the dielectric function calculated in the random-phase approximation (RPA).


2. Ground-state calculation + GW "single-shot" run

Create a working directory tutorial-GW and move inside it.

$ mkdir tutorial-gw
$ cd tutorial-gw

The typical G$_{0}$W$_{0}$ run consists of two steps:

  1. First, one needs to perform a ground-state calculation to obtain the self-consistent Kohn-Sham KS orbitals and potential.
  2. Second, using this data as an input, one performs a G$_{0}$W$_{0}$ cycle to calculate quasiparticle QP energies.

Below is the input file (input.xml) for bulk Si, which we consider in this tutorial.

<input>
 
   <title>Silicon</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal>
         <basevect>5.13 5.13 0.00</basevect>
         <basevect>5.13 0.00 5.13</basevect>
         <basevect>0.00 5.13 5.13</basevect>
      </crystal>
      <species speciesfile="Si.xml" rmt="2.1">
         <atom coord="0.00 0.00 0.00"></atom>
         <atom coord="0.25 0.25 0.25"></atom>
      </species>
   </structure>
 
   <groundstate
      do="fromscratch"
      rgkmax="7.0"
      ngridk="4 4 4"
      xctype="LDA_PW"
      >
   </groundstate>
 
   <gw
      taskname="g0w0"
      ngridq="2 2 2"
      nempty="22"
      ibgw="1" 
      nbgw="10"
      >
   </gw>
 
</input>

Let us look closer at this input file. Please note that in production calculation one should always check the convergence of results with respect to the groundstate parameters (such as rgkmax, ngridk, etc.). The new element gw initializes G$_{0}$W$_{0}$ calculations. Other specific attributes and elements used in this example are

Parameter Description
gw/taskname Specifies the task name.
gw/ ngridq Specifies the number of both the k and the q-points used specifically in GW calculations. Notice that these values are independent from the one specified by groundstate/ ngridq.
gw/nempty The number of empty states used as an input in GW calculations. Note that this value is independent from the one specified inside the groundstate element.
gw/ibgw, gw/nbgw Give the range of resulting quasiparticle states to be calculated. ibgw (nbgw) is the lowest (highest) band index for which the GW correction is calculated. Notice that the numbering of the band indexes is such that the index 1 is the one for the lowest occupied valence band. Furthermore, nbgw must be chosen to be larger than the number of occupied bands.

For all available parameters and options and their description, one should check the GW Section of Input Reference.

To perform the actual calculation, copy and paste the text from above into a file called input.xml within a directory of your choice. Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

Now, you can start the calculation by invoking the exciting executable (serial version)

$ time excitingser

After the run, one can inspect the INFO.OUT file for checking the convergence behavior of the SCF cycle and for other information related to the ground-state calculation. The main output of G$_{0}$W$_{0}$ run is stored in GW_INFO.OUT and EVALQP.DAT. The former file contains information on the GW program work-flow. In particular, energy band-gap values for KS and QP spectra can be found in the end of this file.

The text file EVALQP.DAT contains all important quantities that are calculated during G$_{0}$W$_{0}$ run. Below, we describe the structure of this file. The line "k-point # 1:" starts the block of data and contains the lattice coordinates of the irreducible k-point and the corresponding Brillouin-zone-integration weight. Then, for each electronic state (within the range determined by the ibgw and nbgw input attributes) one can find the following quantities:

  • the starting KS eigenenergies (E_KS);
  • the quasiparticle eneergies as obtained by considering the exchange part of the self-energy only (E_HF);
  • the G$_{0}$W$_{0}$ quasiparticle energies (E_GW);
  • the diagonal matrix elements of the exchange self-energy (Sx);
  • the diagonal matrix elements of the correlation self-energy (Sc);
  • the diagonal matrix elements of the underlying exchange-correlation potential (Vxc);
  • the energy difference E_HF$-$E_KS (DE_HF);
  • the energy difference E_GW$-$E_KS (DE_GW);
  • the linearization prefactor $Z_{n\mathbf{k}}$ (Znk).

It is important to note that the zero of the E_HF and E_GW quasiparticle energies listed in EVALQP.DAT is the KS Fermi energy (that can be found, e.g., in GW_INFO.OUT). Therefore, in order to obtain absolute QP eigenenergies one has to add to the data listed in EVALQP.DAT the corresponding KS Fermi energy. An example of the content of the first lines of EVALQP.DAT is given in the following.

k-point #     1:    0.000000    0.000000    0.000000    0.125000
 state    E_KS      E_HF       E_GW       Sx         Sc         Vxc         DE_HF        DE_GW       Znk
   1    -0.24974   -0.48679   -0.46467   -0.69023    0.21961   -0.45318   -0.23704   -0.21493    0.78986
   2     0.19048    0.09502   -0.00819   -0.59347    0.09863   -0.49801   -0.09546   -0.19867    0.78709
   3     0.19048    0.09502   -0.00856   -0.59347    0.09817   -0.49801   -0.09546   -0.19904    0.78439
   4     0.19048    0.09502   -0.00861   -0.59347    0.09811   -0.49801   -0.09546   -0.19909    0.78417
   5     0.28272    0.48340    0.11444   -0.23178   -0.15909   -0.43246    0.20068   -0.16828    0.79071
   6     0.28272    0.48340    0.11444   -0.23178   -0.15912   -0.43246    0.20068   -0.16828    0.79120
   7     0.28272    0.48340    0.11493   -0.23178   -0.15857   -0.43246    0.20068   -0.16779    0.79254
   8     0.30752    0.54184    0.14817   -0.32549   -0.18035   -0.55982    0.23433   -0.15934    0.77472
   9     0.47138    0.68975    0.29754   -0.12737   -0.18501   -0.34573    0.21836   -0.17385    0.81888
  10     0.47467    0.69897    0.31363   -0.16335   -0.17396   -0.38765    0.22430   -0.16104    0.79701

3. Quasiparticle band structure

i) Bandstructure plot

Now, we will visualize the results to compare the Kohn-Sham (KS) and quasiparticle (QP) bandstructures. To calculate them, one should introduce the following changes into the input.xml file.

  • In the groundstate element change do = "fromscratch" to do = "skip" to avoid running the ground-state calculations from scratch.
...
   <groundstate
      do="skip"
      ...>
   </groundstate> 
...
  • After the groundstate block, add the new element block properties as indicated below.
...  
   <properties>
      <bandstructure>
         <plot1d>
            <path steps="100">
               <point coord=" 0.750   0.500   0.250" label="W"/>
               <point coord=" 0.500   0.500   0.500" label="L"/>
               <point coord=" 0.000   0.000   0.000" label="GAMMA"/>
               <point coord=" 0.500   0.500   0.000" label="X"/>
               <point coord=" 0.750   0.500   0.250" label="W"/>
               <point coord=" 0.750   0.375   0.375" label="K"/>
            </path>
         </plot1d>
      </bandstructure>
   </properties>
...
  • In the gw element replace taskname = "g0w0" by taskname = "band" to calculate the QP bandstructure on the basis of results obtained from the previous G$_{0}$W$_{0}$ run. In this case, Fourier interpolation is used to evaluate the band energies for the k-point path specified in the properties element.
...
   <gw
      taskname="band"
      ... >
   </gw>
...

Finally, run exciting one more time.

$ time excitingser &

The output files BAND.OUT and BAND-QP.OUT can now be used to plot the electronic bandstructure. This can be easily done using the script PLOT-gwbands.py. To execute it, type

$ PLOT-gwbands.py
The script produces a file PLOT-band.png which you can visualize on the screen with standard tools. The result should look like the following.
PLOT-band.png
ii) Density of states

Now, we will compute and compare the Kohn-Sham (KS) and quasiparticle (QP) density of states. To calculate them, one should introduce the following changes into the input.xml file.

  • Modify the properties block in order to include the dos sub-element as indicated below.
...  
   <properties>
      <!--bandstructure>
         <plot1d>
            <path steps="100">
               <point coord=" 0.750   0.500   0.250" label="W"/>
               <point coord=" 0.500   0.500   0.500" label="L"/>
               <point coord=" 0.000   0.000   0.000" label="GAMMA"/>
               <point coord=" 0.500   0.500   0.000" label="X"/>
               <point coord=" 0.750   0.500   0.250" label="W"/>
               <point coord=" 0.750   0.375   0.375" label="K"/>
            </path>
         </plot1d>
      </bandstructure-->
 
      <dos 
         winddos="-0.5 0.5" 
         nwdos="200"/>
 
   </properties>
...
  • In the gw element replace taskname = "band" by taskname = "dos" to calculate the QP density of states on the basis of results obtained from the previous G$_{0}$W$_{0}$ run.
...
   <gw
      taskname="dos"
      ... >
   </gw>
...

Now, you can run exciting:

$ time excitingser &

The results for the KS and G$_{0}$W$_{0}$ density of states are stored in TDOS.OUT and TDOS-QP.OUT correspondingly. To visualize them one can simply execute the script PLOT-gwdos.py.

$ PLOT-gwdos.py
The script produces a file PLOT-dos.png which you can visualize on the screen with standard tools. The result should look like the following.
PLOT-dos.png
Converging results

The quasiparticle energies obtained from G$_{0}$W$_{0}$ depend on many parameters (on some of them in a crucial way). The most important parameters are

Parameter Description Info
gw/ngridq k-mesh very crucial and should always be converged
gw/nempty number of unoccupied states very crucial and sometimes difficult to converge

In the general case, additional attention should be payed to the proper choice of the system-dependent parameters used in the mixbasis, freqgrid, etc., subelements of gw to obtain physically relevant results (see Input Reference).

Exercise
  • Once you have successfully run the previous GW calculation, you can play around with computational parameters in order to achieve the desired convergence for the value of the electronic band gap.
    1. Try to perform a convergence test with respect to the size of the k and q-meshs that appear in the both the groundstate and the gw elements.
    2. Try to run a similar test for the number of empty states. Be careful when increasing this parameter, since the calculation might become quickly very time-consuming.
  • Repeat the complete GW calculation for a different system, e.g., diamond or boron nitride.

Literature

  • Hybertsen1986: M. S. Hybertsen and S. G. Louie, Phys. Rev. B 34, 5390 (1986).
  • Hedin1965: L. Hedin, Phys. Rev. 139, A796 (1965).
  • Tutorial talk (PDF) at the HoW exciting! 2016 workshop in Berlin .
  • More details on the implementation of the GW formalism within the LAPW method: H. Jiang, R. I. Gomez-Abal, X. Li, C. Meisenbichler, C. Ambrosch-Draxl, and M. Scheffler, "FHI-gap: A GW code Based on Augmented Planewaves", Comp. Phys. Commun. 184, 348 (2013).
  • The converged G$_{0}$W$_{0}$ results for Si and other materials, obtained with exciting, can be found in D. Nabok, A. Gulans, and C. Draxl, "Accurate all-electron G$_{0}$W$_{0}$ quasiparticle energies employing the full-potential augmented plane-wave method", Phys. Rev. B 94, 035118 (2016).
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License