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

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

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 

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

bands.in

&bands
  prefix = 'SrVO3',
  outdir = './tmp/',
  filband = 'bands.dat',
/

$ 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などを使用して描きます。

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

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

SrVO3.scf.inをコピーしてSrVO3.nscf.inという名前にして、次のように編集してください。

• calculation='nscf'と変更
• occupationsを削除
• K_POINTS以下を削除

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

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

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

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

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

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

また、次の変数を指定します。

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

$ wannier90.x -pp SrVO3

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

SrVO3. pw2wan.in

 &INPUTPP
   outdir    =  './tmp/',
   prefix    =  'SrVO3',
   seedname  =  'SrVO3',
   write_unk = .true.
 /

ワニエ関数をプロットする場合（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を元に行うのが便利です。

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

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

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

とします。

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