前文回顾

前面说到，使用python构建了自动生成Slater-Koster系数以及连接矢的代码，这一节在此基础上，构建一个六角晶格的TB模型，可以跟前文提到的两个案例和对比。这个体系是量子反常霍尔效应，点赞多的话，后面会更新它的贝里曲率以及陈数的计算代码，做成一个相对通用的二维体系TB模型。

前文回顾

1. 自动获取SK系数的python代码
2. 六角晶格的SK系数
3. SK系数以及隐藏的SK系数获取
4. Kagome晶格TB模型

构建紧束缚模型需要知道连接矢以及对应的hopping参数t，而这个t是可以通过SK方法获得。

TB部分代码

晶格和轨道

### 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]
num = 4 #一共四个轨道

给定键能

dd_sigma = 0
dd_pi = 0.4096
dd_delta = 0.0259

自动生成连接矢量和SK的hopping系数t

我们需要给出原子之间的hopping方向，pythtb软件也是这么做的：

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)

哈密顿 本征值

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)) 
        
    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]

k空间布点

#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

#K点路径G-M
kgm = linspace(G,M,100,endpoint=False)
kmk = linspace(M,K,100,endpoint=False)
kkg = linspace(K,G,100)

两个高对称K之间的距离比例

##K点相对距离
def Dist(r1,r2):
    return linalg.norm(r1-r2) 
lgm=Dist(G,M)
lmk=Dist(M,K)
lkg=Dist(K,G)

lk = linspace(0,1,100)
xgm = lgm * lk
xmk = lmk * lk + xgm[-1]
xkg = lkg * lk + xmk[-1]
kpath = concatenate((xgm,xmk,xkg),axis=0)
Node = [0,xgm[-1],xmk[-1],xkg[-1]] 

计算本征值

Eig_gm = array(list(map(eH,kgm)))
Eig_mk = array(list(map(eH,kmk)))
Eig_kg = array(list(map(eH,kkg)))
eig_1 = hstack((Eig_gm[:,0],Eig_mk[:,0],Eig_kg[:,0]))
eig_2 = hstack((Eig_gm[:,1],Eig_mk[:,1],Eig_kg[:,1]))
eig_3 = hstack((Eig_gm[:,2],Eig_mk[:,2],Eig_kg[:,2]))
eig_4 = hstack((Eig_gm[:,3],Eig_mk[:,3],Eig_kg[:,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;])
plt.show()

结果

完整代码

见下一节推文