Wannier90でLCAOの跳び移り積分を計算
PWscfとWannier90というソフトを使って、LACO法(Linear Combination of Atomic Orbitals method)で用いられる跳び移り積分の計算を行う。
LCAO法は強束縛法(Tight-binding method)と等価である。
電子のスピンは無視して考える。
サイトを中心とした軌道の電子のワニエ型波動関数をとする。
このときのハミルトニアンの行列要素は
すなわち、第二量子化の形でハミルトニアンを書くと
ここではサイトにある軌道の電子の生成演算子である。
LCAO法についてはHarrisonのテキストに詳しく、一読の価値がある。
- 作者: Walter A. Harrison
- 出版社/メーカー: Dover Publications
- 発売日: 1989/07/01
- メディア: ペーパーバック
- クリック: 6回
- この商品を含むブログ (1件) を見る
Elementary Electronic Structure
- 作者: Walter A. Harrison
- 出版社/メーカー: World Scientific Pub Co Inc
- 発売日: 2004/03/01
- メディア: ペーパーバック
- クリック: 7回
- この商品を含むブログを見る
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.xchromium.scf.out
$ pw.xchromium.nscf.out
$ wannier90.x -pp cr
$ pw2wannier90.xchromium.pw2wan.out
$ wannier90.x cr
pw.xとpw2wannier90.xの実行には時間がかかります。
フェルミ面を見たければ
$ xcrysden --bxsf cr.bxsf
で、XcrySDenが起動します。
Band numberの全部にチェックしてやれば、すべてのバンドのフェルミ面が見えます。
(全部のバンドを同時に見るには、タブにある[Merged Bands]を選んでください。)
バンド構造は、GNUPLOTを使って
とするか
で見ることができます。
飛び移り積分はcr_hr.datに出力されています。
例えば
-6 2 -4 1 1 0.003331 0.000000
であれば、格子点から格子点の飛び移りのベクトルを としたときの への軌道1から軌道1への飛び移り積分の大きさの実部が 0.003331 eV、虚部が 0.000000 eV という意味です。