# Calculations with SPEX

*Note*: This tutorial of the FLEUR Workshop 2023 is based on a special version (06.00pre6) of the SPEX code.

## 1. *GW* calculations

The *GW* self-energy is an approximation to the electronic self-energy used in many-body perturbation theory. This theory allows excited state properties
to be calculated on the basis of a mean-field system, which is usually taken from a previous DFT calculation, e.g., using the local-density approximation (LDA).
One says that LDA is the *starting point* for the *GW* calculation.

Thus, in order to carry out a *GW* calculation, we first have to perform a self-consistent DFT calculation with the FLEUR code. For the latter, you will use
your compiled FLEUR program.

### 1.1 Installation

As a first step, we are going to install the SPEX code. Copy the file "spex06.00pre6.tgz" to your working directory and unpack it:

cp spex06.00pre6.tgz work

tar xvzf spex06.00pre6.tgz

Enter the directory `spex06.00pre6`

:

cd spex06.00pre6

In order to configure the compilation for your system (choice and test of compiler, search for libraries, etc.), you have to run the configure script. The script allows for customization through command line arguments:

./configure --with-wan --with-dft=fleurR5.1 --prefix=$HOME --with-hdf5=/usr/lib/x86_64-linux-gnu/hdf5/serial/

In your docker image `judft/fleur:2021`

, the HDF5 library is not in a standard directory. Therefore, we have to specify the HDF5 root directory with the option `--with-hdf5`

.

A short explanation of the options:

`--with-wan`

: Compile with Wannier support (searches for and links in the Wannier90 library),`--with-hdf5=PATH`

: Specify HDF5 root directory,`--prefix=PREF`

: Install executables (and shell scripts) into PREF/bin (defaults to PREF=/usr/local; we choose PREF to be your home directory).`--enable-mpi`

: Compile parallel version using MPI 3.1 (not needed for the tutorial),`--with-dft=fleurR5.1`

: Compile the interface for the MaX Release 5.1 of FLEUR (still compatible with 6.1),

We do not have ScaLAPACK in the prepared Docker image. Therefore, we cannot compile an MPI-parallelized version of SPEX.

Compile the source! This will take a few minutes. At the end, you will be asked to install the binaries.

make

For faster compilation in parallel, use `make -j`

.

You have just compiled the executables:

: For systems with inversion symmetry (and without SOC)`spex.inv`

: For systems without inversion symmetry (or with SOC)`spex.noinv`

: Extraction utility`spex.extr`

These executables can be installed by

make install

The reason for having two executables, one for systems with inversion symmetry (`spex.inv`

) and the other for systems without (`spex.noinv`

),
is that many quantities can be declared as real-valued instead of complex-valued arrays in the case of inversion symmetry,
which speeds up the computation considerably.

Together with the aforementioned executables, we have also installed the "SPEX launcher script" called `spex`

.
The executables (`spex.inv`

or `spex.noinv`

) can be called directly. Alternatively, the launcher script allows simply running `spex`

in both cases.
You can choose one of the two options. In the following, we will use `spex`

for simplicity.

### 1.2 Silicon

#### 1.2.1 DFT as the mean-field starting point

Your first *GW* calculation will be on bulk Silicon. As a first step, we must perform a DFT calculation that creates the necessary input for *GW*.
Silicon crystallizes in a diamond lattice (fcc lattice with a two-atom basis) with a lattice constant of 5.43 Å (10.26 Bohr).

You can create your own FLEUR input file "inp_Si" as you have already learned it or use the prepared one "Si/inp_Si". If you create it yourself, it should contain the line `&expert spex=1 /`

, which tells FLEUR to create some additional files (containing structural information etc.) for the subsequent *GW* calculation.

*Important*: We need to use version MaX R6.1 of Fleur. The executables are called "fleur_MPI-MaX-R6.1" and "inpgen-MaX-R6.1".

We create the Fleur input file "inp.xml" by running "inpgen-MaX-R6.1":

```
inpgen-MaX-R6.1 -f inp_Si
```

Then, perform a DFT calculation for Silicon with "fleur_MPI-MaX-R6.1", which produces a self-consistent density and effective potential. Make sure that the calculation has converged.

#### 1.2.2 Special kpoint set

As the *GW* self-energy is nonlocal, its evaluation involves summations over two kpoints, **k** and **q**, and their sum **k+q**.
Thus, we must make sure that the set of kpoints forms an equidistant regular mesh, which contains **k+q** as an element if **k** and **q** are elements of the set.
Furthermore, the Γ point (k=0) must be included. (In other words, the kpoint set must form a group under vector addition.)
The explicit inclusion of the Γ point is necessary because of the divergence of the Coulomb interaction in the long-wavelength limit.

