Easy ways to generate "teaching" band structures in Python?

Question

I'd like to introduce band structure to a class of undergraduate chemists, along the lines of Roald Hoffmann's Solids and Surfaces.

That is, I'd like to start with a s-band in 1D, which is easy because you can do it with hydrogen:

But then I'd like to look at different 2p-bands in 1D, 3d-bands in 1D, and build up a 2D and 3D band structure.

For example, here's a 2p band in 1D which naturally introduces nodes (from Solids and Surfaces by Roald Hoffmann):

And here's how this builds the 2D case with both s and p orbitals (also from Solids and Surfaces by Hoffmann):

Ideally all of these would make it fairly easy to change the lattice parameters to show the effects, delocalization vs. localization, etc.

There are plenty of codes and notes for real materials. Are there already some tight-binding packages / notebooks for teaching purposes?

Answer 1

I'd suggest having a look at PythTB by Sinisa Coh and David Vanderbilt, which I found very useful and easy to use.

Here's a stab at recreating the band structure of the 2d square lattice shown above

#!/usr/bin/env python
"""s and p orbitals on 2d square lattice.

 - nearest neighbor hopping only
 - no hopping between orbitals of different type (e.g. s <> pz)
"""

from pythtb import tb_model
import matplotlib.pyplot as plt

# 2d square lattice
lat=[[1.0, 0.0],[0.0, 1.0]]
# one s and 3 p orbitals (at the same site)
orb=[[0.,0.] for _ in range(1+3)]

# 2-dimensional k-space, 2-dimensional real space
my_model=tb_model(dim_r=2,dim_k=2,lat=lat,orb=orb)

# p-orbitals are at higher energy
my_model.set_onsite([0,3,3,3])

# set hoppings (one for each connected pair of orbitals)
# (amplitude, i, j, [lattice vector to cell containing j])

# isotropic nearest-neighbor hopping for s and pz
for neighbor in [ [1,0], [0,1]]:
    my_model.set_hop(-0.5, 0, 0, neighbor) # s
    my_model.set_hop(-0.25, 3, 3, neighbor) # pz

# directional nearest-neighbor hopping for px and py
my_model.set_hop(+0.25, 1, 1, [1,0])
my_model.set_hop(-0.25, 1, 1, [0,1])
my_model.set_hop(-0.25, 2, 2, [1,0])
my_model.set_hop(+0.25, 2, 2, [0,1])

# generate list of k-points following a segmented path in the BZ
k_label=[r"$\Gamma$",r"$X$", r"$M$", r"$\Gamma$"]
path=[[0.0,0.0], [0.5,0.0], [0.5, 0.5], [0.0,0.0]]
(k_vec,k_dist,k_node)=my_model.k_path(path,100)

# solve model
evals=my_model.solve_all(k_vec)

# plot band structure
fig, ax = plt.subplots()
for band in evals:
    ax.plot(k_dist, band)
ax.set_title("s and p orbitals on 2d square lattice")
ax.set_xlabel("Path in k-space")
ax.set_ylabel("Band energy")
ax.set_xticks(k_node)
ax.set_xticklabels(k_label)
ax.set_xlim(k_node[0],k_node[-1])
for n in range(len(k_node)):
  ax.axvline(x=k_node[n], linewidth=0.5, color='k')
fig.tight_layout()
fig.savefig("s_and_p.pdf")

Note that matplotlib colors the lines by band energy at each k-point, not by the orbital symmetry of the band.

This may be suitable as a first simple model to understand how the hopping matrix elements influence the band structure.

From here, there are a number of possible extensions to learn more / make it more realistic:

• the nearest-neighbor hopping matrix elements between different orbitals are not entirely independent parameters and can be expressed in terms of bond integrals using the Slater-Koster tables.
• double the size of the unit cell, look at band folding and/or introduce site-dependent onsite energies (e.g. to simulate symmetry breaking when adsorbing a 2d material on a substrate)
• go 3d
• look at the effect of introducing 2nd- and 3rd-nearest neighbor hopping (e.g. in graphene, 2nd-nearest neighbor hopping affects the electronic structure near the Fermi level only by a rigid shift, while 3rd-nearest neighbor hopping plays a significant role)

The extensions above are still within the tight-binding framework, treating electrons as non-interacting Fermions, and can be modeled with pythtb.

P.S. Another conceptually very interesting extension would be to introduce an on-site Coulomb repulsion term between electrons, i.e. to move from tight binding to the Hubbard model. The Hubbard model is no longer solvable exactly (except for 1d), but a solution to its mean-field approximation can be found by iteration to self-consistency (similar to what DFT codes do, but on a lattice). This is not very difficult to do, but it goes beyond the scope of pythtb and I am not aware of a similar, easy-to-use software package that implements this functionality. If someone does, please let me know, I think it would make for a great addition to this toolbox.

P.P.S. There have been at least two attempts (one by me) to convince the authors to host the pythtb source code on GitHub to make it easy to contribute back (e.g. I made some minor modifications for my Ph.D. thesis). So far, they have been reluctant due to possible maintenance/supervision involved, but it may be worth trying again.

Answer 2

Another package (sisl) can do the same thing, but a bit more verbose and dynamic. (disclaimer, I am the author).
One of sisl's main goal is to interact with large TB models, and thus the complexity of the TB models are a bit more verbose, but it allows greater flexibility when one can query stuff on the fly. It is also intrinsically 3D so one has to deal with the extra dimension not used. It may for instance also be used to read in the Hamiltonian from other DFT packages, such as Siesta and Wannier90.

