Wannier90でLCAOの跳び移り積分を計算

PWscfとWannier90というソフトを使って、LACO法(Linear Combination of Atomic Orbitals method)で用いられる跳び移り積分の計算を行う。
LCAO法は強束縛法(Tight-binding method)と等価である。

電子のスピンは無視して考える。
サイトiを中心とした軌道\alphaの電子のワニエ型波動関数\phi_{\alpha}({\bf r}-{\bf R}_i)とする。
このときのハミルトニアン\mathcal{H}の行列要素は


t_{i,j}^{\alpha,\beta} = \int d^3 r \, \phi_\alpha^{\ast}({\bf r}-{\bf R}_i) \mathcal{H} \phi_\beta({\bf r}-{\bf R}_j)
で表される。
すなわち、第二量子化の形でハミルトニアンを書くと

\mathcal{H} = \sum_{i,j} \sum_{\alpha,\beta} t_{ij}^{\alpha,\beta} c_{i,\alpha}^\dagger c_{j,\beta}
となる。
ここでc_{i,\alpha}^\daggerはサイトiにある軌道\alphaの電子の生成演算子である。

LCAO法についてはHarrisonのテキストに詳しく、一読の価値がある。

Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond (Dover Books on Physics)

Electronic Structure and the Properties of Solids: The Physics of the Chemical Bond (Dover Books on Physics)

Elementary Electronic Structure

Elementary Electronic Structure


Wannier90のインストール

Wannier90を使うにはPWscfが必要なので、先にインストールしてください。
Ubuntu11.10 Quantum ESPRESSOのインストール - cmphysの日記
以下では、上の手順に従ってインストールしたものとします。
また、Intel Fortran CompilerとMKLはすでにインストールされているものとします。

次のサイトから Wannier90-1.2.tar.gz をダウンロードします。
http://www.wannier.org/download.html
ダウンロードしたファイルをホーム・フォルダにあるespresso-5.0の中に入れて展開し、展開したフォルダ名を'wannier90'とでもしておきます。
フォルダ内の'make.sys'というファイルをテキストエディタで開き、下から二行を自分のMKLの入っているフォルダに合わせて変更します。
ちなみに自分の場合は次の様に変更しました。

LIBDIR = /opt/intel/composer_xe_2011_sp1.7.256/mkl
LIBS = -L$(LIBDIR) -lmkl_lapack95 -lmkl_intel -lmkl_intel_thread -lmkl_core -openmp -lpthread

できたらmakeしましょう。
端末から

$ cd espresso-5.0/wannier90/
$ make all

でしばらく待つとコンパイルが終わります。
ホーム・フォルダにある.bashrcを開いて

export PATH=$PATH:/home/koudai/espresso-5.0/wannier90/

を付け加えます。

Crの跳び移り積分

ここでは、サンプルプログラムの解説をしてもよいのですが、spin-density waveの代表的な物質であるクロム(元素記号:Cr)の跳び移り積分を求めてみましょう。
最初から作るのは面倒なので、 /espresso-5.0/wannier90/examples/example08/ にある同じ遷移金属である鉄のインプットファイルを流用します。
ここから

iron.scf
iron.nscf
iron.pw2wan
iron.win

をコピーして、適当なフォルダ(ここではexamplesの中に'chromium'という名前のディレクトリを作りましょう)に貼り付けます。
ファイル名はironの部分をchromiumに代えてしまいましょう。
ただし、.winのファイルだけは下に出てくるseednameに合わせます。
これとさらに偽ポテンシャルのデータが必要です。
偽ポテンシャルのファイルはQunatum ESPRESSOのウェブサイトからダウンロードできます。
以下のサイトから'Cr.pw91-sp-van.UPF'を手に入れてchromiumディレクトリに一緒に入れてください。
QUANTUMESPRESSO - QUANTUMESPRESSO
これで、chromiumディレクトリに次の5つのファイルが入っているか確認してください。

chromium.scf
chromium.nscf
chromium.pw2wan
cr.win
Cr.pw91-sp-van.UP

以下、各ファイルの設定の説明をします。
詳しい説明はプログラムに添付のマニュアル(~/espresso-5.0/PW/Doc/INPUT_PW.html, ~/espresso-5.0/Wannier90/doc/user_guid/user_guide.tex)などを参照してください。

chromium.scf

電子密度を求めるのに必要なファイルです。

&control
    calculation='scf'
    restart_mode='from_scratch',
    pseudo_dir = './',
    outdir='./'
    prefix='cr'
 /

pseudo_dirは偽ポテンシャルが記述されているファイルのあるディレクトリを指定します。
あとはprefixをcrにしておきます。

 &system
    ibrav = 3, celldm(1) =5.453, nat= 1, ntyp= 1,
    ecutwfc =120.0, nspin = 1, nbnd= 16, ecutrho=800.0
    occupations='smearing', smearing='cold', degauss=0.02
 /

ibravはブラベー格子の種類で、体心立方格子の場合は3です。
他の種類についてはマニュアルを参照してください。
celldm(1)格子定数aです。
ブラベー格子が立方晶系でない場合は、celldm(2)やcelldm(3)も指定する必要があります。
ボーア半径を1とした場合の長さを入れてください。
natは、単位胞(ブラベー格子でない)に含まれる原子の数、ntypは原子の種類の数を指定します。
ecutwfc波動関数を計算する際の運動エネルギーのカットオフをリュードベリで指定します。
ecutofrhoは電子密度を計算する際のカットオフで、ecutwfcの8〜12倍程度にします。
nspinはスピンがオーダーしているかどうかで、常磁性の場合は1です。
nbndはバンドの数です。多めにとっておきます。
occupationsは電子がどのように局在しているかを指定し、smearingであればガウシアンに分布しています。
このガウシアンの拡がり方には種類があり、smearingで指定します。
coldはMarzari-Vanderbilt cold smearing(PRL 82, 3296 (1999))のものを使います。
ガウシアンの拡がりはdegaussで指定します。

 &electrons
    conv_thr = 1.0e-8
 /

