帯電した点欠陥の形成エネルギー図：FNV法（sxdefectalign） による有限サイズ補正

絶縁体・半導体デバイスの特性は、母材中に存在する点欠陥によって大きく左右されます。どの欠陥がどの電荷状態でどれだけ生成するかは、形成エネルギーを Fermi 準位の関数として描いた図から読み取ることができます。ところが第一原理計算でこの図を描こうとすると、帯電した欠陥は周期境界条件のもとで自身の鏡像や中和用の背景電荷と、実在しない静電相互作用をもってしまいます。この誤差はセルの大きさ $L$ に対して $1/L$ でしか減衰せず、実用的なスーパーセルでは欠陥準位をギャップの外に押し出すほどの大きさになります。本解析事例では、第一原理計算ソフトウェア Advance/PHASE と、Freysoldt–Neugebauer–Van de Walle (FNV) 法 のオリジナル実装である sxdefectalign を連携させてこの誤差を補正し、MgO の酸素空孔とマグネシウム空孔について形成エネルギー図を構築しました。補正に必要な誘電率も PHASE で自己完結的に求めています。

Keywords: 第一原理計算 (DFT), 点欠陥, 帯電欠陥, 形成エネルギー, 電荷遷移準位, 有限サイズ補正, FNV法, sxdefectalign, 誘電率, 化学ポテンシャル, 平衡欠陥濃度, F中心, MgO

1. 形成エネルギーと有限サイズ補正

母結晶に電荷 $q$ の欠陥 $X$ を導入したときの形成エネルギーは、次式で定義されます [1]。

$$E^{\rm f}[X^{q}](E_{\rm F}) = E_{\rm tot}[X^{q}] - E_{\rm tot}[\text{bulk}] - \sum_{i} n_{i}\,\mu_{i} + q\,(E_{\rm VBM} + E_{\rm F}) + E_{\rm corr}[q]$$

ここで $n_i$ は欠陥生成の際に加えた原子数（取り除いた場合は負）、$\mu_i$ はその原子種の化学ポテンシャル、$E_{\rm VBM}$ は母結晶の価電子帯上端、$E_{\rm F}$ は VBM を基準にした Fermi 準位です。最後の $E_{\rm corr}[q]$ が有限サイズ補正で、中性欠陥ではゼロ、帯電欠陥では必須になります。

FNV 法では補正を 2 つの寄与に分けます。

$$E_{\rm corr}[q] = E^{\rm PC}(\varepsilon, L) - q\,\Delta V_{q/{\rm b}}$$

• 点電荷項 $E^{\rm PC}$：欠陥の電荷を模したモデル電荷を、孤立系に置いたときと周期系に置いたときの静電エネルギーの差。電荷がセル内で十分局在していれば Madelung 型の表式 $q^{2}\alpha_{\rm M}/(2\varepsilon L)$ に一致します。
• アラインメント項 $-q\,\Delta V_{q/{\rm b}}$：欠陥セルと完全結晶セルとで静電ポテンシャルの基準がずれる分の補正。欠陥から十分離れた領域で、DFT のポテンシャル差からモデル長距離ポテンシャルを差し引いた「短距離成分」が到達する一定値から読み取ります。

本事例では、FNV 法 [2] との連携で、Advance/PHASE を用いた SCF 計算の後処理として出力する静電ポテンシャルを sxdefectalign [3] に読み込ませ、長距離モデルポテンシャルと短距離成分の分離から点電荷項とアラインメント項を求めました。補正に必要なのは SCF 計算の副産物であるこのポテンシャルだけで、追加の第一原理計算は必要ありません。この連携は、PHASE の GUI である PHASE-Viewer から、 FNV法が実装されたpymatgen-analysis-defects [4] を呼び出して実行することもできます [5]。

2. 計算モデルと計算条件

対象は、ワイドギャップ絶縁体の代表である MgO（岩塩型）です。酸素空孔 $V_{\rm O}$（F 中心）とマグネシウム空孔 $V_{\rm Mg}$ を、立方晶 2×2×2 スーパーセル（64 原子）の 1 サイトに導入しました。

図1. 計算モデル。MgO 岩塩型 2×2×2 立方スーパーセル（32 Mg + 32 O）。緑丸の Mg サイトに $V_{\rm Mg}$ を、赤丸の O サイトに空孔 $V_{\rm O}$（F 中心）を導入します。

