====== 概要 ======

  * 遷移金属酸化物SrVO3を例に、d電子系の強束縛模型を作成します。
  * フェルミ面付近の3本のバンドはt2g軌道 (xy,yz,zx) に由来します。
  * 結晶構造は以下のサイトのものを使いました。
    * http://www.crystallography.net/cod/1512117.html


====== バンド構造 ======

まずはバンド計算を行い、強束縛模型で使う軌道がどのエネルギーにあるかを把握します。

<file - SrVO3.scf.in>
&CONTROL
  calculation = 'scf',
  wf_collect  = .true.,
  prefix      = 'SrVO3',
  verbosity   = 'high'
  outdir      = './tmp/',
  pseudo_dir  = './pseudo/',
/
&SYSTEM
  ibrav       = 1
  A           = 3.8409
  nat         = 5
  ntyp        = 3
  ecutwfc     = 100
  occupations = "tetrahedra_opt"
/
&ELECTRONS
  conv_thr    =  1d-8,
/
ATOMIC_SPECIES 
  Sr   87.6200    Sr.pbesol-spn-kjpaw_psl.1.0.0.UPF 
  V    50.9415     V.pbesol-spnl-kjpaw_psl.1.0.0.UPF 
  O    15.9994     O.pbesol-n-kjpaw_psl.1.0.0.UPF 
ATOMIC_POSITIONS alat
  Sr   0.000000000000000   0.000000000000000   0.000000000000000
  V    0.500000000000000   0.500000000000000   0.500000000000000
  O    0.000000000000000   0.500000000000000   0.500000000000000
  O    0.500000000000000   0.000000000000000   0.500000000000000
  O    0.500000000000000   0.500000000000000   0.000000000000000
K_POINTS automatic
  6 6 6 0 0 0 
</file>

<file - SrVO3.bands.in>
&CONTROL
  calculation = 'bands',
  prefix      = 'SrVO3',
  wf_collect  = .true.,
  verbosity   = 'high'
  outdir      = './tmp/',
  pseudo_dir  = './pseudo/',
/
&SYSTEM
  ibrav       = 1
  A           = 3.8409
  nat         = 5
  ntyp        = 3
  ecutwfc     = 100
  occupations = "tetrahedra_opt"
/
&ELECTRONS
  conv_thr    =  1d-8,
/
ATOMIC_SPECIES 
  Sr   87.6200    Sr.pbesol-spn-kjpaw_psl.1.0.0.UPF 
  V    50.9415     V.pbesol-spnl-kjpaw_psl.1.0.0.UPF 
  O    15.9994     O.pbesol-n-kjpaw_psl.1.0.0.UPF 
ATOMIC_POSITIONS alat
  Sr   0.000000000000000   0.000000000000000   0.000000000000000
  V    0.500000000000000   0.500000000000000   0.500000000000000
  O    0.000000000000000   0.500000000000000   0.500000000000000
  O    0.500000000000000   0.000000000000000   0.500000000000000
  O    0.500000000000000   0.500000000000000   0.000000000000000
K_POINTS crystal_b
6
  0.5  0.5  0.0  20
  0.0  0.0  0.0  20
  0.5  0.0  0.0  20
  0.5  0.5  0.0  20
  0.5  0.5  0.5  20
  0.0  0.0  0.0  0
</file>

<file - bands.in>
&bands
  prefix = 'SrVO3',
  outdir = './tmp/',
  filband = 'bands.dat',
/
</file>

  $ mpirun -n 8 pw.x < SrVO3.scf.in > SrVO3.scf.out
  $ mpirun -n 8 pw.x < SrVO3.bands.in > SrVO3.bands.out
  $ mpirun -n 8 bands.x < bands.in > bands.out

バンドはplotband.xなどを使用して描きます。

{{:quantumespresso:wannier90:srvo3_bands.png?400|}}

フェルミ準位 (12.4501 eV) 近傍に3本のt2g軌道に由来するバンドが、その上に2本のeg軌道に由来するバンドがあることがわかります。
下の方にある9本のバンドは酸素のp軌道です（酸素原子は単位胞に3つ含まれています）。

====== ワニエ化 ======


===== NSCF計算 =====

最局在ワニエ関数を使ったフィッティングを行うk点を指定します。
今回は4*4*4で区切ったk点を使いましょう。
kmesh.plは /(QEをインストールしたディレクトリ)/wannier90-2.1.0/utility に入っていますので、それをコピーしてきて使用します。

SrVO3.scf.inをコピーしてSrVO3.nscf.inという名前にして、次のように編集してください。
  * calculation='nscf'と変更
  * occupationsを削除
  * K_POINTS以下を削除

