前文回顧
之前介紹過 Slater Koster系數的自動擷取，以及M2X3的QAHE的哈密頓分析，
這一節在之前基礎上，将M2X3這個體系，使用一個通用的TB模型代碼，将其能帶以及量子反常霍爾效應展示出來。
這個相對通用的TB模型，還差一點點東西，後面會把SOC代碼塊加進來。
結構
MX結構：M為金屬，X為O族元素或者有機配體
代碼功能
給出晶格，軌道，建能 代碼會自動給出Slater Koster系數，給出哈密頓，給出本征值、能帶，後續會給出貝裡曲率、陳數。
現在這個代碼還沒寫全，還差一個通用的貝裡曲率計算代碼塊，soc矩陣代碼塊以及陳數計算代碼塊。這一節先手動補上，後面有時間會寫成通用的代碼塊。
代碼展示
from numpy import *
import numpy.matlib as mtl
import matplotlib.pyplot as plt
import matplotlib
from matplotlib import cm
### lattice and orbit
# s py pz px dxy dyz dz2 dxz d(x2-y2)
# 0 1 2 3 4 5 6 7 8
lat = [[1,0],[-1/2, sqrt(3)/2]]
lat=array(lat)
o1 = [2/3,1/3, 5]
o2 = [2/3,1/3, 7]
o3 = [1/3,2/3, 5]
o4 = [1/3,2/3, 7]
## bond energy
dd_sigma = 0
dd_pi = 0.4096
dd_delta = 0.0259
num = 4 ##orbit number
## connect vector and direction cosin
def vector(o1,o2,x,y):
r = (o2[0:2]+array([x,y])-o1[0:2]).dot(lat)
l = r[0]/sqrt(r[0]**2+r[1]**2)
m = r[1]/sqrt(r[0]**2+r[1]**2)
return r,l, m
## SK parameters
# s py pz px dxy dyz dz2 dxz d(x2-y2)
# 0 1 2 3 4 5 6 7 8
def sk(l,m,d1,d2,dd_sigma,dd_pi,dd_delta):
sk = zeros((9,9))
sk[4,4] = 3*l**2 *m**2 * dd_sigma + (l**2 + m**2 - 4 * l**2 * m**2) * dd_pi + l**2 * m**2 * dd_delta
sk[4,5] = 0
sk[4,6] = -0.5*sqrt(3)*l*m*(l**2+m**2)*dd_sigma + 0.5*l*m*dd_delta
sk[4,7] = 0
sk[4,8] = 3/2 * l * m *(l**2-m**2) * dd_sigma + 2*l*m*(m**2-l**2) * dd_pi + 0.5*l*m*(l**2-m**2)*dd_delta
sk[5,4] = 0
sk[5,5] = m**2*dd_pi+ l**2*dd_delta
sk[5,6] = 0
sk[5,7] = l*m*dd_pi -l*m*dd_delta
sk[5,8] = 0
sk[6,4] = sk[4,6]
sk[6,5] = 0
sk[6,6] = 1/4 * (l**2+m**2)**2 * dd_sigma + 3/4 * (l**2+m**2)**2 * dd_delta
sk[6,7] = 0
sk[6,8] = sqrt(3)/4 * (l**2-m**2) * dd_delta - sqrt(3)/4 * (l**2-m**2) * (l**2+m**2) * dd_sigma
sk[7,4] = 0
sk[7,5] = sk[5,7]
sk[7,6] = 0
sk[7,7] = l**2 * dd_pi + m**2 * dd_delta
sk[7,8] = 0
sk[8,4] = sk[4,8]
sk[8,5] = 0
sk[8,6] = sk[6,8]
sk[8,7] = sk[7,8]
sk[8,8] = 3/4 * (l**2-m**2) * dd_sigma + (l**2 + m**2 - (l**2-m**2)**2) * dd_pi + \
1/4 * (l**2-m**2)**2 * dd_delta
return sk[d1,d2]
### give connect vector and hopping from SK results
def hop_vec(o1,o2,x,y):
r,l,m =vector(o1,o2,x,y)
d1,d2 = o1[2],o2[2]
t=sk(l,m,d1,d2,dd_sigma,dd_pi,dd_delta)
return t,r
## give hopping t and connect vector for the matrix element
t131,r131 = hop_vec(o1,o3,0,0)
t141,r141 = hop_vec(o1,o4,0,0)
t132,r132 = hop_vec(o1,o3,1,0)
t142,r142 = hop_vec(o1,o4,1,0)
t133,r133 = hop_vec(o1,o3,0,-1)
t143,r143 = hop_vec(o1,o4,0,-1)
t231,r231 = hop_vec(o2,o3,0,0)
t241,r241 = hop_vec(o2,o4,0,0)
t232,r232 = hop_vec(o2,o3,1,0)
t242,r242 = hop_vec(o2,o4,1,0)
t233,r233 = hop_vec(o2,o3,0,-1)
t243,r243 = hop_vec(o2,o4,0,-1)
## eigenvalue of Hamiltonian
def H(K):
H=zeros((num,num),dtype=complex)
H[0,2] = t131 * exp(1.j*K.dot(r131)) + t132 * exp(1.j*K.dot(r132)) + t133 * exp(1.j*K.dot(r133))
H[0,3] = t141 * exp(1.j*K.dot(r141)) + t142 * exp(1.j*K.dot(r142)) + t143 * exp(1.j*K.dot(r143))
H[1,2] = t231 * exp(1.j*K.dot(r231)) + t232 * exp(1.j*K.dot(r232)) + t233 * exp(1.j*K.dot(r233))
H[1,3] = t241 * exp(1.j*K.dot(r241)) + t242 * exp(1.j*K.dot(r242)) + t243 * exp(1.j*K.dot(r243))
H[0,0]= 0
H[1,1]= 0
H[2,2]= 0
H[3,3]= 0
for i in range(num):
for j in range(num):
if j > i:
H[j,i] = conj(H[i,j])
return H
def eH(K):
return linalg.eigh(H(K))[0]
#reciprocal lattice
b1=array([1,1/sqrt(3)])*pi*2
b2=array([0,2/sqrt(3)])*pi*2
#高對稱點
G=array([0,0])
M=0.5*b1
K= 1/3 * b1 + 1/3 * b2
K2= -1/3 * b1 - 1/3 * b2
#K點路徑G-M
kgm = linspace(G,M,100,endpoint=False)
kmk = linspace(M,K,100,endpoint=False)
kkg = linspace(K,G,100)
kgk2 = linspace(G,K2,100)
##K點相對距離
def Dist(r1,r2):
return linalg.norm(r1-r2)
lgm=Dist(G,M)
lmk=Dist(M,K)
lkg=Dist(K,G)
lgk2=Dist(G,K2)
lk = linspace(0,1,100)
xgm = lgm * lk
xmk = lmk * lk + xgm[-1]
xkg = lkg * lk + xmk[-1]
xgk2 = lgk2 * lk + xkg[-1]
kpath = concatenate((xgm,xmk,xkg,xgk2),axis=0)
Node = [0,xgm[-1],xmk[-1],xkg[-1],xgk2[-1]]
## 按照高對稱點求本征值
Eig_gm = array(list(map(eH,kgm)))
Eig_mk = array(list(map(eH,kmk)))
Eig_kg = array(list(map(eH,kkg)))
Eig_gk2 = array(list(map(eH,kgk2)))
eig_1 = hstack((Eig_gm[:,0],Eig_mk[:,0],Eig_kg[:,0],Eig_gk2[:,0]))
eig_2 = hstack((Eig_gm[:,1],Eig_mk[:,1],Eig_kg[:,1],Eig_gk2[:,1]))
eig_3 = hstack((Eig_gm[:,2],Eig_mk[:,2],Eig_kg[:,2],Eig_gk2[:,2]))
eig_4 = hstack((Eig_gm[:,3],Eig_mk[:,3],Eig_kg[:,3],Eig_gk2[:,3]))
##繪圖
plt.plot(kpath,eig_1)
plt.plot(kpath,eig_2)
plt.plot(kpath,eig_3)
plt.plot(kpath,eig_4)
#plt.plot(kpath,eig_5)
#plt.plot(kpath,eig_6)
plt.xticks(Node,[r'$\Gamma#39;,'M','K',r'$\Gamma#39;,'K\''])
plt.show()
不含SOC結果
四帶模型哈密頓
SOC
加入SOC後, 補充到哈密頓中，求本征值
#Hsoc: lambda = 0.02
lam = 0.02
Hsoc = array([[0,1.j*lam,0,0],
[-1.j*lam,0,0,0],
[0,0,0,1.j*lam],
[0,0,-1.j*lam,0]])
def esocH(K):
return linalg.eigh(H(K)+Hsoc)[0]
含soc的能帶
四帶模型的貝裡曲率
dkp =1e-6
numk=201
def dHx(k):
k2 = k - np.array([dkp,0])
return (H(k) - H(k2))/dkp
def dHy(k):
k2 = k - np.array([0,dkp])
return (H(k) - H(k2))/dkp
#按順序比對對應的本征值和本征矢
def ewH(k):
e,w=np.linalg.eigh(H(k)+Hsoc)
w0 = w[:, np.argsort(np.real(e))[0]]
w1 = w[:, np.argsort(np.real(e))[1]]
w2 = w[:, np.argsort(np.real(e))[2]]
w3 = w[:, np.argsort(np.real(e))[3]]
e = np.sort(np.real(e))
return [w0,w1,w2,w3],e
#<v|dH/dk|c> v,c 取0對應價帶波函數，1對應導帶波函數
def vcdHx(v,k,c):
dhc = np.dot(dHx(k), ewH(k)[0][c])
vdhc = np.dot(ewH(k)[0][v].conj(),dhc)
return vdhc
def vcdHy(v,k,c):
dhc = np.dot(dHy(k), ewH(k)[0][c])
vdhc = np.dot(ewH(k)[0][v].conj(),dhc)
return vdhc
def Omega(k,v,c):
return 1.j * (vcdHx(v,k,c) * vcdHy(c,k,v) - vcdHy(v,k,c) * vcdHx(c,k,v))/((ewH(k)[1][v]-ewH(k)[1][c])**2)
xx = np.linspace(-1.6*np.pi,1.6*np.pi,numk)
yy = np.linspace(-1.6*np.pi,1.6*np.pi,numk)
Z= np.zeros((numk,numk))
X,Y = np.meshgrid(xx,yy)
for i in range(numk):
for j in range(numk):
k = np.array([xx[i],yy[j]])
Z[i][j]=np.real(Omega(k,1,2))+np.real(Omega(k,1,3))+ np.real(Omega(k,1,0)) + np.real(Omega(k,0,1)) + np.real(Omega(k,0,2)) + np.real(Omega(k,0,3))
貝裡曲率繪圖
from mpl_toolkits.mplot3d.axes3d import Axes3D
fig = plt.figure()
ax = fig.add_subplot(1, 1, 1, projection='3d')
surf = ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.coolwarm,
linewidth=0, antialiased=False)
#ax.set_xlim3d(0.0, num)
#ax.set_ylim3d(0.0, num)
#ax.set_zlim3d(-0.90, 0.90)
ax.view_init(elev=50,azim=50)
plt.show()
陳數計算
這裡積分區間如上圖，取了三個FBZ的面積是以最後除以3
dk=3.2*np.pi/(numk-1)
C=sum(Z)*dk*dk/2/np.pi/3
print("integral of Berry curvature is %f" %C)
print("Chern number is %d" %round(C))