學習python已經(jīng)1個多月了，最近需要畫一下石墨烯的能帶圖和態(tài)密度，接下來利用python對角化石墨烯中C原子的Hamilton量盼理，取本征值，畫一下能帶圖吧俄删。
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from numpy import *
from math import *
def Hamiltonian3(pos,a1,a2,k,t):
H=mat(zeros((2,2)))
eps=0.01
for i in range(0,2):
for j in range(0,2):
for m in range(-1,2):
for n in range(-1,2):
d=pos[i]+n*a1+m*a2-pos[j]
# 最近鄰判斷條件
if abs(linalg.norm(d)-1)<eps:
# 滿足最近鄰條件的原子間跳躍積分求和
H[i,j]=H[i,j]+t * e ** (1j*dot(k,d))
return H
# 元胞內(nèi)的兩個C原子坐標
pos1=array([0,0])
pos2=array([0,1])
pos=[pos1,pos2]
# 取晶格基矢
a1=array([sqrt(3),0])
a2=array([sqrt(3)/2,3/2])
#新建一個名叫 fig的畫圖窗口
fig = plt.figure()
# 準備繪制三維圖形
ax = Axes3D(fig)
# k點的取值范圍
Kx = arange(-pi, pi, 0.01)
Ky = arange(-pi, pi, 0.01)
n=len(Kx)
# 初始化能量E為n×n的0矩陣
E1=mat(zeros((n,n)))
t=-1
for i in range(n):
for j in range(n):
# 取k點
k=array([Kx[i],Ky[j]])
# 計算Hamilton
H=Hamiltonian3(pos,a1,a2,k,t)
# 取本征值
Es=linalg.eig(H)[0]
# 存儲本征值
E1[i,j]=Es[0]
E2=-E1
#畫圖
ax.plot_surface(Kx,Ky,E1, rstride=1, cstride=1, cmap='rainbow')
ax.plot_surface(Kx,Ky,E2, rstride=1, cstride=1, cmap='rainbow')
plt.show()
畫出來的如如下宏怔。
image.png
