# Tutorial 2: the band structure of bulk silicon (calculated via a supercell)

In this tutorial, we will calculate the KI bandstructure of bulk silicon. The input file used for this calculation can be downloaded `here`

.

## Wannierization

While this calculation will bear a lot of similarity to the previous tutorial, there are several differences between performing a Koopmans calculation on a bulk system vs. a molecule. One major difference is that we use Wannier functions as our variational orbitals.

### What are Wannier functions?

Most electronic-structure codes try to calculate the Bloch states \(\psi_{n\mathbf{k}}\) of periodic systems (where \(n\) is the band index and \(\mathbf{k}\) the crystal momentum). However, other representations are equally valid. One such representation is the Wannier function (WF) basis. In contrast to the very delocalised Bloch states, WFs are spatially localised and as such represent a very convenient basis to work with. In our case, the fact that they are localised means that they are suitable for use as variational orbitals.

Wannier functions \(w_{n\mathbf{R}}(\mathbf{r})\) can be written in terms of a transformation of the Bloch states:

where our Wannier functions belong to a particular lattice site \(\mathbf{R}\), \(V\) is the unit cell volume, the integral is over the Brillouin zone (BZ), and \(U^{(\mathbf{k})}_{mn}\) defines a unitary rotation that mixes the Bloch states with crystal momentum \(\mathbf{k}\). Crucially, this matrix \(U^{(\mathbf{k})}_{mn}\) is not uniquely defined — indeed, it represents a freedom of the transformation that we can exploit.

We choose our \(U^{(\mathbf{k})}_{mn}\) that gives rise to WFs that are “maximially localised”. We quantify the spread \(\Omega\) of a WF as

The Wannier functions that minimise this spread are called maximally localised Wannier functions (MLWFs). For more details, see Ref. [20].

### How do I calculate Wannier functions?

MLWFs can be calculated with Wannier90, an open-source code that is distributed with `Quantum ESPRESSO`

. Performing a Wannierization with Wannier90 requires a series of calculations to be performed with `pw.x`

, `wannier90.x`

, and `pw2wannier90.x`

. This workflow is automated within `koopmans`

, as we will see in this tutorial.

Note

This tutorial will not discuss in detail how perform a Wannierization with Wannier90. The Wannier90 website contains lots of excellent tutorials.

One important distinction to make for Koopmans calculations — as opposed to many of the Wannier90 tutorials — is that we need to separately Wannierize the occupied and empty manifolds.

Let’s now inspect this tutorial’s `input file`

. At the top you will see that

```
2 "workflow": {
3 "task": "wannierize",
4 "functional": "ki",
```

which tells the code to perform a standalone Wannierization calculation. Meanwhile, near the bottom of the file there are some Wannier90-specific parameters provided in the `w90`

block

```
38 "w90": {
39 "bands_plot": true,
40 "projections": [[{"fsite": [ 0.25, 0.25, 0.25 ], "ang_mtm": "sp3"}],
41 [{"fsite": [ 0.25, 0.25, 0.25 ], "ang_mtm": "sp3"}]],
42 "dis_froz_max": 10.6,
43 "dis_win_max": 16.9
44 },
```

The `w90`

block format is explained more fully here.

We run this calculation as per usual:

```
koopmans si.json | tee si_wannierize.out
```

After the usual header, you should see something like the following:

```
Wannierization
==============
Running wannier/scf... done
Running wannier/nscf... done
Running wannier/block_1/wann_preproc... done
Running wannier/block_1/pw2wan... done
Running wannier/block_1/wann... done
Running wannier/block_2/wann_preproc... done
Running wannier/block_2/pw2wan... done
Running wannier/block_2/wann... done
Running wannier/bands... done
Running pdos/projwfc... done
Workflow complete
```

These various calculations that are required to obtain the MLWFs of bulk silicon. You can inspect the and various output files will have been generated in a new `wannier/`

directory.

## Click here for detailed descriptions of each calculation

- scf
a

`pw.x`

self-consistent DFT calculation performed with no empty bands. This obtains the ground-state electronic density- nscf
a

`pw.x`

non-self-consistent DFT calculation that determines the Hamiltonian, now including some empty bands- block_1/wann_preproc
a preprocessing

`wannier90.x`

calculation that generates some files required by`pw2wannier90.x`

- block_1/pw2wan
a

