補足３：じゅりあ「そのままで…いいよ？」ぱいそん「コンプレックス、あるんだ…」
In [1]:
#5. グラフェンと呼ばれるタイトバインディング模型：アームチェア版
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
#アームチェア型グラフェンの結晶構造
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-3,4)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-4,4)]
x0= [2*j*1.5*a for j in range(-2,3)]
x1= [(2*j+1)*1.5*a for j in range(-2,2)]
c="red"
lw=0.5
ls="--"
for y in y0:
for x in x0:
plt.plot([x,x+a ],[y,y ],color=c, lw=lw, ls=ls)
plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
for y in y1:
for x in x1:
plt.plot([x,x+a ],[y,y ],color=c, lw=lw, ls=ls)
plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*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]:
#2種類の原子の色分け
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-2,3)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-2,3)]
x0= [2*j*1.5*a for j in range(-1,3)]
x1= [(2*j+1)*1.5*a for j in range(-1,2)]
c="red"
b="blue"
lw=0.5
ls="--"
for y in y0:
for x in x0:
plt.plot([x,x+a ],[y,y ],color=c, lw=lw)
plt.plot([x],[y],color=b,marker="o")
plt.plot([x+a],[y],color=c,marker="o")
plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
for y in y1:
for x in x1:
plt.plot([x,x+a ],[y,y ],color=c, lw=lw)
plt.plot([x],[y],color=b,marker="o")
plt.plot([x+a],[y],color=c, marker="o")
plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw, ls=ls)
plt.gca().set_aspect("equal")
plt.xlim([-1.8,3])
plt.ylim([-2.2,2])
Out[2]:
In [3]:
#ユニットセルのプロット
a = 0.5
y0= [i*np.sqrt(3)*a for i in range(-2,3)]
y1= [i*np.sqrt(3)*a + np.sqrt(3)*a/2 for i in range(-2,3)]
x0= [2*j*1.5*a for j in range(-1,3)]
x1= [(2*j+1)*1.5*a for j in range(-1,2)]
c="red"
b="blue"
lw=0.5
ls="--"
x=0
y=0
plt.plot([x,x+a ],[y,y ],color=c, lw=lw)
plt.plot([x],[y],color=b,marker="o")
plt.plot([x+a],[y],color=c,marker="o")
plt.plot([x,x-a/2],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x,x-a/2],[y,y+np.sqrt(3)*a/2],color=c, lw=lw)
plt.plot([x-a/2],[y-np.sqrt(3)*a/2],color=c,marker="o")
plt.plot([x-a/2],[y+np.sqrt(3)*a/2],color=c,marker="o")
plt.plot([x+a,x-a/2+2*a ],[y,y-np.sqrt(3)*a/2],color=c, lw=lw,ls=ls)
plt.plot([x+a,x-a/2+2*a ],[y,y+np.sqrt(3)*a/2],color=c, lw=lw)
plt.plot([x-a/2+2*a],[y-np.sqrt(3)*a/2],color=b,marker="o")
plt.plot([x-a/2+2*a],[y+np.sqrt(3)*a/2],color=b,marker="o")
plt.gca().set_aspect("equal")
plt.xlim([-0.5,1])
plt.ylim([-1,1])
Out[3]:
In [4]:
#アームチェア型グラフェンのハミルトニアン
def calc_HGrapheneA(Nx,Ny,μ):
N=Nx*Ny*4
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,5):
for b in range(1,5):
hop=False
ii=((iy-1)*Nx+ix-1)*4+a
if a==1:
if dx==-1 and dy==0 and b==4:
hop = True
elif dx==0 and dy==1 and b==2:
hop=True
elif dx==0 and dy==0 and b==2:
hop=True
elif a==2:
if dx==0 and dy==0 and b==1:
hop = True
elif dx==0 and dy==0 and b==3:
hop=True
elif dx==0 and dy==-1 and b==1:
hop=True
elif a==3:
if dx==0 and dy==0 and b==2:
hop = True
elif dx==0 and dy==0 and b==4:
hop=True
elif dx==0 and dy==-1 and b==4:
hop=True
elif a==4:
if dx==0 and dy==0 and b==3:
hop = True
elif dx==1 and dy==0 and b==1:
hop=True
elif dx==0 and dy==1 and b==3:
hop=True
if hop and 1<=jx<=Nx and 1<=jy<=Ny:
jj=((jy-1)*Nx+jx-1)*4+b
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [5]:
#アームチェア型グラフェンのエネルギー固有値
Nx=20
Ny=40
μ=0
mat_H=calc_HGrapheneA(Nx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
plt.hist(energy,bins=100,label="y1",ec="black")
Out[5]:
In [6]:
#片側フーリエ変換したハミルトニアン
def calc_HGrapheneAkx(kx,Ny,μ):
N=Ny*4
mat_Htb=np.zeros([N,N],dtype=np.complex )#複素数の行列
mat_Htb += (-μ)*np.eye(N)
t0=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
t=0.0
for a in range(1,5):
for b in range(1,5):
hop=False
t=0.0
ii=(iy-1)*4+a
if a==1:
if dy==0 and b==4:
hop = True
t+=t0*np.exp(-1j*kx)
elif dy==1 and b==2:
hop=True
t+=t0
elif dy==0 and b==2:
hop=True
t+=t0
elif a==2:
if dy==0 and b==1:
hop = True
t+=t0
elif dy==0 and b==3:
hop=True
t+=t0
elif dy==-1 and b==1:
hop=True
t+=t0
elif a==3:
if dy==0 and b==2:
hop = True
t+=t0
elif dy==0 and b==4:
hop=True
t+=t0
elif dy==-1 and b==4:
hop=True
t+=t0
elif a==4:
if dy==0 and b==3:
hop = True
t+=t0
elif dy==0 and b==1:
hop=True
t+=t0*np.exp(1j*kx)
elif dy==1 and b==3:
hop=True
t+=t0
if hop and 1<=jy<=Ny:
jj=(jy-1)*4+b
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [7]:
#kxを横軸としてエネルギー固有値をプロット
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*4])
cnt=0
for kx in vkx:
cnt+=1
mat_H=calc_HGrapheneAkx(kx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
for i in range(1,Ny*4+1):
#print(energy[i-1])
ep[cnt-1,i-1]=energy[i-1]
plt.plot(vkx,ep)
In [8]:
#端状態を出すためにy方向の周期条件を外したハミルトニアン
def calc_HGrapheneAkxw(kx,Ny,μ):
N=Ny*4
mat_Htb=np.zeros([N,N],dtype=np.complex )
mat_Htb += (-μ)*np.eye(N)
t0=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
t=0.0
for a in range(1,5):
for b in range(1,5):
hop=False
t=0.0
ii=(iy-1)*4+a
if a==1:
if dy==0 and b==4:
hop = True
t+=t0*np.exp(-1j*kx)
elif dy==1 and b==2:
hop=True
t+=t0
elif dy==0 and b==2:
hop=True
t+=t0
elif a==2:
if dy==0 and b==1:
hop = True
t+=t0
elif dy==0 and b==3:
hop=True
t+=t0
elif dy==-1 and b==1:
hop=True
t+=t0
elif a==3:
if dy==0 and b==2:
hop = True
t+=t0
elif dy==0 and b==4:
hop=True
t+=t0
elif dy==-1 and b==4:
hop=True
t+=t0
elif a==4:
if dy==0 and b==3:
hop = True
t+=t0
elif dy==0 and b==1:
hop=True
t+=t0*np.exp(1j*kx)
elif dy==1 and b==3:
hop=True
t+=t0
if hop and 1<=jy<=Ny:
jj=(jy-1)*4+b
mat_Htb[ii-1,jj-1]=t
return mat_Htb
In [9]:
#バンド構造。端状態は出ない。
μ=0.0
Ny=50
nkx=100
vkx=np.linspace(-np.pi,np.pi,nkx)
ep=np.zeros([nkx,Ny*4])
cnt=0
for kx in vkx:
cnt+=1
mat_H=calc_HGrapheneAkx(kx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
for i in range(1,Ny*4+1):
#print(energy[i-1])
ep[cnt-1,i-1]=energy[i-1]
plt.plot(vkx,ep)
コメント
コメントを投稿