表1. 計算条件

項目	設定
結晶構造	岩塩型、格子定数 4.2645 Å（PBE 平衡値。実験値 4.2123 Å）
スーパーセル	立方 2×2×2（64 原子）
構造緩和	全原子を緩和（セル格子定数は固定）、力の最大値\(|F|_{\rm max}\) < 5×10−4 hartree/bohr
スピン	スピン分極あり
擬ポテンシャル	PAW
交換相関汎関数	GGA-PBE
カットオフ（波動関数 / 電荷密度）	30 / 270 Ry
k 点	2×2×2 メッシュ
欠陥種と電荷状態	\(V_{\rm O}\) (0, +1, +2)、\(V_{\rm Mg}\) (0, −1, −2)
化学ポテンシャルの参照相	hcp-Mg、O2 分子（三重項）
価電子帯上端 \(E_{\rm VBM}\)	5.37 eV（完全結晶セル。固有値の基準は本文参照）
バンドギャップ	4.26 eV（実験値 7.78 eV）

構造緩和は全電荷状態で行っています。格子定数は実験値ではなく PBE を用いた計算の平衡値を用いました。実験値のまま固定すると、PBE にとっては圧縮された状態なので 6.0 GPa 程度の残留圧力が残り、欠陥の形成体積との積が形成エネルギーに系統誤差として乗るためです。

形成エネルギーの基準となる価電子帯上端 $E_{\rm VBM}$ は、欠陥セルと同じ 64 原子の完全結晶セルの SCF 計算から求めました。MgO の VBM は Γ 点にあり、2×2×2 メッシュには Γ 点が含まれるので、Γ 点の最高占有 Kohn–Sham 準位をそのまま使えます。ただし PHASE の固有値は、局所擬ポテンシャルの $\mathbf{G}=0$ 成分を Kohn–Sham ポテンシャルから除いた基準で出力され、全エネルギーにはこの成分が別途含まれています。電子数の異なるセルの全エネルギー差と整合させるため、固有値にこの定数を足し戻した値を $E_{\rm VBM}$ としました（本計算では 0.27 eV + 5.10 eV = 5.37 eV）。$E_{\rm VBM}$ は母結晶の定数で、欠陥の種類や電荷状態にはよらず、形成エネルギー図の横軸 $E_{\rm F}$ の原点を与えます。

3. 誘電率の決定

補正項の大きさは誘電率 $\varepsilon$ に反比例するため、$\varepsilon$ の選択が結果を大きく左右します。PHASE では、SCF 計算の後処理として電子系の誘電テンソル $\varepsilon_{\infty}$ を、Berry 位相にもとづく計算から Born 有効電荷 $Z^{*}$ を、さらに Γ 点のフォノン振動数を求めることができ、これらから Lyddane–Sachs–Teller の関係で静的誘電率 $\varepsilon_{0}$ を組み立てられます。外部から実験値を持ち込む必要はありません。

$$\varepsilon_{0} = \varepsilon_{\infty} + \dfrac{4\pi (Z^{*})^{2}}{\Omega\,\mu\,\omega_{\rm TO}^{2}} = \varepsilon_{\infty}\left(\dfrac{\omega_{\rm LO}}{\omega_{\rm TO}}\right)^{2}$$

表2. PHASE で算出した MgO の誘電応答

量	PHASE	実験値
\(\varepsilon_{\infty}\)（電子系）	3.10	2.96
\(Z^{*}\)（陽イオン）	1.971	1.97
\(\omega_{\rm TO}\) (cm−1)	369.7	401
\(\omega_{\rm LO}\) (cm−1)	682.9	718
\(\varepsilon_{0}\)（LST, 静的）	10.57	9.83

どちらの誘電率を使うかは、欠陥電荷を遮蔽しているのが何かで決まります。原子位置を固定した計算では遮蔽の担い手は電子系だけなので $\varepsilon_{\infty}$、原子緩和を許すと格子分極も加わるので $\varepsilon_{0}$ を使います [1]。本事例は構造緩和を行っているので $\varepsilon_{0}$ = 10.57 を用いました。MgO では両者が 3 倍以上違うため、この区別を誤ると補正量が 3 倍ずれることになります。