We can now proceed in two ways:

- SPEX creates the new kpoint set. FLEUR runs a
*one-shot*calculation on that new kpoint set. Finally run the*GW*calculation with SPEX. - Directly run the
*GW*calculation with SPEX with an additional preprocessing step (`JOB DFT0`

) in the input file "spex.inp".

Obviously, the second procedure is much simpler as it spares us two runs (with SPEX and FLEUR). The first procedure is more flexible, though. For simplicity, we will use the second procedure in this tutorial.

#### 1.2.3 Perform the *GW* calculation

Just like FLEUR, SPEX uses an input file containing all kinds of parameters for the SPEX calculation.
This is a simple input file for a *GW* calculation of silicon:

```
BZ 4 4 4
KPT X=1/2*(0,1,1) L=1/2*(0,0,1)
JOB DFT0 GW 1:(1,2,5) X:(1,3,5) L:(1-3,5)
NBAND 80
```

A few remarks about the input file: The order of keywords (`BZ`

, `KPT`

, ...) is not important. A keyword must not be given more than once. Empty lines are allowed. Everything after `#`

is interpreted as a comment (and ignored). All keywords can also be given in lowercase letters (e.g., `kpt`

instead of `KPT`

).

The second line of the above input file defines some kpoint labels, here for the X and L point. These labels are used in the job definition (third line):
`GW`

stands for the *GW* approximation; then the states (on which the quasiparticle correction is to be applied) are defined
with the syntax "kpoint:(band indices)". The kpoint can be specified by its index (see list of kpoints in the output, e.g, "1" for the Γ point)
or by a kpoint label.

We have chosen to define the band indices in the job definition in such a way that, in the case of degeneracy, only one state out of a group of degenerate states is
calculated. For example, the second, third, and fourth state at Γ (k-point index 1) are degenerate. You can also include all states with `1:(1-5)`

.

The keyword `NBAND`

defines the number of bands (occupied and unoccupied) used in the *GW* calculation. Obviously, this number must not exceed the number of LAPW basis functions.
If you omit the `NBAND`

definition, SPEX takes into account all available bands.

Now we can perform the *GW* calculation with SPEX.
When SPEX is run, it writes information to the two standard streams of the shell: standard output (stdout) and standard error (stderr).
The first stream is used for the normal output, the second for warning, info, and error messages. Normally, both streams are directed to the screen, but they can also
be directed to a file using suitable shell commands. Here are some examples:

: Both output and error streams are written to the screen (not recommended).`spex`

: The output stream is written to the file "spex.out", which you can later open with your favourite editor. The error stream is directed to the screen.`spex > spex.out`

: Both streams are written to the file "spex.out".`spex &> spex.out`

: The output stream is written to the screen but also to the file "spex.out" (recommended). This allows you to monitor the calculation while it is running. The error stream is written to the screen.`spex | tee spex.out`

: The output stream is written to the terminal, but also to the file "spex.out". The error stream is written to the file "spex.err".`spex 2> spex.err | tee spex.out`

#### 1.2.4 A longer input file

Before we investigate the output file, let us have a look at the input file again.
The one used in the previous section was a minimal one. A typical *GW* calculation requires a quite large number of parameters. When these parameters
are not given in the input file (as above), SPEX uses default parameters or determines reasonable ones from the DFT data. Of course, all parameters can also be
defined explicitly in the input file. A longer input file for the above calculation could look like this:

```
# silicon bulk
BZ 4 4 4
KPT X=1/2*(0,1,1) L=1/2*(0,0,1)
JOB DFT0 GW 1:(1,2,5) X:(1,3,5) L:(1-3,5)
NBAND 80
SECTION MBASIS
LCUT 4
GCUT 3.6
SELECT 2;3
OPTIMIZE MB 4.0
END
SECTION SUSCEP
HILBERT 50 30
END
SECTION SENERGY
MESH 10+2 10.0
CONTINUE
END
```

*Note*: The empty lines and the two-character indentation in the sections are used only for clarity. Everything after `#`

is interpreted as a comment.

Some explanations:

The section `MBASIS`

contains the definition of the mixed product basis. The keyword `HILBERT`

defines the real-frequency mesh for the Hilbert transformation
of the susceptibility (polarization function).
The section `SENERGY`

defines parameters for the evaluation of the self-energy. The keyword `MESH`

defines the imaginary-frequency mesh.
The keyword `CONTINUE`

tells SPEX to use analytic continuation (PadÃ© extrapolation) of the self-energy from the imaginary to the real-frequency axis.
(This is the default.)