<file - SrVO3.nscf.in>
&CONTROL
  calculation = 'nscf',
  wf_collect  = .true.,
  prefix      = 'SrVO3',
  verbosity   = 'high'
  outdir      = './tmp/',
  pseudo_dir  = './pseudo/',
/
&SYSTEM
  ibrav       = 1
  A           = 3.8409
  nat         = 5
  ntyp        = 3
  ecutwfc     = 100
/
&ELECTRONS
  conv_thr    =  1d-8,
/
ATOMIC_SPECIES 
  Sr   87.6200    Sr.pbesol-spn-kjpaw_psl.1.0.0.UPF 
  V    50.9415     V.pbesol-spnl-kjpaw_psl.1.0.0.UPF 
  O    15.9994     O.pbesol-n-kjpaw_psl.1.0.0.UPF 
ATOMIC_POSITIONS crystal
  Sr   0.000000000000000   0.000000000000000   0.000000000000000
  V    0.500000000000000   0.500000000000000   0.500000000000000
  O    0.000000000000000   0.500000000000000   0.500000000000000
  O    0.500000000000000   0.000000000000000   0.500000000000000
  O    0.500000000000000   0.500000000000000   0.000000000000000
</file>

  $ kmesh.pl 4 4 4 >> SrVO3.nscf.in
  $ mpirun -n 8 pw.x < SrVO3.nscf.in > SrVO3.nscf.out

これで準備が完了しました。


===== Wannier90による計算 =====

次にWannier90を使って、実際に強束縛模型を作ります。
入力ファイル名は//seedname//.winの形式にします（seednameとprefixは異なっていても構いません）。

<file - SrVO3.win>
 num_wann = 3
 num_bands = 25
 mp_grid = 4 4 4 
 dis_win_min = 11.0
 dis_win_max = 14.5
 dis_froz_min = 11.0
 dis_froz_max = 13.5
 write_hr     = .true. 
 bands_plot   = .true.
 bands_num_points = 100
 wannier_plot = .true.

Begin Unit_Cell_Cart
Ang
   3.8409     0.0000     0.0000
   0.0000     3.8409     0.0000
   0.0000     0.0000     3.8409 
End Unit_Cell_Cart

Begin Atoms_Cart
Ang
Sr      0.00000   0.00000   0.00000
V       1.92045   1.92045   1.92045
O       0.00000   1.92045   1.92045
O       1.92045   0.00000   1.92045
O       1.92045   1.92045   0.00000
End Atoms_Cart

Begin Kpoint_path
R  0.5  0.5  0.0    G  0.0  0.0  0.0
G  0.0  0.0  0.0    X  0.5  0.0  0.0
X  0.5  0.0  0.0    R  0.5  0.5  0.0
R  0.5  0.5  0.0    M  0.5  0.5  0.5
M  0.5  0.5  0.5    G  0.0  0.0  0.0
End Kpoint_path

Begin Projections
 V:l=2,mr=2,3,5
End Projections

Begin Kpoints
0.00000000	0.00000000	0.00000000
0.00000000	0.00000000	0.25000000
（中略）
0.75000000	0.75000000	0.75000000
End Kpoints
</file>

変数の意味は次のとおりです。

^変数^説明^
|num_wann|強束縛模型で取り出す軌道の数。今回はt2g軌道の3つの軌道（xy,yz,zx）です。|
|num_bands|密度汎関数理論計算で使ったすべての軌道の数。SrVO3.nscf.out内のnumber of Kohn-Sham statesがそれに対応します。|
|mp_grid|k点のメッシュ。kmesh.plで指定したものをそのまま入力します。| 
|dis_win_min|ワニエ化するバンドが含まれているエネルギーの下限。|
|dis_win_max|ワニエ化するバンドが含まれているエネルギーの上限。|
|dis_froz_min|ワニエ化するバンドが含まれ、かつ他のバンドが含まれないエネルギーの下限。|
|dis_froz_max|ワニエ化するバンドが含まれ、かつ他のバンドが含まれないエネルギーの上限。|
|write_hr|飛び移り積分を出力するかどうか。.true.にするとSrVO3_hr.datに出力されます| 
|bands_plot|ワニエ関数によりできたハミルトニアンのバンドを出力するかどうか。|
|bands_num_points|バンドを出力する場合、そのk点数|
|wannier_plot|最局在ワニエ関数を出力するかどうか。|