4. 静電ポテンシャルのアラインメント

補正の要は、欠陥セルと完全結晶セルの静電ポテンシャル差からモデル長距離成分を差し引いた短距離ポテンシャル $V^{\rm sr}$ です。モデルが妥当であれば、$V^{\rm sr}$ は欠陥近傍でのみ有限で、離れた領域では一定値 $C$ に収束します。なお、以下のポテンシャルはいずれも電子のポテンシャルエネルギーとして表示します。正の欠陥電荷は電子を引きつけるので、$V_{\rm O}^{2+}$ のモデルポテンシャルは欠陥位置で谷になります。

この手順がどれだけうまく働くかを、まず原子を完全結晶の位置に固定した計算で確かめました（図2 上）。この場合の遮蔽は電子系だけなので $\varepsilon_{\infty}$ = 3.10 を用います。欠陥から離れた領域では、DFT のポテンシャル差（青）がモデル（橙）に一定値を足しただけの形で重なり、短距離成分（緑）は $C$ = −0.30 eV で平坦になります。読み取り領域でのばらつき（標準偏差）は 0.01 eV です。

ただし構造緩和を行った場合は、動いたイオンの芯ポテンシャルが欠陥セルと完全結晶セルで打ち消し合わないため、平面平均ポテンシャルが原子層ごとに大きく振動します。64 原子セルでは欠陥から最も遠い面でも原子がわずかに動いているため、そのままでは平坦な領域が現れません。そこで sxdefectalign のマニュアル [3] が同じ MgO の空孔を例に推奨しているとおり、原子層の間隔（[100] 方向に 4.03 bohr）で窓平均をとって振動を平滑化し、欠陥から離れた領域での平均値として $C$ を求めました（図2 下）。

図2. MgO $V_{\rm O}^{2+}$（64 原子セル）の [100] 方向平面平均ポテンシャル。上：原子を完全結晶の位置に固定した場合（$\varepsilon_{\infty}$ = 3.10）、下：構造緩和後（$\varepsilon_{0}$ = 10.57、原子層間隔 4.03 bohr で窓平均）。青：DFT のポテンシャル差、橙：モデル長距離ポテンシャル、緑：両者の差である短距離成分。灰色の帯が $C$ の読み取り領域（欠陥中心から 0.30$L$ 以上）、破線がそこでの平均値 $C$ です。ポテンシャルは電子のポテンシャルエネルギーで表示しているため、正の欠陥電荷では欠陥位置が谷になります。

窓平均をとっても、図2 下のように ±0.1 eV 程度の起伏は残ります。この起伏は DFT のポテンシャル差（青）にそのまま現れている原子層ごとの振動で、モデル（橙）は読み取り領域で滑らかに変化するだけです。読み取り領域での短距離成分のばらつき（標準偏差）は 0.03〜0.07 eV で、窓平均をとらない場合（0.09〜0.20 eV）の約 3 分の 1 です。一方、$C$ の値そのものは窓平均の有無で 0.03 eV 未満しか変わらず、読み取り領域の境界を 0.25$L$〜0.35$L$ の間で変えても電荷遷移準位の変化は 0.01 eV 以内です。窓平均は値を動かすのではなく読み取りの不確かさを減らすものといえます。残る起伏は 64 原子という小さなセルの限界で、セルを大きくするほど平坦な領域が広がります。なお sxdefectalign は欠陥の電荷を「過剰電子数」$-q$ で数える規約なので、プラトー $C$ からアラインメント項は $+q\,C$ と書けます。表3 にはこの符号でまとめています。

5. 補正はどれだけ効くか

表3. 64 原子セルにおける有限サイズ補正

欠陥	\(q\)	点電荷項 (eV)	プラトー \(C\) (eV)	アラインメント項 \(+qC\) (eV)	合計 \(E_{\rm corr}\) (eV)
VO+	+1	0.227	−0.119	−0.119	0.108
VO2+	+2	0.906	−0.124	−0.249	0.657
VMg−	−1	0.227	−0.094	+0.094	0.320
VMg2−	−2	0.906	−0.042	+0.084	0.991

