Complex bandstructure of Si(100)¶
Version: 2015.1
In this tutorial you will calculate the complex bandstructure of a silicon crystal along the (100) direction.
In particular you will:
create the Si(100) surface;
set up and run the calculation;
plot the complex bandstructure (3D and 2D).
Background¶
In a periodic solid the eigenstates \(\psi_{n{\bf k}}\) of the Schrödinger equation \(H \psi_{n{\bf k}} = E_{n{\bf k}} S \psi_{n{\bf k}}\), (\(S\) is the overlap matrix) can be written as \(\psi_{n{\bf k}}({\bf r}) = e^{-i {\bf k}\cdot {\bf r}} U_{n{\bf k}}({\bf r})\), where \(U_{n{\bf k}}({\bf r})\) is a periodic function with the same periodicity as the crystal itself. In a usual bandstructure calculation, the wave vector \({\bf k}\) is a real number, and by solving the Schrödinger equation above for various fixed values of \({\bf k}\) (typically located along different symmetry lines in the first Brillouin zone) an eigenproblem is defined, from which the eigenenergies \(E_{n{\bf k}}\) (i.e. the bandstructure) can be determined.
When calculating the complex bandstructure another approach is taken [1]. Instead the energy \(E\) is fixed, and the values of \({\bf k}\) which solve the Schrödinger equation are sought. Such solutions will be found with both real and complex \({\bf k}\); the solutions with a real \({\bf k}\) are the usual Bloch states, while the solutions with an imaginary part are exponentially decaying in one direction and increasing in the other. Such solutions cannot exist in a bulk material, and so they are normally ignored in a bandstructure calculation. They may however exist at a surface or interface, and they give information about how electronic states decay in the material, for instance through a thin tunneling barrier.
Note
You will primarily use the graphical user interface QuantumATK for setting up and analyzing the results. If you are new to QuantumATK, it is recommended to go through the Basic QuantumATK Tutorial.
The calculations in this tutorial will be performed using the semi-empirical
models of QuantumATK. A complete description of all the
parameters, and in many cases a longer discussion about their physical
relevance, can be found in the ATK Reference Manual. In particular,
the Reference Manual entry for complex bandstructure calculations is
of relevance: ComplexBandstructure
.
Si(100) surface¶
Start QuantumATK and create a new project named “ComplexBandstructure”. Then click Open and launch the Builder via the icon on the toolbar.
In the Builder, click icon in the lower right-hand corner.
and locate the “Silicon (alpha)” structure. Double-click the line to add the structure to the Stash, or click theUnfold Builders from the plugin panel and open the Surface (Cleave) tool, and follow the next steps to set up the Si(100) surface:
Keep the default Miller indices (100) and click Next.
Keep the default lattice vectors and click Next.
Select Periodic (bulk-like) for the out-of-plane cell vector.
Set the thickness of the surface to 1 layer.
Click Finish to add the Si(100) surface structure to the Stash.
With these settings we use a surface representation with a minimum number of layers in order to reduce the calculation time and avoid zone folding.
Note
The cleave plane – in the present case (100) – is always spanned by the two first unit cell vectors, \({\bf}\) and \({\bf B}\). In the “electrode” mode, the normal to this plane coincides with \({\bf C}\), the third unit cell vector, but in the “bulk-like” mode this is not the case. In QuantumATK the complex bandstructure is always computed along the third reciprocal vector, \({\bf g}_C\), which is of course parallel to \({\bf A} \times {\bf B}\). With the present geometry you will therefore obtain the complex bandstructure along (100).
Complex bandstructure calculation¶
Send the structure to the Script Generator using the icon in the lower right-hand corner of the Builder and follow the next steps:
Next, open the inserted New Calculator block (double-click it), and set up a tight-binding calculation:
Select the ATK-SE: Extended Hückel calculator.
Set the k-point sampling to 5x5x5.
Under Huckel Basis set, select the Cerda.Silicon [diamond GW] basis set. This basis set has been fitted to GW calculations, and gives an excellent description of the bandstructure, including the size of the band gap.
Close the dialogue by clicking OK.
Next, open the ComplexBandstructure analysis block, and edit it:
Set the energy range to -15 to 10 eV with 501 points.
Note
The projection of \({\bf k}\) on the cleave plane is kept constant in the complex andstructure calculation, and the value of the projection is specified by the \((k_A,k_B)\) parameters, which are given in units of the two first reciprocal lattice vectors, \({\bf g}_A\) and \({\bf g}_B\). Therefore, the obtained solutions will lie along the third reciprocal lattice vector, \({\bf g}_C\), which is parallel to the cleave plane normal. A new set of solutions \(k_C+i\kappa_C\) is obtained for each value of \((k_A,k_B)\).
Please note that \(k_C\) is therefore the real part of the complex bandstructure, while \(kappa_C\) is the complex part.
You have now finished the Python script. Save it as si_100_cbs.py
and send it to the Job Manager for execution.
It should only take a few minutes for this small system. If needed,
you can also download the script here: si_100_cbs.py
.
Tip
This type of calculation parallelizes extremely well, so for larger structures it is warmly recommended to execute the script in parallel.
Analysing the results¶
The file si_100_cbs.hdf5
should now have appeared on the QuantumATK LabFloor.
It contains the saved Hückel calculation and the ComplexBandstructure
analysis object:
Select that analysis object and click Show 2D Plot tool from the right-hand side plugins panel. A window with a plot of the complex bandstructure pops up. Zoom in a bit to filter out the solutions with large values of the complex part \(\kappa_C\), since these are not really relevant:
Note
In the plot above, the solutions \(\kappa_C\) are given in
reciprocal Cartesian coordinates,
to make it easier to compare different structures. For the real part
of the bandstructure, on the other hand, the solutions \(k_C\)
are normalized by \(L\), the layer separation perpendicular to
the cleave plane. In the case of a fcc crystal cleaved along (100),
the layer separation is \(L=a/2\), where \(a\) is the lattice
constant. The layer separation can be extracted from the
ComplexBandstructure object using the layerSeparation() method,
see ComplexBandstructure
.
3D and 2D visualizations¶
In the complex part of the bandstructure, the solutions actually have
a real part too. Thus, they could be visualized in a three-dimensional
plot, \((x,y,z)=(k_C, \kappa_C, E)\). This is possible with the
following script. Open the QuantumATK Editor and copy-paste
the script and save it as 3D_plot.py
. Make sure that the indentation
is correct. Alternatively, you can simply download it here: 3D_plot.py
.
1from QuantumATK import *
2from mpl_toolkits.mplot3d import Axes3D
3import matplotlib.pyplot as plt
4import math
5
6# Read the complex bandstructure object from the NC file
7cbs = nlread('si_100_cbs.nc', ComplexBandstructure)[-1]
8energies = cbs.energies().inUnitsOf(eV)
9k_real, k_complex = cbs.evaluate()
10L = cbs.layerSeparation()
11
12fig = plt.figure()
13ax = fig.add_subplot(111, projection='3d')
14
15# First plot the real bands
16kvr = numpy.array([])
17e = numpy.array([])
18
19# Loop over energies, and pick those where we have solutions
20for (j, energy) in enumerate(energies):
21 k = k_real[j]*L/math.pi
22 if len(k)>0:
23 e = numpy.append(e,[energy,]*len(k))
24 kvr = numpy.append(kvr,k)
25
26# Plot the bands with red
27ax.scatter([0.0,]*len(kvr), kvr, e, c='r', marker='o', linewidths=0, s=10)
28
29# Next plot the complex bands
30kvr = []
31kvi = []
32e = []
33
34# Again loop over energies and pick solutions
35for (j, energy) in enumerate(energies):
36 if len(k_complex[j])>0:
37 for x in numpy.array(k_complex[j]*L/math.pi):
38 kr = numpy.abs(x.real)
39 ki = -numpy.abs(x.imag)
40 # Discard rapidly decaying modes which clutter the plot
41 if ki>-0.3:
42 e += [energy]
43 kvr += [kr]
44 kvi += [ki]
45
46# Plot the complex bands with blue
47ax.scatter(kvi, kvr, e, c='b', marker='o', linewidths=0, s=10)
48
49# Put on labels
50ax.set_xlabel('$\kappa$ (1/Ang)')
51ax.set_ylabel('$kL/\pi$')
52ax.set_zlabel('Energy / eV')
53
54plt.show()
Execute the script using the Job Manager or from a Terminal. The plot will resemble the figure below, but your points will be more scattered, since the figure below was produced from a Hückel calculation with 10,001 energy points instead of 501.
Another way to visualize the real part of the complex bandstructure
is using colors. The script 2D_plot.py
does this (you can see the script
below). It is a bit complicated towards the end, but that part is just
for placing the colorbar, and could be omitted.
The “forest” of complex bands with a rather large value of \(\kappa\) are not usually seen in tight-binding simulations. However, for DFT and Hückel, the basis sets are larger, so there are more complex bands as well, connecting unoccupied levels with various valence bands.
1from QuantumATK import *
2import matplotlib.pyplot as plt
3import math
4
5# Read the complex bandstructure object from the NC file
6cbs = nlread('si_100_cbs.nc', ComplexBandstructure)[-1]
7energies = cbs.energies().inUnitsOf(eV)
8k_real, k_complex = cbs.evaluate()
9
10ax = plt.axes()
11cmap="Spectral"
12
13# First plot the real bands
14kvr = numpy.array([])
15e = numpy.array([])
16for (j, energy) in enumerate(energies):
17 k = k_real[j]*cbs.layerSeparation()/math.pi
18 if len(k)>0:
19 e = numpy.append(e,[energy,]*len(k))
20 kvr = numpy.append(kvr,k)
21
22# Plot
23ax.scatter(kvr, e,
24 c=numpy.abs(kvr),
25 cmap=cmap, marker='o', linewidths=0, s=10)
26
27# Next plot the complex bands
28kvr = numpy.array([])
29kvi = numpy.array([])
30e = numpy.array([])
31
32for (j, energy) in enumerate(energies):
33 if len(k_complex[j])>0:
34 kr = [numpy.abs(x.real) for x in k_complex[j]]
35 ki = [numpy.abs(x.imag) for x in k_complex[j]]
36 e = numpy.append(e,[energy,]*len(kr))
37 kvr = numpy.append(kvr,kr)
38 kvi = numpy.append(kvi,ki)
39
40# Plot with color depending on the imaginary part (corresponding to real k-points)
41sc = ax.scatter(-kvi, e,
42 c=kvr,
43 cmap=cmap, marker='o', linewidths=0, s=10)
44
45# Put on labels and decorations
46ax.axvline(0,color='b')
47ax.grid(True, which='major')
48ax.set_xlim(-1, 1)
49ax.set_ylim(-15, 10)
50ax.annotate('$\kappa$ (1/Ang)', xy=(0.25,-0.07), xycoords="axes fraction", ha="center")
51ax.annotate('$kL / \pi$', xy=(0.75,-0.07), xycoords="axes fraction", ha="center")
52ax.set_ylabel('Energy / eV')
53
54# Add a colorbar
55fig = plt.gcf()
56x1, x2, y1, y2 = 0., 1, ax.get_ylim()[0], ax.get_ylim()[0]+1
57trans = ax.transData + fig.transFigure.inverted()
58ax_x1, ax_y1 = trans.transform_point([x1, y1])
59ax_x2, ax_y2 = trans.transform_point([x2, y2])
60ax_dx, ax_dy = ax_x2 - ax_x1, ax_y2 - ax_y1
61cmap_axes = plt.axes([ax_x1, ax_y1, ax_dx, ax_dy])
62a = numpy.outer(numpy.arange(0,1,0.01),numpy.ones(10)).transpose()
63cmap_plt = plt.imshow(a,aspect='auto',cmap=plt.get_cmap(cmap),origin=(0,0))
64ax = plt.gca()
65ax.set_xticks([])
66ax.set_yticks([])
67ax.set_xticklabels([])
68ax.set_yticklabels([])
69
70plt.show()
Hint
You may note that there are “gaps” in the visualizations. The reason is that, unlike a normal bandstructure plot, the data is not plotted with lines, but with dots. In a standard bandstructure plot you can – with some level of confidence – define “bands”, which continuously run between the symmetry points (only band crossings can create some small problems). In the present case, the number of solutions is first of all different at each energy (particularly on the complex side), and depending on the density of the energy sampling, you may not hit a particular band close to points where the bands are very flat (heavy effective mass). This can to some extent be alleviated by increasing the number of energy points.
References¶