補足2:じゅりあ「そのままで…いいよ?」ぱいそん「コンプレックス、あるんだ…」
In [1]:
#4. グラフェンと呼ばれるタイトバインディング模型:ジグザグ版
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
#ジグザグ型グラフェンの結晶構造
a = 0.5
#単位ベクトル
x0= [i*np.sqrt(3)*a for i in range(-3,4)]
x1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-4,4)]
y0= [2*j*1.5*a for j in range(-2,3)]
y1= [(2*j+1)*1.5*a for j in range(-3,4)]
c="red"
lw=0.5
ls="--"
for y in y0:
for x in x0:
plt.plot([x,x ],[y,y+a ],color=c, lw=lw, ls=ls)
plt.plot([x,x-np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
plt.plot([x,x+np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
for y in y1:
for x in x1:
plt.plot([x,x ],[y,y+a ],color=c, lw=lw, ls=ls)
plt.plot([x,x-np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
plt.plot([x,x+np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
plt.gca().set_aspect("equal")
plt.xlim([-3,3])
plt.ylim([-2.5,3])
Out[1]:
In [2]:
a = 0.5
#2種類の格子点を色づけてプロット
x0= [i*np.sqrt(3)*a for i in range(-3,4)]
x1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-4,4)]
y0= [2*j*1.5*a for j in range(-2,3)]
y1= [(2*j+1)*1.5*a for j in range(-3,4)]
c="red"
b="blue"
lw=0.5
ls="--"
for y in y0:
for x in x0:
plt.plot([x,x ],[y,y+a ],color=c, lw=lw, ls=ls)
plt.plot([x],[y],color=c,marker="o")
plt.plot([x-np.sqrt(3)*a/2],[y-a/2],color=b,marker="o")
plt.plot([x,x-np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw)
plt.plot([x,x+np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
for y in y1:
for x in x1:
plt.plot([x,x ],[y,y+a ],color=c, lw=lw, ls=ls)
plt.plot([x,x-np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw)
plt.plot([x],[y],color=c,marker="o")
plt.plot([x-np.sqrt(3)*a/2],[y-a/2],color=b, marker="o")
plt.plot([x,x+np.sqrt(3)*a/2],[y,y-a/2],color=c, lw=lw, ls=ls)
plt.gca().set_aspect("equal")
plt.xlim([-3,3])
plt.ylim([-2.5,3])
Out[2]:
In [3]:
#グラフェンのハミルトニアン
def calc_HGraphene(Nx,Ny,μ):
N=Nx*Ny*2
mat_Htb = np.zeros([N,N])
mat_Htb += (-μ)*np.eye(N)
t=1.0
for ix in range(1,Nx+1):
for iy in range(1,Ny+1):
for dx in range(-1,2):
for dy in range(-1,2):
jx=ix+dx
jx+= -Nx if jx>Nx else 0
jx+= Nx if jx<1 else 0
jy=iy+dy
jy+= -Ny if jy>Ny else 0
jy+= Ny if jy<1 else 0
for a in range(1,3):
b=2 if a==1 else 1
ii=((iy-1)*Nx+ix-1)*2+a
jj=((jy-1)*Nx+jx-1)*2+b
if dx==0 and dy==0:
mat_Htb[ii-1,jj-1]=t
elif dx==+1 and dy==0 and a==1:
mat_Htb[ii-1,jj-1]=t
elif dx==0 and dy==1 and a==1:
mat_Htb[ii-1,jj-1]=t
elif dx==-1 and dy==0 and a==2:
mat_Htb[ii-1,jj-1]=t
elif dx==0 and dy==-1 and a==2:
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [8]:
#エネルギー固有値のヒストグラム
Nx=40
Ny=40
μ=0
mat_H=calc_HGraphene(Nx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
plt.hist(energy,bins=100,label="y1",ec="black")
Out[8]:
In [4]:
#x方向だけフーリエ変換したハミルトニアン
def calc_HGraphenekx(kx,Ny,μ):
N=Ny*2
mat_Htb=np.zeros([N,N], dtype=np.complex)#複素数を含むN×N行列はおまじないが必要
mat_Htb+=(-μ)*np.eye(N)#np.eye()でN×N単位行列
t=1.0
for iy in range(1,Ny+1):
for dy in range(-1,2):
jy = iy + dy
jy+= -Ny if jy>Ny else 0
jy+= Ny if jy<1 else 0
for a in range(-1,3):
b=2 if a==1 else 1
ii=(iy-1)*2+a
jj=(jy-1)*2+b
if dy==0 and a==1:
mat_Htb[ii-1,jj-1]=t+t*np.exp(1j*kx)#複素数は"1j"みたいに数字とj
elif dy==1 and a==1:
mat_Htb[ii-1,jj-1]=t
elif dy==0 and a==2:
mat_Htb[ii-1,jj-1]=t+t*np.exp(-1j*kx)
elif dy==-1 and a==2:
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [ ]:
#バンド構造のプロット
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*2])
cnt=0
for kx in vkx:
cnt+=1
mat_H=calc_HGraphenekx(kx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)#複素数でもnp.linalg.eigh()で固有値、固有ベクトル
for i in range(1,Ny*2+1):
#print(energy[i-1])
ep[cnt-1,i-1]=energy[i-1]
plt.plot(vkx,ep)
In [5]:
#xy方向ともにフーリエ変換、ky方向の長さを調整したハミルトニアン
def calc_HGraphenekxky(kx,ky,μ):
N=2
mat_Htb=np.zeros([N,N], dtype=np.complex)
mat_Htb+=(-μ)*np.eye(N)
t=1.0
c=np.sqrt(3)/2
# for iy in range(1,Ny+1):
# for dy in range(-1,2):
# jy = iy + dy
# jy+= -Ny if jy>Ny else 0
# jy+= Ny if jy<1 else 0
for a in range(-1,3):
b=2 if a==1 else 1
ii=a
jj=b
if a==1:
mat_Htb[ii-1,jj-1]=t+t*np.exp(1j*kx)+t*np.exp(1j*ky*c)
# elif dy==1 and a==1:
# mat_Htb[ii-1,jj-1]=t
elif a==2:
mat_Htb[ii-1,jj-1]=t+t*np.exp(-1j*kx)+t*np.exp(-1j*ky*c)
# elif dy==-1 and a==2:
# mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [7]:
#対角要素がないハミルトニアンであることがわかる。
μ=0.0
calc_HGraphenekxky(1.0,2.2,μ)
Out[7]:
In [8]:
#固有値が解析的に求まるのでプロット。
def f(kx,ky):
f=np.sqrt(3 + 2*np.cos(kx)+2*np.cos(ky*np.sqrt(3)/2)+2*np.cos(kx-ky*np.sqrt(3)/2))
return f
x=np.array([i*2*np.pi/100 for i in range(-100,101)])
y=np.array([i*2*np.pi/100 for i in range(-100,101)])
z=[[f(i,j) for i in x] for j in y]
plt.contour(x,y,z)
plt.gca().set_aspect("equal")
In [9]:
#xy軸の60度の傾きを調整して再計算。
def f(kx,ky):
f=np.sqrt(3 + 2*np.cos(kx+ky/2)+2*np.cos(ky)+2*np.cos(kx-ky/2))
return f
x=np.array([i*2*np.pi/100 for i in range(-100,101)])
y=np.array([i*2*np.pi/100 for i in range(-100,101)])
z=[[f(i,j) for i in x] for j in y]
plt.contour(x,y,z)
plt.gca().set_aspect("equal")
In [10]:
#周期的境界条件を外した片側フーリエ変換ハミルトニアン。端状態を示す。
def calc_HGraphenekx_w(kx,Ny,μ):
N=Ny*2
mat_Htb=np.zeros([N,N], dtype=np.complex)
mat_Htb+=(-μ)*np.eye(N)
t=1.0
for iy in range(1,Ny+1):
for dy in range(-1,2):
jy = iy + dy
#jy+= -Ny if jy>Ny else 0 #周期的境界条件を外す
#jy+= Ny if jy<1 else 0
for a in range(-1,3):
b=2 if a==1 else 1
ii=(iy-1)*2+a
jj=(jy-1)*2+b
if 1<=jy<=Ny:
if dy==0 and a==1:
mat_Htb[ii-1,jj-1]=t+t*np.exp(1j*kx)
elif dy==1 and a==1:
mat_Htb[ii-1,jj-1]=t
elif dy==0 and a==2:
mat_Htb[ii-1,jj-1]=t+t*np.exp(-1j*kx)
elif dy==-1 and a==2:
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [ ]:
#端状態を含むバンド構造のプロット
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*2])
cnt=0
for kx in vkx:
cnt+=1
mat_H=calc_HGraphenekx_w(kx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
for i in range(1,Ny*2+1):
#print(energy[i-1])
ep[cnt-1,i-1]=energy[i-1]
plt.plot(vkx,ep)
In [72]:
#kx=π付近の固有関数のプロット。
kx=np.pi-0.1
Ny=50
mat_H=calc_HGraphenekx_w(kx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
yv=[]
ψ1=[]
ψ2=[]
for i in range(1,Ny+1): #各リストに固有関数の各成分を収納
yv.append(i-1)
ψ1.append(mat_v[2*i-1,Ny-1])
ψ2.append(mat_v[2*i-2,Ny])
#for i in range(1,Ny+1):
# np.append(yv,i-1)
# np.append(ψ1,mat_v[2*(i-1),Ny])
# np.append(ψ2,mat_v[2*(i-1)-1,Ny+1])
plt.plot(np.array(yv),np.array(np.real(ψ1)), label="1:Real")
plt.plot(np.array(yv),np.array(np.imag(ψ1)), label="1:Imag")
plt.plot(np.array(yv),np.array(np.real(ψ2)), label="2:Real")
plt.plot(np.array(yv),np.array(np.imag(ψ2)), label="2:Imag")
plt.legend()
#plt.ylim([-0.6,0.6])
Out[72]:
コメント
コメントを投稿