図3. 補正項の内訳。点電荷項が支配的で、電荷の 2 乗に比例するため多価の電荷状態ほど強く効きます。アラインメント項は小さいものの、無視できる大きさではありません。

補正量は最大 0.99 eV です。数値としては控えめに見えますが、この 1 eV 弱が欠陥準位の位置を決定的に左右します。補正を入れない場合と入れた場合で電荷遷移準位を比べると、次のようになります。

表4. 有限サイズ補正の有無による電荷遷移準位の違い（VBM からの eV）

欠陥	遷移	補正なし	補正あり	効果
VO	ε(2+/+)	1.91	1.36	0.55 eV 下がる
VO	ε(+/0)	2.58	2.47	0.11 eV 下がる
VMg	ε(0/−)	−0.27	0.05	価電子帯の下からギャップ内へ
VMg	ε(−/2−)	−0.10	0.57	価電子帯の下からギャップ内へ

補正を入れないと、$V_{\rm Mg}$ の 2 つの遷移準位はどちらも価電子帯上端より下に沈みます。この場合 $V_{\rm Mg}$ はギャップ内に準位をもたず、Fermi 準位によらず常に 2 価の負イオンとして存在することになり、電荷状態が変わるという描像そのものが失われます。補正を入れて初めて両準位がギャップ内に現れ、Fermi 準位に応じて 0 / −1 / −2 と電荷状態を変える二重アクセプタとしての振る舞いが得られます。

$V_{\rm O}$ でも、$\varepsilon$(2+/+) が補正によって 0.55 eV 下がり、ギャップに対する相対位置が 45 % から 32 % に移ります。8 節で見るように、この補正後の位置がハイブリッド汎関数による文献値とよく一致します。

欠陥濃度で見ると、補正の影響はさらにはっきりします。濃度は形成エネルギーに対して指数関数的に変わり、1500 K では 1 eV の違いが約 3.4 桁に相当するためです。9 節で述べる電荷中性条件を補正なしで解くと、O 過剰条件では Fermi 準位が 0.64 eV から 0.32 eV に下がり、支配的な欠陥である Mg 空孔の濃度は 1.0 桁、少数欠陥である酸素空孔の濃度は 4.4 桁、それぞれ過大に見積もられます。

6. 化学ポテンシャル

形成エネルギーは化学ポテンシャル $\mu_i$ に依存します。母結晶との平衡条件と、各元素の単体相が析出しない条件から、$\mu_i$ の可動範囲は生成エンタルピー $\Delta H_{\rm f}$ の幅に限られます。本事例では両端（Mg 過剰／O 過剰）の 2 条件で図を描きました。実験的には酸素分圧を変えることに対応します。

本計算の $\Delta H_{\rm f}$ は −6.08 eV/f.u. で、実験値 −6.23 eV/f.u. [6] との差は 0.15 eV に収まりました。ここでは O 側の基準となる O₂ 分子のエネルギーに補正を加えています。GGA では O₂ 分子の結合エネルギーの誤差などにより、酸化物の生成エンタルピーの大きさが系統的に過小評価されます。Wang ら [7] は、非遷移金属酸化物の生成エンタルピーへのフィットから O₂ 分子あたり 1.36 eV（$\mu_{\rm O}$ に換算して +0.68 eV/原子）の定数補正を提唱しており、補正前の −5.40 eV/f.u.（実験値との差 0.83 eV）がこれによって改善されます。

7. 形成エネルギー図と電荷遷移準位

形成エネルギーは電荷状態 $q$ ごとに $E_{\rm F}$ の一次関数になり、その傾きが $q$ に等しくなります。ある $E_{\rm F}$ で実際に観測されるのは、そのうち最もエネルギーの低い電荷状態です。したがって形成エネルギー図とは、各電荷状態の直線群の下側包絡線（最小値）を取ったものにほかなりません。図4 では、細い線で各電荷状態の直線を、太い線でその最小値を取った包絡線を示しています。包絡線の折れ曲がり点（●）が電荷遷移準位に対応し、そこで安定な電荷状態が入れ替わります。

図4. MgO の欠陥形成エネルギー図（64 原子セル、構造緩和あり、$\varepsilon_{0}$ = 10.57）。上：Mg 過剰（O 欠乏）条件、下：O 過剰条件。横軸は PBE のバンドギャップ 4.26 eV の範囲。破線は形成エネルギーがゼロの位置で、これより下は欠陥が自発的に生成することを意味します。

