GW eigenvalues and Theoretical Spectroscopy
In this tutorial we will assess the suitability of densityfunctional theory (DFT), HartreeFock (HF) and manybody perturbation theory (MBPT) in the GW approach for the calculation of electronic excitations. Most of the calculations will be performed on the ethylene molecule C\(_2\)H\(_4\), with the purpose of comparing the performance of different theoretical approaches with experimental photoemission spectra. Hence, we invite you to organize the results of each exercise in a text file, or in the table reported at the end of this page. Please note that all energies you obtain are relative to the vacuum level.
Reproduced from M. E. Casida^{1}.
(In)adequacy of DFT eigenvalues for the description of charged electronic excitations
Let us first investigate the meaningfulness of the KohnSham eigenvalues by comparison with the experimentally measured ionization potentials (IPs). The first IP and the highestoccupied KohnSham level (HOMO) can be related as follows: $$ IP =  e_\text{KS}(N) $$ where \(N\) refers to the number of electrons in the system and \(e_\text{KS}(N)\) is the \(N\)th KS eigenvalue.
First step is to perform a KohnSham DFT calculation with the PBE exchangecorrelation functional and Tier 2 basis set (tight settings) for ethylene C\(_2\)H\(_4\). You can proceed as follows:
 Generate the geometry file for ethylene from the experimental data (Herzberg 1966^{7}; you can find it from this website https://cccbdb.nist.gov/expgeom1x.asp by searching for C2H4):
atom 0.0000 0.0000 0.6695 C
atom 0.0000 0.0000 0.6695 C
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H
 Generate the template for
control.in
to set up a spinunpolarized DFT calculation using the PBE functional:
xc pbe
 Copy the required basis sets from tight into the
control.in
file  Start a parallel calculation by typing:
mpirun n 4 aims.210716_2.scalapack.mpi.x  tee output

Compare the three KS eigenvalues  corresponding to the highest occupied molecular orbital (HOMO) HOMO1, HOMO2, and HOMO3  with the first three experimental ionization potentials (IP1, IP2, and IP3) given with this table^{2}:
C\(_2\)H\(_4\) Experimental IP1 10.68 eV IP2 12.80 eV IP3 14.80 eV
Optional: If you like, you can repeat the above steps for a different functional, e.g. you could use xc hf
(the HartreeFock Method) and compare the results for this setting.
Electron removal energies from deltaSCF
The ionization potentials (\(IP\)) of C\(_2\)H\(_4\) can be evaluated with the deltaselfconsistentfield (\(\Delta\)SCF) approach^{2}. Following the definition of the IP,
the total energy difference between the neutral (\(E^{PBE}_{tot}(N)\) ) and positively (\(E^{PBE}_{tot}(N1)\)) charged species is computed from two separate DFT PBE total energy calculations where \(N\) is the number of electrons of the neutral molecule.
To evaluate IP, \(E^{PBE}_{tot}(N)\) can be extracted from the output file of Section 1. In addition we need to compute \(E^{PBE}_{tot}(N1)\), which requires a second DFT calculation. Note that the charged molecule requires a spinpolarized treatment. For more details regarding spinpolarized simulations in FHIaims, cf. Part 2 of Basics of running FHIaims tutorial. You can proceed as follows:
 Copy the input files from the first exercise.
 Modify the
control.in
file and set the necessary flags for performing a spinpolarized calculation of the charged molecule:
xc pbe
spin collinear
charge +1
 Also modify the
geometry.in
for \(C_2H_4\) to account for the spin initialization (for more details, cf. Part 2 of Basics of running FHIaims tutorial. Since we have a charge +1 in our system, the difference between the number of electrons in the spin up and the spin down channel will be one: \(N_\uparrowN_\downarrow=1\). Based on a basic chemical intuition, we place an initial moment of 0.5 on each of the carbon atoms:
atom 0.0000 0.0000 0.6695 C
initial_moment 0.5
atom 0.0000 0.0000 0.6695 C
initial_moment 0.5
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H
atom 0.0000 0.9289 1.2321 H

Start the calculation for the charged system. For the neutral system, you can reuse the output file of the previous section

Compute the ionization potential of C\(_2\)H\(_4\) using the formula for the IP. How do these values compare to the bare PBE eigenvalue and experiment?
Perturbative G\(_0\)W\(_0\) and quasiparticle corrections
An improved description of charged electronic excitations is obtained by the perturbative inclusion of manybody effects through the selfenergy \(\Sigma\). In the GW approximation^{3} the selfenergy is calculated as:
where \(G(r, r^\prime, \omega)\) is the oneparticle Green's function and \(W(r, r^\prime, \omega)\) is the screened Coulomb interaction^{3} ^{4} (also see a recent educational review article^{5} for details). The GW selfenergy can be used to perturbatively correct the DFT or HF eigenvalues by means of the quasiparticle equation:
where \(\Sigma_{\mathrm x}\) is the exactexchange operator, and \(\Sigma_{\mathrm c}^{GW}\) is the correlation part of the GW selfenergy. \(V_{xc}^{KS}\) is the exchangecorrelation potential of the preceding DFT/HF calculation, \(\epsilon_i^{KS}\) and \(\psi^{KS}_i\) are the corresponding eigenvalues and eigenvectors. This approximation is known as G\(_0\)W\(_0\) or oneshot GW, because the selfenergy is calculated only once, whereas a more rigorous approach would require a fully selfconsistent evaluation of \(\Sigma\). Since the quasiparticle energies are evaluated perturbatively on top of a preceding singleparticle calculation (generally DFT or HartreeFock), the G\(_0\)W\(_0\) approach depends on the initial reference calculation. To refer to PBE as starting often the term G\(_0\)W\(_0\)@PBE is used.
G\(_0\)W\(_0\)@PBE quasiparticle energies
Lets actually perform a G\(_0\)W\(_0\) calculation to obtain the quasiparticle energies of ethylene. For this exercise, proceed along the following steps:
 Copy the
control.in
andgeometry.in
files from the first section.  Modify the
control.in
file including the following flags:xc pbe qpe_calc gw anacon_type 1
The keyword anacon_type
specifies the analytical continuation fo the selfenergy. The type can be either 0
(the twopole fitting) or 1
(Padé Approximation). In principle, the twopole fitting is numerically more robust, however, both methods have the downsides. Usually, the Padé Approximation performs well for the highestoccupied states of the system. A better solution for lower lying states (core levels) are the contour deformation methods, which will be discussed later in the core level chapter. A detailed study of the keyword anacon_type
can be found in the GW Compendium^{5}.
Now start the FHIaims calculation (mpirun n ...
). The output file will contain a table  similar to this figure  with the quasiparticle corrections to the singleparticle eigenvalues. Extract the quasiparticle energies for the HOMO, HOMO1, HOMO2, and HOMO3 levels and compare them with the results
from the previous calculations.
G\(_0\)W\(_0\) basis set convergence
All beyondDFT methods share a severe issue: Observables converge very slowly with the number of basis functions.
In this exercise, we will plot the convergence of the first G\(_0\)W\(_0\)@PBE quasiparticle energy (i.e. the G\(_0\)W\(_0\)@PBE HOMO level) for ethylene using the Tier 1, Tier 2, and Tier 3 basis sets. Proceed as follows:
For each tier, create a separate folder. Copy control.in
and geometry.in
of the previous exercise to each folder. For tier 1, please only uncomment the tier1 basis sets for all elements. For tier 2, additionally uncomment the tier2 basis sets, and so on. Once ready, proceed with executing the calculations.
As a result you should find two plots similar to the following (basis set tier vs highest occupied state). Look at the change of the highest occupied state from tier to tier. We find that the HOMO level calculated with the GW method converges much slower than the corresponding DFT HOMO eigenvalue. Further, even with tight species defaults the GW HOMO level is not converged. (What is the origin of the qualitative differences between the convergence behavior in PBE and G\(_0\)W\(_0\)@PBE?)
The FHIaims standard NAO species defaults are in general not suited to systematically obtain converged observables with beyondDFT methods. To achieve this a method called infinitebasisset extrapolation along with the evenly tempered Dunning Gaussian basis sets is used. We will introduce this method in the following subsection.
Infinitebasisset extrapolation
(Sometimes this method is also referred to as "extrapolation to the complete basis set limit")
Here, we will shortly introduce the concept of the infinitebasisset extrapolation. We will use the Dunning basis sets ccpVNZ, which you can find in the species_defaults
directory of the source code:
species_defaults/nonstandard/gaussian_tight_770/ccpVNZ
They are ordered as N\(\zeta\) (or in the folder as NZ), where N refers to the highest main quantum number used in the basis set. N is to be replaced by D (double), T (triple), Q (quadruple), 5.
For each of the available ccpVNZ
basis sets create a subfolder and copy the C2H4 geometry.in
into it. Also for all of the folders prepare a control.in
file with the following content:
xc pbe
qpe_calc gw
anacon_type 1
and attach the corresponding species defaults from the ccpVNZ
folder for the TZ, QZ, and 5Z basis sets, e.g. for TZ attach the species defaults for H and C from the folder nonstandard/gaussian_tight_770/ccpVDZ
to the control.in
and similar for the other N\(\zeta\) species defaults.
Run all of the calculations and extract the GW HOMO quasiparticle level from the output file and the number of basis functions. For the number of basis functions \(1/N_\text{basis}\) look in the main output for the following line starting with:
 Total number of basis functions :
Plot the GW levels over \(1/N_\text{basis}\) and perform a linear fit for the data points. You should be able to reproduce a plot as shown below. We found for the G0W0@PBE HOMO level an infinitebasisset level of 10.42 eV.
Visualization of the G\(_0\)W\(_0\) spectra
The quasiparticle energies calculated in the previous task are the peak positions of the molecule’s electronic excitation spectrum. Now, you will visualize the spectra for the G0W0@PBE calculation for C2H4.
 Use the script
create_spectrum.py
, located in the folderPart1/solutionsGW/G0W0atPBE
, to extract the quasiparticle energies from the output files and transform them into a spectrum where the energies were broadened by 0.05 eV to facilitate the comparison with experimental data. Supply the name of your FHIaims output file as the first argument and call the script withpython3 create_spectrum.py output 1> spectrum.dat
where the output spectrum was redirected from the terminal to the file spectrum.dat
.

Visualize the generated spectra (e.g. with
climsxyplot
,qtiplot
, orxmgrace
, ...) together with the experimental photoemission spectroscopy data provided in the filePart1/solutionsGW/G0W0atPBE/C2H4PES.dat
.climsxyplot spectrum.dat C2H4PES.dat

How large is the deviation from the experimental HOMO levels?
Selfconsistent GW
In this exercise, you will perform a fully selfconsistent GW calculation. Differently from G\(_0\)W\(_0\), the Green's function is calculated by solving the Dyson's equation selfconsistently. The Dyson equation relates the input Green's function \(G_0\) to the GW Green's function \(G\) via the selfenergy \(\Sigma\)
or in inverted form
where we used the shorthand notation \(1\equiv \lbrace {\mathbf r}_1,t_1, \sigma_1\rbrace\)  see e.g. the Hedin paper^{4} for an introduction. \(v_{\mathrm H}\) is the Hartree potential. We refer to the publication by Caruso et al.^{6} for details of the scGW implementation in FHIaims.
Spectral function from the selfconsistent Green's function
To perform a selfconsistent GW calculation for C\(_2\)H\(_4\), create a new directory and copy the input files from the first section. Modify the first part of the control.in
file:
xc pbe
sc_self_energy scgw
spin none
and choose Tier 1 settings for all species at the bottom of the control.in
file.
After running FHIaims, the file spectrum_sc.dat
will be created.
The file spectrum_sc.dat
contains the spectral function calculated from the selfconsistent GW Green's function
where \(G\) has been determined selfconsistently from the Dyson equation.
You can visualize the spectral function using xmgrace
, qtiplot
or another available plotting tool. Also, clims
can be used as:
climsxyplot spectrum_sc.dat
The first three IPs of C\(_2\)H\(_4\) must be extracted directly from the spectral function (like in experiment) by reading of the peak positions.

M. E. Casida. In: Recent Advances in Density Functional Methods (1995). ↩

G. Bieri and L. Asbrink. In: J. of Electron Spectrosc. and Rel. Phenom. 20 (1980), p. 149. ↩↩

D. Golze, M. Dvorak, and P. Rinke. “The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy”. In: Front. Chem. 7 (2019), p. 377. ↩↩

F. Caruso et al. In: Phys. Rev. B 88 (2013), p. 75105. ↩

Herzberg, G., Electronic spectra and electronic structure of polyatomic molecules,Van Nostrand, New York, 1966 ↩