帯電欠損状態の評価 (バージョン2021.01以降)
概要
PHASE/0 では、電荷非中性状態を計算する際、周期的境界条件によるエネルギー発散を抑えるために一様な背景電荷を考慮して処理しています。この処理により、全エネルギーは周期的境界条件による背景電荷間の相互作用を含んだ値となります。このため、例えば、帯電した原子欠損などを含む半導体の生成エネルギーは、実際の値とは異なるものとなります。これを補正する手法として、近年、FNV 法 [Freysoldt09] が注目されており、対応するプログラム( sxdefectalign , CoFFEE など) がいくつか公開されています。PHASE/0 2021年版以降、これらのプログラムに対するインターフェースが配備され、CUBE ファイル形式で出力することができるようになりました。また、これらのプログラムを利用せずに補正量を評価する仕組みも実装されています。さらに、異なる荷電状態間の生成エネルギーを比較した図を描画するためのPythonスクリプトが付属します。
理論
電荷qをもつ欠陥 (D) の生成エネルギーは
で書けます。 \(E_{\rm DFT}^{\rm D, q}\) および \(E_{\rm DFT}^{\rm bulk}\) はDFT計算における全エネルギー値で、前者は帯電欠陥、後者は非帯電・無欠陥の完全結晶 (Pristine) に対応します。 \(n_i\) は完全結晶からの元素iの原子数の増減、 \(\mu_i\) はその化学ポテンシャルです。 \(E_{\rm VBM}\) は完全結晶の価電子帯上端のエネルギー、 \(E_{\rm corr}^q\) は帯電量に依存する補正項です。
FNV法
FNV法では、 \(E_{\rm corr}^{q}\) を2つの寄与に分けて考えます [Freysoldt09] 。
第1項は、孤立した点電荷qの静電エネルギーと、一様な背景電荷を加えて周期的境界条件を課した際の静電エネルギーの差です。
この値は、スーパーセルのサイズ・誘電率から、DFT計算とは無関係に決定されます。一方、第2項は、potential alignment 補正項で、静電ポテンシャルの差
を、欠陥から十分に離れた位置で評価した値です。ここで、 \(V_{\rm PC}^q\) は、静電エネルギー に対応する静電ポテンシャルです。
外部解析プログラム
sxdefectalign、CoFFEEなどの解析プログラムは、以下の入力を要求します。
結果の位置情報
誘電率の値
\(V_{\rm DFT}^{D,q} \left( \mathbf{r} \right)\) や \(V_{\rm DFT}^{\rm bulk}\) の空間分布。たとえばCUBE形式。
出力されるのは以下の情報です。
\(E_{\rm PC}^q\) の値
\(\Delta V \left( \mathbf{r} \right)\) の値。各格子ベクトル方向に対するプロット。他の格子ベクトル方向については面平均をする。欠陥からの距離に対して、原子位置近傍で平均した値を出力することもある。
既知の問題として、以下が挙げられます。
sxdefectalign では、欠陥からの距離に対するプロットが正常に出力されないことがある。
CoFFEE では、各格子ベクトル方向に対するプロットも、ユーザーが python スクリプトを、都度書く必要がある。また、欠陥からの距離に対するプロットは得られない。
CoFFEE では、 \(E_{\rm PC}^q\) の値は、セルサイズを変えて複数計算し、fitting の結果として得る。このため、ユーザーの手間が多い。
extended FNV 法
Kumagaiら [Kumagai14] のextended FNV法では、 \(E_{\rm PC}^q\) は
で表されます。ここで、 \(\gamma\) は適当に選んだ収束パラメータです。また、 \(\varepsilon\) は誘電率テンソル、 \(\Omega\) は系の体積です。一方、ポテンシャル \(V_{\rm PC}^q\) は、
で表されます。ここで、 \(\mathbf{R}_d\) は欠陥の位置です。
\(V_{\rm far}\) の評価に関しては、以下の手順をとります。まず、各原子サイトにおける \(\Delta V\left( \mathbf{r} \right)\) を計算します。このうち、欠陥からの距離 \(R_{WS}\) より遠方にあるデータを平均し、 \(\Delta V_{\rm far}\) とします。なお、 \(R_{WS}\) は、欠陥間の最短距離の半分とします。 (5.85) , (5.86) による計算を用いることによってsxdefectalign, CoFFEEなどを用いずとも(自動的に) \(E_{\rm corr}^q\) の評価が可能です。
使い方
外部解析プログラムの使用
sxdefectalign や CoFFEE に受け渡す静電ポテンシャルを CUBE 形式で出力するには、入力パラメーターファイルを以下のように設定します。
postprocessing{
electrostatic_potential{
sw_write_electrostatic_pot = on
unit = Rydberg ! { Rydberg | Hartree | eV }
}
}
単位系は unit で指定し、Rydberg、Hartree、eV の選択が可能です。default 値は eV です。
結果はelectrostatic_pot.cubeファイルに出力されます。ヘッダー部に単位系が表示され、以降、各軸方向の分割数、メッシュの刻み幅、各メッシュ点における値と続きます。
Calculated by phase
Local+Hartree potential in Rydberg
215 0.0000 0.0000 0.0000
180 0.180912 0.000000 0.000000
180 0.000000 0.180912 0.000000
180 0.000000 0.000000 0.180912
33 33.000000 2.396007 2.396052 2.396080
(後略)
帯電欠陥、及び完全結晶の静電ポテンシャルを出力しておくようにします(sxdefectalignを利用する場合には、unitは Rydberg にしておく必要があります。)。以下に、これらのCUBEファイルを用いたsxdefectalign の実行例を示します。
sxdefectalign --charge 3 --eps 12.88
--vdef ../electrostatic_pot.cube --vref ../bulk/electrostatic_pot.cube
--center 0.0,0.0,0.0 --relative
--ecut 270 --qe > Log
系の帯電量を q としたとき、charge には -q を渡します。eps には誘電率の値を指定します。vdef及びvrefには、帯電欠陥及びバルクにおける静電ポテンシャルファイルを、それぞれ指定します。center には欠陥位置の内部座標を与えます。ecut には nfinp.data で指定した cutoff_cd を Ry 単位で渡します。
sxdefectalignのLog ファイルの末尾には、以下のような出力がされます。
Defect correction (eV): 0.828347 (incl. screening & alignment)
ここで表示されているのは、 値 ( 式 (5.83) 参照 )であり、potential alignment の寄与は含まれていません。
Logファイル以外には、 vline-eV-an.dat ( n=0, 1, 2 )
及び vAtoms.dat
が出力されます。前者は、式 (5.84) の各結晶軸方向に対する、ポテンシャルの1次元プロファイルを出力したものです。一方、後者は、各原子の欠陥位置からの距離、及び原子上でのポテンシャルの値を出力したものです。ただし、 vAtoms.dat
は正常に出力されないことが多々あるようです。
extended FNV 法
まず帯電なし・欠陥なしの bulk の計算を行います。入力パラメーターファイルは以下のように設定します。
postprocessing{
electrostatic_potential{
sw_write_electrostatic_pot = on
}
}
結果はelecpot_bin.data ファイルが出力されます。静電ポテンシャルのG空間での値をバイナリ形式で出力しているため、外部ソフトと連携することはできません。ファイル名は、F_ELECPOT_BINにより変更できます。
&fnames
F_ELECPOT_BIN = './elecpot_bin.data'
/
ついで帯電欠陥の計算を行います。入力パラメーターファイルは以下のように設定します。
postprocessing{
electrostatic_potential{
sw_write_electrostatic_pot = on
}
charged_defect{
correction{
sw_calc_extfnv_correction = on ! default : off
dielectric_constant{
exx = 12.88
eyy = 12.88
ezz = 12.88
}
position{
x = 0.0, y = 0.0, z = 0.0
}
}
}
}
charged_defect ブロック内の correction ブロックで、sw_calc_extfnv_correction = on を設定します。また、dielectric_constat ブロック内で誘電率テンソルの値、position ブロック内で欠陥の内部座標を指定します。欠損以外の場合には、内部座標ではなく、原子番号( atom_id )による指定も可能です (例:atom_id= 2 )。 さらに、参照する bulk の静電ポテンシャルファイル名をF_ELECPOT_BIN_REF で指定します。
&fnames
F_ELECPOT_BIN_REF = '../bulk/elecpot_bin.data'
F_ELECPOT_BIN = './elecpot_bin.data'
/
結果はdefect_pot_correction.direction_n ( n=1,2,3, )、defect_pot_correction.atomsに出力されます。前者は、式 (5.84) の各結晶軸方向に対する、ポテンシャルの1次元プロファイルを出力したものです。一方、後者は、各原子の欠陥位置からの距離、及び原子上でのポテンシャルの値を出力したものです。
以下に、defect_pot_correction.direction_1 の出力例を示します。
# dist. (Ang), pot_diff, Vpc, pot_diff -Vpc (eV)
0.00000 -0.13205 -0.03907 -0.09299
0.09573 -0.13743 -0.03902 -0.09841
(後略)
第2列及び第3列は、それぞれ式 (5.84) の \(V_{\rm DFT}^{D,q} \left( \mathbf{r} \right) - V_{\rm DFT}^{\rm bulk} \left( \mathbf{r} \right)\) 及び \(V_{\rm PC}^q \left( \mathbf{r} \right)\) に対応します。第4列は \(\Delta V \left( \mathbf{r} \right)\) に対応します。なお、符号については、電子の電荷が負であることを考慮しています。
以下に、defect_pot_correction.atoms の出力例を示します。
# no., dist. (Ang), pot_diff, Vpc, pot_diff -Vpc (eV)
1 2.1961 -0.0550653932 -0.2936706923 0.2386052991
2 3.9534 -0.1647261198 -0.1693351641 0.0046090443
(中略)
# Correction energy (eV): 0.7827480523
なお、最終行に表示されている “Correction energy” は、生成エネルギー図作成に必要な \(E_{\rm corr}^q\) です。
生成エネルギー図の作成1
欠陥ごとの生成エネルギー図の作成方法を説明します。生成エネルギー図作成にはcalc_defect_formation_energy.pyスクリプトを用います。
平衡状態のとき、化学ポテンシャルは
を満たします。ここで、 \(\mu_{\rm Ga}^0\) 及び \(\mu_{\rm As}^0\) は、それぞれGa及びAs単体の1原子あたりのエネルギーです。また、 \(\mu_{\rm GaAs}^0\) は、GaAs 結晶の2原子あたりのエネルギーです。特に、Ga-rich 極限の場合は、
となり、As-rich 極限では
となります。これら化学ポテンシャルの情報は、以下で説明するcalc_defect_formation_energy.pyスクリプトの入力に記述する必要があるものです。
例えば、GaAs のスーパーセルの中にGa欠損を導入した計算を行い、帯電量 q =-3, -2, -1, 0 の計算結果が得られているとします。生成エネルギー図作成のために、以下のようなファイルを作成します。ファイル名に制限はありません。以下の例ではtmp1.in とします。
&VBM #(eV)
6.00887
&band_gap #(eV)
1.424
&Chemical_potentials #(Ha)
Ga -138.9838873703 # mu_GaAs_bulk -mu_As
As -87.9848825305
&Defects #elements, number ( negative==vacancy, positive==impurity )
Ga -1
&Host_supercell_energy #(Ha)
-7263.0006368270
&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0 -7123.9292947570
&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0
以下、用語の説明をします。
キーワード |
単位 |
意味 |
&VBM |
eV |
\(E_{\rm VBM}\) |
&band_gap |
eV |
バンドギャップ値 |
&Chemical_potentials |
Hartree |
化学ポテンシャル \(\mu_i\) 。系を構成する全元素について記述する。 |
&Defects |
導入した欠陥 |
|
&Host_supercell_energy |
Hartree |
|
&Defective_supercell_energy |
Hartree |
|
&Correction |
eV |
|
以下の要領でcalc_defect_formation_energy.pyスクリプトを実行します。
python3 calc_defect_formation_energy.py input
[-o OUTFILE]
[--emin EMIN] [--emax EMAX] [--de DE]
[--vmin VMIN] [--vmax VMAX]
[--image_format IMAGE_FORMAT]
括弧内は省略可能なオプションで、その意味は以下のとおりです。
引数 |
意味 |
デフォルト値 |
-o |
出力するファイルの名称 |
result |
--emin |
エネルギーEF の最小値 |
-1.0 |
--emax |
エネルギーEF の最大値 |
6.0 |
--de |
エネルギーEFの刻み幅 |
0.01 |
--vmin |
生成エネルギーの表示範囲の最小値 |
-5.0 |
--ymax |
生成エネルギーの表示範囲の最大値 |
5.0 |
--image_format |
可視化画像の形式 (png/eps) |
png |
実行例を以下に示します。
python3 calc_defect_formation_energy.py tmp1.in -o result1
上述のコマンドを実行した結果、生成されるファイルは result1.qdep、result1.min、result1.gnu、result1.png です。以下に、result1.qdep及びresult1.min の一部を示します。前者には、各帯電状態における生成エネルギーの、フェルミエネルギー依存性が出力されています。後者は、各フェルミエネルギー値における最小の生成エネルギー値が出力され、ファイル末尾に Charge Transition Level が追記されています。
result1.qdepファイルの内容
# Formation energy
# Ef (eV) q=-3 q=-2 q=-1 q=0
-1.00000 6.48430 4.90715 3.53961 2.37978
(後略)
result1.minファイルの内容
# Formation energy
# Ef (eV) min
-1.00000 2.37978
(中略)
6.00000 -14.51570
# Charge Transtion level [eV]
#-2/-3 0.57715
#-1/-2 0.36754
# 0/-1 0.15983
result1.gnu は可視化のための gnuplot 用ファイル、result1.png は gnuplot の出力です。
生成エネルギー図の作成2
興味ある全ての欠陥構造について、 生成エネルギー図の作成1 の作業が終了しているものとします。ここでは、これらをまとめた図の作成を行うため、以下のようなファイルを作成します。ファイル名に制限はありません。以下の例ではgather1.in とします。
#
# title filename ( excluding ".min" )
#
&List
Vac_Ga Vacancy_Ga/result1
Vac_As Vacancy_As/result1
Ga_As Ga_for_As/resultaa
As_Ga As_for_Ga/resultaa
&band_gap #(eV)
1.424
以下、用語の説明をします。
ワード |
単位 |
意味 |
&List |
欠陥の名称、及び計算の出力( 生成エネルギー図の作成1 指定した outfile 名) |
|
&band_gap |
eV |
バンドギャップ値 |
以下の要領でplot_multiple_defect_formation_energy.pyスクリプトを実行します。
python3 plot_multiple_defect_formation_energy.py input
[-o OUTFILE]
[--emin EMIN] [--emax EMAX]
[--vmin VMIN] [--vmax VMAX]
[--image_format IMAGE_FORMAT]
[--keypos_h KEYPOS_H] [--keypos_v KEYPOS_V]
括弧内は省略可能なオプションで、その意味は以下のとおりです。なお、EMIN、EMAX値は、 生成エネルギー図の作成1 の指定と揃えた方がよいです。
引数 |
意味 |
デフォルト値 |
-o |
出力するファイルの名称 |
result_all |
--emin |
エネルギーEF の最小値 |
なし |
--emax |
エネルギーEF の最大値 |
なし |
--vmin |
生成エネルギーの表示範囲の最小値 |
なし |
--ymax |
生成エネルギーの表示範囲の最大値 |
なし |
--image_format |
可視化画像の形式 (png/eps) |
png |
--keypos_h |
凡例の水平位置 (left/center/right) |
right |
--keypos_v |
凡例の垂直位置 (top/center/bottom) |
top |
以下に実行例を示します。
python3 plot_multiple_defect_formation_energy.py gather1.in
実行するとresult_all.gnu 及び result_all.png が生成されます。前者は gnuplot 用ファイルで、後者はこれを可視化したものです。
例題
概要
GaAs 64 原子をホストとして、帯電欠陥の計算を行った例を紹介します。計算条件は以下の通りです。なお、GaAs の格子定数は、基本格子で最適化を行いました。補正エネルギーは、PHASE/0に実装されている extended FNV 法により評価しました。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
Monkhorst-Pack (2×2×2) |
交換相関相互作用 |
GGAPBE, PAW |
単位胞の1辺 [Å] |
11.48882 |
SCF 収束条件 |
[Ha/atom] 1.0E-8 |
力の収束条件 |
[Ha/bohr] 2.0E-4 |
擬ポテンシャル |
Ga_ggapbe_paw_02.pp, As_ggapbe_paw_02.pp |
(補正項の計算で使用する) |
誘電率 ※12.88 |
※ https://www.microwaves101.com/encyclopedias/gallium-arsenide
例題のディレクトリー構成および計算のながれ
本例題は複数の計算を行い、その結果をとりまとめ、スクリプトで処理することによって結果が得られる仕組みになっています。ここではディレクトリーの構成と計算のながれについて説明します。
必要な入力ファイルはサンプルディレクトリーの下の samples/defectq
以下のサブディレクトリーに配置されています。以下のようなディレクトリー構成になっています。
defectq
以下には Preparation
ディレクトリーと GaAs_64_lattice_opt
ディレクトリーが存在します。
Preparationディレクトリー
GaAs_64_lattice_optディレクトリー
各種生成エネルギー計算の入力ファイルが格納されています。ベースとなる結晶は格子最適化によってもとまった格子定数から作成した64原子系です。GaAs_64_lattice_optディレクトリーは、さらに以下のようなサブディレクトリー群が存在します。
GaAs_64_lattice_opt以下のディレクトリー構成 ディレクトリー名
説明
Pristine
欠損のない結晶の入力ファイルが納められたディレクトリー。
Ga_for_As
AsをGaで置換した欠陥構造の入力ファイルが納められたディレクトリー。
As_for_Ga
GaをAsで置換した欠陥構造の入力ファイルが納められたディレクトリー。
As_interstitial_As4
Asが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。
As_interstitial_Ga4
Gaが最近接原子となる位置に入り込んだAs interstitialの入力ファイルが納められたディレクトリー。
Ga_interstitial_As4
Asが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。
Ga_interstitial_Ga4
Gaが最近接原子となる位置に入り込んだGa interstitialの入力ファイルが納められたディレクトリー。
それぞれのディレクトリーにはさらにq_ q ディレクトリーが存在します。ここで q は電荷をあらわす数値です。
Pristineの場合はq_0とq_0.2ディレクトリーが存在します。 \(E_{\rm VBM} = \frac{E\left( N \right) - E\left( N-\Delta N \right)}{\Delta N}\) (今の場合 \(\Delta N = 0.2\) )という関係から \(E_{\rm VBM}\) をもとめ、スクリプト入力の&VBMに記述します。またq_0のエネルギーをスクリプト入力の&Host_supercell_energyに記述します。
ほかの欠損に対応するディレクトリーでは、エネルギーの計算結果をスクリプト入力の&Defective_supercell_energyに記述します。また、q_0以外の計算ではextended FNV法による補正エネルギー \(E_{\rm corr}^q\) の計算がなされます。結果はdefect_pot_correction.atomsファイルの末尾に記録されるので、その計算値をスクリプト入力の&Correction に記述します (q=0の項には0を記述します)。
各欠損ディレクトリーにおいて calc_defect_formation_energy.py スクリプトの入力を作成し、実行することによって各欠損の生成エネルギー図を作成することができます。またすべての欠陥のデータを集約し 生成エネルギー図の作成2 の手続きをふむことによってすべての欠陥の結果をまとめた生成エネルギー図を作成することができます。
補正エネルギー比較
以下に、Ga欠損 (q=-3) における補正エネルギーについて、sxdefectalign とPHASE/0実装のextentend FNV 法による評価の比較を示します。両者がおおよそ一致していることが分かります。
sxdefectalign |
PHASE/0 |
|
Epc |
1.24251 |
1.24272 |
dV |
-0.00731 |
-0.02779 |
Ecorr |
1.22058 |
1.15935 |
生成エネルギー図の作成 1.
以下では、Ga欠損に対して、2つの極限における入力を作成し可視化します。以降、Ga-rich 極限及びAs-rich極限の入力を示します。
calc_defect_formation_energy.py 用の入力 ( Ga-rich 極限におけるGa欠損; samples/defectq/GaAs_64_lattice_opt/Vacancy_Ga/cond_Ga_rich.in
)
&VBM #(eV)
6.00887
&band_gap #(eV)
1.424
&Chemical_potentials #(Ha)
Ga -138.9586691142 #mu_Ga
As -88.0101007866 #mu_GaAs -mu_Ga
&Defects #elements, number ( negative==vacancy, positive==impurity )
Ga -1
&Host_supercell_energy #(Ha)
-7263.0006368270
&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0 -7123.9292947570
&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0
calc_defect_formation_energy.py 用の入力 ( As-rich 極限におけるGa欠損; samples/defectq/GaAs_64_lattice_opt/Vacancy_Ga/cond_As_rich.in
)
&VBM #(eV)
6.00887
&band_gap #(eV)
1.424
&Chemical_potentials #(Ha)
Ga -138.9838873703 # mu_GaAs_bulk -mu_As
As -87.9848825305 #mu_As
&Defects #elements, number ( negative==vacancy, positive==impurity )
Ga -1
&Host_supercell_energy #(Ha)
-7263.0006368270
&Defective_supercell_energy # charge_state(q) and energy (Ha)
-3 -7123.2688494462
-2 -7123.4909799361
-1 -7123.7111092854
0 -7123.9292947570
&Correction energy # charge_state(q) and energy (eV)
-3 1.1593520601
-2 0.6178614391
-1 0.2315198082
0 0.0
つぎのコマンドを実行すると、生成エネルギーの図 (result-Ga-rich.png及びresult-As-rich.png ) が出来ます。結果を 図 5.97, 図 5.98 に示します。
calc_defect_formation_energy.pyの実行 (Ga欠損)
python3 calc_defect_formation_energy.py
cond_Ga_rich.in -o result_Ga_rich --emin -0.5 --emax 2.0
python3 calc_defect_formation_energy.py
cond_As_rich.in -o result_As_rich --emin -0.5 --emax 2.0