2 つの条件で主役となる欠陥が入れ替わるのが見て取れます。Mg 過剰条件では $V_{\rm O}$（ドナー）が支配的で、$E_{\rm F}$ が 1.32 eV より下では形成エネルギーが負になり、自発的に生成します。O 過剰条件では逆に $V_{\rm Mg}$（アクセプタ）が支配的になり、$E_{\rm F}$ が 1.29 eV より上で形成エネルギーが負になります。Fermi 準位をバンド端に近づけようとすると、それを打ち消す向きの欠陥が自発的に生じるわけで、これが酸化物における自己補償の典型的な構図です。

表5. 電荷遷移準位（VBM からの eV）

欠陥	遷移	本計算	ギャップに対する位置
VO	ε(2+/+)	1.36	32 %
VO	ε(+/0)	2.47	58 %
VMg	ε(0/−)	0.05	1 %
VMg	ε(−/2−)	0.57	13 %

8. 文献値との比較

得られた結果が妥当かどうかを、学術誌に公表された計算値・実験値と照らし合わせて確認しました。

誘電応答。$Z^{*}$ = 1.971 は実験値 1.97 とほぼ完全に一致します。$\varepsilon_{\infty}$ と $\varepsilon_{0}$ も実験値を 5〜8 % で再現しており、補正の入力として十分な精度です。PBE 平衡体積ではフォノンがやや軟化するため $\varepsilon_{0}$ は実験値より大きめに出ますが、これは PBE の既知の性質です。

生成エンタルピー。−6.08 eV/f.u. は実験値 −6.23 eV/f.u. [6] と 0.15 eV の差です。形成エネルギーの基準を決める化学ポテンシャルが信頼できる精度で得られていることになります。

バンドギャップ。4.26 eV は実験値 7.78 eV の 55 % にとどまります。これは PBE 汎関数に共通する既知の過小評価で、Ertekin らの PBE 計算でも 4.83 eV [8] と、実験値を大きく下回っています。以下の電荷遷移準位の比較では、この点を踏まえてギャップに対する相対位置で比べます。

電荷遷移準位。本事例の中心的な結果です。比較の相手には、本計算と同じ 64 原子セルで構造緩和まで行ったハイブリッド汎関数（HSE06 型、厳密交換の混合率 0.35）による計算 [9] を選びました。図5 に、両者の準位を、それぞれのバンドギャップを同じ高さに揃えて並べています。

図5. MgO の電荷遷移準位。左：本計算（PBE + FNV 補正、構造緩和あり、$E_g$ = 4.26 eV）、右：HSE06 による文献値 [9]（$E_g$ = 7.92 eV）。バンドギャップの幅を揃えて描いているため、準位のギャップ内での相対位置（括弧内の %）が直接比較できます。$V_{\rm Mg}$ は文献 [9] では報告されていません。

$V_{\rm O}$ の 2 つの遷移準位は、ギャップに対する相対位置で 58 % と 32 % となり、HSE06 の 54 % と 35 % [9] とよく一致します。絶対値（2.47 eV と 1.36 eV）が HSE06 の 4.26 eV と 2.74 eV より低いのは、PBE のギャップが実験の 55 % しかないためで、ギャップ幅の違いを除けば描像は同じです。

まとめると、補正の入力となる誘電応答は実験値を 8 % 以内で、化学ポテンシャルの基準となる生成エンタルピーは 0.15 eV 以内で再現しています。電荷遷移準位については、有限サイズ補正を正しく入れることで、PBE でもハイブリッド汎関数と同じ描像に到達できることが示されました。残る差はバンドギャップの過小評価という PBE の既知の性質に帰着し、手法の枠組み自体には問題がありません。

9. 考察：欠陥平衡と Fermi 準位

形成エネルギー図から一歩進めて、電荷中性条件（自由キャリアと帯電欠陥の総和がゼロ）を解けば、ある温度・化学ポテンシャル条件での自己無撞着な Fermi 準位と平衡欠陥濃度が得られます。1500 K（代表的な熱処理温度）として解いた結果を表6 に示します。自由キャリア濃度には放物バンド近似と代表的な有効質量（$m_e^{*}$ = 0.35、$m_h^{*}$ = 2.0）を用いました。

