補足１：じゅりあ「そのままで…いいよ？」ぱいそん「コンプレックス、あるんだ…」
In [15]:
#3. タイトバインディング模型とフーリエ変換
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
a=1.0
m=10
#格子点のプロット
vecx=[i for i in range(0,m+1)] #リスト内包表現、0からmまでのリスト
lattice=[0.0 for i in range(0,m+1)]
plt.plot(vecx,lattice, marker="o", label="y1")
plt.legend()
Out[15]:
In [3]:
ns=100
vecx1=[i/ns for i in range(0,ns*m+1)]
k=2*2*(np.pi)/m
wave1=np.cos(k*np.array(vecx1))
#np.array()の中に入れることで各成分ごとのCos成分リストができる
plt.plot(vecx,lattice, marker="o", label="y1")
plt.plot(vecx1,wave1, label="y2")
plt.legend()
Out[3]:
In [4]:
k += 2*np.pi
wave2=np.cos(k*np.array(vecx1))
plt.plot(vecx,lattice, marker="o", label="y1")
plt.plot(vecx1,wave1, label="y2")
plt.plot(vecx1,wave2)
c='red'
lw=0.5
ls="--"
for x in vecx:
plt.plot([x,x],[-1,1],color=c, ls=ls, linewidth = 0.5)
In [5]:
#ハミルトニアン
def calc_HtbModel(Nx,μ):
mat_Htb = np.zeros([Nx,Nx])
t = 1.0
for i in range(1,Nx+1):
for dx in range(-1,2):
j = i + dx
j += -Nx if j > Nx else 0 #周期境界条件
j += Nx if j < 1 else 0
if dx ==0:
mat_Htb[i-1,j-1] = -μ
elif abs(dx) == 1:
mat_Htb[i-1,j-1] = -t
return mat_Htb
In [6]:
μ=0.0
Nx=100
mat_H=calc_HtbModel(Nx,μ)
energy,mat_v=np.linalg.eigh(mat_H)
#固有値のヒストグラム、状態密度を反映
plt.hist(energy,bins=10,label="y1",ec="black")#ec="black"無しだとヒストグラムに線が出ない
plt.legend()
Out[6]:
In [7]:
μ=0.0
Nx=2000
mat_H=calc_HtbModel(Nx,μ)
energy,mat_v=np.linalg.eigh(mat_H)
#±2の状態密度を見るために粒子数を増やしたヒストグラム
plt.hist(energy,bins=100,label="y1",ec="black")
plt.legend()
Out[7]:
In [8]:
#2次元正方格子のハミルトニアン、最近接ホッピングのみを考える。
def calc_HtbModel2D(Nx,Ny,μ):
N = Nx * Ny
mat_Htb=np.zeros([N,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
ii=(iy-1)*Nx+ix
jj=(jy-1)*Nx+jx
if dx ==0 and dy==0:
mat_Htb[ii-1,jj-1] = -μ
elif abs(dx)==1 and dy==0:
mat_Htb[ii-1,jj-1] = -t
elif abs(dy)==1 and dx==0:
mat_Htb[ii-1,jj-1] = -t
return mat_Htb
In [9]:
μ=0.0
Nx=50
Ny=50
mat_H=calc_HtbModel2D(Nx,Ny,μ)
energy,mat_v=np.linalg.eigh(mat_H)
#2次元正方格子の状態密度
plt.hist(energy,bins=50,label="y1",ec="black")
plt.legend()
Out[9]:
In [10]:
#波数空間でのエネルギーのプロット。
def f(kx,ky):
f=-2*(np.cos(kx)+np.cos(ky))
return f
x=np.array([i*np.pi/100 for i in range(-100,101)])
y=np.array([i*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)#コンター図はxもyもリストじゃなくて配列にしておく必要あり。
plt.gca().set_aspect("equal")#アスペクト比と色使いの設定
In [11]:
#波数空間の勾配のプロット
def grad(kx,ky):
g = np.sqrt((4*(np.sin(kx)**2 + np.sin(ky)**2)))
return g
x=np.array([i*np.pi/100 for i in range(-100,101)])
y=np.array([i*np.pi/100 for i in range(-100,101)])
z=[[grad(i,j) for i in x] for j in y]
plt.contour(x,y,z)
plt.gca().set_aspect("equal")
In [12]:
#片側だけフーリエ変換
def calc_HtbModelX(Nx,ky,μ):
mat_Htb = np.zeros([Nx,Nx])
t = 1.0
for ix in range(1,Nx+1):
for dx in range(-1,2):
jx = ix + dx
jx += -Nx if jx>Nx else 0
jx += Nx if jx < 1 else 0
if dx == 0:
mat_Htb[ix-1,jx-1] = -μ -2*t*np.cos(ky)
elif abs(dx) == 1:
mat_Htb[ix-1,jx-1] = -t
return mat_Htb
In [14]:
#エネルギーのky依存性
μ=0.0
Nx = 100
nky = 100
vky = np.linspace(-np.pi,np.pi,nky)
#print(vky)
ep = np.zeros([nky,Nx])
cnt = 0
for ky in vky:
cnt += 1
mat_H = calc_HtbModelX(Nx,ky,μ)
energy,mat_v = np.linalg.eigh(mat_H)
for i in range(1,Nx+1):
#println(energy[i])
ep[cnt-1,i-1] = energy[i-1]
#バンド構造のプロット
plt.plot(vky,ep)
コメント
コメントを投稿