by Dmitrii Nabok, Keith Gilmore, & Sven Lubeck for exciting neon
Purpose: In this tutorial you will learn how to perform a basic G$_{0}$W$_{0}$ calculation. As an example, the electronic band structure of bulk Si is calculated. Notice that selfconsistent 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.
 PLOTbandstructure.py: Python visualization script for plotting and comparing energy bands.
 PLOTdos.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: G_{0}W_{0} approximation.
The quasiparticle energies in the G$_{0}$W$_{0}$ approximation are given as a solution of the linearized quasiparticle (QP) equation [Hybertsen1986] as
(1)where $\epsilon_{n\mathbf{k}}$ are the KohnSham (KS) eigenvalues, $\Sigma_{n\mathbf{k}}$ and $V_{n\mathbf{k}}^{xc}$ are, respectively, the diagonal matrix elements of the selfenergy and the exchangecorrelation (xc) potential that is employed in the singleparticle KS Hamiltonian. The QP renormalization factor $Z_{n\mathbf{k}}$ accounts for the energydependence of the selfenergy. The selfenergy $\Sigma(\mathbf{r},\mathbf{r}';\omega)$ is given by [Hedin1965]
(2)where $G_0$ is the noninteracting singleparticle Green function obtained from the KohnSham states and $W_0$ is the dynamically screened Coulomb potential
(3)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 randomphase approximation (RPA).
2. Groundstate calculation
Create a working directory SiGW and move inside it.
$ mkdir SiGW
$ cd SiGW
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="3 3 3" xctype="LDA_PW" > </groundstate> </input>
To perform the actual calculation, copy and paste the text from above into a file called input.xml. Make sure to set $EXCITINGROOT to the correct exciting root directory in the speciespath attribute using the command
$ SETUPexcitingroot.sh
Now, you can start the groundstate calculation by invoking the exciting executable (shared memory version)
$ time exciting_smp
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. Please note that the quality of the groundstate calculations (the starting point) plays an important role when calculating G$_{0}$W$_{0}$ QP corrections. Therefore one should always check the convergence of results with respect to the groundstate parameters (such as rgkmax, ngridk, etc.) in production calculations.
3. GW "singleshot" run
Using the data from the previous run as an input, one performs a G$_{0}$W$_{0}$ calculation to obtain the QP energies.
Disable execution of the groundstate run by changing the attribute do = "fromscratch" to do = "skip" in the groundstate element of the input.xml file.
... <groundstate do="skip" ...> </groundstate> ...
To initialize G$_{0}$W$_{0}$ calculations, copy and paste the following block should into your input file:
... <gw taskname="g0w0" ngridq="3 3 3" nempty="22" ibgw="1" nbgw="10" coreflag="xal" > <mixbasis lmaxmb="3" epsmb="1.d4" gmb="1.0" ></mixbasis> <freqgrid nomeg="12" ></freqgrid> </gw> ...
Now, you can start the calculation by invoking the exciting executable (serial version)
$ time exciting_smp
While the calculation is running, let us look closer at the gw element. 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. It is interesting to notice that the computational setup given for the mixedproduct basis and the frequency grid is rather universal and could be applied for most materials without modification.
Once the calculation has finished, the main output of the G$_{0}$W$_{0}$ run will be stored in GW_INFO.OUT and EVALQP.DAT. The former file contains information on the GW program workflow. In particular, energy bandgap values for KS and QP spectra can be found at the end of this file.
The text file EVALQP.DAT contains all important quantities that are calculated during the G$_{0}$W$_{0}$ run. Below, we describe the structure of this file. The line "kpoint # 1:" starts the block of data and contains the lattice coordinates of the irreducible kpoint and the corresponding Brillouinzoneintegration 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 energies as obtained by considering the exchange part of the selfenergy only (E_HF);
 the G$_{0}$W$_{0}$ quasiparticle energies (E_GW);
 the diagonal matrix elements of the exchange selfenergy (Sx);
 the diagonal matrix elements of the correlation selfenergy (Sc);
 the diagonal matrix elements of the underlying exchangecorrelation 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.
4. Quasiparticle band structure
i) Band structure plot
Now, we will visualize the results to compare the KohnSham (KS) and quasiparticle (QP, G$_{0}$W$_{0}$) band structures. To calculate them, one should introduce the following changes into the input.xml file.
 After the groundstate block, add the new element block properties as indicated below.
... <properties> <bandstructure> <plot1d> <path steps="100"> <point coord="1.0 0.0 0.0" label="Gamma"/> <point coord="0.625 0.375 0.0" label="K"/> <point coord="0.5 0.5 0.0" label="X"/> <point coord="0.0 0.0 0.0" label="Gamma"/> <point coord="0.5 0.0 0.0" label="L"/> </path> </plot1d> </bandstructure> <dos inttype="trilin" winddos="0.5 0.5" nwdos="200" ></dos> </properties> ...
 In the gw element replace taskname = "g0w0" by taskname = "band" to calculate the QP band structure 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 kpoint path specified in the properties element.
... <gw taskname="band" ... > </gw> ...
Finally, run exciting one more time.
$ time exciting_smp
The output files BAND.OUT and BANDQP.OUT can now be used to plot the electronic band structure. This can be easily done using the script PLOTbandstructure.py, which is described in details in The python script "PLOTbandstructure.py". To execute it, type
$ PLOTbandstructure.py a KS GW e 15 10
The script produces a file PLOT.png which you can visualize on the screen with standard tools. The result should look like the following.
An interesting result is obtained if the KS and G$_{0}$W$_{0}$ electronic band structures are compared by setting the energy zero at the valence band maximum (VBM):
$ PLOTbandstructure.py a KS GW e 15 10 z vbM
The result should look like the following.
This plot shows that the main difference between the two band structures is due to the conduction bands which seem to be uniformly shifted upwards by a constant value. This is the theoretical justification of the socalled scissor approximation.
ii) Density of states
Now, we will compute and compare the KohnSham (KS) and quasiparticle (QP, G$_{0}$W$_{0}$) density of states. To calculate them, one should introduce the following changes into the input.xml file.
 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 exciting_smp
The results for the KS and G$_{0}$W$_{0}$ density of states are stored in TDOS.OUT and TDOSQP.OUT correspondingly. To visualize them one can simply execute the script PLOTdos.py (discussed in The python script "PLOTdos.py" ) .
$ PLOTdos.py a KS GW
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 (some of them in a crucial way). In order to keep the runtimes of the above calculations short, certain parameters were reduced from their converged values. The most important parameters are

In the general case, additional attention should be payed to the proper choice of the systemdependent 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.
 Try to perform a convergence test with respect to the size of the k and qmeshs that appear in the both the groundstate and the gw elements.
 Try to run a similar test for the number of empty states. Be careful when increasing this parameter, since the calculation might quickly become very timeconsuming.
 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. GomezAbal, X. Li, C. Meisenbichler, C. AmbroschDraxl, and M. Scheffler, "FHIgap: 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 allelectron G$_{0}$W$_{0}$ quasiparticle energies employing the fullpotential augmented planewave method", Phys. Rev. B 94, 035118 (2016).