Band structures#
To compute a band structure we do two things. First, we compute the self-consistent band structure. Then we compute the band structure at the desired \(k\)-points. We will use Si as an example (adapted from http://bbs.sciencenet.cn/bbs/upload/20083418325986.pdf).
First, we get the self-consistent electron density in a calculation.
from vasp import Vasp
from ase import Atom, Atoms
from ase.visualize import view
a = 5.38936
atoms = Atoms([Atom('Si', [0, 0, 0]),
Atom('Si', [0.25, 0.25, 0.25])])
atoms.set_cell([[a / 2., a / 2., 0.0],
[0.0, a / 2., a / 2.],
[a / 2., 0.0, a / 2.]], scale_atoms=True)
calc = Vasp(label='bulk/Si-selfconsistent',
xc='PBE',
prec='Medium',
lcharg=True,
lwave=True,
kpts=(4, 4, 4),
atoms=atoms)
calc.calculate()
Now, we run a new calculation along the k-point path desired. The standard VASP way of doing this is to modify the INCAR and KPOINTS file and rerun VASP. We will not do that. Doing that results in some lost information if you overwrite the old files. We will copy the old directory to a new directory, using code to ensure this only happens one time.
from vasp import Vasp
from vasp.exceptions import VaspNotConverged
wd = 'bulk/Si-bandstructure'
calc = Vasp(label='bulk/Si-selfconsistent').clone(wd)
kpts = [[0.5, 0.5, 0.0], # L
[0, 0, 0], # Gamma
[0, 0, 0],
[0.5, 0.5, 0.5]] # X
calc.set(kpts=kpts,
reciprocal=True,
kpts_nintersections=10,
icharg=11)
try:
calc.calculate()
except VaspNotConverged as e:
print(f'VaspNotConverged: {e}')
print('Note: Run the self-consistent cell above first to generate CHGCAR.')
-3.62224484
We will learn how to manually parse the EIGENVAL file here to generate the band structure. The structure of the EIGENVAL file looks like this:
head -n 20 bulk/Si-bandstructure/EIGENVAL
2 2 1 1
0.1956688E+02 0.3810853E-09 0.3810853E-09 0.3810853E-09 0.5000000E-15
1.000000000000000E-004
CAR
unknown system
8 20 8
0.5000000E+00 0.5000000E+00 0.0000000E+00 0.5000000E-01
1 -1.826747
2 -1.826743
3 3.153321
4 3.153347
5 6.743989
6 6.744017
7 16.392596
8 16.393943
0.4444444E+00 0.4444444E+00 0.0000000E+00 0.5000000E-01
1 -2.669487
2 -0.918463
2 2 1 1
0.1956688E+02 0.3810853E-09 0.3810853E-09 0.3810853E-09 0.5000000E-15
1.000000000000000E-004
CAR
unknown system
8 20 8
0.5000000E+00 0.5000000E+00 0.0000000E+00 0.5000000E-01
1 -1.826747
2 -1.826743
3 3.153321
4 3.153347
5 6.743989
6 6.744017
7 16.392596
8 16.393943
0.4444444E+00 0.4444444E+00 0.0000000E+00 0.5000000E-01
1 -2.669487
2 -0.918463
We can ignore the first five lines.
import os
import matplotlib.pyplot as plt
eigenval_path = 'bulk/Si-bandstructure/EIGENVAL'
if not os.path.exists(eigenval_path):
print(f'EIGENVAL file not found at {eigenval_path}')
print('Run the band structure calculation cell above first.')
else:
with open(eigenval_path, 'r') as f:
line1 = f.readline()
line2 = f.readline()
line3 = f.readline()
line4 = f.readline()
comment = f.readline()
unknown, nkpoints, nbands = [int(x) for x in f.readline().split()]
blankline = f.readline()
band_energies = [[] for i in range(nbands)]
for i in range(nkpoints):
x, y, z, weight = [float(x) for x in f.readline().split()]
for j in range(nbands):
fields = f.readline().split()
id, energy = int(fields[0]), float(fields[1])
band_energies[id - 1].append(energy)
blankline = f.readline()
for i in range(nbands):
plt.plot(range(nkpoints), band_energies[i])
ax = plt.gca()
ax.set_xticks([]) # no tick marks
plt.xlabel('k-vector')
plt.ylabel('Energy (eV)')
ax.set_xticks([0, 10, 19])
ax.set_xticklabels(['$L$', '$\\Gamma$', '$X$'])
plt.savefig('images/Si-bandstructure.png')
Next we will examine the connection between band structures and density of states. In this example, we will compute the band structure of TiO2 using a function built into mod:vasp to do the analysis described above.
from vasp import Vasp
calc = Vasp(label='bulk/tio2/step3')
print(calc.get_fermi_level())
calc.abort()
n, bands, p = calc.get_bandstructure(kpts_path=[('$\Gamma$', [0.0, 0.0, 0.0]),
('X', [0.5, 0.5, 0.0]),
('X', [0.5, 0.5, 0.0]),
('M', [0.0, 0.5, 0.5]),
('M', [0.0, 0.5, 0.5]),
('$\Gamma$', [0.0, 0.0, 0.0])])
p.savefig('images/tio2-bandstructure-dos.png')
create example showing band dispersion with change in lattice constant#
In this section, we examine the effect of the lattice constant on the band structure. Since the lattice constant affects the overlap of neighboring atoms, we expect that smaller lattice constants will show more dispersion, i.e. broader bands. Larger lattice constants, in contrast, should show narrower bands. We examine this in silicon.
from vasp import Vasp
from ase import Atom, Atoms
calcs = []
for i, a in enumerate([4.7, 5.38936, 6.0]):
atoms = Atoms([Atom('Si', [0, 0, 0]),
Atom('Si', [0.25, 0.25, 0.25])])
atoms.set_cell([[a/2., a/2., 0.0],
[0.0, a/2., a/2.],
[a/2., 0.0, a/2.]], scale_atoms=True)
calc = Vasp(label=f'bulk/Si-bs-{i}',
xc='PBE',
lcharg=True,
lwave=True,
kpts=(4, 4, 4),
atoms=atoms)
print(calc.run())
calcs += [calc]
Vasp.wait(abort=True)
for i, calc in enumerate(calcs):
n, bands, p = calc.get_bandstructure(kpts_path=[('L', [0.5,0.5,0.0]),
('$\Gamma$', [0, 0, 0]),
('$\Gamma$', [0, 0, 0]),
('X', [0.5, 0.5, 0.5])],
kpts_nintersections=10)
if p is not None:
png = f'images/Si-bs-{i}.png'
p.savefig(png)
-7.55662509
-10.80024435
-10.13735105
You can see the band structure for a=6.0 is notably sharper than the band structure for a=4.0.