# GW eigenvalues and Theoretical Spectroscopy

In this tutorial we will assess the suitability of density-functional theory (DFT), Hartree-Fock (HF) and many-body 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. Casida1.

## (In)adequacy of DFT eigenvalues for the description of charged electronic excitations

Let us first investigate the meaningfulness of the Kohn-Sham eigenvalues by comparison with the experimentally measured ionization potentials (IPs). The first IP and the highest-occupied Kohn-Sham 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 Kohn-Sham DFT calculation with the PBE exchange-correlation functional and Tier 2 basis set (tight settings) for ethylene C$$_2$$H$$_4$$. You can proceed as follows:

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 spin-unpolarized 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) HOMO-1, HOMO-2, and HOMO-3 -- with the first three experimental ionization potentials (IP-1, IP-2, and IP-3) given with this table2:

C$$_2$$H$$_4$$ Experimental
IP-1 10.68 eV
IP-2 12.80 eV
IP-3 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 Hartree-Fock Method) and compare the results for this setting.

## Electron removal energies from delta-SCF

The ionization potentials ($$IP$$) of C$$_2$$H$$_4$$ can be evaluated with the delta-self-consistent-field ($$\Delta$$SCF) approach2. Following the definition of the IP,

$I=E^{PBE}_{tot}(N-1) - E^{PBE}_{tot}(N)\quad,$

the total energy difference between the neutral ($$E^{PBE}_{tot}(N)$$ ) and positively ($$E^{PBE}_{tot}(N-1)$$) 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}(N-1)$$, which requires a second DFT calculation. Note that the charged molecule requires a spin-polarized treatment. For more details regarding spin-polarized simulations in FHI-aims, cf. Part 2 of Basics of running FHI-aims 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 spin-polarized 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 FHI-aims 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_\uparrow-N_\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 re-use 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 quasi-particle corrections

An improved description of charged electronic excitations is obtained by the perturbative inclusion of many-body effects through the self-energy $$\Sigma$$. In the GW approximation3 the self-energy is calculated as:

$\Sigma^{GW}(r, r^\prime, \omega) = \frac{i}{2\pi}\int d\omega' G(r, r^\prime, \omega') W(r, r^\prime, \omega'+\omega)\quad,$

where $$G(r, r^\prime, \omega)$$ is the one-particle Green's function and $$W(r, r^\prime, \omega)$$ is the screened Coulomb interaction3 4 (also see a recent educational review article5 for details). The GW self-energy can be used to perturbatively correct the DFT or HF eigenvalues by means of the quasi-particle equation:

$\epsilon_i^{QP} = \epsilon_i^{KS} - \left\langle \psi_i^{KS} \middle| \hat{V}_{xc}^{KS} - \hat{\Sigma}_{\mathrm c}^{GW}(\epsilon_i^{QP}) - \hat \Sigma_{\mathrm x} \middle| \psi^{KS}_i \right\rangle \quad,$

where $$\Sigma_{\mathrm x}$$ is the exact-exchange operator, and $$\Sigma_{\mathrm c}^{GW}$$ is the correlation part of the GW self-energy. $$V_{xc}^{KS}$$ is the exchange-correlation 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 one-shot GW, because the self-energy is calculated only once, whereas a more rigorous approach would require a fully self-consistent evaluation of $$\Sigma$$. Since the quasi-particle energies are evaluated perturbatively on top of a preceding single-particle calculation (generally DFT or Hartree-Fock), 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 quasi-particle energies of ethylene. For this exercise, proceed along the following steps:

• Copy the control.in and geometry.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 self-energy. The type can be either 0 (the two-pole fitting) or 1 (Padé Approximation). In principle, the two-pole fitting is numerically more robust, however, both methods have the downsides. Usually, the Padé Approximation performs well for the highest-occupied 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 Compendium5.

Now start the FHI-aims calculation (mpirun -n ...). The output file will contain a table -- similar to this figure -- with the quasi-particle corrections to the single-particle eigenvalues. Extract the quasi-particle energies for the HOMO, HOMO-1, HOMO-2, and HOMO-3 levels and compare them with the results from the previous calculations.

### G$$_0$$W$$_0$$ basis set convergence

All beyond-DFT 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 quasi-particle 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 tier-1 basis sets for all elements. For tier 2, additionally uncomment the tier-2 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 FHI-aims standard NAO species defaults are in general not suited to systematically obtain converged observables with beyond-DFT methods. To achieve this a method called infinite-basis-set extrapolation along with the evenly tempered Dunning Gaussian basis sets is used. We will introduce this method in the following subsection.

### Infinite-basis-set 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 infinite-basis-set extrapolation. We will use the Dunning basis sets cc-pVNZ, which you can find in the species_defaults directory of the source code:

species_defaults/non-standard/gaussian_tight_770/cc-pVNZ


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 cc-pVNZ 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 cc-pVNZ folder for the TZ, QZ, and 5Z basis sets, e.g. for TZ attach the species defaults for H and C from the folder non-standard/gaussian_tight_770/cc-pVDZ to the control.in and similar for the other N$$\zeta$$ species defaults.

Run all of the calculations and extract the GW HOMO quasi-particle 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 infinite-basis-set level of -10.42 eV.

### Visualization of the G$$_0$$W$$_0$$ spectra

The quasi-particle 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 folder Part-1/solutions-GW/G0W0-at-PBE, to extract the quasi-particle 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 FHI-aims output file as the first argument and call the script with
python3 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 clims-xyplot, qtiplot, orxmgrace, ...) together with the experimental photo-emission spectroscopy data provided in the file Part-1/solutions-GW/G0W0-at-PBE/C2H4-PES.dat.

clims-xyplot spectrum.dat C2H4-PES.dat

• How large is the deviation from the experimental HOMO levels?

## Self-consistent GW

In this exercise, you will perform a fully self-consistent GW calculation. Differently from G$$_0$$W$$_0$$, the Green's function is calculated by solving the Dyson's equation self-consistently. The Dyson equation relates the input Green's function $$G_0$$ to the GW Green's function $$G$$ via the self-energy $$\Sigma$$

$G(1,2) = G_0 (1,2) + \int d34\; G_0(1,3)\left[v_{\mathrm H}(3)\delta(3,4)+\Sigma(3,4)\right]G(4,2) \quad,$

or in inverted form

$G^{-1}(1,2) = G_0^{-1} (1,2) - v_{\mathrm H}(1) \delta(1,2) - \Sigma(1,2)\quad,$

where we used the shorthand notation $$1\equiv \lbrace {\mathbf r}_1,t_1, \sigma_1\rbrace$$ -- see e.g. the Hedin paper4 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 FHI-aims.

### Spectral function from the self-consistent Green's function

To perform a self-consistent 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 FHI-aims, the file spectrum_sc.dat will be created. The file spectrum_sc.dat contains the spectral function calculated from the self-consistent GW Green's function

$A(\omega) = -\frac{1}{\pi}\int \lim_{{r^\prime} \rightarrow {r}} Im G({r}, {r^\prime},\omega)d{r}\quad$

where $$G$$ has been determined self-consistently from the Dyson equation. You can visualize the spectral function using xmgrace, qtiplot or another available plotting tool. Also, clims can be used as:

clims-xyplot 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.

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

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

3. Bagus. In: Phys. Rev. A 139 (1965), p. 619.

4. Hedin. In: Phys. Rev. A 139 (1965), p. 796.

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

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

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