#### 1.2.5 Investigating the results in "spex.out"

Have a look at the output file "spex.out". After reading the parameters (routine "getinput"), running a one-shot DFT calculation (routine "iterate"),
and setting up the mixed product basis (routine "mixedbasis"), SPEX calculates the Coulomb matrix (routines "structureconstant"
and "coulombmatrix"). Then, the correlation self-energy is constructed
in a kpoint loop, which consumes most of the CPU time. For each kpoint the susceptibility (or polarization function *P*) is calculated
(routine "susceptibility"), from which the screened interaction *W*=(1-*vP*)^{-1} is obtained, and the correlation self-energy
Σ=*iG*(*W-v*) is updated. Finally, the self-energy is analytically continued to the real-frequency axis and the quasiparticle equation is solved.
The results are written in a table:

```
Bd vxc sigmax sigmac Z KS HF GW lin/dir
1 -12.11172 -19.10561 6.57060 0.64929 -6.20997 -13.20387 -6.37436 0.72185
1.04269 -0.10593 -6.37825 0.71237
2 -13.21577 -14.83300 1.14256 0.76185 5.60804 3.99081 5.24641 0.00147
0.00088 -0.00169 5.24744 0.00162
5 -11.42949 -6.88796 -4.16444 0.75863 8.16243 12.70396 8.44847 -0.00690
-0.00667 -0.00489 8.44752 -0.00710
```

The columns are

: the band index,`Bd`

: expectation value of v`vxc`

_{xc}, e.g., -12.11172 eV,: expectation value of the exchange self-energy Σ`sigmax`

_{x}, e.g., -19.10561 eV,: expectation value of the correlation self-energy Σ`sigmac`

_{c}, e.g., (6.57060+1.04269i) eV (note that the self-energy is complex),: renormalization factor`Z`

*Z*, e.g., (0.64929-0.10593i) eV,: Kohn-Sham value (from the previous DFT calculation), e.g., -6.20997 eV,`KS`

: Hartree-Fock value (one-shot), e.g., -13.20387 eV,`HF`

:`GW`

*GW*quasiparticle values, e.g., (-6.37436+0.72185i) eV (solution of a linearized, i.e., approximated, quasiparticle equation) and (-6.37852+0.71237i) eV (exact solution of the non-linear quasiparticle equation). Usually, there is only little difference between the two values.

Look up the Kohn-Sham, HF, and *GW* direct and indirect gaps (i.e., Γ→Γ and Γ→X) and compare with the experimental values,
3.4 eV and 1.17 eV, respectively. Note that the conduction band minimum is not exactly at the X point, though (see below).
The fourth band is the highest valence band: There are two silicon atoms per unit cell. Each contributes four electrons. So, there are eight electrons in the
unit cell. (You can find this line in the output: `Number of valence electrons: 8`

.)
Then, dividing by two (two electron spins) gives four.

#### 1.2.6 Contour integration

Instead of analytic continuation, you can use contour integration for evaluating the self-energy. This is more accurate, but takes more CPU time.

Simply define `CONTOUR 0.01`

in the file "spex.inp" inside a section `SENERGY`

:

```
SECTION SENERGY
CONTOUR 0.01
END
```

(or replace `CONTINUE`

by `CONTOUR 0.01`

if you are using the input file of Section 1.2.4). Then rerun the calculation.
(You might want to use a different output file.) The argument "0.01" tells SPEX to calculate the self-energy at two energies symmetrically around
each KS energy ε (ε-0.01 htr and ε+0.01 htr), which enables a linearization of the self-energy.
Therefore, we only get the "linearized" quasiparticle result.
(We can also obtain both solutions again by specifying a range of values instead of the single argument, e.g., `CONTOUR +-{-0.05:0.05,0.01}`

.)

Compare the *GW* results with the previous calculation.

#### 1.2.7 Band convergence

You can play around with the parameters, e.g., increase the number of bands (keyword `NBAND`

) to 100 and 120 and check the band-gap convergence.
You can also experiment with other parameters, e.g., those of the mixed product basis.

#### 1.2.8 Denser kpoint sets

If you have time, you can try denser kpoint sets, e.g., 6×6×6.

#### 1.2.9 Band structure calculation

The calculation of a *GW* band structure is not as straightforward as in the case of DFT, because we cannot arbitrarily add kpoints to the kpoint set.
On the other hand, if we restrict ourselves to the set of kpoints defined by the keyword `BZ`

, the band structure will have very low k resolution.
Anyway, this will be our first calculation.

##### 1.2.9.1 Keyword `PATH`