`pw2wannier90.x`

calculation that extracts from the eariler`pw.x`

calculations several key quantities required for generating the Wannier orbitals for the occupied manifold: namely, the overlap matrix of the cell-periodic part of the Block states (this is the`wann.mmn`

file) and the projection of the Bloch states onto some trial localised orbitals (`wann.amn`

)- block_1/wann
the

`wannier90.x`

calculation that obtains the MLWFs for the occupied manifold- block_2/…
the analogous calculations as those in

`occ/`

, but for the empty manifold- bands
a

`pw.x`

calculation that calculates the band structure of silicon explicitly, used for verification of the Wannierization (see the next section)- pdos/projwfc
a

`projwfc.x`

calculation that calculates the projected density of states, also used for checking the Wannierization

The main output files of interest in `wannier/`

are files `block_1/wann.wout`

and `block_2/wann.wout`

, which contain the output of `wannier90.x`

for the Wannierization of the occupied and empty manifolds. If you inspect either of these files you will be able to see a lot of starting information, and then under a heading like

```
*------------------------------- WANNIERISE ---------------------------------*
+--------------------------------------------------------------------+<-- CONV
| Iter Delta Spread RMS Gradient Spread (Ang^2) Time |<-- CONV
+--------------------------------------------------------------------+<-- CONV
```

you will then see a series of steps where you can see the Wannier functions being optimised and the spread (labelled `SPRD`

) decreasing from one step to the next. Scrolling down further you should see a statement that the Wannierization procedure has converged, alongside with a summary of the final state of the WFs.

### How do I know if the Wannier functions I have calculated are “good”?

Performing a Wannierization calculation is not a straightforward procedure, and requires the tweaking of the Wannier90 input parameters in order to obtain a “good” set of Wannier functions.

One check is to see if an interpolated bandstructure generated by the MLWFs resembles an explicitly-calculated band structure. (For an explanation of how one can use Wannier functions to interpolate band structures, refer to Ref. [20].) You might have noticed that when we ran `koopmans`

earlier it also generated a file called `si_wannierize_bandstructure.png`

. It should look something like this:

Clearly, we can see that the interpolation is no good! (The interpolated band structure ought to lie on top of the explicitly calculated band structure.) The reason for this is that the Brillouin zone is undersampled by our \(2\times2\times2\) \(k\)-point grid. Try increasing the size of the k-point grid and see if the interpolated bandstructure improves.

Tip

Trying different grid sizes can be very easily automated within `python`

. Here is a simple script that will run the Wannierization for three different grid sizes:

```
from koopmans.io import read
from koopmans.utils import chdir
for grid_size in [2, 4, 8]:
# Read in the input file
wf = read('si.json')
# Modify the kgrid
wf.kpoints.grid = [grid_size, grid_size, grid_size]
# Run the workflow in a subdirectory
with chdir('{0}x{0}x{0}'.format(grid_size)):
wf.run()
```

Tip

You may also notice a file called `si_wannierize_bandstructure.fig.pkl`

. This is a version of the figure that you can load in python and modify as you see fit. e.g. here is a script that changes the y-axis limits and label:

```
import pickle
import matplotlib.pyplot as plt
with open('2x2x2/si_wannierize_bandstructure.fig.pkl', 'rb') as fd:
pickle.load(fd)
ax = plt.gca()
ax.set_ylim([-5, 5])
ax.set_ylabel(r'$\omega$ (eV)')
# plt.show() uncomment to view the figure interactively
plt.savefig('2x2x2/si_wannierize_bandstructure_rescaled.png')
```

## The KI calculation

Having obtained a Wannierization of silicon that we are happy with, we can proceed with the KI calculation. In order to do this simply change the `task`

keyword in the input file from `wannierize`

to `singlepoint`

.

Tip

Although we just discovered that a \(2\times2\times2\) \(k\)-point grid was inadequate for producing good Wannier functions, this next calculation is a lot more computationally intensive and will take a long time on most desktop computers. We therefore suggest that for the purposes of this tutorial you switch back to the small \(k\)-point grid. (But for any proper calculations, always use high-quality Wannier functions!)

### Initialization

If you run this new input the output will be remarkably similar to that from the previous tutorial, with a couple of exceptions. At the start of the workflow you will see there is a Wannierization procedure, much like we had earlier when we running with the `wannierize`

