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

tight-binding4
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]:
(-2.5, 3)
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]:
(-2.5, 3)
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]:
(array([  31.,   24.,   30.,   24.,   30.,   24.,   36.,   12.,   42.,
          30.,   30.,   36.,   18.,   30.,   36.,   30.,   36.,   18.,
          42.,   36.,   42.,   36.,   18.,   54.,   36.,   54.,   36.,
          18.,   66.,   36.,   54.,   66.,   30.,  123.,   42.,   72.,
          18.,   24.,   42.,   18.,   24.,   24.,   12.,   24.,   12.,
          12.,    0.,    6.,    6.,    0.,    0.,    6.,    6.,    0.,
          12.,   12.,   24.,   12.,   24.,   24.,   18.,   42.,   24.,
          18.,   72.,   42.,  123.,   30.,   66.,   54.,   36.,   66.,
          18.,   36.,   54.,   36.,   54.,   18.,   36.,   42.,   36.,
          42.,   18.,   36.,   30.,   36.,   30.,   18.,   36.,   30.,
          30.,   42.,   12.,   36.,   24.,   30.,   24.,   30.,   24.,   31.]),
 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 [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]:
array([[ 0.00000000+0.j        ,  1.21204358+1.78605884j],
       [ 1.21204358-1.78605884j,  0.00000000+0.j        ]])
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]:
<matplotlib.legend.Legend at 0x29c1a99c470>

コメント

このブログの人気の投稿

有機物量子スピン液体の熱伝導論争の流れを振り返る

【改題】室温超伝導ふたたび!~大丈夫じゃなかった、Natureの論文だもん!~

2023年7月の気になった論文(完全版)