Electronic Bandstructure from GW

for lithium version of Exciting

Purpose: In this tutorial you will learn how to perform a basic GW calculation. As an example, the electronic bandstructure of bulk Si is calculated.

0. Define relevant shell variables and download scripts

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

Before starting, be sure that relevant shell variables are already defined and that the excitingscripts directory has already been downloaded, as specified in Tutorial scripts and environment variables. Here is 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.

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. Ground state calculation + GW "single-shot" run

The typical G$_{0}$W$_{0}$ run consists of two parts. First, one needs to perform a groundstate calculation to obtain the self-consistent Kohn-Sham orbitals and potential. Second, using this data as an input, one performs a G$_{0}$W$_{0}$ cycle to calculate quasiparticle energies.

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

<input>
 
  <title>Bulk Silicon: G0W0 example</title>
 
  <structure speciespath="$EXCITINGROOT/species" autormt="false">
    <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_LAPW.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"
    gmaxvr="14"
    ngridk="2 2 2"
    xctype="LSDAPerdew-Wang"
    tetra="true"
    nempty="17"
    lradstep="4"
    >
  </groundstate>
 
  <!--properties>
    <bandstructure character="false">
      <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-->
 
  <gw
    taskname="gw"
    ibgw="1" nbgw="15"
    >
    <mixbasis
      lmaxmb="3"
      epsmb="1.0d-4"
      gmb="1.0">
    </mixbasis>
  </gw>
 
</input>

Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command

$ SETUP-excitingroot.sh

Let us look closer at this input file. First, notice that the species file we are using is $EXCITINGROOT/species/Si_LAPW.xml.

Secondly, rather low values for the groundstate parameters (such as rgkmax, ngridk, etc.) are chosen simply for time-saving reasons. Please note, the following important attribute which must be specified explicitly to make the GW code run: tetra="true". This option instructs exciting to use the tetrahedral integration method (implemented into the library LIBBZINT) for performing the Brillouin-zone integration. The new element gw initializes G$_{0}$W$_{0}$ calculations. Other specific attributes and elements used in this example are

  • taskname, which specifies the task name;
  • ibgw and nbgw, which give the range of resulting quasiparticle states to be calculated;
  • mixbasis, which defines the parameters of the mixed basis.

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

As you can see above, the properties block of the input file has been commened for the calculation performed in this first part. To perform the actual calculation, copy and paste the text from above into a file called input.xml within a directory of your choice and start the calculation by invoking the exciting executable (serial version).

$ excitingser

During the run one can inspect the INFO.OUT file for convergence and other information related to the groundstate calculations. The main output of G$_{0}$W$_{0}$ run is stored in GWINFO.OUT and QPENE-eV.OUT. The former file contains information on the GW program work-flow. In particular, energy band-gap values for KS and GW spectra can be found in the end of this file. The latter file contains quasiparticle energies and other important quantities such as exchange and correlation self-energies, diagonal matrix elements of the exchange-correlation potential, etc.

2. Quasiparticle band structure

Now we will visualize the results to compare the Kohn-Sham (KS) and quasiparticle bandstructures (QP). 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 groundstate calculations from scratch;
  • uncomment the properties element (by removing !-- and --) to obtain the KS bandstructure;
  • in the gw element replace taskname="gw" by taskname="band" to calculate QP bandstructure based on the results obtained from the previous G$_{0}$W$_{0}$ run. In this case the Fourier interpolation is used to evaluate the band energies for the k-point path specified in the properties element.

Run exciting one more time.

$ excitingser

The files BANDS.OUT and BANDS-QP.OUT can now be used to plot the electronic bandstructure "spaghetti" plot. In addition, the corresponding density of states are also calculated and stored in TDOS-KS-QP.OUT. The results can be easily visualized with PLOT-gwbands.py. To execute it, type in your working directory

$ PLOT-gwbands.py

The script produces a file Figure.png which you can visualize on the screen with standard tools. The result should look like the following.

Si_KS-QP.png
Converging the results

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

Description Parameter Info
basis-set size rgkmax defines the quality of eigenvalues and wavefunctions
k-mesh ngridk very crucial and should always be converged
energy cutoff (number of empty states) nempty very crucial and sometimes difficult to converge
mixed-basis size mixbasis defines the quality of the mixed-basis functions
Exercise
  • Once you have successfully run the previous GW calculation, you can play around with computational parameters in order to achieve the desidered convergence for the value of the electronic band gap.
    1. Try to perform a convergence test with respect to the k-mesh size.
    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.
    3. A convergence study for the mixed basis set size is required, especially for materials which are more complex than Si. Typically, the convergence with respect to the mixed-basis parameters is slower than for the previous two parameters.

Literature

  • Tutorial talk (PDF)
  • 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" (submitted to Computer Physics Communication)
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License