task:

```
20 Wannierization
21 ==============
22 Running wannier/scf... done
23 Running wannier/nscf... done
24 Running wannier/block_1/wann_preproc... done
25 Running wannier/block_1/pw2wan... done
26 Running wannier/block_1/wann... done
27 Running wannier/block_2/wann_preproc... done
28 Running wannier/block_2/pw2wan... done
29 Running wannier/block_2/wann... done
```

which replaces the previous series of semi-local and PZ calculations that we used to initialize the variational orbitals for a molecule.

There is then an new “folding to supercell” subsection:

```
32 Folding to supercell
33 --------------------
34 Running block_1/w2kcp... done
35 Running block_2/w2kcp... done
```

In order to understand what these calculations are doing, we must think ahead to the next step in our calculation, where we will calculate the screening parameters using the ΔSCF method. These calculations, where we remove/add an electron from/to the system, require us to work in a supercell. This means that we must transform the \(k\)-dependent primitive cell results from previous calculations into equivalent \(\Gamma\)-only supercell quantities that can be read by `kcp`

. This is precisely what the above `wan2odd`

calculations do.

### Calculating the screening parameters

Having transformed into a supercell, the calculation of the screening parameters proceeds as usual. The one difference to tutorial 1 that you might notice at this step is that we are skipping the calculation of screening parameters for some of the orbitals:

```
42 Calculating screening parameters
43 ================================
44 Running calc_alpha/ki... done
45
46 Orbitals 1-31
47 -------------
48 Skipping; will use the screening parameter of an equivalent orbital
49
50 Orbital 32
51 ----------
```

The code is doing this because of what we provided for the `orbital_groups`

in the input file:

```
9 "alpha_guess": 0.077,
10 "orbital_groups": [0, 0, 0, 0, 1, 1, 1, 1],
11 "pseudo_library": "pseudo_dojo_standard",
```

which tells the code to use the same parameter for orbitals belonging to the same group. In this instance we are calculating a single screening parameter for all four filled orbitals, and a single screening parameter for the empty orbitals.

### The final calculation and postprocessing

The final difference for the solids calculation is that there is an additional preprocessing step at the very end:

```
89 Postprocessing
90 ===============
91
92 Wannierization
93 ==============
94 Running wannier/scf... done
95 Running wannier/nscf... done
96 Running wannier/block_1/wann_preproc... done
97 Running wannier/block_1/pw2wan... done
98 Running wannier/block_1/wann... done
99 Running wannier/block_2/wann_preproc... done
100 Running wannier/block_2/pw2wan... done
101 Running wannier/block_2/wann... done
102 Running wannier/bands... done
103 Running pdos/projwfc... done
104 Running occ/ki... done
105 Running emp/ki... done
106
107 Workflow complete
```

Here, we transform back our results from the supercell sampled at \(\Gamma\) to the primitive cell with \(k\)-space sampling. This allows us to obtain a bandstructure. The extra Wannierization step that is being performed is to assist the interpolation of the band structure in the primitive cell, and has been performed because in the input file we specified

```
45 "ui": {
46 "smooth_int_factor": 4
47 }
```

For more details on the “unfold and interpolate” procedure see here and Ref. [11].

## Extracting the KI bandstructure and the bandgap of Si

The bandstructure can be found in `postproc/bands_interpolated.dat`

as a raw data file, but there is a more flexible way for plotting the final bandstructure using the python machinery of `koopmans`

:

```
from koopmans.io import read
# Load the workflow object
wf = read('si.kwf')
# Access the band structure from the very last calculation
results = wf.calculations[-1].results
bs = results['band structure']
# Print the band strucutre to file
bs.plot(filename='si_bandstructure.png')
# Extract the band gap
n_occ = wf.number_of_electrons() // 2
gap = bs.energies[:, :, n_occ:].min() - bs.energies[:, :, :n_occ].max()
print(f'KI band gap = {gap:.2f} eV')
```

Running this script will generate a plot of the bandstructure (`si_bandstructure.png`

) as well as printing out the band gap. You should get a result around 1.35 eV. Compare this to the PBE result of 0.68 eV and the experimental value of 1.22 eV. If we were more careful with the Wannier function generation, our result would be even closer (indeed in Ref. [24] the KI band gap was found to be 1.22 eV!)