Let us make a band structure calculation for the path from L over Γ to X. In principle, we would have to identify all kpoint indices along this path and define separate jobs in the job definition, but SPEX can do this for you. To define the kpoint path, add the line

```
KPTPATH (L,1,X)
```

to "spex.inp". The *GW* quasiparticle energies are to be calculated on each kpoint along this path.
To this end, use the keyword `PATH`

in the job definition, replacing the previous definition by

```
JOB DFT0 GW PATH:(1-8)
```

(To get a somewhat finer band structure, you can also change `BZ 4 4 4`

to `BZ 6 6 6`

.)

When you now run SPEX, it will automatically choose the kpoints that lie on the path defined by `KPTPATH`

and calculate the respective quasiparticle energies.
SPEX writes the file "bandsp" containing the band-structure data as an *xy* dataset. The band structure can then be plotted with any plot program. In the Docker image, we have two possibilities:

**Python/Matplotlib**: File → New → Console (or Notebook) → Select kernel "Python 3"**Gnuplot**: File → New → Console (or Notebook) → Select kernel "Gnuplot"

*Note*: We will only suggest Gnuplot plot commands in the following. Using the Gnuplot kernel, you can pretty much work as in "classic Gnuplot". For example, enter `plot "bandsp" with linespoints`

to plot the band structure. To submit a command, press SHIFT+ENTER. Make sure that you are in the correct directory. `pwd`

(same as in BASH) tells you the current directory, and `cd "<path>"`

(note the quotation marks!) lets you change the directory.

As expected, the band structure is not smooth. One solution would be to use denser kpoint sets, but this would entail very expensive calculations.

A simple solution is to interpolate between the explicit energies using splines. This can be done with the extraction utility "spex.extr":

```
spex.extr g -b -i spex.out > band
```

The extraction utility takes the values from the output file (here named "spex.out"). The interpolated band structure is written to the file "band". Of course, the interpolation is purely mathematical and does not always result in a physically correct band structure.

*Note*: The absolute quasiparticle energies (*y* axis) are plotted. The Fermi energy is *not* subtracted. The valence-band maximum is at around 5.2 eV.

##### 1.2.9.2 Additional kpoints (`KPT +...`

)

Another possibility is to use a feature that allows us to add arbitrary kpoints, one at a time.
Each additional kpoint **q** must be accompanied with all shifted kpoints **q+k**.
So, for each single kpoint **q** we must run FLEUR to generate wavefunctions
and energies at the shifted kpoint set and then run SPEX. Let us do this for a single kpoint first.

We define the conduction-band minimum by

```
KPT +=1/2*(0,0.75,0.75)
```

*Note*: A keyword (such as `KPT`

) is not allowed to be used more than once in the input file. Therefore, if you already have a `KPT`

definition in your input file, either replace it with the one above or simply add the argument `+=1/2*(0,0.75,0.75)`

to it (e.g., `KPT X=1/2*(0,1,1) L=1/2*(0,0,1) +=1/2*(0,0.75,0.75)`

).

The label `+`

is a special kpoint label for additional kpoints.
Before running SPEX, the job definition line should be replaced by

```
JOB GW 1:(1-5) +:(1-5)
```

The bands 1-5 are calculated for the Γ point and the additional kpoint. After the calculation is done, you can read the fundamental band gap from the output file. It should be smaller by about 0.1 eV than the result for the Γ→X transition calculated above.

*Note*: It is recommendable to reduce the kpoint set to `BZ 4 4 4`

for the following excercises.

Let us now make a full band structure calculation for the same path in section 1.2.9.1. Specify eight bands with

```
JOB GW +:(1-8) CYCLE
```

We have added `CYCLE`

to the end of the job line. This will tell SPEX to run *GW* multiple times, each time with a new kpoint from the kpoint list defined in `KPTPATH`

. For this, we have to change the `KPT`

definition in the following way

```
KPT X=1/2*(0,1,1) L=1/2*(0,0,1) +=PATH
```

Also add

```
RESTART
```

to the input file. SPEX will then store the screened interaction to harddisc, when it calculates the first kpoint, and reuse it for all other kpoints.

The file "bandsp" again contains the data in the form of an *xy* dataset. The band structure can be plotted with one of the three possibilities described above. You can also compare to the spline-interpolated band structure from section 1.2.9.1.

Note that the conduction band minimum is, in fact, between the Γ and X point at about k=1/2*(0,0.75,0.75).

If you inspect the file "bandsp", you will see that it contains three columns: (1) the kpoint path, (2) the real part of the quasiparticle energies, and (3) the third column contains the imaginary parts of the quasiparticle energies. Remember that these imaginary parts are proportional to the line width and inversely proportional to the excitation lifetime. The line widths can be plotted together with the band structure in gnuplot with

