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):
h=np.zeros((len(orb),len(orb)),dtype='complex')
for i in range (len(hop)):
x=np.array(hop[i][1])
y=np.array(hop[i][2])
latr=np.dot((np.array(hop[i][3:5])+np.array(orb[y])-np.array(orb[x])),lat)
newterm=hop[i][0]*np.exp(1j*(np.dot(latr,k)))
h[x,y]=h[x,y]+newterm
h=h+np.conjugate(h).transpose()
return h
def hamx(lat,orb,hop,k):
h=np.zeros((len(orb),len(orb)),dtype='complex')
for i in range (len(hop)):
x=np.array(hop[i][1])
y=np.array(hop[i][2])
latr=np.dot((np.array(hop[i][3:5])+np.array(orb[y])-np.array(orb[x])),lat)
newterm=hop[i][0]*np.exp(1j*(np.dot(latr,k)))*1j*latr[0]
h[x,y]=h[x,y]+newterm
h=h+np.conjugate(h).transpose()
return h
def hamy(lat,orb,hop,k):
h=np.zeros((len(orb),len(orb)),dtype='complex')
for i in range (len(hop)):
x=np.array(hop[i][1])
y=np.array(hop[i][2])
latr=np.dot((np.array(hop[i][3:5])+np.array(orb[y])-np.array(orb[x])),lat)
newterm=hop[i][0]*np.exp(1j*(np.dot(latr,k)))*1j*latr[1]
h[x,y]=h[x,y]+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
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=[0]*len(orb)
nx=100
ny=100
kxlim=2*np.pi
kylim=2*np.pi
ds=(4*kxlim*kylim)/(nx)/ny
origin=[0,0]
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))
kxlist=np.zeros(nx+1)
kylist=np.zeros(ny+1)
print("开始计算K空间数据...")
for x in range(nx+1):
kx=(x-nx/2)/(nx/2)*kxlim+origin[0]
kxlist[x]=kx
for y in range(ny+1):
ky=(y-ny/2)/(ny/2)*kylim+origin[1]
kylist[y]=ky
H=ham(lat,orb,hop,[kx,ky])
for i in range(len(H)):
H[i,i]= H[i,i]+onsite[i]
Hx=hamx(lat,orb,hop,[kx,ky])
Hy=hamy(lat,orb,hop,[kx,ky])
e,v=np.linalg.eigh(H)
Q=np.zeros((2,2),dtype='complex')
for i in range(1,len(orb)):
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[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[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[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)
# print(Hy)
# print(e)
# print(v)
g=Q.real
evals,evec=np.linalg.eigh(g)
c=np.sqrt(np.abs(evals[0]-evals[1]))
if evals[0]>evals[1]:
a=np.sqrt(evals[0])
else:
a=np.sqrt(evals[1])
e=c/a
F=-2*Q.imag
B=F[0][1]
eccent[x,y]=e
area[x,y]=np.sqrt(np.linalg.det(g))*np.pi
curva[x,y]=B*ds
differ[x,y]=(np.trace(g)-np.abs(B))*ds
fig, (ax1,ax2) = plt.subplots(1,2,figsize=(12,6))
ax1.imshow(curva.T, origin="lower", cmap=cm.RdBu_r,
extent=[kxlist.min(), kxlist.max(), kylist.min(), kylist.max()])
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)
ax1.set_title(str(round(chen,2))+'Berry Curvature')
fig.colorbar(cmap,ax=ax1,orientation='vertical',shrink=0.945)
ax2.imshow(differ.T,origin="lower",cmap=cm.RdBu_r)
norm = mpl.colors.Normalize(vmin=differ.min(), vmax=differ.max())
cmap = mpl.cm.ScalarMappable(norm=norm, cmap=mpl.cm.RdBu)
ax2.set_title('tr(g)-|B|')
fig.colorbar(cmap,ax=ax2,orientation='vertical',shrink=0.945)
plt.tight_layout()
plt.show()
我这么写有什么问题latr=np.dot((np.array(hop[i][3:5])+np.array(orb[y])-np.array(orb[x])),lat)这么写对不对