岩爆预测的多维云模型M2X3类型二维材料的量子反常霍尔效应

发布时间：2026-07-25阅读(1)

导读前文回顾之前介绍过SlaterKoster系数的自动获取，以及M2X3的QAHE的哈密顿分析，这一节在之前基础上，将M2X3这个体系，使用一个通用的TB模型....

前文回顾

之前介绍过 Slater Koster系数的自动获取，以及M2X3的QAHE的哈密顿分析，

这一节在之前基础上，将M2X3这个体系，使用一个通用的TB模型代码，将其能带以及量子反常霍尔效应展示出来。

这个相对通用的TB模型，还差一点点东西，后面会把SOC代码块加进来。

结构

MX结构：M为金属，X为O族元素或者有机配体

代码功能

给出晶格，轨道，建能代码会自动给出Slater Koster系数，给出哈密顿，给出本征值、能带，后续会给出贝里曲率、陈数。

现在这个代码还没写全，还差一个通用的贝里曲率计算代码块，soc矩阵代码块以及陈数计算代码块。这一节先手动补上，后面有时间会写成通用的代码块。

代码展示

from numpy import *import numpy.matlib as mtlimport matplotlib.pyplot as pltimport matplotlibfrom 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 energydd_sigma = 0dd_pi = 0.4096dd_delta = 0.0259num = 4 ##orbit number## connect vector and direction cosindef 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 8def 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 resultsdef 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 elementt131,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 Hamiltoniandef 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 Hdef eH(K): return linalg.eigh(H(K))[0] #reciprocal latticeb1=array([1,1/sqrt(3)])*pi*2b2=array([0,2/sqrt(3)])*pi*2#高对称点G=array([0,0])M=0.5*b1K= 1/3 * b1 1/3 * b2K2= -1/3 * b1 - 1/3 * b2#K点路径G-Mkgm = 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 * lkxmk = 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$,M,K,r$\Gamma$,K\])plt.show()

不含SOC结果

四带模型哈密顿

SOC

加入SOC后, 补充到哈密顿中，求本征值

#Hsoc: lambda = 0.02lam = 0.02Hsoc = 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-6numk=201def dHx(k): k2 = k - np.array([dkp,0]) return (H(k) - H(k2))/dkpdef 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 vdhcdef 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 Axes3Dfig = 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()