```
plot "bandsp" using 1:2:(6*abs($3)) with points ps var
```

You will see that the bands get wider the further they are away from the Fermi energy. (We will come back to this point in Section 1.2.10.)

##### 1.2.9.3 Wannier interpolation

As a third alternative, we can use Wannier interpolation. We have to modify "spex.inp" in the following way

```
BZ 4 4 4
JOB DFT0 GW IBZ:(1-4)
KPT X=1/2*(0,1,1) L=1/2*(0,0,1)
KPTPATH (L,1,X)
SECTION WANNIER
ORBITALS 1 4 [sp3]
MAXIMIZE
INTERPOL
END
```

The job definition lets SPEX calculate quasiparticle energies in the whole irreducible Brillouin zone (IBZ).
The section `WANNIER`

defines a set of maximally localized Wannier functions (MLWFs) of bonding sp^{3} hybrid orbitals.
These four hybrid orbitals are generated from the four lowest valence bands (first two arguments after `ORBITALS`

).
Running SPEX will construct the MLWFs and use them to interpolate the band structure.
The band structure data is written to the files "bands0" for the Kohn-Sham energies and "bands1" for the *GW* quasiparticle energies.
Please compare the interpolated band structure stored in "bands1" and the one that we have calculated
explicitly before (file "bandp"). (In Gnuplot: `plot "bandsp","bands1" with lines`

)
When you open the files with an editor, you will see that "bands1" has a column more than "bands0".
This additional column contains the interpolated imaginary parts of the quasiparticle energies.
You can plot the band structure together with the variable point size as described in the previous section.
Since we have defined four Wannier functions, we only get the four valence bands of silicon and no conduction bands.

We will need Wannier functions again when we calculate Hubbard *U* parameters in Section 2.

*Note*: If interested, you can now shortcut to other topics: (1.3) dielectric functions, (1.4) *GW* including spin-orbit coupling, (1.5) Hartree-Fock and *GW* for metals, (2) Hubbard *U* parameters.

#### 1.2.10 Spectral functions

The central quantity in the *GW* method is the renormalized Green function *G*, which is related to the self-energy Σ
through the Dyson equation *G*=*G*_{0}+*G*_{0}Σ*G*.
Every *GW* program calculates Σ and, then, solves the Dyson equation.
The last step is often simplified in such a way that an effective single-particle equation
(the so-called quasiparticle equation) is solved instead,
yielding the quasiparticle energies and lifetimes instead of the Green function *G*, as we have done in this tutorial
up to now. One can say that the quasiparticle energies and lifetimes *parameterize* *G*.
This works well for the quasiparticle peaks but not for many-body states such as plasmon satellites.
SPEX can also be used to calculate *G* directly. The imaginary part of its trace is called the *spectral function*.
To this end, we modify "spex.inp" in the following way:

```
SECTION SENERGY
CONTOUR {-2:1,0.01}
SPECTRAL
END
```

Here, we recommend to use the contour deformation method (`CONTOUR`

).
If you now run a *GW* calculation, you will get a file "spectral" containing spectral functions for each kpoint defined in the job definition. For example,
with

```
JOB GW 1:(1-8)
```

you would get the spectral function at the Γ point. Of course, you can choose any kpoint (also more than one kpoint).
Plot the spectral function.
The *x* axis is the energy axis (in eV).
One can imagine the spectral function as a vertical cut through the band structure at the respective kpoint.
If you compare with the band structure of silicon,
you will see that there is a (quasiparticle) peak in the spectral function for each band at the respective kpoint.
The peaks broaden the further they are away from the Fermi energy.
One can also observe a broad peak at around -30 eV. (You will have to reduce the y range, e.g., `plot [] [:20] "spectral" with lines`

to see the shallow peaks.) This is the plasmon satellite. You should see another
broad feature above the Fermi energy, at 20-30 eV, which is the plasmon satellite for the inverse photoemission process.

It is also possible to calculate the spectral functions along a path in the Brillouin zone,
to obtain a *renormalized* band structure, in a similar way as above.
In particular, we would follow the steps of Section 1.2.9.2 or Section 1.2.9.3.
(With the procedure of Section 1.2.9.1 one would not get sufficient k resolution.)

Using the latter approach (Section 1.2.9.3), SPEX produces a file named "spectralw", which contains a k-resolved spectral function (renormalized band structure). Do not forget to set `JOB GW IBZ:(1-4)`

and include the section `WANNIER`

