"Julia - Julia のパッケージ作り"
http://goropikari.hatenablog.com/entry/julia_make_package
を参考にして、パッケージを作ってみました。
"TightBinding.jl"
https://github.com/cometscome/TightBinding.jl
強束縛模型を作って、バンド図などを描いたりするパッケージです。
サンプルとして、グラフェン、グラフェンナノリボン、鉄系超伝導体有効模型、を作ってみました。
バージョン情報
Julia 1.0.0
インストール方法
インストールは、]でパッケージモードにしてから
add https://github.com/cometscome/TightBinding.jl
でできます。
これを使う場合には、
using TightBinding
とします。ここで、
WARNING: eval from module Plots to ...
みたいなエラーメッセージが出てしまっていますが、パッケージ自体は問題なく使えます。エラーを出さなくする方法は今のところわかりません。調査中です。
模型の構築
強束縛模型（TightBinding模型）では、原子を並べて、その原子を電子が飛び移る量子力学的な模型です。飛び移り方を設定することで、実際の固体中の電子の挙動を調べることができます。
量子力学については
"Juliaで学ぶ量子力学"
https://github.com/cometscome/QM
強束縛模型については
"In Japanese. Juliaで学ぶタイトバインディング模型とトポロジカル物質"
https://github.com/cometscome/TightBinding
これらが役に立つと思います。
グラフェンの場合
炭素が蜂の巣状に並びそれが積み重なったものがグラファイトですが、蜂の巣状のシートが1枚だけのときをグラフェン（あるいはグラフィーン）と呼び、原子層1層の物質ということで盛んに研究がなされています。作成方法ですが、グラファイトに粘着テープをくっつけて、それをシリコンにはりつけて剥がすとできる
http://www.c.u-tokyo.ac.jp/info/about/booklet-gazette/bulletin/535/open/B-2-2.html
ということが知られています。
その強束縛模型は
http://www.f-denshi.com/000okite/500metal/nanotube01.html
こちらを参考にするとわかりやすいでしょう。
基本並進ベクトル
まず、グラフェンの基本並進ベクトルを設定します。2次元なので2本のベクトルでよくて、
a1 = [sqrt(3)/2,1/2]
a2= [0,1]
です。これを、格子としてセットします。
la = set_Lattice(2,[a1,a2])
格子の追加
次に、この基本並進ベクトルで囲まれた平行四辺形の中に原子を置きます。グラフェンの場合、平行四辺形の中には
add_atoms!(la,[1/3,1/3])
add_atoms!(la,[2/3,2/3])
の二つがあります。
ホッピングの追加
ここまでで、平行四辺形タイルを無限に敷き詰めることで、グラフェンの原子の位置がきまったことになります。次に、電子がどう飛び移るかを指定します。
実際に設定する前に、どのような飛び移りがありえるか、みてみましょう。
show_neighbors(la)
を実行すると、
Possible hoppings
(1,1), x:-1//1, y:-1//1
(1,2), x:-2//3, y:-2//3
(2,2), x:-1//1, y:-1//1
(1,1), x:-1//1, y:0//1
(1,2), x:-2//3, y:1//3
(2,2), x:-1//1, y:0//1
(1,1), x:-1//1, y:1//1
(1,2), x:-2//3, y:4//3
(2,2), x:-1//1, y:1//1
(1,1), x:0//1, y:-1//1
(1,2), x:1//3, y:-2//3
(2,2), x:0//1, y:-1//1
(1,1), x:0//1, y:0//1
(1,2), x:1//3, y:1//3
(2,2), x:0//1, y:0//1
(1,1), x:0//1, y:1//1
(1,2), x:1//3, y:4//3
(2,2), x:0//1, y:1//1
(1,1), x:1//1, y:-1//1
(1,2), x:4//3, y:-2//3
(2,2), x:1//1, y:-1//1
(1,1), x:1//1, y:0//1
(1,2), x:4//3, y:1//3
(2,2), x:1//1, y:0//1
(1,1), x:1//1, y:1//1
(2,2), x:1//1, y:1//1
という出力が得られます。ここでの1は原子1番目、2は原子2番目をあらわし、x,yは原子を結ぶベクトルにどのようなものがあるかを記したものです。一番単純なグラフェンの模型の場合には、電子の飛び移りは原子1から原子2に飛び移るものだけです。つまり、(1,2)を見ればよくて、その中で距離が短そうなものは、
(1,2), x:1//3, y:1//3
(1,2), x:-2//3, y:1//3
(1,2), x:1//3, y:-2//3
の三つです。ですので、この三つをホッピングとして加えます。
t = 1.0
add_hoppings!(la,-t,1,2,[1/3,1/3])
add_hoppings!(la,-t,1,2,[-2/3,1/3])
add_hoppings!(la,-t,1,2,[1/3,-2/3])
2から1への飛び移りもありますが、この飛び移りは1から2と等価なので、追加しません。
ここまで作った結晶の構造を見るには、
pls = plot_lattice_2d(la)
using Plots
savefig("structure.png")
とすればよく、出力結果は
となります。
ちゃんと蜂の巣格子になっているのがわかります。
状態密度のプロット
この模型の電子状態の状態密度を計算するには、
pdos = plot_DOS(la, 100)
savefig("dos.png")
とします。ここで、100はk空間のメッシュ数となっており、100であれば、100x100にk空間を区切っています。
得られた状態密度の図は
となります。
バンド図のプロット
次に、バンド図のプロットをします。バンド図とは、波数空間におけるある線上でのエネルギー分散をプロットしてもので、結晶の電子状態をみるときによく使われます。
グラフェンの場合
グラフェンは六角格子ですので、波数空間におけるブリルアンゾーンも六角形になっています。
ブリルアンゾーン中の$\Gamma$点-$K$点のライン、$K$点から$M$点のライン、そして、$M$点から$\Gamma$点のライン上でのエネルギー固有値を見てみましょう。
ラインの設定
プロットするラインを設定します。まず、ラインを初期化：
klines = set_Klines()
して、3本のラインを追加します。
kmin = [0,0]
kmax = [2π/sqrt(3),0]
add_Kpoints!(klines,kmin,kmax,"G","K")
kmin = [2π/sqrt(3),0]
kmax = [2π/sqrt(3),2π/3]
add_Kpoints!(klines,kmin,kmax,"K","M")
kmin = [2π/sqrt(3),2π/3]
kmax = [0,0]
add_Kpoints!(klines,kmin,kmax,"M","G")
これで、klinesにはプロットしたいラインの情報が入りました。バンド図をプロットするには、
pls2 = calc_band_plot(klines,la)
とします。ファイルを保存するには、
savefig("band.png")
とすればよくて、
という図が得られます。$M$点に"ディラックコーン"があるのが見てとれます。
有限系の計算
これまでは、波数が良い量子数でした。次は、一つの方向に壁あるいは周期境界条件を課してみましょう。
グラフェンの場合
片方の方向を実空間で有限サイズにすることにより、波数は一つだけになります。
その時のバンド図を描いてみましょう。
まず、k点の設定を
klines = set_Klines()
kmin = [-π]
kmax = [π]
add_Kpoints!(klines,kmin,kmax,"-pi","pi")
としました。
モデルはすでに作ってあるので、あとはバンドをプロットするだけです。
どちらの方向を有限サイズにするか、ということをdirectionで指定して、
direction = 1
pls3= calc_band_plot_finite(klines,la,direction,periodic=true)
savefig("PBC.png")
とすると、有限系のバンド図を描くことができます。ここでは、periodic=trueなので、周期系です。
得られたバンド図は
となります。バンドが複数あるのは、有限系になったためその方向の波数が飛び飛びの値しか取れなくなったためです。
periodic=falseとすると
pls4= calc_band_plot_finite(klines,la,direction,periodic=false)
savefig("OBC.png")
壁を入れた系のバンド図を描くことができて、図は
となります。ゼロエネルギー付近にフラットなバンドが見えていると思います。これがグラフェンにおける端状態です。
鉄系超伝導体の有効模型
もう一つサンプルとして、鉄系超伝導体の模型を作ってみましょう。
2軌道有効模型
一番簡単な例として、dxz軌道とdyz軌道からなる2軌道有効模型[S. Raghu et al. S. Rachu et al. Phys. Rev. B 77, 220503(R) (2008)]
を作ってみます。
ハミルトニアンは
H =
\sum_k
\left(
\begin{matrix}
c^{\dagger}_{1k} & c^{\dagger}_{2k}
\end{matrix}
\right)
\left( \begin{matrix}
\epsilon_x(k)-\mu & \epsilon_{xy}(k) \\
\epsilon_{xy}(k) & \epsilon_y(k)-\mu
\end{matrix} \right)
\left(
\begin{matrix}
c_{1k} \\
c_{2k}
\end{matrix}
\right)
となっており、
\epsilon_x(k) = -2 t_1 \cos k_x - 2 t_2 \cos k_y - 4 t_3 \cos k_x \cos k_y \\
\epsilon_y(k) = -2 t_2 \cos k_x - 2 t_1 \cos k_y - 4 t_3 \cos k_x \cos k_y\\
\epsilon_{xy}(k) = -4 t_4 \sin k_x \sin k_y
です。$\cos k_x = (\exp (i k_x) + \exp (- ik_x))/2$を使えば、どのようなホッピングがあるかはわかります。$\epsilon_x(k)$は、第一項が[1,0]方向、第二項が[0,1]方向、第三項が[1,1]と[1,-1]方向のホッピングです。
$\epsilon_{xy}(k) $は[1,1]方向には$-t_4$、[1,-1]方向には$t_4$でホッピングしています。
2軌道模型なので原子は2つ同じ場所にあるとして、
la = set_Lattice(2,[[1,0],[0,1]]) #正方格子
add_atoms!(la,[0,0]) #dxz軌道
add_atoms!(la,[0,0]) #dyz軌道
とします。
ホッピングは、
#hoppings
t1 = -1.0
t2 = 1.3
t3 = -0.85
t4 = t3
μ = 1.45
#dxz軌道
add_hoppings!(la,-t1,1,1,[1,0])
add_hoppings!(la,-t2,1,1,[0,1])
add_hoppings!(la,-t3,1,1,[1,1])
add_hoppings!(la,-t3,1,1,[1,-1])
#dyz軌道
add_hoppings!(la,-t2,2,2,[1,0])
add_hoppings!(la,-t1,2,2,[0,1])
add_hoppings!(la,-t3,2,2,[1,1])
add_hoppings!(la,-t3,2,2,[1,-1])
#dxzとdyzの飛び移り
add_hoppings!(la,-t4,1,2,[1,1])
add_hoppings!(la,-t4,1,2,[-1,-1])
add_hoppings!(la,t4,1,2,[1,-1])
add_hoppings!(la,t4,1,2,[-1,1])
#化学ポテンシャル
set_μ!(la,μ)
となります。
この模型のバンドをプロットするには、
klines = set_Klines()
kmin = [0,0]
kmax = [π,0]
add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)")
kmin = [π,0]
kmax = [π,π]
add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)")
kmin = [π,π]
kmax = [0,0]
add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(0,0)")
pls = calc_band_plot(klines,la)
とすればよくて、得られるバンドは
となり、論文の図とよく一致していることがわかります。
なお、ハミルトニアンが欲しい時には、
ham = hamiltonian_k(la) #ハミルトニアンを返す関数を生成
kx = 0.1
ky = 0.2
hamk = ham([kx,ky]) #hamはk=[kx,ky]の関数
println(hamk)
とすればよく、以後はhamをkの関数として呼び出せば好きなように使うことができます。
5軌道有効模型
最後に、鉄系超伝導体の5軌道模型[K. Kuroki et al., Phys. Rev. Lett. 101, 087004 (2008)]を作ってみましょう。
この論文にあるテーブルに従ってホッピングを設定すればよいです。
コードは、
la = set_Lattice(2,[[1,0],[0,1]])
add_atoms!(la,[0,0])
add_atoms!(la,[0,0])
add_atoms!(la,[0,0])
add_atoms!(la,[0,0])
add_atoms!(la,[0,0])
tmat = [
-0.7 0 -0.4 0.2 -0.1
-0.8 0 0 0 0
0.8 -1.5 0 0 -0.3
0 1.7 0 0 -0.1
-3.0 0 0 -0.2 0
-2.1 1.5 0 0 0
1.3 0 0.2 -0.2 0
1.7 0 0 0.2 0
-2.5 1.4 0 0 0
-2.1 3.3 0 -0.3 0.7
1.7 0.2 0 0.2 0
2.5 0 0 0.3 0
1.6 1.2 -0.3 -0.3 -0.3
0 0 0 -0.1 0
3.1 -0.7 -0.2 0 0
]
tmat = 0.1.*tmat
imap = zeros(Int64,5,5)
count = 0
for μ=1:5
for ν=μ:5
count += 1
imap[μ,ν] = count
end
end
Is = [1,-1,-1,1,1,1,1,-1,-1,1,-1,-1,1,1,1]
σds = [1,-1,1,1,-1,1,-1,-1,1,1,1,-1,1,-1,1]
tmat_σy = tmat[:,:]
tmat_σy[imap[1,2],:] = -tmat[imap[1,3],:]
tmat_σy[imap[1,3],:] = -tmat[imap[1,2],:]
tmat_σy[imap[1,4],:] = -tmat[imap[1,4],:]
tmat_σy[imap[2,2],:] = tmat[imap[3,3],:]
tmat_σy[imap[2,4],:] = tmat[imap[3,4],:]
tmat_σy[imap[2,5],:] = -tmat[imap[3,5],:]
tmat_σy[imap[3,3],:] = tmat[imap[2,2],:]
tmat_σy[imap[3,4],:] = tmat[imap[2,4],:]
tmat_σy[imap[3,5],:] = -tmat[imap[2,5],:]
tmat_σy[imap[4,5],:] = -tmat[imap[4,5],:]
hoppingmatrix = zeros(Float64,5,5,5,5)
hops = [-2,-1,0,1,2]
hopelements = [[1,0],[1,1],[2,0],[2,1],[2,2]]
for μ = 1:5
for ν=μ:5
for ii=1:5
ihop = hopelements[ii][1]
jhop = hopelements[ii][2]
#[a,b],[a,-b],[-a,-b],[-a,b],[b,a],[b,-a],[-b,a],[-b,-a]
#[a,b]
i = ihop +3
j = jhop +3
hoppingmatrix[μ,ν,i,j]=tmat[imap[μ,ν],ii]
#[a,-b] = σy*[a,b] [1,1] -> [1,-1]
if jhop != 0
i = ihop +3
j = -jhop +3
hoppingmatrix[μ,ν,i,j]=tmat_σy[imap[μ,ν],ii]
end
if μ != ν
#[-a,-b] = I*[a,b] [1,1] -> [-1,-1],[1,0]->[-1,0]
i = -ihop +3
j = -jhop +3
hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
#[-a,b] = I*[a,-b] = I*σy*[a,b] #[2,0]->[-2,0]
if jhop != 0
i = -ihop +3
j = jhop +3
hoppingmatrix[μ,ν,i,j]=Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
end
end
#[b,a],[b,-a],[-b,a],[-b,-a]
if jhop != ihop
#[b,a] = σd*[a,b]
i = jhop +3
j = ihop +3
hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat[imap[μ,ν],ii]
#[-b,a] = σd*σy*[a,b]
if jhop != 0
i = -jhop +3
j = ihop +3
hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
end
if μ != ν
#[-b,-a] = σd*[-a,-b] = σd*I*[a,b]
i = -jhop +3
j = -ihop +3
hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat[imap[μ,ν],ii]
#[b,-a] = σd*[-a,b] = σd*I*[a,-b] = σd*I*σy*[a,b] #[2,0]->[-2,0]
if jhop != 0
i = jhop +3
j = -ihop +3
hoppingmatrix[μ,ν,i,j]=σds[imap[μ,ν]]*Is[imap[μ,ν]]*tmat_σy[imap[μ,ν],ii]
end
end
end
end
end
end
for μ=1:5
for ν=μ:5
for i = 1:5
ih = hops[i]
for j = 1:5
jh = hops[j]
if hoppingmatrix[μ,ν,i,j] != 0.0
add_hoppings!(la,hoppingmatrix[μ,ν,i,j],μ,ν,[ih,jh])
end
end
end
end
end
onsite = [10.75,10.96,10.96,11.12,10.62]
set_onsite!(la,onsite)
set_μ!(la,10.96) #set the chemical potential
となります。
これでハミルトニアンが作れたので、バンドは、
nk = 100
klines = set_Klines()
kmin = [0,0]
kmax = [π,0]
add_Kpoints!(klines,kmin,kmax,"(0,0)","(pi,0)",nk=nk)
kmin = [π,0]
kmax = [π,π]
add_Kpoints!(klines,kmin,kmax,"(pi,0)","(pi,pi)",nk=nk)
kmin = [π,π]
kmax = [0,0]
add_Kpoints!(klines,kmin,kmax,"(pi,pi)","(0,0)",nk=nk)
pls = calc_band_plot(klines,la)
savefig("Fe5band.png")
で描くことができます。
その結果、
というバンドが得られます。
このバンド図はPRL論文の図と少し違いますが、これは、PRL論文がtableにない微小なホッピングを考慮しているためであることが知られていて、実際、tableから作った場合のバンド図は[T. Nomura, J. Phys. Soc. Jpn. 78, 034716 (2009)]のFig.2にありまして、これと完全に一致していることがわかります。
フェルミ面の描画
フェルミ面を描く関数を追加しました。
これは、
pls = plot_fermisurface_2D(la)
でプロットすることができて、
鉄系超伝導体の場合には、
となります。