import numpy as np
import matplotlib.pyplot as plt
import matplotlib.cm as cm
import matplotlib.colors as mcolors
import matplotlib as mpl
# 确保所有必需的函数都被定义
# --- 哈密顿量及其导数计算函数 ---
def ham(lat,orb,hop,k):
"""
计算给定k点的哈密顿量H(k)。
"""
h=np.zeros((len(orb),len(orb)),dtype='complex')
lat_np = np.array(lat) # 确保lat是numpy数组以便进行点积
for i in range (len(hop)):
t_val = hop[i][0] # 跳跃项的值
from_orb_idx = hop[i][1] # 起始轨道索引
to_orb_idx = hop[i][2] # 结束轨道索引
Rvec_coords = np.array(hop[i][3:5]) # 晶格位移矢量 (rx, ry)
# 计算真实空间中的位移矢量 R + r_to - r_from
# (Rvec_x, Rvec_y) + (orb_to_x - orb_from_x, orb_to_y - orb_from_y)
orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx])
latr = np.dot((Rvec_coords + orbital_diff), lat_np) # 转换为笛卡尔坐标系下的真实空间位移
newterm = t_val * np.exp(1j * (np.dot(latr, k)))
h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm
# 哈密顿量是厄米的,所以添加共轭转置部分
h = h + np.conjugate(h).transpose()
return h
def hamx(lat,orb,hop,k):
"""
计算哈密顿量对kx的导数 dH/dkx。
"""
h=np.zeros((len(orb),len(orb)),dtype='complex')
lat_np = np.array(lat)
for i in range (len(hop)):
t_val = hop[i][0]
from_orb_idx = hop[i][1]
to_orb_idx = hop[i][2]
Rvec_coords = np.array(hop[i][3:5])
orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx])
latr = np.dot((Rvec_coords + orbital_diff), lat_np)
newterm = t_val * np.exp(1j * (np.dot(latr, k))) * (1j * latr[0]) # 乘以 i * (R_x + r_y_x - r_x_x)
h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm
h = h + np.conjugate(h).transpose()
return h
def hamy(lat,orb,hop,k):
"""
计算哈密顿量对ky的导数 dH/dky。
"""
h=np.zeros((len(orb),len(orb)),dtype='complex')
lat_np = np.array(lat)
for i in range (len(hop)):
t_val = hop[i][0]
from_orb_idx = hop[i][1]
to_orb_idx = hop[i][2]
Rvec_coords = np.array(hop[i][3:5])
orbital_diff = np.array(orb[to_orb_idx]) - np.array(orb[from_orb_idx])
latr = np.dot((Rvec_coords + orbital_diff), lat_np)
newterm = t_val * np.exp(1j * (np.dot(latr, k))) * (1j * latr[1]) # 乘以 i * (R_y + r_y_y - r_x_y)
h[from_orb_idx, to_orb_idx] = h[from_orb_idx, to_orb_idx] + newterm
h = h + np.conjugate(h).transpose()
return h
# --- 用户提供的模型参数 ---
lat=np.array([[1.0,0.0],[0.0,1.0]])
orb=np.array([[0.0,0.0],[0.5,0.0],[0.0,0.5]])
t1=1.0
t2=0.4*np.exp((1.j)*np.pi*0.37)
t3=0.13
# 将pythtb的set_hop定义转换为我们ham函数所需的hop列表格式
# 格式: [hoppping_term, orbital_from, orbital_to, lattice_vector_x, lattice_vector_y]
hop=[
[t1, 0, 1, 0, 0],
[t1, 0, 2, 0, 0],
[t1, 0, 1, -1, 0],
[t1, 0, 2, 0, -1],
[t2, 1, 2, 0, 0],
[t2, 2, 1, -1, 0],
[t2, 1, 2, 1, -1],
[t2, 2, 1, 0, 1],
[t3, 0, 0, 1, 0],
[t3, 0, 0, 0, 1],
[t3, 1, 1, 1, 0],
[t3, 1, 1, 0, 1],
[t3, 2, 2, 1, 0],
[t3, 2, 2, 0, 1],
]
# 在位能 (onsite energies),根据原始代码设置为零
onsite=[0]*len(orb)
# --- K空间网格设置 ---
nx=200 # x方向K点数
ny=230 # y方向K点数
kxlim=2*np.pi # Kx方向的半范围 (从-kxlim到+kxlim)
kylim=2*np.pi # Ky方向的半范围 (从-kylim到+kylim)
# K空间总面积为 (2*kxlim) * (2*kylim)
# ds 为每个K点代表的K空间面积元
ds=(4*kxlim*kylim)/(nx*ny)
origin=[0,0] # K空间原点
# 初始化存储计算结果的数组
eccent=np.zeros((nx+1,ny+1)) # 特征椭球离心率
area=np.zeros((nx+1,ny+1)) # 特征椭球面积
curva=np.zeros((nx+1,ny+1)) # 贝里曲率 (未积分)
differ=np.zeros((nx+1,ny+1)) # tr(g)-|B|
kxlist=np.zeros(nx+1)
kylist=np.zeros(ny+1)
# --- K空间循环与计算 ---
print("🚀 开始计算K空间数据...")
for x in range(nx+1):
kx=(x-nx/2)/(nx/2)*kxlim+origin[0] # 线性映射kx到 [-kxlim, kxlim]
kxlist[x]=kx
for y in range(ny+1):
ky=(y-ny/2)/(ny/2)*kylim+origin[1] # 线性映射ky到 [-kylim, kylim]
kylist[y]=ky
# 1. 计算哈密顿量H(k)
H=ham(lat,orb,hop,[kx,ky])
for i in range(len(H)):
H[i,i]= H[i,i]+onsite[i] # 添加在位能
# 2. 计算哈密顿量对k的导数Hx和Hy
Hx=hamx(lat,orb,hop,[kx,ky])
Hy=hamy(lat,orb,hop,[kx,ky])
# 3. 对角化H(k)得到本征值e和本征矢量v
e,v=np.linalg.eigh(H) # e: eigenvalues, v: eigenvectors (columns are eigenvectors)
# 4. 计算量子度规张量Q
Q=np.zeros((2,2),dtype='complex')
# 根据量子度规张量的定义 Q_ab = sum_{n!=m} <psi_m|d_a H|psi_n><psi_n|d_b H|psi_m> / (E_m - E_n)^2
# 这里计算的是基态 (m=0) 的量子度规张量
for i in range(1,len(orb)): # 遍历所有激发态 (n=1, 2, ...)
# Q_xx
Q[0,0]=Q[0,0]+np.dot(v[:,0].conjugate(),np.dot(Hx,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hx,v[:,0]))/((e[0]-e[i])**2)
# Q_xy
Q[0,1]=Q[0,1]+np.dot(v[:,0].conjugate(),np.dot(Hx,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hy,v[:,0]))/((e[0]-e[i])**2)
# Q_yx
Q[1,0]=Q[1,0]+np.dot(v[:,0].conjugate(),np.dot(Hy,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hx,v[:,0]))/((e[0]-e[i])**2)
# Q_yy
Q[1,1]=Q[1,1]+np.dot(v[:,0].conjugate(),np.dot(Hy,v[:,i]))*np.dot(v[:,i].conjugate(),np.dot(Hy,v[:,0]))/((e[0]-e[i])**2)
# 5. 从Q中提取量子度规张量g和贝里曲率B
g=Q.real # g 是量子度规张量的实部 (g_ab = Re(Q_ab))
F=-2*Q.imag # 贝里曲率张量 (F_ab = -2 Im(Q_ab))
B=F[0][1] # 贝里曲率分量 (F_xy)
# 6. 计算特征椭球的参数
evals,evec=np.linalg.eigh(g) # g 的本征值 (代表半轴长的平方)
c=np.sqrt(np.abs(evals[0]-evals[1])) # 焦点距离
a=np.sqrt(np.max(evals)) # 最长半轴
e=c/a if a != 0 else 0.0 # 离心率,避免除以零
# 7. 存储结果
eccent[x,y]=e
area[x,y]=np.sqrt(np.linalg.det(g))*np.pi # 椭球面积 = pi * a * b = pi * sqrt(g1*g2) = pi * sqrt(det(g))
curva[x,y]=B*ds # 贝里曲率密度
differ[x,y]=(np.trace(g)-np.abs(B))*ds # tr(g)-|B| 密度
print("✅ K空间数据计算完成!")
# --- 绘制图像 ---
print("📊 正在生成图表...")
fig, ([ax1,ax2],[ax3,ax4]) = plt.subplots(2,2,figsize=(12,12))
# 1. 贝里曲率图 (Berry Curvature)
ax1.imshow(curva.T, origin="lower", cmap=cm.RdBu_r,
extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()])
# 计算 Chern Number (假设K空间范围是4个布里渊区)
chen = np.sum(curva) / (2 * np.pi) / 4
norm = mpl.colors.Normalize(vmin=curva.min(), vmax=curva.max())
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r) # 使用RdBu_r以使正值显示红色,负值显示蓝色
ax1.set_title(f'Berry Curvature (C = {round(chen, 3)})') # 标题显示切恩数
fig.colorbar(cmap, ax=ax1, orientation='vertical', shrink=0.945)
ax1.set_xlabel('$k_x$')
ax1.set_ylabel('$k_y$')
ax1.set_aspect('equal', adjustable='box') # 保持宽高比
# 2. tr(g)-|B| 图
ax2.imshow(differ.T, origin="lower", cmap=cm.RdBu_r,
extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()])
norm = mpl.colors.Normalize(vmin=differ.min(), vmax=differ.max())
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r)
ax2.set_title('$\\mathrm{tr}(g)-|B|$') # LaTeX格式的标题
fig.colorbar(cmap, ax=ax2, orientation='vertical', shrink=0.945)
ax2.set_xlabel('$k_x$')
ax2.set_ylabel('$k_y$')
ax2.set_aspect('equal', adjustable='box')
# 3. 特征椭球面积图 (Characteristic ellipsoid area)
ax3.imshow(area.T, origin="lower", cmap=cm.RdBu_r,
extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()])
norm = mpl.colors.Normalize(vmin=area.min(), vmax=area.max())
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r)
ax3.set_title('Characteristic Ellipsoid Area')
fig.colorbar(cmap, ax=ax3, orientation='vertical', shrink=0.945)
ax3.set_xlabel('$k_x$')
ax3.set_ylabel('$k_y$')
ax3.set_aspect('equal', adjustable='box')
# 4. 特征椭球离心率图 (Characteristic ellipsoid eccentricity index)
ax4.imshow(eccent.T, origin="lower", cmap=cm.RdBu_r,
extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()])
norm = mpl.colors.Normalize(vmin=eccent.min(), vmax=eccent.max())
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu_r)
ax4.set_title('Characteristic Ellipsoid Eccentricity Index')
fig.colorbar(cmap, ax=ax4, orientation='vertical', shrink=0.945)
ax4.set_xlabel('$k_x$')
ax4.set_ylabel('$k_y$')
ax4.set_aspect('equal', adjustable='box')
plt.tight_layout() # 自动调整子图参数,使之填充整个图像区域
plt.show()
print("🎉 图表已成功显示!")
这有什么问题吗
最新发布