. You can use "gnuplot" to plot the resulting renormalized band structure as a color-coded plot by

```
set cbrange [0:50]
plot [] [-45:5] "spectralw" palette with lines
```

Check if you can see the plasmon satellites in the plots! You might have to change the color range, e.g., to `set cbrange [0:5]`

The first approach (Section 1.2.9.2) is computationally more expensive (perhaps too expensive for the limited time of the hands-on session). In this case, an explicit frequency range has to be defined, e.g., `SPECTRAL {-1.5:0.25,0.001}`

.

The corresponding data is written to "spectralp" and can be plotted in the same way as "spectralw" before.

Of course, the k resolution is not very high. One can produce a finer kmesh by `KPTPATH (L,1,X) 50`

.
The additional parameter sets the number of kpoints along a unit length. The default is 20.

### 1.3 Gallium nitride

GaN crystallizes in the Wurtzite (ZnS) structure. Its hexagonal unit cell ("hP" in the input for "inpgen") has the lattice parameters a=3.19 Å and c=5.189 Å with the space group 186 (P63mc). The representative atoms are at (internal coordinates):

atom species | x | y | z | Wyckoff symbol |
---|---|---|---|---|

Ga | 1/3 | 2/3 | 0.0 | 2b |

N | 1/3 | 2/3 | 0.623 | 2b |

#### 1.3.1 *GW* band gap and 3d levels

Set up the input for the material (or use the prepared one "GaN/inp_GaN"), calculate the *GW* band gap.
It is at the Γ point. You also need to know the band index of the valence band maximum.
The number of electrons in the unit cell is given in the output (`Number of valence electrons: 36`

).
Dividing this number by 2 (for the two spin directions) yields the band index of the valence-band maximum.
So, the band gap is between band 18 and band 19. We are also interested in the 3d levels of Ga. These bands are the first
ten bands (each Ga atom contributes five). Just add the first band to the job definition.

GaN has twice as many atoms as silicon. So, we can reduce the BZ sampling:

```
BZ 3 3 3
```

Compare the *GW* band gap to the DFT one and to experiment (3.39 eV). What wavelength does a GaN laser diode have?

Compare the resulting 3d levels with the Kohn-Sham and HF values as well as with experiment (-17.1 eV). Note that the semicore level is measured with respect to the valence-band maximum (Band 18 at Γ).

#### 1.3.2 Dielectric function

It is also possible to calculate the dielectric function. For this, change the job definition to

```
JOB DIELEC 1:{0:30eV,0.1eV}
```

Here, we request a plot of the dielectric function between 0 and 5 eV with a step size of 0.01 eV. (Generally, energy values without a unit are interpreted as values in "Hartree". Alternatively, the values can be given in eV as shown.) This is done for kpoint label "1", meaning for a displacement vector of (0,0,0), corresponding to optical transitions. (The momentum of the photon is much smaller than that of the electrons. The long-wavelength limit k=0 is, thus, a good approximation). A SPEX run will produce the files "dielec" and "dielecR", corresponding to the bare and the renormalized dielectric function. The optical absorption spectrum (excluding excitonic effects) can be plotted as a function of frequency ω with ("gnuplot")

```
plot "dielecR" using 1:3 with lines
```

For an electron-energy loss spectrum (EELS), we have to plot the imaginary part of the reciprocal value:

```
plot "dielecR" using 1:(-imag(1/($2+{0,1}*$3))) with lines
```

The literature value of the dielectric constant (ω=0) is about 5.8. The refractive index *n* is the square root of the dielectric constant.
The literature value is *n*=2.41.

#### 1.3.3 Self-consistent QS*GW* calculation

The one-shot *GW* band gap calculated in section 1.3.1 still underestimates the experimental band gap of 3.39 eV. Perhaps a self-consistent *GW* calculation brings us closer to experiment. You can try such a QS*GW* calculation with the following modifications of the input file

```
JOB DFT GW FULL IBZ:(11-24) CYCLE 4
SECTION DFT
EXC NONLOCAL
END
RESTART 0
```

This will run four QS*GW* iterations. Compare the resulting QS*GW* band gap with the experimental value.

### 1.5 Sodium

If you still have time, you can run a *GW* calculation for metallic Na. We provide "inp_Na", an input file for "inpgen", in the directory "Na/".

Try to calculate a *GW* band structure from P over Γ to H as explained in Section 1.2.9.2. You can use a 6×6×6 kpoint mesh and the kpoint definition

```
KPT P=1/4*(1,1,1) H=1/2*(-1,1,1)
```

