where the periodic function $u_{\vec{k}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$ has been expressed as a Fourier series. The reciprocal lattice vectors have been relabeled as running over $\vec{G}'$ instead of $\vec{G}$. It does not matter how they are labeled since we sum over all of the reciprocal lattice vectors, $\sum\limits_{\vec{G}}C_{\vec{G}}e^{i\vec{G}\cdot\vec{r}}=\sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i\vec{G}'\cdot\vec{r}}$. This form for the wavefunction is substituted into the Schrödinger equation,
$$\sum\limits_{\vec{G}'}\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m} C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}} + \sum_{\vec{G}}\sum\limits_{\vec{G}''} U_{\vec{G}} C_{\vec{G}''}e^{i(\vec{k}+\vec{G}+\vec{G}'')\cdot\vec{r}} = E \sum\limits_{\vec{G}'}C_{\vec{G}'}e^{i(\vec{k}+\vec{G}')\cdot\vec{r}}. $$
In the middle term with the double sum, the sum over $\vec{G}'$ has been relabeled as a sum over $\vec{G}''$. Again it does not matter that the label has changed since the sum is over all of the states but we need a way to keep track of the product terms in the double sum. Next we collect like terms. The exponential factors can be written as $e^{i \vec{k} \cdot \vec{r}}= \cos(\vec{k} \cdot \vec{r})+ i\sin(\vec{k} \cdot \vec{r})$. Only terms that have the same wavelength can be equal to each other so only the terms where $\vec{k} +\vec{G}'= \vec{k} + \vec{G}+\vec{G}''$ can be equal to each other. This results in the condition $\vec{G}'' = \vec{G}' - \vec{G}$ and for each $\vec{G}'$ there is an equation,
$$\dfrac{\hbar^2(\vec{k}+\vec{G}')^2}{2m}C_{\vec{G}'} + \sum_{\vec{G}} U_{\vec{G}} C_{\vec{G}'-\vec{G}} = E C_{\vec{G}'}. $$
This set of equations are called the central equations. The Schrödinger equation, a differential equation for $\psi_{\vec{k}}$, has been replaced with the central equations which are algebraic equations for the coefficients $C_{\vec{G}'}$. The algebraic equations can be put in the form of an eigenvalue problem. These equations written out in the one-dimensional case for $-G_0, 0, G_0$ are,
If the terms indicated by $\cdots$ are negelected, this can be written in matrix form as an eigenvalue problem that can be solved for the energies.
$$\textbf{M} \vec{C} = E \vec{C}.$$
When more equations for $\vec{G}'$ are included, the matrix gets bigger and since the coefficients $U_{\vec{G}}$ typically get smaller as $\vec{G}$ gets larger, it becomes more justified to neglect the terms indicated by $\cdots$. An $N\times N$ matrix will have $N$ solutions for the energy $E$ at each value of $\vec{k}$. The different solutions correspond to different bands.
In three dimensions the reciprocal lattice vectors are typically indexed by three integers $h$, $k$, and $l$. However to construct the matrix, it is necessary to index the reciprocal lattice vectors by a single integer. This can be done by making a list of $[hkl]$ triples and assigning each triple to an integer $i$. The python code below defines $h$, $k$, and $l$ vectors and assigns them to an integer $i$ like this:
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
The periodic potential could be a periodic Coulomb potential, a muffin tin potential, or a pseudopotential. To calculate the dispersion relation for a simple cube, bcc, or fcc lattice with a periodic potential not listed below, use the code for the Coulomb potential and it will only be necessary to change the line that specifies the off-diagonal elements of the matrix $\textbf{M}$ for that potential.
Comparisons with DFT calculations: Li, Na, Mg, Al, Cu.
1-D zero potential
1-D comb potential
If all of the components of a Fourier series are real and have the same amplitude, there is constructive interference at the Bravais lattice points and this produces a comb function with one high peak per Bravais lattice point. (see: Fourier synthesis).
The individual Fourier components.
The interference generates the negative peaks.
Below is some python code that calculates the band structure of a one-dimensional comb potential. The results of the numerical 1-D band structure solver are shown for comparison.
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
import math
pi = math.pi
hbar = 1.054571817E-34
me = 9.1093837015E-31
e = 1.602176634E-19
a = 3E-10
alpha = hbar*hbar/(2*me*a*a*e)
kvec = np.empty(101) # wavevector
energy = np.empty((9,101))
M = -1*np.ones((9,9)) # all of the off-diagonal elements are the same (-1)
#calculate the diagonal matrix elements for every k
for j in range(101):
k = pi*j/100
for i in range(9):
G = 2*pi*(i-4)
M[i,i] = alpha*(k -G)*(k-G)
kvec[j] = k
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,9):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
plt.ylim([-1,30]) #plot only the low energies, the high energies are unreliable
plt.show()
1-D square-wave potential
It is possible to find analytic solutions to the problem of an electron in one dimension moving in a square wave potential. This is called the Kronig-Penney model. The plane wave method can be used to approximate a square wave using the first few Fourier components. Below is some python code that calculates the band structure of a one-dimensional square-wave potential. The results of the Kronig-Penney model are shown for comparison.
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
import math
pi = math.pi
hbar = 1.054571817E-34
me = 9.1093837015E-31
e = 1.602176634E-19
a = 3E-10
alpha = hbar*hbar/(2*me*a*a*e)
kvec = np.empty(101) # wavevector
energy = np.empty((9,101))
M = np.zeros((9,9))
# calculate the off-diagonal matrix elements
# the formula for a square wave is given at
# http://lampz.tugraz.at/~hadley/ss1/crystaldiffraction/fourier/fourier1.php
for i in range(9):
for j in range(i+1,9):
n = j-i
M[i,j]=4*math.sin(pi*n/2)/(pi*n)
M[j,i] = M[i,j]
print(M)
#calculate the diagonal matrix elements for every k
for j in range(101):
k = pi*j/100
for i in range(9):
G = 2*pi*(i-4)
M[i,i] = alpha*(k -G)*(k-G)
kvec[j] = k
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,9):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
plt.ylim([-1,30]) #plot only the low energies, the high energies are unreliable
plt.show()
1-D cosine potential
The plane wave method was used to calculate the band structure of an electron moving in a cosine potential. The results of the numerical 1-D band structure solver are shown for comparison.
import numpy as np
import scipy.linalg as la
import matplotlib.pyplot as plt
import math
pi = math.pi
hbar = 1.054571817E-34
me = 9.1093837015E-31
e = 1.602176634E-19
a = 3E-10
alpha = hbar*hbar/(2*me*a*a*e)
kvec = np.empty(101) # wavevector
energy = np.empty((9,101))
# only the subdiagonal and superdiagonal elements are non zero
# 2 cos(x) = exp(ix) + exp(-ix)
M = np.zeros((9,9))
for i in range(1,9):
M[i,i-1]= 1
M[i-1,i] = M[i,i-1]
print(M)
#calculate the diagonal matrix elements for every k
for j in range(101):
k = pi*j/100
for i in range(9):
G = 2*pi*(i-4)
M[i,i] = alpha*(k -G)*(k-G)
kvec[j] = k
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,9):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
plt.ylim([-1,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Body Centered Cubic with a zero potential
A good way to check if your band structure code is working is to calculate the band structure for a zero potential. This should reproduce the results of the empty lattice approximation. Below is some python code that uses the plane wave method to calculate the band structure of a three-dimensional bcc lattice where the periodic potential has zero amplitude.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
#The off-diagonal elements are zero for a zero potential
M = np.zeros((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
#Gamma 0, 0, 0
#H 0, 0, 2*pi
#P pi, pi, pi
#N 0, pi, pi
#Gamma - H
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#H - N
k1x = 0 #initial k point
k1y = 0
k1z = 2*pi
k2x = 0 #final k point
k2y = pi
k2z = pi
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#N - Gamma
k1x = 0 #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma -P
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#P - H
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','H','N','$\Gamma$','P','H']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([0,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Body Centered Cubic with a Comb Potential
If all of the components of a Fourier series are real and have the same amplitude, there is constructive interference at the Bravais lattice points and this produces a comb function with one large amplitude peak per Bravais lattice point. Below is python code to calculate the band structure of a bcc comb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
# Define the off-diagonal elements of the matrix M
# For a Comb potential, all of the off-diagonal matrix elements are the same
for i in range(125):
for j in range(125):
M[i,j] = -5; #Comb potential
#Gamma 0, 0, 0
#H 0, 0, 2*pi
#P pi, pi, pi
#N 0, pi, pi
#Gamma - H
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#H - N
k1x = 0 #initial k point
k1y = 0
k1z = 2*pi
k2x = 0 #final k point
k2y = pi
k2z = pi
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#N - Gamma
k1x = 0 #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma -P
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#P - H
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','H','N','$\Gamma$','P','H']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([5,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Body Centered Cubic with a Coulomb Potential
Consider a periodic function consisting of a Coulomb potential at every Bravais lattice point. This potential is,
Here $Z$ is the effective nuclear charge, $\vec{r}_j$ is the position of the Bravais lattice point, and $\epsilon_0 = 8.854187817\times 10^{-12}\text{ F/m}$ is the permittivity constant. The Fourier series for this potential is,
where $V_{u.c.}$ is the volume of the primitive unit cell. For bcc with a conventional lattice constant $a$, there are two primitive unit cells in a conventional unit cell so $V_{u.c.}= \frac{a^3}{2}$
Below is python code to calculate the band structure of a bcc Coulomb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
# Define the off-diagonal elements of the matrix M
# The Fourier series for a Columb potential is example 4 at
# http://lampz.tugraz.at/~hadley/ss1/crystaldiffraction/fourier/fourier_examples.php
for i in range(125):
Gxi =2*pi*(h[i]+l[i])
Gyi =2*pi*(h[i]+k[i])
Gzi =2*pi*(k[i]+l[i])
for j in range(i+1,125):
Gxj =2*pi*(h[j]+l[j])
Gyj =2*pi*(h[j]+k[j])
Gzj =2*pi*(k[j]+l[j])
Gx = Gxj - Gxi
Gy = Gyj - Gyi
Gz = Gzj - Gzi
G2 = Gx*Gx + Gy*Gy + Gz*Gz #(Gi -Gj)^2
#G is in units of 1/a
# M is in units of eV, Z = 1, Vuc = 0.5*a^2
M[i,j]=-2*e/(a*eps0*G2)
M[j,i] = M[i,j]
#Gamma 0, 0, 0
#H 0, 0, 2*pi
#P pi, pi, pi
#N 0, pi, pi
#Gamma - H
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#H - N
k1x = 0 #initial k point
k1y = 0
k1z = 2*pi
k2x = 0 #final k point
k2y = pi
k2z = pi
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#N - Gamma
k1x = 0 #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma -P
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#P - H
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = 0
k2z = 2*pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+l[i])
Gy =2*pi*(h[i]+k[i])
Gz =2*pi*(k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','H','N','$\Gamma$','P','H']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([-5,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Face Centered Cubic with a zero potential
A good way to check if your band structure code is working is to calculate the band structure for a zero potential. This should reproduce the results of the empty lattice approximation. Below is some python code that uses the plane wave method to calculate the band structure of a three-dimensional fcc lattice where the periodic potential has zero amplitude.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
#The off-diagonal elements are zero for a zero potential
M = np.zeros((125,125)) #plane wave matrix
kvec = np.empty(909) # wavelector
energy = np.empty((125,909)) #energies
#Gamma 0, 0, 0
#X 0, 2*pi, 0
#L pi, pi, pi
#W pi, 2*pi, 0
#U pi/2, 2*pi, pi/2
#K 3*pi/2, 3*pi/2, 0
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 2*pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - W
k1x = 0 #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#W - K
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#K - Gamma
k1x = 3*pi/2 #initial k point
k1y = 3*pi/2
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - L
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
#L - U
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = pi/2 #final k point
k2y = 2*pi
k2z = pi/2
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+505] = kl+k0
energy[:,j+505] = np.sort(np.real(la.eigvals(M))) # get energies
#U - W
k1x = pi/2 #initial k point
k1y = 2*pi
k1z = pi/2
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+606] = kl+k0
energy[:,j+606] = np.sort(np.real(la.eigvals(M))) # get energies
#W - L
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+707] = kl+k0
energy[:,j+707] = np.sort(np.real(la.eigvals(M))) # get energies
#L - K
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+808] = kl+k0
energy[:,j+808] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[505],kvec[606],kvec[707],kvec[807],kvec[908])
my_xticks=['$\Gamma$','X','W', 'K','$\Gamma$','L','U','W','L','K']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([0,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Face Centered Cubic with a Comb Potential
If all of the components of a Fourier series are real and have the same amplitude, there is constructive interference at the Bravais lattice points and this produces a comb function with one large amplitude peak per Bravais lattice point. Below is python code to calculate the band structure of a fcc comb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(909) # wavelector
energy = np.empty((125,909)) #energies
# Define the off-diagonal elements of the matrix M
# For a Comb potential, all of the off-diagonal matrix elements are the same
for i in range(125):
for j in range(125):
M[i,j] = -5; #Comb potential
#Gamma 0, 0, 0
#X 0, 2*pi, 0
#L pi, pi, pi
#W pi, 2*pi, 0
#U pi/2, 2*pi, pi/2
#K 3*pi/2, 3*pi/2, 0
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 2*pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - W
k1x = 0 #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#W - K
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#K - Gamma
k1x = 3*pi/2 #initial k point
k1y = 3*pi/2
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - L
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
#L - U
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = pi/2 #final k point
k2y = 2*pi
k2z = pi/2
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+505] = kl+k0
energy[:,j+505] = np.sort(np.real(la.eigvals(M))) # get energies
#U - W
k1x = pi/2 #initial k point
k1y = 2*pi
k1z = pi/2
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+606] = kl+k0
energy[:,j+606] = np.sort(np.real(la.eigvals(M))) # get energies
#W - L
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+707] = kl+k0
energy[:,j+707] = np.sort(np.real(la.eigvals(M))) # get energies
#L - K
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+808] = kl+k0
energy[:,j+808] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[505],kvec[606],kvec[707],kvec[807],kvec[908])
my_xticks=['$\Gamma$','X','W', 'K','$\Gamma$','L','U','W','L','K']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([5,50]) #plot only the low energies, the high energies are unreliable
plt.show()
Face Centered Cubic with a Coulomb Potential
Consider a periodic function consisting of a Coulomb potential at every Bravais lattice point. This potential is,
Here $Z$ is the effective nuclear charge, $\vec{r}_j$ is the position of the Bravais lattice point, and $\epsilon_0 = 8.854187817\times 10^{-12}\text{ F/m}$ is the permittivity constant. The Fourier series for this potential is,
where $V_{u.c.}$ is the volume of the primitive unit cell. For fcc with a conventional lattice constant $a$, there are four primitive unit cells in a conventional unit cell so $V_{u.c.}= \frac{a^3}{4}$
Below is python code to calculate the band structure of a fcc Coulomb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
# Define the off-diagonal elements of the matrix M
# The Fourier series for a Columb potential is example 4 at
# http://lampz.tugraz.at/~hadley/ss1/crystaldiffraction/fourier/fourier_examples.php
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(909) # wavelector
energy = np.empty((125,909)) #energies
for i in range(125):
Gxi =2*pi*(h[i]+k[i]-l[i])
Gyi =2*pi*(-h[i]+k[i]+l[i])
Gzi =2*pi*(h[i]-k[i]+l[i])
for j in range(i+1,125):
Gxj =2*pi*(h[j]+k[j]-l[j])
Gyj =2*pi*(-h[j]+k[j]+l[j])
Gzj =2*pi*(h[j]-k[j]+l[j])
Gx = Gxj - Gxi
Gy = Gyj - Gyi
Gz = Gzj - Gzi
G2 = Gx*Gx + Gy*Gy + Gz*Gz #(Gi -Gj)^2
#G is in units of 1/a
# M is in units of eV, Z = 1, Vuc = 0.25*a^3
M[i,j]=-4*e/(a*eps0*G2)
M[j,i] = M[i,j]
#Gamma 0, 0, 0
#X 0, 2*pi, 0
#L pi, pi, pi
#W pi, 2*pi, 0
#U pi/2, 2*pi, pi/2
#K 3*pi/2, 3*pi/2, 0
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = 2*pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - W
k1x = 0 #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#W - K
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#K - Gamma
k1x = 3*pi/2 #initial k point
k1y = 3*pi/2
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - L
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
#L - U
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = pi/2 #final k point
k2y = 2*pi
k2z = pi/2
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+505] = kl+k0
energy[:,j+505] = np.sort(np.real(la.eigvals(M))) # get energies
#U - W
k1x = pi/2 #initial k point
k1y = 2*pi
k1z = pi/2
k2x = pi #final k point
k2y = 2*pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+606] = kl+k0
energy[:,j+606] = np.sort(np.real(la.eigvals(M))) # get energies
#W - L
k1x = pi #initial k point
k1y = 2*pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+707] = kl+k0
energy[:,j+707] = np.sort(np.real(la.eigvals(M))) # get energies
#L - K
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 3*pi/2 #final k point
k2y = 3*pi/2
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*(h[i]+k[i]-l[i])
Gy =2*pi*(-h[i]+k[i]+l[i])
Gz =2*pi*(h[i]-k[i]+l[i])
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+808] = kl+k0
energy[:,j+808] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[505],kvec[606],kvec[707],kvec[807],kvec[908])
my_xticks=['$\Gamma$','X','W', 'K','$\Gamma$','L','U','W','L','K']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([-2,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Simple Cubic with a zero potential
A good way to check if your band structure code is working is to calculate the band structure for a zero potential. This should reproduce the results of the empty lattice approximation. Below is some python code that uses the plane wave method to calculate the band structure of a three-dimensional simple cubic lattice where the periodic potential has zero amplitude.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
#The off-diagonal elements are zero for a zero potential
M = np.zeros((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
#Gamma 0, 0, 0
#X 0, pi, 0
#M pi, pi, 0
#R pi, pi, pi
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - M
k1x = 0 #initial k point
k1y = pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#M - Gamma
k1x = pi #initial k point
k1y = pi
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - R
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#R - X
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','X','M','$\Gamma$','R','X']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([0,20]) #plot only the low energies, the high energies are unreliable
plt.show()
Simple Cubic with a Comb Potential
If all of the components of a Fourier series are real and have the same amplitude, there is constructive interference at the Bravais lattice points and this produces a comb function with one large amplitude peak per Bravais lattice point. Below is python code to calculate the band structure of a simple cubic comb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
# Define the off-diagonal elements of the matrix M
# For a Comb potential, all of the off-diagonal matrix elements are the same
for i in range(125):
for j in range(125):
M[i,j] = -5; #Comb potential
#Gamma 0, 0, 0
#X 0, pi, 0
#M pi, pi, 0
#R pi, pi, pi
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - M
k1x = 0 #initial k point
k1y = pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#M - Gamma
k1x = pi #initial k point
k1y = pi
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - R
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#R - X
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','X','M','$\Gamma$','R','X']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([5,30]) #plot only the low energies, the high energies are unreliable
plt.show()
Simple Cubic with a Coulomb Potential
Consider a periodic function consisting of a Coulomb potential at every Bravais lattice point. This potential is,
Here $Z$ is the effective nuclear charge, $\vec{r}_j$ is the position of the Bravais lattice point, and $\epsilon_0 = 8.854187817\times 10^{-12}\text{ F/m}$ is the permittivity constant. The Fourier series for this potential is,
where $V_{u.c.}$ is the volume of the primitive unit cell. For simple cubic with a conventional lattice constant $a$, there is one primitive unit cell in a conventional unit cell so $V_{u.c.}= a^3$
Below is python code to calculate the band structure of a simple cubic Coulomb potential. One hundred twenty five reciprocal lattice points are included in the calculation where $h$, $k$, and $l$ are in the range $-2,\,-2,\,0,\,1,\,2$.
Every band is drawn in a different color and there must be an energy for that band for every $k$. If a color is not visible for some range of $k$ that band must be degenerate in energy with another band.
import numpy as np
import scipy.linalg as la
import math
import matplotlib.pyplot as plt
pi = math.pi
m = 9.1093837015E-31 # [kg] electron mass
e = 1.602176634E-19 # [C] elementary charge
hbar = 1.054571817E-34 # [J s] reduced Planck's constant
eps0 = 8.8541878128E-12 # [F/m] permittivity of vacuum
a = 3E-10 # lattice constant
alpha = hbar*hbar/(2*m*a*a*e) # prefactor for kinetic energy term
#define the integers that label the G vectors
h = np.empty(125)
k = np.empty(125)
l = np.empty(125)
i = 0
for nx in range(5):
for ny in range(5):
for nz in range(5):
h[i] = nx-2 # h = -2,-1,0,1,2
k[i] = ny-2 # k = -2,-1,0,1,2
l[i] = nz-2 # l = -2,-1,0,1,2
i= i+1
M = np.empty((125,125)) #plane wave matrix
kvec = np.empty(505) # wavelector
energy = np.empty((125,505)) #energies
# Define the off-diagonal elements of the matrix M
# The Fourier series for a Columb potential is example 4 at
# http://lampz.tugraz.at/~hadley/ss1/crystaldiffraction/fourier/fourier_examples.php
for i in range(125):
Gxi = 2*pi*h[i]
Gyi = 2*pi*k[i]
Gzi = 2*pi*l[i]
for j in range(i+1,125):
Gxj = 2*pi*h[j]
Gyj = 2*pi*k[j]
Gzj = 2*pi*l[j]
Gx = Gxj - Gxi
Gy = Gyj - Gyi
Gz = Gzj - Gzi
G2 = Gx*Gx + Gy*Gy + Gz*Gz #(Gi -Gj)^2
#G is in units of 1/a
# M is in units of eV, Z = 1, Vuc = a^3
M[i,j]=-e/(a*eps0*G2)
M[j,i] = M[i,j]
#Gamma 0, 0, 0
#X 0, pi, 0
#M pi, pi, 0
#R pi, pi, pi
#Gamma - X
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = 0 #final k point
k2y = pi
k2z = 0
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
# define the diagonal elements of the matrix M
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j] = kl
energy[:,j] = np.sort(np.real(la.eigvals(M))) # get energies
#X - M
k1x = 0 #initial k point
k1y = pi
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = 0
k0 = klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+101] = kl+k0
energy[:,j+101] = np.sort(np.real(la.eigvals(M))) # get energies
#M - Gamma
k1x = pi #initial k point
k1y = pi
k1z = 0
k2x = 0 #final k point
k2y = 0
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+202] = kl+k0
energy[:,j+202] = np.sort(np.real(la.eigvals(M))) # get energies
#Gamma - R
k1x = 0 #initial k point
k1y = 0
k1z = 0
k2x = pi #final k point
k2y = pi
k2z = pi
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+303] = kl+k0
energy[:,j+303] = np.sort(np.real(la.eigvals(M))) # get energies
#R - X
k1x = pi #initial k point
k1y = pi
k1z = pi
k2x = 0 #final k point
k2y = pi
k2z = 0
k0 = k0 + klength
klength = math.sqrt((k2x - k1x)*(k2x - k1x)+(k2y - k1y)*(k2y - k1y)+(k2z - k1z)*(k2z - k1z))
for j in range(101):
kx = k1x + (k2x -k1x)*j/100
ky = k1y + (k2y -k1y)*j/100
kz = k1z + (k2z -k1z)*j/100
kl = klength*j/100
for i in range(125):
Gx =2*pi*h[i]
Gy =2*pi*k[i]
Gz =2*pi*l[i]
M[i,i] = alpha*((kx-Gx)*(kx-Gx)+(ky-Gy)*(ky-Gy)+(kz-Gz)*(kz-Gz))
kvec[j+404] = kl+k0
energy[:,j+404] = np.sort(np.real(la.eigvals(M))) # get energies
plt.figure()
for pp in range(0,125):
plt.plot(kvec,energy[pp,:])
plt.xlabel('$ka$')
plt.ylabel('$E$ [eV]')
positions = (kvec[0],kvec[101],kvec[202],kvec[303],kvec[404],kvec[504])
my_xticks=['$\Gamma$','X','M','$\Gamma$','R','X']
plt.xticks(positions, my_xticks)
plt.grid(axis = 'x')
plt.ylim([-5,20]) #plot only the low energies, the high energies are unreliable
plt.show()