JuliaでTightBinding模型を作ってみた

"Julia - Julia のパッケージ作り"
http://goropikari.hatenablog.com/entry/julia_make_package
を参考にして、パッケージを作ってみました。

"TightBinding.jl"
https://github.com/cometscome/TightBinding.jl
強束縛模型を作って、バンド図などを描いたりするパッケージです。

サンプルとして、グラフェン、グラフェンナノリボン、鉄系超伝導体有効模型、を作ってみました。

#バージョン情報
Julia 1.0.0

#インストール方法
インストールは、]でパッケージモードにしてから

```julia
add https://github.com/cometscome/TightBinding.jl
```
でできます。
これを使う場合には、

```julia
using TightBinding
```
とします。ここで、

```bash
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本のベクトルでよくて、

```julia
a1 = [sqrt(3)/2,1/2]
a2= [0,1]
```
です。これを、格子としてセットします。

```julia
la = set_Lattice(2,[a1,a2])
```
### 格子の追加
次に、この基本並進ベクトルで囲まれた平行四辺形の中に原子を置きます。グラフェンの場合、平行四辺形の中には

```julia
add_atoms!(la,[1/3,1/3])
add_atoms!(la,[2/3,2/3])
```
の二つがあります。

### ホッピングの追加
ここまでで、平行四辺形タイルを無限に敷き詰めることで、グラフェンの原子の位置がきまったことになります。次に、電子がどう飛び移るかを指定します。
実際に設定する前に、どのような飛び移りがありえるか、みてみましょう。

```julia
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
```
の三つです。ですので、この三つをホッピングとして加えます。

```julia
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と等価なので、追加しません。
ここまで作った結晶の構造を見るには、

```julia
pls = plot_lattice_2d(la)
using Plots
savefig("structure.png")
```
とすればよく、出力結果は
![structure.png](https://qiita-image-store.s3.amazonaws.com/0/246113/74f30ecf-5a7c-d376-925d-eaf5c4647cbc.png)
となります。
ちゃんと蜂の巣格子になっているのがわかります。

### 状態密度のプロット
この模型の電子状態の状態密度を計算するには、

```julia
pdos = plot_DOS(la, 100)
savefig("dos.png")
```
とします。ここで、100はk空間のメッシュ数となっており、100であれば、100x100にk空間を区切っています。
得られた状態密度の図は
![dos.png](https://qiita-image-store.s3.amazonaws.com/0/246113/9465d2c1-d4f3-3647-606e-86bbc14b150b.png)
となります。

# バンド図のプロット
次に、バンド図のプロットをします。バンド図とは、波数空間におけるある線上でのエネルギー分散をプロットしてもので、結晶の電子状態をみるときによく使われます。
## グラフェンの場合
グラフェンは六角格子ですので、波数空間におけるブリルアンゾーンも六角形になっています。
ブリルアンゾーン中の$\Gamma$点-$K$点のライン、$K$点から$M$点のライン、そして、$M$点から$\Gamma$点のライン上でのエネルギー固有値を見てみましょう。

### ラインの設定
プロットするラインを設定します。まず、ラインを初期化：

```julia
klines = set_Klines()
```
して、3本のラインを追加します。

```julia
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にはプロットしたいラインの情報が入りました。バンド図をプロットするには、

```julia
pls2 = calc_band_plot(klines,la)
```
とします。ファイルを保存するには、

```julia
savefig("band.png")
```
とすればよくて、

![band.png](https://qiita-image-store.s3.amazonaws.com/0/246113/2ae0e839-c29d-31f3-3e33-1641d7d1cb08.png)

という図が得られます。$M$点に"ディラックコーン"があるのが見てとれます。

# 有限系の計算
これまでは、波数が良い量子数でした。次は、一つの方向に壁あるいは周期境界条件を課してみましょう。
## グラフェンの場合
片方の方向を実空間で有限サイズにすることにより、波数は一つだけになります。
その時のバンド図を描いてみましょう。
まず、k点の設定を

```julia
klines = set_Klines()
kmin = [-π]
kmax = [π]
add_Kpoints!(klines,kmin,kmax,"-pi","pi")
```
としました。

モデルはすでに作ってあるので、あとはバンドをプロットするだけです。
どちらの方向を有限サイズにするか、ということをdirectionで指定して、

```julia
direction = 1
pls3= calc_band_plot_finite(klines,la,direction,periodic=true)
savefig("PBC.png")
```
とすると、有限系のバンド図を描くことができます。ここでは、periodic=trueなので、周期系です。
得られたバンド図は

![PBC.png](https://qiita-image-store.s3.amazonaws.com/0/246113/f203ecec-9852-091a-3c28-6bc4c8a85ff1.png)

となります。バンドが複数あるのは、有限系になったためその方向の波数が飛び飛びの値しか取れなくなったためです。

periodic=falseとすると

```julia
pls4= calc_band_plot_finite(klines,la,direction,periodic=false)
savefig("OBC.png")
```
壁を入れた系のバンド図を描くことができて、図は

![OBC.png](https://qiita-image-store.s3.amazonaws.com/0/246113/610861bc-1e80-43dc-00d2-05d27f8e845e.png)
となります。ゼロエネルギー付近にフラットなバンドが見えていると思います。これがグラフェンにおける端状態です。

# 鉄系超伝導体の有効模型
もう一つサンプルとして、鉄系超伝導体の模型を作ってみましょう。
## 2軌道有効模型
一番簡単な例として、dxz軌道とdyz軌道からなる2軌道有効模型[S. Raghu et al. S. Rachu et al. Phys. Rev. B 77, 220503(R) (2008)]
を作ってみます。
ハミルトニアンは

```math
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)
```
となっており、

```math
\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つ同じ場所にあるとして、

```julia
la = set_Lattice(2,[[1,0],[0,1]]) #正方格子
add_atoms!(la,[0,0]) #dxz軌道
add_atoms!(la,[0,0]) #dyz軌道
```
とします。

ホッピングは、

```julia
#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,μ) 
```
となります。
この模型のバンドをプロットするには、

```julia
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)
```
とすればよくて、得られるバンドは

![Fe.png](https://qiita-image-store.s3.amazonaws.com/0/246113/ac41fc7d-a1a1-b00e-1f1b-e6d8fd739efa.png)
となり、論文の図とよく一致していることがわかります。

なお、ハミルトニアンが欲しい時には、

```julia
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)]を作ってみましょう。
この論文にあるテーブルに従ってホッピングを設定すればよいです。
コードは、

```julia
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
```
となります。
これでハミルトニアンが作れたので、バンドは、

```julia
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")
```
で描くことができます。
その結果、

![Fe5band.png](https://qiita-image-store.s3.amazonaws.com/0/246113/c5b666bd-fb44-1515-a855-5abc24511cef.png)
というバンドが得られます。
このバンド図はPRL論文の図と少し違いますが、これは、PRL論文がtableにない微小なホッピングを考慮しているためであることが知られていて、実際、tableから作った場合のバンド図は[T. Nomura, J. Phys. Soc. Jpn. 78, 034716 (2009)]のFig.2にありまして、これと完全に一致していることがわかります。

### フェルミ面の描画
フェルミ面を描く関数を追加しました。
これは、

```
pls = plot_fermisurface_2D(la)
```
でプロットすることができて、
鉄系超伝導体の場合には、

![FeFermi.png](https://qiita-image-store.s3.amazonaws.com/0/246113/515d7634-f9f0-a7d2-f5c4-f368f6c714ff.png)

となります。