Sodium has four semicore states: 2s and 2p. So, the band indices in the job definition should start with 5. We intend to compare with a Hartree-Fock calculation. Please include the keyword `LOGDIV`

into the section `SENERGY`

(see explanation below). The band structure data is again written to the file "bandsp".

Compare this band structure to the corresponding Hartree-Fock band structure. To get a HF band structure, we do not have to run an extra HF calculation because the HF self-energy is contained in the *GW* self-energy, and every *GW* calculation also creates HF results. Do extract the HF energies from the output of the *GW* calculation, type

```
spex.extr h -b -O spex.out > bands_hf
```

(assuming that the output was piped into the file "spex.out").

While the band gradient of HF shows a logarithmic divergence at the Fermi energy,
there is no divergent band anomaly in the *GW* band structure because the correlation self-energy cancels the divergence of the HF term exactly by a term of opposite sign. The HF band anomaly in metals is unphysical. Therefore, SPEX calculates it explicitly only if the keyword `LOGDIV`

is specified in the input file.

Note that the band width is reduced in the *GW* band structure compared to the Kohn-Sham band structure. (Can be extracted with `k`

after `spex.extr`

.)
This reduction is caused by the increased effective mass of the quasiparticles with respect to the bare electrons.
As the electron moves, it has to drag its polarization cloud (Coulomb hole) along, which makes the quasiparticle effectively heavier.

## 2 Hubbard *U* parameters

### 2.1 Strontium vanadate (SrVO_{3})

As a first example, we will calculate the Hubbard *U* (partially screened Coulomb interaction) parameter for the V t_{2g} orbitals of
SrVO_{3} using the cRPA method. You find an input file ("inp_SrVO3") for "inpgen" in the directory "SrVO3/".

First, we calculate a DFT band structure with SPEX for the kpoint path R-Γ-X-M-Γ.
The special kpoints are at R=(0.5,0.5,0.5), Γ=(0,0,0), X=(0.5,0,0), and M=(0.5,0.5,0).
You have to specify the keywords `KPT`

and `KPTPATH`

accordingly as described in Sec. 1.2.11. To let SPEX calculate a DFT band structure, you specify `JOB BANDS`

. SPEX will stop afterwards. The file containing the band-structure data is called "bands". It usually contains so many bands that, in order to see separate bands, you need to decrease the plotted energy range:

`plot [] [2:20] "bands" with lines`

Before we calculate the Hubbard *U* parameter, we have to find out which correlated subspace to define. We expect it to be a vanadium *d* subspace, but will the subspace be composed of t_{2g} or e_{g} or all *d* orbitals?
To find this out, we calculate the band structure again, this time with orbital projections. This is achieved with
the keyword `PROJECT`

(section `ANALYZE`

). By default, we would get a decomposition into *s*, *p*, *d*, etc. orbitals of the first atom, then the second atom and so on. In the present case, we are only interested in the *d* states of vanadium. Therefore, we write

```
SECTION ANALYZE
PROJECT () (t2g,eg)
END
```

(the first atom is strontium, which is skipped). The next SPEX run will add orbital information to "bands", which can be plotted with

`set cbrange [0:1]`

`plot [] [2:20] "bands" u 1:2:3 palette with lines`

(for the t_{2g} states)

`plot [] [2:20] "bands" u 1:2:4 palette with lines`

(for the e_{g} states)

You will see that the t_{2g} bands form an isolated set of bands at the Fermi energy (which is at about 12.7 eV). Therefore, we use this
set of bands as the correlated subspace. The corresponding band indices are 21-23. A basis set for the correlated subspace will be formed by Wannier orbitals.

*Note*: You can get the Fermi energy from the "out" file of FLEUR (`grep fermi out`

) and also from SPEX (if you run with `JOB DFT0`

). Look for "Fermi energy"! The Fermi energies can be different due to the usage of different algorithms in the two codes.

Now, we calculate only the bands forming the t_{2g} subspace using
Wannier interpolation, as described in 1.2.11.3.
You can use the following input for "spex.inp".

```
BZ 4 4 4
JOB DFT0
KPT R=1/2*(1,1,1) G=(0,0,0) X=1/2*(1,0,0) M=1/2*(1,1,0)
KPTPATH (R,G,X,M,G)
SECTION WANNIER
ORBITALS 21 23 () (t2g)
INTERPOL
END
```

The Wannier-interpolated band structure is written to the file "bands0". (If you still have the keyword `PROJECT`

in your input file, you will see that "bands0" also contains Wannier-interpolated orbital projections. You can plot them in the same way as "bands", and you will see that the orbital projections are practically identical to the previous result.)

You can compare the two band structures in "gnuplot":