conv_thrは自己無撞着に解く際の収束の大きさです。

ATOMIC_SPECIES
 Cr 51.9961 Cr.pw91-sp-van.UPF
ATOMIC_POSITIONS
 Cr 0.0 0.0 0.0
K_POINTS (automatic)
16 16 16 0 0 0

ATOMIC_SPECIESでは原子量と偽ポテンシャルのファイルを指定し、それに'Cr'というインデックスをつけます。
ATOMIC POSITIONSでは単位胞内での各原子の位置を指定します(今は単体のクロムを考えているので楽ですね)。

chromium.nscf

固有エネルギーを求めます。
calculation='nscf'を除いて、K_POINTS {crystal}より上の部分はchuromium.scfと同じように書き換えます。
それと、&systemにnosym=.true.を追加しないと、うまくバンドが表現できないみたいです。

 &system
    ibrav = 3, celldm(1) =5.453, nat= 1, ntyp= 1,
    ecutwfc =120.0, nspin = 1, nbnd= 16, ecutrho=800.0 
    nosym=.true.
    occupations='smearing', smearing='cold', degauss=0.02
 /
chromium.pw2wan

次の様に書き直してください。

&inputpp 
   outdir = './'
   prefix = 'cr'
   seedname = 'cr'
   write_mmn = .true.
   write_amn = .true.
   write_unk = .true.
/
cr.win
num_bands         =   16    
num_wann          =   9

num_bandsはバンドの数、num_wannは軌道の数です。
今考えているクロムは5つの3d軌道と3つの4p軌道と1つの4s軌道の計9つです。

bands_plot         = true
hr_plot            = true
wannier_plot       = true
fermi_surface_plot = true
dis_win_max        = 70.0d0
dis_win_min        = 0.0d0
dis_froz_max       = 25.0d0
dis_num_iter       = 1000
dis_mix_ratio      = 1.0
dis_conv_tol       = 1.0E-10
dis_conv_window    = 3

num_iter           = 1000
conv_tol           = 1.0E-10
conv_window        = 3
mp_grid            = 8 8 8

guiding_centres    = true

ワニエ関数を計算する際の設定です。
多いので詳しくはマニュアルを見てください。
mp_gridは.nscfのファイルに合わせます。

begin atoms_frac
Cr  0.000  0.000  0.000
end atoms_frac

begin projections
Cr:sp3d2;dxy;dxz,dyz
end projections

Fermi_energy     18.4914

projectionは考えるWannier関数の軌道です。
Fermi energyはPWscfでchromium.nscfを実行した際にできるchromium.scf.outの中に書かれている値を使います。

begin kpoint_path  
G 0.0000 0.0000 0.0000	 H 0.500 -0.5000 -0.5000	 
H 0.500 -0.5000 -0.5000	 P 0.7500 0.2500 -0.2500	 
P 0.7500 0.2500 -0.2500	 N 0.5000 0.0000 -0.5000	 
N 0.5000 0.0000 -0.5000	 G 0.0000 0.0000 0.000	 
G 0.0000 0.0000 0.000	 H 0.5 0.5 0.5		 
H 0.5 0.5 0.5		 N 0.5 0.0 0.0		 
N 0.5 0.0 0.0		 G 0.0 0.0 0.0		 
G 0.0 0.0 0.0		 P 0.75 0.25 -0.25	 
P 0.75 0.25 -0.25	 N 0.5 0.0 0.0            
end kpoint_path

バンド図をプロットする際の、k空間の経路です。

begin unit_cell_cart
bohr
 2.7265  2.7265 2.7265
-2.7265  2.7265 2.7265
-2.7265 -2.7265 2.7265
end unit_cell_cart

begin unit_cell_cartには、ボーア半径を単位として、単位胞の基本並進ベクトルを入れます。

以下はそのままにしておきましょう。

プログラムの実行

$ cd chromium
$ pw.x chromium.scf.out
$ pw.x chromium.nscf.out
$ wannier90.x -pp cr
$ pw2wannier90.x chromium.pw2wan.out
$ wannier90.x cr

pw.xとpw2wannier90.xの実行には時間がかかります。

フェルミ面を見たければ

$ xcrysden --bxsf cr.bxsf

で、XcrySDenが起動します。
Band numberの全部にチェックしてやれば、すべてのバンドのフェルミ面が見えます。
(全部のバンドを同時に見るには、タブにある[Merged Bands]を選んでください。)
バンド構造は、GNUPLOTを使って

$ gnuplot cr_band.gnu

とするか

$ gnuplot
> load 'cr_band.gnu'

で見ることができます。

飛び移り積分はcr_hr.datに出力されています。
例えば

   -6    2   -4    1    1    0.003331    0.000000

であれば、格子点から格子点の飛び移りのベクトルを {\bf R} = n_1 {\bf a}_1 + n_2 {\bf a}_2 + n_3 {\bf a}_3 としたときの n_1 = -6, n_2  = 2, n_3 = -4 への軌道1から軌道1への飛び移り積分の大きさの実部が 0.003331 eV、虚部が 0.000000 eV という意味です。