Ga-rich limit (ファイル名:result-Ga-rich.png )

As-rich limit (ファイル名:result-As-rich.png )
生成エネルギー図の作成 2.
以下で、種々の欠陥構造に対する生成エネルギーをまとめた図を作成します。Ga-rich 及びAs-rich 極限のファイルが、それぞれ、 samples/defectq/GaAs_64_lattice_opt
の下の gather_Ga_rich.in
及び gather_As_rich.in
にあります。後者は、つぎに紹介する入力のresult-Ga-rich を result-As-rich に置換したものです。 なお、欠陥の種類としては、文献 [Broberg18] 記載のものを選択しました。
gather_Ga_rich.inの内容
#
# title filename ( excluding ".min" )
#
&List
Vac_Ga Vacancy_Ga/result-Ga-rich
Vac_As Vacancy_As/result-Ga-rich
Ga_As Ga_for_As/result-Ga-rich
As_Ga As_for_Ga/result-Ga-rich
Ga_i_As4 Ga_interstitial_As4/result-Ga-rich
As_i_As4 As_interstitial_As4/result-Ga-rich
Ga_i_Ga4 Ga_interstitial_Ga4/result-Ga-rich
As_i_Ga4 As_interstitial_Ga4/result-Ga-rich
&band_gap #(eV)
1.424
つぎのコマンドを実行すると、生成エネルギーの図 (results-Ga-rich.png及びresults-As-rich.png) が出来ます。結果を 図 5.99 および 図 5.100 に示します。
python3 plot_multiple_defect_formation_energy.py
gather_Ga_rich.in -o results_Ga_rich --emin -0.2 --emax 2.0 --vmin -4.0 --vmax 5.0
--keypos_h left --keypos_v bottom
python3 plot_multiple_defect_formation_energy.py
gather_As_rich.in -o results_As_rich --emin -0.2 --emax 2.0 --vmin -4.0 --vmax 5.0
--keypos_h left --keypos_v bottom

