じゅりあ「そのままで…いいよ？」ぱいそん「コンプレックス、あるんだ…」
[イントロ]
時間は残酷です。気づいたらもう一年経ってしまいました。
明日やろう、明日やろうと先延ばししたことは結局できず、日々をこなすだけで時間だけが過ぎてしまいました。
もう少し計画的に時間を使っていきたいです。明日から頑張ろう。
とはいえ少しずつ前に進んでいる事もあったりします。
例えば写経とか。
[方法]
前回のブログに引き続き、
(http://buhin-blog.blogspot.jp/2018/03/blog-post_18.html)
@cometscomeさんが公開されている「Juliaで学ぶタイトバインディング模型とトポロジカル物質」 (https://github.com/cometscome/TightBinding/blob/master/TB01.ipynb)
のプログラムをPythonで写経しました。Numpyの使い方に少し詳しくなれました。
今回は「TB02」～「TB05」までの内容を写経しました。
TB02： この記事
TB03： http://buhin-blog.blogspot.jp/2018/03/blog-post_64.html
TB04： http://buhin-blog.blogspot.jp/2018/03/blog-post_16.html
TB05： http://buhin-blog.blogspot.jp/2018/03/blog-post_71.html
ついにグラフェンのバンド構造を描くことができました。すごい。
今回躓いた点は下記のとおりです。
・Pythonの場合、「A if B else C」 で「BならばA、それ以外ならC」の計算を実行
・例えば、「vecx=[i for i in range(0,m+1)] 」はリスト内包表現で0からmまでのリストを作れる。
・例えば、「wave1=np.cos(np.array(vecx1))」のようにnp.array()の中に入れることで、
各成分ごとにCos化した配列ができる。リストのままだとリストが帰ってくる。
・ヒストグラムを書く際は、「ec="black"」無しだとヒストグラムを区切る線が出ない。なんでやねん。
・「plt.contour(x,y,z)」でコンター図を書く際はxもyもリストじゃなくて、配列にしておく必要あり。
・複素数を含む行列は「A=np.zeros([N,N], dtype=np.complex)」のようにおまじないが必要。
Juliaだと要らないみたいでなくてもそのままかける。
Pythonはコンプレックスがあるみたい。かわいい。
・「A=np.eye(N)」でN×N単位行列が作れる。
・複素数でもnp.linalg.eigh()で固有値、固有ベクトルが出せる。
引き続き頑張るぞい。
[参考文献]
１，「Juliaで学ぶタイトバインディング模型とトポロジカル物質」 https://github.com/cometscome/TightBinding
#2. タイトバインディング模型と固体中の電子
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
def calc_HtbModel(Nx):
mat_Htb=np.zeros([Nx,Nx])#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 # A if B else C で「BならばA、それ以外ならC」の計算
if dx==0:
mat_Htb[i-1,j-1]==0.0
elif abs(dx)==1:
mat_Htb[i-1,j-1]=-t
return mat_Htb
#固有値の計算
Nx=100
mat_H=calc_HtbModel(Nx)
energy,mat_v=np.linalg.eigh(mat_H) #固有値と固有ベクトルの計算
plt.plot(energy, label="y1")
plt.legend()
#ハーフフィリングのときの一番大きなエネルギー(ゼロになる)
print(energy[int(Nx/2-1)])
#化学ポテンシャルを入れたタイトバインディングモデルハミルトニアン
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
Nx=100
μ=-1.5
mat_H=calc_HtbModel(Nx,μ)
energy, mat_v=np.linalg.eigh(mat_H)
count=0
#付のエネルギーをもつ状態の個数
for ene in energy:
if ene>0:
break
else:
count+=1
print("Filling =", count/Nx)
コメント
コメントを投稿