また、次の変数を指定します。
^変数^説明^
|Unit_Cell_Cart|デカルト座標による結晶の基本並進ベクトル。SrVO3.nscf.outのcrystal axesを見ながら入力します。単位はオングストローム(Ang)あるいはボーア半径(Bohr)のどちらかを指定します。|
|Atoms_Cart|デカルト座標による結晶の原子位置。こちらもSrVO3.nscf.outを見ながら入力します。単位はオングストローム(Ang)あるいはボーア半径(Bohr)のどちらかを指定します。結晶の基本並進ベクトルの分率座標を使いたければ<nowiki>ATOMS_FRAC</nowiki>を使用します。|
|Kpoint_path|ワニエ化後に出力するバンドの経路。bands_plot=.true.としたときに使用します。逆格子ベクトルの分率座標で指定します。|
|Projections|どの軌道でフィッテイングするか指定します。今回はt2g軌道の3つの軌道（xz,yz,xy）です（指定の仕方はWannier90のマニュアルを参照してください）。|
|Kpoints|ワニエ化に使うブロッホ関数のk点。SrVO3.nscf.inのものをそのまま使用します。|

できたらWannier90を実行し、計算に必要なファイルを生成します。
指定するのは//seedname//のみです。
MPIが使えないことに注意します。

  $ wannier90.x -pp SrVO3

次に、pw.xで生成したファイルをWanner90が読み込める形式に変換します。

<file - SrVO3. pw2wan.in>
 &INPUTPP
   outdir    =  './tmp/',
   prefix    =  'SrVO3',
   seedname  =  'SrVO3',
   write_unk = .true.
 /
</file>

ワニエ関数をプロットする場合（SrVO3.winでwannier_plot = .true.とした場合）は write_unk=.true. とする必要があります。
ただし計算時間が非常にかかります。
今回の例でも普通のパソコンだと半日から丸一日くらいかかると思いますので注意してください。

  $ mpirun -n 8 pw2wannier90.x < SrVO3.pw2wan.in > SrVO3.pw2wan.out

これで準備が完了したので、実際にワニエ化を行います。

  $ wannier90.x SrVO3


強束縛模型によるバンドはSrVO3_band.datに、飛び移り積分はSrVO3_hr.datに出力されます。

バンドのプロットは、GNUPLOTを使ってSrVO3_band.gnuを元に行うのが便利です。

{{:quantumespresso:wannier90:srvo3_bands_wan.png?400|}}

紫の線がDFTによるバンド、緑の点がワニエフィットにより得られたバンドです。
bands.xで作成したバンドの出力ファイルは波数が2π/aを単位としているので、横軸のとり方に注意してください。


====== SrVO3_hr.dat ======

SrVO3_hr.datには、ファイルの最初にコメントやワニエ軌道の数、飛び移り先のサイトの数などが出力されたあと、実際の飛び移り積分の値が並んでいます。

単位胞の位置を$\mathbf{r}$、単位胞内の軌道のインデックスを$\alpha$, $\beta$として、強束縛模型のハミルトニアンが
\begin{equation}
 H = \sum_{\mathbf{r}, \mathbf{R}} t_{\alpha,\beta} (\mathbf{R}) c^\dagger_{\mathbf{r}, \alpha} c_{\mathbf{r} + \mathbf{R}, \beta}
\end{equation}
と書けるとします。
ここで$\mathbf{R}$は結晶の基本並進ベクトルを使って
\begin{equation}
 \mathbf{R} = m_1 \mathbf{a}_1 + m_2 \mathbf{a}_2 + m_3 \mathbf{a}_3
\end{equation}
（$m_1$, $m_2$, $m_3$は整数）です。
$m_1=m_2=m_3=0$かつ$\alpha=\beta$のときはオンサイトのエネルギーを表します。

SrVO3_hr.dat内の飛び移り積分の読み方は、例えば

  1   -2   -1    3    2   -0.001905    0.000000

の場合は$m_1=1$, $m_2=-2$, $m_3=-1$, $\beta=3$, $\alpha=2$, で、$t_{\alpha,\beta}(\mathbf{R})=-0.001905+i0.000000$という意味です（$i$は虚数単位）。

軌道の種類とインデックスの関係は、ワニエ関数をプロットすることで確かめてください。


====== 最局在ワニエ関数のプロット ======

//seedname//.win で wannier_plot = .true. とした場合は最局在ワニエ関数が //seedname//_//n//.xsf に出力されます（//n//は軌道の指標）。
プロットするにはXCrysDenを使用します。

例えば上で計算したSrVO3の軌道1のワニエ軌道が見たい場合は

  $ xcrysden --xsf SrVO3_00001.xsf

  - [Modify] -> [Atomic Radius] -> SpaceFill/Ball factor を0.5に設定 -> [OK]
  - [Tools] -> [Data Grid] -> [Isosurface] -> "Isovalue"を1に設定、"Render +/- isovalue"にチェック -> [Submit]

とします。

{{:quantumespresso:wannier90:srvo3_00001.png?400|}}

赤が正、青が負の部分のワニエ軌道です。
l=2, mr=2の軌道（dxz軌道）の形になっていますが、酸素のp軌道が少し混成して酸素サイトにもワニエ軌道が広がっています。
これを防ぐには酸素軌道も含めてワニエ化しなければなりません。