表6. 電荷中性条件から求めた平衡 Fermi 準位と欠陥濃度

化学ポテンシャル条件	T (K)	EF (eV, VBM 基準)	平衡欠陥濃度 (cm−3)
Mg 過剰（O 欠乏）	1500	3.20	VO 4.9×1018 VMg 1.2×1015
O 過剰	1500	0.64	VO 6.6×106 VMg 3.4×1018

Mg 過剰条件では酸素空孔が 10¹⁸ cm⁻³ 台に達し、Fermi 準位はギャップ中央より上の 3.20 eV に押し上げられます。この位置は $V_{\rm O}$ の $\varepsilon$(+/0) より上なので、酸素空孔のほとんどは電子 2 個を捕獲した中性の F 中心として存在します。O 過剰条件では逆に Mg 空孔が支配的になり、その負電荷を正孔が補償する形で Fermi 準位は 0.64 eV まで下がります。すなわち酸素分圧を変えるだけで Fermi 準位を 2.56 eV 動かせることになり、支配的な欠陥の種類も電荷状態も入れ替わります。

10. まとめ

第一原理計算ソフトウェア Advance/PHASE と FNV 法のオリジナル実装である sxdefectalign を連携させ、MgO の点欠陥（空孔V）について有限サイズ補正込みの形成エネルギー図を構築しました。補正は PHASE が後処理で出力する静電ポテンシャルだけから得られ、追加の第一原理計算を必要としません。補正に必要な誘電率も、PHASE の誘電応答計算と Born 有効電荷・Γ 点フォノンから自己完結的に決定できます。64 原子セルにおける補正量は最大 0.99 eV で、$V_{\rm Mg}$ の 2 つのアクセプタ準位が価電子帯の下に沈むか、ギャップ内に正しく現れるかを分けるほどの大きさでした。補正を入れた結果、$V_{\rm O}$ の電荷遷移準位はギャップに対する相対位置でハイブリッド汎関数による文献値とよく一致し、$V_{\rm O}$ が二重ドナー、$V_{\rm Mg}$ が二重アクセプタという正しい描像が得られました。本手法は、パワーデバイス材料のキャリア補償、酸化物のイオン伝導や色中心、半絶縁性基板の Fermi 準位制御など、点欠陥が性能を支配する材料設計に広く適用できます。

本解析の詳細や、研究への適用可能性に関するご相談はこちら

お問い合わせ

参考文献

1. C. Freysoldt, B. Grabowski, T. Hickel, J. Neugebauer, G. Kresse, A. Janotti, and C. G. Van de Walle, "First-principles calculations for point defects in solids", Rev. Mod. Phys. 86, 253 (2014).
2. C. Freysoldt, J. Neugebauer, and C. G. Van de Walle, "Fully Ab Initio Finite-Size Corrections for Charged-Defect Supercell Calculations", Phys. Rev. Lett. 102, 016402 (2009).
3. sxdefectalign (S/PHI/nX add-on), C. Freysoldt, https://sxrepo.mpie.de/.
4. https://github.com/materialsproject/pymatgen-analysis-defects
5. 帯電した点欠陥の形成エネルギー補正計算の解析事例 (FNV法 & eFNV法)
6. D. R. Lide (ed.), CRC Handbook of Chemistry and Physics, 84th ed. (CRC Press, 2003).
7. L. Wang, T. Maxisch, and G. Ceder, "Oxidation energies of transition metal oxides within the GGA+U framework", Phys. Rev. B 73, 195107 (2006).
8. E. Ertekin, L. K. Wagner, and J. C. Grossman, "Point-defect optical transitions and thermal ionization energies from quantum Monte Carlo methods: Application to the F-center defect in MgO", Phys. Rev. B 87, 155210 (2013).
9. A. G. Marinopoulos, "First-Principles Study of the Formation and Stability of the Interstitial and Substitutional Hydrogen Impurity in Magnesium Oxide", Condens. Matter 11, 2 (2026).

関連ページ

• 第一原理計算ソフトウェア Advance/PHASE
• 解析分野：ナノ・バイオ
• 産業分野：材料・化学