Here is the same TB models shown for the simple cubic, but also for 2x1 and 2x2 times the basic system, for checking band-folding as mentioned by @leopold.talirz

import numpy as np
import sisl
from matplotlib import pyplot as plt
"""s and p orbitals on 2d square lattice.

 - nearest neighbor hopping only
 - no hopping between orbitals of different type (e.g. s <> pz)
"""
# lattice constant
alat = 1.0 # Ang
# Create a square lattice with lattice constant as defined above
# Each atom have s and p orbitals and their orbital ranges are alat + .01 Ang
# (just slightly above the lattice constant)
s = sisl.AtomicOrbital("s", R=alat + 0.01)
px = sisl.AtomicOrbital("px", R=alat + 0.01)
py = sisl.AtomicOrbital("py", R=alat + 0.01)
pz = sisl.AtomicOrbital("pz", R=alat + 0.01)
# Since each orbital has a specific orbital range, one can also
# use different ranges per orbital (if needed)

# Create the atom, in sisl atoms needs to have a specie associated.
# They aren't internally used in the electronic structure calculation,
# so it doesn't matter which specie you choose here.
atom = sisl.Atom("H", [s, px, py, pz])
geometry = sisl.geom.sc(alat, atom)
# the sc is a simple cubic lattice, but we will truncate one dimension
# to make it 2D, the below code will remove periodicity along the 3rd
# lattice vector
geometry.set_nsc(c=1)

# Now define how the tight-binding model is denoted
# Since there are 2 orbitals per atom, you need to do define a
# simple method that will create the hoppings
def construct(H, ia, atoms, atoms_xyz=None):
    # ia: is the atom that we are currently assigning hopping elements to
    # atoms: is a subset of all atoms that is ensured to be reachable by ia
    #        in this small example it isn't needed, but when very large
    #        TB models are made, it can greatly speed up the construction
    #        since it reduces the search space.

    # convert to the first orbital on the atom
    io = H.geometry.a2o(ia)
    # this will give a list of:
    #   jas_onsite: R <= 0.1
    #   jas_nn: 0.1 < R <= s.R
    jas_onsite, jas_nn = H.geometry.close(ia, R=[0.1, s.R], atoms=atoms, atoms_xyz=atoms_xyz)

    # we know that there are two orbitals per atom,
    # the 1st is s (as defined in the sisl.Atom(...)
    # the 2nd is pz
    jos_onsite = H.geometry.a2o(jas_onsite)
    H[io, jos_onsite] = 0 # s
    for i in (1, 2, 3):
        H[io+i, jos_onsite+i] = 3 # p

    jos_nn = H.geometry.a2o(jas_nn)
    H[io, jos_nn] = -0.5 # s
    H[io+3, jos_nn+3] = -0.25 # pz

    # Now the anisotropic hoppings for px, py
    # first we need to get the supercell hoppings
    jas_nn_isc = H.geometry.a2isc(jas_nn)
    for jo, isc in zip(jos_nn, jas_nn_isc):
        if isc[0] != 0:
            # this *must* be hopping along x
            H[io+1, jo+1] = +0.25 # px
            H[io+2, jo+2] = -0.25 # py
        elif isc[1] != 0:
            # this *must* be hopping along y
            H[io+1, jo+1] = -0.25 # px
            H[io+2, jo+2] = +0.25 # py

# Define the Hamiltonian
H = sisl.Hamiltonian(geometry)
# Now construct the TB model
H.construct(construct)


def plot_bs(H):
    bs = sisl.physics.BandStructure(H,
                                    [[0, 0, 0], [0.5, 0, 0],
                                     [0.5, 0.5, 0], [0, 0, 0]],
                                    100, names=["Gamma", "X", "M", "Gamma"])

    # Calculate the eigenvalues
    evals = bs.apply.ndarray.eigh()

    # plot band structure
    fig, ax = plt.subplots()
    # retrieve the linear spacings of the k-points
    lk, kt, kl = bs.lineark(True)

    ax.set_xticks(kt)
    for band in evals.T:
        ax.plot(lk, band)
    ax.set_title("s and p orbitals on 2d square lattice")
    ax.set_xlabel("Path in k-space")
    ax.set_ylabel("Band energy")
    ax.set_xticks(kt)
    ax.set_xticklabels(kl)
    fig.tight_layout()

plot_bs(H)
# the same bandstructure, twice the size along x
plot_bs(H.tile(2, axis=0))
# the same bandstructure, twice the size along x and y
plot_bs(H.tile(2, axis=0).tile(2, axis=1))
plt.show()

The same plot:

The 2x1 system:

The 2x2 system:

See here for additional details: https://zerothi.github.io/sisl/tutorials/tutorial_es_1.html There one can also see how one can expand the eigenstates into real-space quantities.

Please remark that this TB model is not efficiently solved in one Hamiltonian. Since all orbitals are disjointed from the each other, there is no need to have 1 Hamiltonian, instead for better efficiency one could do 4 Hamiltonians. In sisl extracting sub orbitals is a breeze:

H_s = H.sub_orbital(atom, s)
# or via index
H_py = H.sub_orbital(atom, 2)

if you diagonalize and plot all 4 individually you'll see exactly the same band-structure.

Note that Hubbard is based on sisl and allows the Hubbard model quite easily.