Skip to content

Instantly share code, notes, and snippets.

@marekyggdrasil

marekyggdrasil/run.py

Last active May 7, 2020
Embed
What would you like to do?
Simulation of the 1D Tight-Binding Model, full tutorial https://mareknarozniak.com/2020/05/07/tight-binding/
import matplotlib.pyplot as plt
import numpy as np
from qutip import basis
def _ket(n, N):
return basis(N, n)
def _bra(n, N):
return basis(N, n).dag()
# spacing between atoms
a = 2
# number of atoms
N = 23
# total length of the 1D chain
L = N*a
# onsite energy
eps = 0.25
# hopping matrix element
t = 0.5
# construct the Hamiltonian
H = sum([eps*_ket(n, L)*_bra(n, L) for n in range(0, L)])
H -= sum([t*_ket(n, L)*_bra(n + 1, L) for n in range(0, L - 1)])
H -= sum([t*_ket(n, L)*_bra(n - 1, L) for n in range(1, L)])
if N <= 8:
print('Hamiltonian')
print(H)
evals, ekets = H.eigenstates()
if N <= 8:
print()
print('Eigenvalues')
print(evals)
# satisfy periodic boundary conditions we need E(-k) = E(k)
numerical = np.concatenate((np.flip(evals, 0), evals), axis=0)
k = np.linspace(-np.pi/a, np.pi/a, 2*L)
exact = eps - 2.*t*np.cos(k*a)
xticks = np.linspace(-np.pi/a, np.pi/a, 9)
xlabels = ['' for k in xticks]
xlabels[0] = '$-\\frac{\pi}{a}$'
xlabels[-1] = '$\\frac{\pi}{a}$'
fig, axs = plt.subplots()
axs.set_xlim(-np.pi/a, np.pi/a)
axs.set_title('Tight-Binding Model, $N='+str(N)+', a='+str(a)+'$')
axs.set_ylabel('$E$')
axs.set_xlabel('$ka$')
axs.axvline(x=0., color='k')
axs.plot(k, numerical, 'ro', label='Eigenvalues of H')
axs.plot(k, exact, label='$E(k) = \epsilon_0 - 2 t \cos(ka)$')
axs.set_yticks([eps-2.*t, eps-t, eps, eps+t, eps+2.*t])
axs.set_yticklabels(['$\epsilon_0-2t$', '$\epsilon_0-t$', '$\epsilon_0$', '$\epsilon_0+t$', '$\epsilon_0+2t$'])
axs.set_xticks(xticks)
axs.set_xticklabels(xlabels)
axs.legend()
axs.grid(True)
# fig.savefig('tight_binding_theory.png')
fig.savefig('tight_binding_diag.png')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.