`plot "bands" with lines,"bands0" with lines`

.

You might be surprised that we do not use the keyword `MAXIMIZE`

here. Indeed, the maximal localization procedure is often unnecessary in the case of an isolated set of bands (presently the t_{2g} bands). So, Wannier90 is not employed here.

A minimal SPEX input file for Hubbard *U* calculations looks like this:

```
BZ 3 3 3
JOB DFT0 SCREENW {0}
SECTION SUSCEP
HUBBARD
END
SECTION WANNIER
ORBITALS 21 23 () (t2g)
END
```

We have reduced the Brillouin-zone sampling to 3×3×3 to speed up the calculation.
The keyword `SCREENW`

defines the job for calculating the partially screened Coulomb potential *U*. The expression in the brackets {...} can be
a general frequency range (in the same form as for the dielectric function, Section 1.3). Here, {0} means to calculate the Hubbard *U* parameter for
the static limit (ω=0). We only need the static *U* for model Hamiltonians or LDA+U calculations.

In the section `SUSCEP`

we provide the keyword `HUBBARD`

to exclude transitions in the correlated subspace, otherwise SPEX calculates the fully screened Coulomb interaction *W*.

Wannier functions are constructed internally in the SPEX code.
For this purpose, we have to provide the keyword `ORBITALS`

and specify the states from which the Wannier functions are to be constructed.
For SrVO_{3}, we construct three V t_{2g} orbitals from the states 21, 22, 23. The V atom is the second type as specified in the FLEUR input.
Therefore, we use an "empty" definition `()`

for the first atom (Sr).
In this calculation, we simply use the *first-guess* Wannier functions, thus the keyword `MAXIMIZE`

is excluded.

When SPEX is run, it will calculate the Coulomb matrix *U* in the mixed product basis
and then projects it onto the Wannier basis. You will get some effective parameters at the end of the program run,
Kanamori parameters for e_{g} and t_{2g} orbitals (and Hubbard-Hund parameters if full shells are specified).
The full interaction matrix is written to the file "spex.cou".

You can run the same calculation without the keyword `HUBBARD`

. Do the interaction parameters differ from the previous ones? If so, why?

Let us now also include the oxygen *p* orbitals in the Wannier set:

```
SECTION WANNIER
ORBITALS 12 23 () (t2g) (p)
END
```

The oxygen *p* bands are just below the vanadium t_{2g} bands. There are nine *p* bands (three atoms, three p bands per atom).

Calculate the Hubbard *U* parameters (do not forget the keyword `HUBBARD`

!) and compare the new *U* parameter for the t_{2g} states with the old one! Do they differ? What do you conclude from the difference?

Also try running a Wannier interpolation for the Wannier set including the oxygen *p* states and perhaps one including the e_{g} states in addition. Compare them with the standard DFT band structure ("bands").

### 2.2 Nickel

If you still have time, you can calculate the Hubbard *U* parameter for Nickel. To be more precise, we want to calculate the *U* parameter for
the 3d states of Ni.

We can mainly proceed as in the previous section. The Brillouin-zone sampling (`BZ`

) should be a bit finer:

```
BZ 6 6 6
```

And we need to define the Wannier orbitals (`ORBITALS`

) differently:

```
SECTION WANNIER
ORBITALS 4 9 (d)
END
```

Why do the band indices start with 4 in this case? (Have a look at the file "inp.xml"!)

As bulk nickel is ferromagnetic and Wannier functions have a formal spin dependence (since they are constructed from wavefunctions), we have to specify the spin indices for the Wannier products explicitly in the job definition. (Usually, the spin dependence is numerically small.)

```
JOB DFT0 SCREENW uu{0}
```

Since we have defined a full shell, the 3d shell, we get Hubbard-Hund parameters (instead of Kanamori parameters).

As the Ni *d* bands are entangled with the *s* band, the Wannier interpolation is a bit more complicated. We employ the Wannier90 library for the disentanglement procedure (`FROZEN`

) and for the maximal localization (`MAXIMIZE`

). The Wannier section could look like this:

```
SECTION WANNIER
ORBITALS 4 * (s,d)
MAXIMIZE
FROZEN
INTERPOL
END
```

Here, we have used a wildcard (`*`

) for the upper band index (of the Wannier window), which tells SPEX to choose a reasonable index. Also, the "frozen window" for disentanglement is chosen automatically. (It can also be set manually for fine tuning.) You can use the L-Γ-X path `KPTPATH (4,1,13)`

. (Here, instead of kpoint labels, we select the kpoint indices directly. The numbers are for 6×6×6). Run a Wannier interpolation and compare with the explicit DFT band structure.