補足3:じゅりあ「そのままで…いいよ?」ぱいそん「コンプレックス、あるんだ…」

tight-binding5
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]:
(-2.5, 3)
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]:
(-2.2, 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]:
(-1, 1)
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]:
(array([  27.,   24.,   28.,   32.,   28.,   24.,   36.,   14.,   38.,
          34.,   26.,   30.,   34.,   26.,   32.,   36.,   24.,   34.,
          34.,   36.,   44.,   32.,   30.,   38.,   40.,   48.,   32.,
          42.,   46.,   44.,   60.,   50.,   50.,  105.,   52.,   56.,
          30.,   32.,   32.,   26.,   24.,   12.,   24.,   16.,   12.,
          10.,    2.,   10.,    2.,    2.,    2.,    2.,   10.,    2.,
          10.,   12.,   16.,   24.,   12.,   24.,   26.,   32.,   32.,
          30.,   56.,   52.,  105.,   50.,   50.,   60.,   44.,   46.,
          42.,   32.,   48.,   40.,   38.,   30.,   32.,   44.,   36.,
          34.,   34.,   24.,   36.,   32.,   26.,   34.,   30.,   26.,
          34.,   38.,   14.,   36.,   24.,   28.,   32.,   28.,   24.,   27.]),
 array([-3.  , -2.94, -2.88, -2.82, -2.76, -2.7 , -2.64, -2.58, -2.52,
        -2.46, -2.4 , -2.34, -2.28, -2.22, -2.16, -2.1 , -2.04, -1.98,
        -1.92, -1.86, -1.8 , -1.74, -1.68, -1.62, -1.56, -1.5 , -1.44,
        -1.38, -1.32, -1.26, -1.2 , -1.14, -1.08, -1.02, -0.96, -0.9 ,
        -0.84, -0.78, -0.72, -0.66, -0.6 , -0.54, -0.48, -0.42, -0.36,
        -0.3 , -0.24, -0.18, -0.12, -0.06,  0.  ,  0.06,  0.12,  0.18,
         0.24,  0.3 ,  0.36,  0.42,  0.48,  0.54,  0.6 ,  0.66,  0.72,
         0.78,  0.84,  0.9 ,  0.96,  1.02,  1.08,  1.14,  1.2 ,  1.26,
         1.32,  1.38,  1.44,  1.5 ,  1.56,  1.62,  1.68,  1.74,  1.8 ,
         1.86,  1.92,  1.98,  2.04,  2.1 ,  2.16,  2.22,  2.28,  2.34,
         2.4 ,  2.46,  2.52,  2.58,  2.64,  2.7 ,  2.76,  2.82,  2.88,
         2.94,  3.  ]),
 <a list of 100 Patch objects>)
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)

コメント

このブログの人気の投稿

学振採用者はどこへ消えた?

物理系研究関係者、ツイッターやりすぎランキング(ぶひん調べ)

オレ達はあと何本論文を書けば東大教授になれるんだ?