じゅりあ「えっちなしじゃだめ?」 ぱいそん「えっちないとだめ」
[イントロ]
最近暖かくなってきて春が近づいてきた感じがします。人生の春は遠そうですが。
春といえば入学式や入社式など新しい生活が始まる季節です。
そんな春という季節は新しいことを始めるにはピッタリの季節とも言えます。
そこで、今回写経をはじめました。
[方法]
@cometscomeさんが公開されている「Juliaで学ぶタイトバインディング模型とトポロジカル物質」 (https://github.com/cometscome/TightBinding/blob/master/TB01.ipynb)
というJuliaでタイトバインディング模型を学べるすばらしいブログがあります。
せっかくなのでこのJuliaで書かれたプログラムをPythonを使って写経することで、
物理とプログラミングを同時に勉強することにしました。
物理初心者なので。
今回は記事の1つ目、「TB01」を写経してみました。
基本的にそのままPython化しただけなのでJuliaでの書き方は上記ブログをご参照ください。
[結果]
写経の結果は下記に示しています。つかれた(^o^)。
基本的には同じような記述ですが、以下の点で何度か躓きました。
・エルミート行列Hの固有値ε、固有関数ψを計算する場合、
Juliaであれば"ε,ψ=eig(H)"でかけますが、
Pythonの場合Numpyを使用して"ε,ψ=np.linalg.eigh(H)"と"eigh()"のように
エルミート行列専用のhがないとだめみたいです。えっちがないと。
ただし、これは1000x1000行列の場合で、100x100の場合はえっちがなくてもJuliaと
一致しました。なぜだ。。。
・固有値と固有関数を計算すると、当然JuliaとPythonで固有値は一致しますが,
固有関数は符号±だけ、一致しない場合がありました。JuliaとPythonで
行列計算の実装が微妙に違うのが原因かな?
・Juliaで"ifelse(条件式,A,B)"と書くところが、Pythonだと"A if 条件式 else B"と
書き方が違うことが学べました。
・Juliaだと多変数関数V(x,hoge)(hogeは変数)を関数f(hoge,V)のように書いて、
V(x)の形で関数の引数にできるが、Pythonの場合は、1変数関数V(x)だとJuliaと同じように
かけますが、多変数関数V(x,hoge)だとうまく行かず
def V1(x):
V1=V(x,hoge)
return V1
のように1変数に書き直して対処しました。
この辺もう少し上手く書けそうですが、今後の課題です。
[まとめ]
JuliaをPythonで写経することで、数値計算の方法とPythonの書き方を同時に学ぶ事ができました。
この調子で、残りのタイトバインディング模型の記事も写経して勉強したいと思っています。
物理初心者なのでがんばります><
[謝辞]
困ってる点に関して、Twitterで@ceptreeさん、@cometscomephysさん、@FockSpaceさんから
コメントを頂いて解決できました。ありがとうございましたm( _)m。
[参考文献]
1,「Juliaで学ぶタイトバインディング模型とトポロジカル物質」 https://github.com/cometscome/TightBinding/blob/master/TB01.ipynb
2,BloggerにJupyter Notebookを貼り付ける方法
https://qiita.com/Minerva/items/2447130411fb2faee951
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
#1.1 束縛状態を持つ1次元シュレーディンガー方程式
#1.1.1 ポテンシャルを一つ持つ場合の束縛状態
V0=-5.0
def V(x): #ポテンシャル
V=V0*np.exp(-x**2)
return V
N=1000
xmax=10.0
xmin=-10.0
xvec=np.linspace(xmin, xmax, N) #横軸を示すリスト
#ポテンシャルのプロット
plt.plot(xvec, V(xvec), label="Potential")
def make_H1d(N, xmax, xmin, V): # ポテンシャルVのときのNxN行列化
mat_H=np.zeros([N,N]) #空の行列を作る
#V(x)=-1*V0*np.exp(-1*x*x/(ξ*ξ))
a = (xmax-xmin)/(N-1)
for i in range(1,N+1):
x= (i-1)*a + xmin
for dx in range(-1, 2):
j=i+dx
v=0.0
if dx==0:
v=(2/a**2 + V(x))
elif dx==1:
v=-1/a**2
elif dx== -1 :
v=-1/a**2
if 1 <= j <= N:
mat_H[i-1,j-1]=v #ハミルトニアンの[i-1,j-1]成分
return mat_H
mat_H = make_H1d(N,xmax,xmin,V)
energy,mat_v = np.linalg.eigh(mat_H) #エルミート行列の固有値と固有関数、ここ注意
print(energy[0:10])
#ポテンシャル1個の場合の固有関数
plt.plot(xvec,mat_v[:,0])
#1.1.2 ポテンシャルを二つ持つ場合の束縛状態:遠い場合
V0=-5.0
x0=5
N=1000
xmax=10.0
xmin=-10.0
xvec=np.linspace(xmin,xmax,N)
def V2(x): #ポテンシャル2個の場合
V2=V0*(np.exp(-(x-x0)**2)+np.exp(-(x+x0)**2))
return V2
#ポテンシャル2個の場合のプロット
plt.plot(xvec, V2(xvec), label="potential")
plt.legend()
V0 = -5.0
x0 =5
N = 1000
xmax = 10.0
xmin = -10.0
mat_H2=make_H1d(N,xmax,xmin,V2)
energy2p,mat_v2=np.linalg.eigh(mat_H2)
#ポテンシャル2個の場合の固有値
print("ポテンシャル2つのときのエネルギー")
print(energy2p[0:10])
print("ポテンシャル1つのときのエネルギー")
print(energy[0:10])
##ポテンシャル2個の場合のプロット
plt.plot(xvec,mat_v2[:,0], c="blue", label="y1")
plt.plot(xvec,mat_v2[:,1], c="red", label="y2")
plt.legend()
#1.1.3 ポテンシャルを二つ持つ場合の束縛状態:近づけた場合
V0=-5.0
x0=3
N=1000
xmax=10.0
xmin=-10.0
xvec=np.linspace(xmin,xmax,N)
def V2(x):
V2=V0*(np.exp(-(x-x0)**2)+np.exp(-(x+x0)**2))
return V2
#ポテンシャル2個の場合のプロット(前より位置を近づけてる)
plt.plot(xvec,V2(xvec),label="Potential")
plt.legend()
V0 = -5.0
x0 = 3
N = 1000
xmax = 10.0
xmin = -10.0
xvec = np.linspace(xmin,xmax,N)
mat_H2 = make_H1d(N,xmax,xmin,V2)
energy2p,mat_v2 = np.linalg.eigh(mat_H2) #ポテンシャル2個の場合の固有値と固有関数
print("ポテンシャル二つの時のエネルギー")
print(energy2p[0:10])
print("ポテンシャル一つの時のエネルギー")
print(energy[0:10])
#ポテンシャル2個の場合の固有関数のプロット
plt.plot(xvec,mat_v2[:,0],c="blue")
plt.plot(xvec,mat_v2[:,1], c="red")
#1.1.3 周期的なポテンシャルの場合
def Vp(x,V0,xmax,xmin,npr): #周期ポテンシャル
dx=(xmax-xmin)/(npr-1)
V3=0.0
for i in range(1,npr):
x0=(i-1)*dx + xmin+dx/2
V3 += V0*np.exp(-(x-x0)**2)+V0*np.exp(-(x-(xmax-xmin)-x0)**2)+V0*np.exp(-(x+(xmax-xmin)-x0)**2)
return V3
N=2000
xmax=10.0
xmin=-10.0
#周期ポテンシャルのプロット
plt.plot(xvec,Vp(xvec,V0,xmax,xmin,5),label="Potential")
plt.legend()
def make_H1dp(N,xmax,xmin,V): #周期ポテンシャルの場合の行列化
mat_H = np.zeros([N,N])
a = (xmax-xmin)/(N-1)
for i in range(1,N+1):
x = (i-1)*a + xmin
for dx in range(-1,2):
j = i + dx
j += -N if j>N else 0
j += N if j<1 else 0
v = 0.0
if dx == 0:
v = (2/a**2 +V(x))
elif dx == 1:
v = -1/a**2
elif dx == -1:
v = -1/a**2
if 1 <= j <= N:
mat_H[i-1,j-1] = v
return mat_H
N = 2000
xmax = 10.0
xmin = -10.0
V0=-5.0
npr=5
xvec = np.linspace(xmin,xmax,N)
def V3(x):
V3=Vp(x,V0,xmax,xmin,npr)
return V3
mat_H2 = make_H1dp(N,xmax,xmin,V3)
energy4p,mat_v2 = np.linalg.eigh(mat_H2)
#周期ポテンシャルの場合の固有関数のプロット
print("周期ポテンシャルの時のエネルギー")
print(energy4p[0:10])
print("ポテンシャル一つの時のエネルギー")
print(energy[0:10])
plt.plot(xvec,mat_v2[:,0],label="y1")
plt.plot(xvec,mat_v2[:,1],label="y2")
plt.plot(xvec,mat_v2[:,2],label="y3")
plt.plot(xvec,mat_v2[:,3],label="y4")
plt.legend()
#1.2 束縛状態だけを考えたシュレーディンガー方程式
x0=3
N=2000
xmax=10.0
xmin=-10.0
xvec=np.linspace(xmin,xmax,N)
def V1(x):
V1=V0*np.exp(-(x-x0)**2)
return V1
mat_H1=make_H1dp(N,xmax,xmin,V1)
energy1,mat_v1=np.linalg.eigh(mat_H1)
ψ1=mat_v1[:,0]
print(energy1[0])
#ポテンシャルがV1のときの固有関数
plt.plot(xvec,ψ1)
def V2(x):
V1=V0*np.exp(-(x+x0)**2)
return V1
mat_H2=make_H1dp(N,xmax,xmin,V2)
energy2,mat_v2=np.linalg.eigh(mat_H2)
ψ2=mat_v2[:,0]
print(energy2[0])
#ポテンシャルがV2のときの固有関数
plt.plot(xvec,ψ2)
def calc_Htb(energy1,ψ1,ψ2,vecV1):
mat_Htb=np.zeros([2,2])
t=0.0
v0=0.0
N=len(ψ1)
for i in range(1,N+1):
t += (energy[0]+vecV1[i-1])*ψ1[i-1]*ψ2[i-1]
v0 += vecV1[i-1]*ψ2[i-1]*ψ2[i-1]
mat_Htb[0,0]=energy[0]+v0
mat_Htb[0,1]=t
mat_Htb[1,0]=t
mat_Htb[1,1]=energy[0]+v0
return mat_Htb
N = 2000
xmax = 10.0
xmin = -10.0
vecV1 = np.zeros(N)
dx = (xmax-xmin)/(N-1)
for i in range(1,N+1):
x = (i-1)*dx+xmin
vecV1[i-1] = V1(x)
mat_Htb = calc_Htb(energy1,ψ1,ψ2,vecV1)
print(mat_Htb)
ε,c=np.linalg.eigh(mat_Htb)
print(ε)
print(c[0:1])
V0 = -5.0
x0 = 3
xvec = np.linspace(xmin,xmax,N)
def V12(x):
V12 = V0*(np.exp(-(x-x0)**2)+np.exp(-(x+x0)**2))
return V12
mat_H2 = make_H1dp(N,xmax,xmin,V12)
energy2p,mat_v2 = np.linalg.eigh(mat_H2)
print("ポテンシャル二つの時のエネルギー")
print(energy2p[0:10])
#1.3 タイトバインディング模型
def calc_HtbModel(energy1,ψ1,ψ2,vecV1,Nx):
mat_Htb=np.zeros([Nx,Nx])
t=0.0
v0=0.0
N=len(ψ1)
for i in range(1,N+1):
t+=(energy[0]+vecV1[i-1])*ψ1[i-1]*ψ2[i-1]
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] = energy[0]
elif abs(dx) == 1:
mat_Htb[i-1,j-1] = t
return mat_Htb
N = 2000
xmax = 10.0
xmin = -10.0
vecV1 = np.zeros(N)
dx = (xmax-xmin)/(N-1)
for i in range(1,N+1):
x = (i-1)*dx+xmin
vecV1[i-1] = V1(x)
mat_Htb = calc_HtbModel(energy1,ψ1,ψ2,vecV1,4)
print(mat_Htb)
ε4,c4 = np.linalg.eigh(mat_Htb)
print(ε4)
print(c4[:,0])
print(energy4p[0:10])
コメント
コメントを投稿