Ga-rich limit (ファイル名:results-Ga-rich.png)

As-rich limit (ファイル名:results-As-rich.png)
文献 [Broberg18] と同様の生成エネルギー図が得られました (生成エネルギーの絶対値や Charge Transition Level は、完全には一致しません)。
- Freysoldt09(1,2)
C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, Phys. Rev. Lett. 102 (2009) 016402.
- Kumagai14
Y. Kumagai and F. Oba, Phys. Rev. B 89 (2014) 195205.
- Broberg18(1,2)
D. Broberg, B. Medasani, N. Zimmermann, A. Canning, M. Haranczyk, M. Asta, and G. Hautier,, Comput. Phys. Commun. 226 (2018) 165.
高精度な局所状態密度 (バージョン2021.01以降)
概要
原子分割局所状態密度とは,空間を原子ごとの領域に分割し,その領域における電子状態密度をもとめる計算機能です。積算状態密度のフェルミエネルギーにおける値を調べることによって,対象原子の電子数をもとめるという用途に用いることも可能です。
原子分割局所状態密度は,波動関数および電荷密度の実空間におけるFFTメッシュを空間上で最も近い原子に割り当て,その寄与分を加算することによって計算します。このような方法の場合,実空間のFFTメッシュと原子位置の関係によって対称性から等価な原子間でも割り当たるメッシュ数などが異なり,結果が微妙に異なる場合があります。言い換えると,原子配置を単位胞に対してどのように定義するかによって計算結果が変化します。たとえば,対称性から等価なはずの原子の電子数が互いに異なる値となってしまう場合があります。この振る舞いを改善するため,2021年版以降電荷密度の割り当てを実空間FFTメッシュではなく原子中心メッシュに切り替えて評価する方法を用いることができます。このようにするとメッシュは常に同じように各原子に割り当てられるため,原子配置を単位胞に対してどう定義するかによって結果が左右されなくなります。原子中心メッシュは原子ごとに定義され,その値は実空間FFTメッシュ上の値の三次元線形補間によって求められます。
使い方
原子中心メッシュによる局所状態密度計算の指定は,postprocessingブロックにおけるldosブロックにおいて行います。以下のような指定を行います。
postprocessing{
dos{
sw_dos = on
method = t
}
ldos{
sw_aldos = on
aldos{
sw_atom_centered_mesh = on
atom_centered_mesh_factor = 1
}
}
}
postprocessingブロックにdosブロックを作ると,全状態密度の設定を行うことができます。局所状態密度の計算の基本設定は全状態にならう形式になっているので,ここでの設定は局所状態密度にもあてはまります。ldosブロックにおいて局所状態密度計算の設定を,さらにその下のaldosブロックにおいて原子分割局所状態密度計算の設定を行います。aldosブロックにおいて新たに利用できるようになった設定項目は下記の通り。
sw_atom_centered_meshもしくはsw_ac_mesh : onとすることによって,原子分割局所状態密度計算のメッシュが実空間FFTメッシュから原子中心メッシュに切り替わります。デフォルト値はoff.
atom_centered_mesh_factorもしくはac_mesh_factor : 原子中心メッシュはデフォルトの振る舞いではその“濃さ”は実空間FFTメッシュと同じですが,ここで指定する係数分増やすことも可能です。たとえば2とすると,3方向のメッシュ数がそれぞれ2倍となります。その結果計算量は8倍となる点は注意が必要です。デフォルト値は1.
例題
Si結晶およびグラファイトの例題がサンプルディレクトリーの samples/dos_band/aldos_by_acmesh
以下にあります。Si結晶の例題はディレクトリーSi2以下,グラファイトの例題はディレクトリーgraphite以下に配置されており,それぞれ acmesh
と fftmesh
フォルダーが存在し,前者が原子中心メッシュを用いた局所状態密度計算,後者がFFTメッシュを用いた局所状態密度計算の入力ファイルが格納されています。
いずれの結晶も,精度の高い結果が得られるよう比較的大きなカットオフエネルギー(80 Rydberg)を採用しています。また,精密な状態密度を得るため四面体法を用いる設定を施しています。すなわち,入力パラメーターファイルには以下のような設定が施されています。
accuracy{
cutoff_wf = 80 rydberg
...
ksampling{
method = mesh
...
}
smearing{
method = t
}
}
それぞれの最下層のディレクトリーにおいて計算を実行すると,SCF計算のあと状態密度の計算が行われ,結果がdos.dataファイルに記録されます。原子分割局所状態密度の計算結果は ALDOS num_atom = aid
という文字列からはじまる行以降に記録されます。ここで aid
は原子のIDで,各原子の入力ファイルにおける定義順に対応します。各原子の電子数はフェルミエネルギーにおける積算状態密度の値です。
dos.dataファイルは(スピンを考慮していない場合)4カラム目のデータがフェルミエネルギーを原点としたエネルギー,6カラム目が積算状態密度に対応するので,4カラム目が0となる行の6カラム目の値がその原子の電子数に対応することになります。たとえば以下のようなデータがdos.dataに記録されている場合,1番目の原子の電子数はおおよそ3.9803となります。
...
...
ALDOS num_atom = 1
No. E(hr.) dos(hr.) E(eV) dos(eV) sum
6 -0.33730 0.0000000000 -11.955933 0.0000000000 0.0000000000
16 -0.33630 0.0000000000 -11.928722 0.0000000000 0.0000000000
...
...
4366 0.09870 0.6019928210 -0.091764 0.0221228201 3.9797782305
4376 0.09970 0.2103606325 -0.064553 0.0077306078 3.9801773503
4386 0.10070 0.0555466368 -0.037342 0.0020413005 3.9802911964
4396 0.10170 0.0056133904 -0.010130 0.0002062882 3.9803167566
4406 0.10270 0.0000000000 0.017081 0.0000000000 3.9803179066
4416 0.10370 0.0000000000 0.044293 0.0000000000 3.9803179066
...
...
以下,この例題によって得られるSi結晶およびグラファイトのある原子の電子数を報告します。
Si, fftmesh |
Si, acmesh |
C, fftmesh |
C, acmesh |
3.9803 |
4.0000 |
4.0537 |
3.9992 |
いずれの例題もすべての原子は対称性から等価なため,電子数としては4という値が得られるはずです。 表の値から明らかなように,原子中心メッシュともにFFTメッシュの結果よりも正しい解に近い結果が得られています。
高精度なXPS計算 (バージョン2021.01以降)
概要
2020年版までのPHASE/0 では、XPS・XANES のピーク位置が実験値と異なることがありました。原因として、内殻電子の相対論的効果、スピン分裂効果が考えられます。前者については、CIAO を改変し相対論的効果を含む (DIRAC型の) 運動エネルギーを出力可能なようにしました。後者については、PHASE/0 のopencore 法を改変し、その効果を取り込めるようにしました。これによって、SiC、AlN、Li7Ti5O12 結晶の内殻軌道の結合エネルギーについて、従来法に比べて実験値との差が小さくなることなどを確認しています。
理論
相対論的運動エネルギー
Kohn-Sham 方程式は、
で表されます。 \(T\) 及び \(V\) はそれぞれ、運動エネルギー及びポテンシャル演算子、また \(\varepsilon_i\) 及び \(\psi_i\) は固有値・固有波動関数です。ここで、運動エネルギー演算子は、相対論的効果の取り入れ方により、いくつかの形状が提案されています [Lenthe96] 。
\(T^{\rm NR}\) は非相対論、 \(T^{\rm ZORA}\) はZORA (Zero Order Regular Approximation)、 \(T^{\rm DIRAC}\) はDirac相対論によるものです。ZORAの場合は、scalar relativistic 項とスピン軌道相互作用項に分解できます。
Dirac の場合、上式にて \(V\) を \(V-\varepsilon_i\) と読み替えればよく、
となります。 ただし、擬ポテンシャル作成時はscalar-relativisitic で解く必要があり、スピン軌道相互作用による効果は考慮されません。このため、CIAO が擬ポテンシャルファイルに出力する、内殻電子の運動エネルギーは scalar-relativisitic 項のみで評価します。
さて、内殻電子の運動エネルギーは
で表されます。上記は DIRAC の場合ですが、ZORA の場合も同様です。 最後に、CIAOではポテンシャルは球対称として扱うので、波動関数は動径方向成分のみをもちます。このため、運動エネルギー演算子は以下のように展開できます。
内殻電子のスピン
CIAOでは、擬ポテンシャルはスピン中性で作成します。内殻電子についても同様であり、この時の内殻電子密度を \(\rho_{\rm core}\) とします。ここで、ある軌道 (n,l) におけるスピン占有数 \(f_{n,l}^{\uparrow}\) 及び \(f_{n,l}^{\downarrow}\) が変わり、 \(f_{n,l}^{\uparrow} \neq f_{n,l}^{\downarrow}\) になったとします。ただし、 \(f_{n,l}^{\uparrow} + f_{n,l}^{\downarrow}\) は擬ポテンシャル作成時と同じであるとします。このとき。各スピンの内殻電子密度は以下のように表されます。
ここで、
はスピン密度の差で、波動関数 \(\phi_{n,l}^{\rm AE}\) はスピン状態により不変と仮定しています。 このスピン密度変化により影響を受けるエネルギーは、交換相関相互作用 (XC) エネルギー
です。第1項は非PAW項、第2、3項はPAW由来の項です。なお、スピン密度の和は不変であるので、Hartreeポテンシャルは変更を受けません。よって、XC ポテンシャルのみ \(\Delta \rho \left( r \right)\) を考慮します。
XPSのピークエネルギーの計算方法
非金属の固体のXPSにおけるピークエネルギーは、以下の結合エネルギー \(E_{\rm B}\) に対応します [Ozaki17] 。
ここで、 \(N\) は系の電子数です。 \(\mu_0\) は系の化学ポテンシャルで、gapのある系では、価電子帯上端 (VBM) から伝導体下端 (CBM) まで取りえます。特に、真性半導体では、バンドギャップを \(E_g\) として、
となります。なお、VBMのエネルギー \(E_{\rm VBM}\) は
で計算できます。また、帯電状態を周期的境界条件で計算するので、その補正法の1つとして、式 (5.82) の ( q =1 ) を加えることもできます。
一方、系が金属の場合には、
となります。
なお、固体ではなく分子の場合は、
となります。
使い方
相対論的運動エネルギー
相対論的運動エネルギーを出力するためのCIAOの入力は下記の通り。
# PAW
sw_paw 1
# CORE ELECTRON INFO
sw_with_dipole_cor2val 1
method_ekin_core 1 ! default:0 (非相対論), 1:DIRAC, 2: ZORA
method_ekin_core 1 或いは 2 で、内殻電子の相対論的運動エネルギーを出力します。なお、”sw_with_dipole_cor2val 1” は、これまでと同様、内殻電子のXPS・XANES計算には必須です。
“method_ekin_core 1”の場合の3つのスイッチの設定は、以下の入力でも可能です。
sw_write_core_full 1
結果はgncpp2ファイル(擬ポテンシャルファイル)に出力されます。
CORE ENERGY CONTRIB
0.319190836204E+02 kin ! 内殻電子の運動エネルギー
-0.677766126155E+02 ion
0.697542582590E+01 hartr
内殻電子のスピン
ある原子の内殻電子が励起された状態を考える場合、該当する軌道に正孔を入れた擬ポテンシャルを使う必要があります。例えば、1s軌道励起の場合は 1s 軌道に正孔、2p軌道励起の場合は2p軌道に正孔を入れます。 これらの原子は、内殻電子は開殻となるので、スピン自由度をもちます。内殻電子のスピンを扱うには、sw_opencore = on とします。デフォルトでは、式 (5.98) の全項でスピンを考慮しますが、sw_xc_opencore_ae_only = onとすると、第3項 (AE部分)のみ考慮します。 価電子が磁気モーメントをもつ磁性材料では、spin_orientation の指定により、内殻電子と (原子近傍の) 価電子の磁気モーメントの向きを、平行・半平行のいずれかに指定できます。ただし、SCF iteration で価電子の磁気モーメントが反転したりすると、内殻電子も追随するため、収束性が悪くなる場合があります。そこで、sw_fix_core_spin_pol = onとすると、初期磁気モーメントの向きに固定することが出来ます。
内殻電子のスピンを考慮するための入力例。
accuracy{
paw = on
core_electrons{
sw_opencore = on ! default :off
sw_xc_opencore_ae_only = on ! default :off
spin_orientation = anti_parallel ! anti_parallel or parallel (default)
sw_fix_core_spin_pol = on ! default: off
}
}
postprocessing{
corelevels{
sw_calc_core_energy = on ! XPS 計算に必要
}
}
また、s軌道以外から励起する場合には、内殻電子のスピン軌道相互作用により、XPSのピークが2つに分裂します。例えば、2p軌道の場合には、2p1/2と2p3/2 に分裂します。この場合、以下のようにすると、両者の全エネルギー値を計算することが出来ます。
内殻電子準位のスピン軌道分裂を考慮するための入力例。
postprocessing{
corelevels{
sw_calc_core_energy = on
corehole{
atom_id = 1 ! 励起する原子
orbital = 2p ! 励起する軌道
}
}
}
(内殻電子のエネルギーを含む) 系の全エネルギーは、core_energy.data に出力されます。sw_calc_core_energy = onとした場合には、core_energy.data 末尾に、スピン軌道分裂したエネルギー値が追記されます。 (5.99) (5.102) (5.104) などの結合ネルギーの計算における \(E_{\rm final}\) や \(E_{\rm initial}\) はここに記録された値を利用します。
core_energy.dataの出力例
# Etotal (Core+Valence)
-18558.3245482622
# Etotal (Core+Valence+Soc_corehole)
# J = 3/2: -18558.3331659845
# J = 1/2: -18558.3073128174
例題
例題の計算結果を紹介します。XPSの例題は samples/XPS
ディレクトリーに配置されています。用いる擬ポテンシャルは 相対論的運動エネルギー において説明したsw_write_core_full の値を1として作成した擬ポテンシャルで、それぞれの例題のpseudoサブディレクトリーに配置されています。
4H-SiC結晶 (C 1s)
4H-SiC結晶の例は samples/XPS/4H-SiC
以下に配置されています。そのディレクトリー構成は下記の通りです。
第一階層 |
第二階層 |
説明 |
Final |
終状態の入力ファイル |
|
q_0_no_opencore |
中性・内殻スピンなし |
|
q_0_with_opencore |
中性・内殻スピンあり |
|
q_1_no_opencore |
電荷+1・内殻スピンなし |
|
q_1_with_opencore |
電荷+1・内殻スピンあり |
|
Initial |
初期状態の入力ファイル |
|
q_0 |
中性 |
|
q_0.2 |
電荷+0.2 |
|
q_0.5 |
電荷+0.5 |
|
q_1 |
電荷+1 |
計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.108 Å、c= 10.170 Å )を \(3 \times 3 \times 3\) 倍しました。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
Monkhorst-Pack 2×2×2 |
交換相関相互作用 |
GGAPBE, PAW |
SCF 収束条件 [Ha/atom] |
1.0E-8 |
内殻電子のスピン |
sw_xc_opencore_ae_only = on |
CIAOの入力 |
sw_write_core_full 1 |
(補正項計算で使用する)誘電率 |
※9.76 (ab面内), 10.32 (c軸) |
※https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118313534.app3
計算結果は以下の通りです。
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
\(E_{\rm final} (N-1) - E_{\rm initial} (N)\) |
279.258 |
272.309 |
\(E_{\rm corr}^1\) |
0.237 |
0.237 |
\(E_{\rm VBM}\) ( \(\Delta N=0.2\)) |
9.307 |
|
\(E_g\) |
2.435 |
|
\(E_{\rm B}\) |
289.975 |
283.026 |
実験値 |
283 +/- 0.8 a) |
a) http://www.xpsfitting.com/2012/01/silicon.html
内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。 なお、式 (5.103) を用いた場合、 \(E_{\rm B}\) は以下のようになりました。
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
|
290.815 |
283.864 |
w-AlN 結晶 (N 1s)
w-AlN結晶 (N 1s) の例は samples/XPS/w-AlN
以下に配置されています。そのディレクトリー構成は下記の通りです。
第一階層 |
第二階層 |
説明 |
Final_N1s |
終状態の入力ファイル |
|
q_0_no_opencore |
中性・内殻スピンなし |
|
q_0_with_opencore |
中性・内殻スピンあり |
|
q_1_no_opencore |
電荷+1・内殻スピンなし |
|
q_1_with_opencore |
電荷+1・内殻スピンあり |
|
Initial |
初期状態の入力ファイル |
|
q_0 |
中性 |
|
q_0.2 |
電荷+0.2 |
|
q_0.5 |
電荷+0.5 |
|
q_1 |
電荷+1 |
計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.140 Å、c= 5.040 Å ) \(3 \times 3 \times 2\) 倍しました。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
Monkhorst-Pack 2×2×2 |
交換相関相互作用 |
GGAPBE, PAW |
SCF 収束条件 [Ha/atom] |
1.0E-8 |
内殻電子のスピン |
sw_xc_opencore_ae_only = on |
CIAOの入力 |
sw_write_core_full 1 |
(補正項計算で使用する)誘電率 |
※8.23 (ab面内), 9.74 (c軸) |
※https://materialsproject.org/materials/mp-661/
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
\(E_{\rm final} (N-1) - E_{\rm initial} (N)\) |
395.772 |
387.703 |
\(E_{\rm corr}^1\) |
0.310 |
0.310 |
\(E_{\rm VBM}\) ( \(\Delta N=0.2\)) |
6.735 |
|
\(E_g\) |
4.404 |
|
\(E_{\rm B}\) |
405.019 |
396.550 |
実験値 |
397.4 a) |
a) [Mahmood03]
内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。 なお、式 (5.103) を用いた場合、 \(E_{\rm B}\) は以下のようになりました。
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
|
406.976 |
398.889 |
w-AlN 結晶 (Al 2p3/2)
w-AlN結晶 (Al 2p3/2) の例は samples/XPS/w-AlN
以下に配置されています。そのディレクトリー構成は下記の通りです。
第一階層 |
第二階層 |
説明 |
Final_Al2p |
終状態の入力ファイル |
|
q_0_no_opencore |
中性・内殻スピンなし |
|
q_0_with_opencore |
中性・内殻スピンあり |
|
q_1_no_opencore |
電荷+1・内殻スピンなし |
|
q_1_with_opencore |
電荷+1・内殻スピンあり |
|
Initial |
初期状態の入力ファイル |
|
q_0 |
中性 |
|
q_0.2 |
電荷+0.2 |
|
q_0.5 |
電荷+0.5 |
|
q_1 |
電荷+1 |
計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.140 Å、c= 5.040 Å ) \(3 \times 3 \times 2\) 倍しました。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
Monkhorst-Pack 2×2×2 |
交換相関相互作用 |
GGAPBE, PAW |
SCF 収束条件 [Ha/atom] |
1.0E-8 |
内殻電子のスピン |
sw_xc_opencore_ae_only = on |
CIAOの入力 |
sw_write_core_full 1 |
(補正項計算で使用する)誘電率 |
※8.23 (ab面内), 9.74 (c軸) |
※https://materialsproject.org/materials/mp-661/
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
\(E_{\rm final} (N-1) - E_{\rm initial} (N)\) |
64.980 |
64.212 |
\(E_{\rm corr}^1\) |
0.327 |
0.327 |
\(E_{\rm VBM}\) ( \(\Delta N=0.2\)) |
6.735 |
|
\(E_g\) |
4.404 |
|
\(E_{\rm B}\) |
74.244 |
73.476 |
実験値 |
73.3 a) |
a) [Mahmood03]
内殻電子のスピンを考慮することにより、実験値に近い結合エネルギーが得られています。 なお、式 (5.103) を用いた場合、 \(E_{\rm B}\) は以下のようになりました。
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
|
76.221 |
75.446 |
fcc Pt結晶 (Pt 4f)
fcc PtN結晶 (Pt 4f) の例は samples/XPS/Pt
以下に配置されています。そのディレクトリー構成は下記の通りです。
第一階層 |
第二階層 |
説明 |
Final |
終状態の入力ファイル |
|
q_0_no_opencore |
中性・内殻スピンなし |
|
q_0_with_opencore |
中性・内殻スピンあり |
|
Initial |
初期状態の入力ファイル |
|
q_0 |
中性 |
計算条件は以下のとおりです。なお、スーパーセルの格子定数は、自動最適化により得られた値 ( a = 3.963 Å )を \(2 \times 2 \times 2\) 倍しました。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
Monkhorst-Pack 2×2×2 |
交換相関相互作用 |
GGAPBE, PAW |
SCF 収束条件 [Ha/atom] |
1.0E-8 |
内殻電子のスピン |
sw_xc_opencore_ae_only = on |
CIAOの入力 |
sw_write_core_full 1 |
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
\(E_{\rm B}\) |
70.673 |
70.380 |
実験値 |
71 a) |
a) http://techdb.podzone.net/xps/index.cgi?element=Pt
内殻電子のスピン考慮なし |
内殻電子のスピン考慮あり |
|
\(E_{\rm B}\) |
74.333 |
74.040 |
実験値 |
71 a) |
a) http://techdb.podzone.net/xps/index.cgi?element=Pt
なお、CIAOにおいて内殻電子の運動エネルギーを非相対論で評価した擬ポテンシャルを用いた場合、4f 軌道の結合エネルギー ( 4f7/2 と 4f5/2 の重みつき平均 ) は 83.052 eV となり、実験値と離れてしまいます。
O 2 分子 (O 1s)
第一階層 |
説明 |
final_spin_antiparallel |
終状態・スピン反平行 |
final_spin_parallel |
終状態・スピン平行 |
initial |
始状態 |
計算条件は以下のとおりです。
平面波カットオフ [Ry] |
30.0 |
電荷密度カットオフ [Ry] |
270.0 |
k 点サンプリング |
|
交換相関相互作用 |
GGAPBE, PAW |
単位胞の1辺 [Å] |
14.0 |
SCF 収束条件 [Ha/atom] |
1.0E-8 |
内殻電子のスピン |
sw_xc_opencore_ae_only = on |
CIAOの入力 |
sw_write_core_full 1 |
(補正後計算で使用する)誘電率 |
※1.00 |
※真空中
S=1/2 |
S=3/2 |
|
spin orientationの指定 |
anti parallel |
parallel |
|
542.000 |
541.297 |
|
1.436 |
1.439 |
|
543.436 |
542.736 |
実験値 a) |
544.2 |
543.1 |
a) https://t-ozaki.issp.u-tokyo.ac.jp/vps_pao_core2019/O/index.html and references therein.
2つのスピン状態とも、実験値に近い値が得られました。