by Dmitrii Nabok for exciting beryllium

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 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 Tutorial scripts and environment variables. Here, you can find a list of the scripts which are relevant for this tutorial, with a short description.

• EXECUTE-single.sh: (Bash) shell script for running a single exciting calculation.
• 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 consider in this tutorial.

<input>
 
   <title>Bulk Silicon: G0W0 example</title>
 
   <structure speciespath="$EXCITINGROOT/species">
      <crystal scale="10.26">>
         <basevect>0.0   0.5   0.5</basevect>
         <basevect>0.5   0.0   0.5</basevect>
         <basevect>0.5   0.5   0.0</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"
      gmaxvr="14"
      ngridk="2 2 2"
      xctype="LDA_PW">
   </groundstate>
 
   <gw
      taskname="gw"
      ngridq="2 2 2"
      nempty="50"
      ibgw="1" 
      nbgw="15">
      <mixbasis/>
      <freqgrid/>
   </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

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)

$ EXECUTE-single.sh gw-test-1

The output files of this run are stored inside the directory gw-test-1. In order to analyze the results, move to the directory gw-test-1 when the execution has been completed.

$ cd gw-test-1

After the run, one can inspect the INFO.OUT file for checking the convergence behaviour of the SCF cycle and for other information related to the groundstate calculation. The main output of G$_{0}$W$_{0}$ run is stored in GWINFO.OUT and EVALQP.TXT. 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 (QP) bandstructures. To calculate them, one should introduce the following changes into the input.xml file inside the directory gw-test-1. (For being on the safe side, make first a safety copy of the old input.xml to input-gs.xml).

• In the groundstate element change do = "fromscratch" to do = "skip" to avoid running the groundstate 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 = "gw" 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 inside gw-test-1.

$ 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 inside the gw-test-1 directory

$ PLOT-gwbands.py

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

Converging results

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

In the general case, additional attention should be payed to the proper choice of the system-dependent parameters used in gw/mixbasis and gw/freqgrid 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 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.

Literature

• Tutorial talk (PDF) at the HoW exciting! 2012 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).