ハイブリッド汎関数（HSE06）計算の高速化：圧縮交換演算子（ACE）法#

PHASE logo

半導体・絶縁体のバンドギャップは、材料の電気的・光学的性質を決める最も基本的な物性値の一つです。第一原理計算で標準的に用いられる一般化勾配近似（GGA）はバンドギャップを大きく過小評価しますが、Hartree-Fock型の厳密交換を一部取り入れたハイブリッド汎関数HSE06を用いると、実験値に近いバンドギャップが得られます。その一方で、HSE06では非局所的な厳密交換（Fock交換）項の計算が必要になり、計算時間はGGA計算の数百倍以上に達することもあります。この計算コストが、大きな系や多数の材料を対象とした計算にハイブリッド汎関数を適用する際の大きな障害になっていました。第一原理計算ソフトウェアAdvance/PHASE (ver.5.0)では、この厳密交換の計算を高速化するため、圧縮交換演算子（ACE: Adaptively Compressed Exchange）法を実装しました。本解析では、代表的な半導体・絶縁体6種類を対象にHSE06計算を行い、バンドギャップなどの計算結果が従来法と一致したまま、計算時間が約5～10分の1に短縮されることを示します。

Keywords: 第一原理計算 (DFT), ハイブリッド汎関数, HSE06, 厳密交換, 圧縮交換演算子 (ACE), 高速化, バンドギャップ, 半導体

ハイブリッド汎関数計算のボトルネックとACE法#

HSE06汎関数では、電子間の交換相互作用を遮蔽パラメータ $\omega$ で短距離（SR）成分と長距離（LR）成分に分け、短距離部分の交換エネルギーの25%をHartree-Fock型の厳密交換で置き換えます [1,2]。

$E_\text{xc}^\text{HSE06} = \alpha E_\text{x}^\text{HF,SR}(\omega) + (1-\alpha) E_\text{x}^\text{PBE,SR}(\omega) + E_\text{x}^\text{PBE,LR}(\omega) + E_\text{c}^\text{PBE},\\ \quad \alpha = 0.25,\ \omega = 0.106\ \text{bohr}^{-1}$

厳密交換の部分は、各電子の波動関数 $\psi_n$ に次の非局所演算子を作用させる形でハミルトニアンに現れます。

$V_\text{x}\psi_n(\mathbf{r}) = -\alpha \sum_{m}^{\text{occ}} \sum_{\mathbf{q}} \phi_m(\mathbf{r}) \int d\mathbf{r}'\, \phi_m^*(\mathbf{r}')\, \psi_n(\mathbf{r}')\, v_\text{SR}(|\mathbf{r}-\mathbf{r}'|)$

平面波基底の計算では、この作用1回ごとに、すべての占有軌道 $\phi_m$ とすべてのk点（q点）の組について高速フーリエ変換（FFT）を用いた積分を行う必要があります。自己無撞着場（SCF）計算では、この作用を全バンド・ソルバーの各ステップ・各SCF反復で繰り返すため、計算時間の大部分がこの部分に費やされます。

ACE法 [3,4] は、この交換演算子を「圧縮」して再利用する方法です（図1）。まず、占有軌道 $\phi_i$ に対してだけ $W_i = V_\text{x}\phi_i$ を計算し、行列 $M_{ij} = \langle\phi_i|W_j\rangle$ の固有値分解から少数の射影ベクトル $\xi_a$ を作ります。以降の計算では、交換演算子を $V_\text{x}^\text{ACE} = -\sum_a |\xi_a\rangle\langle\xi_a|$ で置き換えます。この演算子の作用は内積と行列積だけで計算でき、FFTを必要としません。ACE演算子は占有軌道に対して元の演算子と厳密に一致する（ $V_\text{x}^\text{ACE}\phi_i = V_\text{x}\phi_i$ ）ため、SCF計算が収束した時点の電子状態（波動関数や電荷密度）は従来法と同じになります。
Advance/PHASEでは、ACE演算子を数回のSCF反復ごとに作り直す二重ループ方式を採用しており、重いFFT計算の回数を大幅に削減しています。

従来法とACE法の比較
図1. 厳密交換項の計算方法の比較。(a) 従来法では、試行波動関数1本ごとに全占有軌道・全q点のFFT計算が必要です。(b) ACE法では、占有軌道に対する交換演算子を少数の射影ベクトルに圧縮し、SCF反復の中では安価な行列積だけで交換項を計算します。下段は、後述するSiの計算で実測した1 SCF反復あたりの平均計算時間（1コアあたり）を、従来法を1とした相対値で示したものです。

計算モデルと計算条件#

ACE法の効果を確認するため、解析事例 [5] でも取り上げた代表的な半導体（Si、ダイヤモンド、3C-SiC、AlP、AlAs）に、ワイドギャップ酸化物のMgOを加えた6物質を対象としました（表1(a)、(b)）。いずれも実験の格子定数を用いた基本単位胞を使用しました。計算は、PBE汎関数でSCF計算を行ったあと、その波動関数と電荷密度を初期値としてHSE06のSCF計算を行う、標準的な2段階の手順で実施しました。HSE06のSCF計算は、同じ入力で従来法とACE法の2通りを実行し、計算結果と計算時間を比較しました。バンドギャップは、ACE法のSCF計算で得られた電荷密度と占有軌道を用いて、L–Γ–X線上のバンド構造計算（固定電荷計算）を行って求めています。

表1(a) 対象物質

物質	結晶構造	格子定数 (Å)	擬ポテンシャル
Si	ダイヤモンド	5.431	NC
C（ダイヤモンド）	ダイヤモンド	3.567	US
3C-SiC	閃亜鉛鉱	4.358	NC（Si）/ US（C）
AlP	閃亜鉛鉱	5.463	NC
AlAs	閃亜鉛鉱	5.661	NC
MgO	岩塩	4.212	NC（Mg）/ US（O）

表1(b) 計算条件

項目	設定
交換相関汎関数	PBE（初期計算）、HSE06
波動関数のカットオフエネルギー	25 Rydberg
k点（交換項のq点）	6×6×6（Γ点を含むメッシュ）
SCF収束条件	全エネルギー変化 1×10^-9 Hartree/原子（3回連続）
計算時間の比較	1コアあたりの計算時間（従来法を1とした相対値）

NC：ノルム保存型擬ポテンシャル、US：ウルトラソフト擬ポテンシャル。

計算結果と考察#

精度：従来法と同じ結果が得られる#

表2に、PBEとHSE06（ACE法）によるバンドギャップと実験値、および従来法とACE法の差を示します。HSE06を用いることで、PBEで大きく過小評価されていたバンドギャップが、いずれの物質でも実験値に近づくことが確認できます（MgOのようなイオン結晶では、HSE06でもなお実験値より小さくなることが知られています）。ACE法の結果から求めた、バンドギャップを決める価電子帯上端と伝導帯下端のエネルギー（Γ点・X点・L点で比較）は、すべての物質で従来法と0.058 meV以内で一致しました。全エネルギーの差は、ノルム保存型擬ポテンシャルのみを用いた物質（Si、AlP、AlAs）では最大でも1原子あたり2×10^-9 Hartreeで、SCF計算の収束判定値（1原子あたり1×10^-9 Hartree）と同程度でした。ウルトラソフト擬ポテンシャルを含む物質（C、3C-SiC、MgO）でも、差は1原子あたり最大3×10^-8 Hartree（約0.8 μeV）と十分に小さい値でした。ACE法は計算精度を損なわないことがわかります。

表2. バンドギャップ（eV）と、従来法とACE法の差

物質	PBE	HSE06 （ACE法）	実験値	全エネルギーの差 (Hartree/単位胞)	バンド端エネルギーの差 (meV)
Si	0.59	1.17	1.17 [6]	2.7×10^-9	0.008
C（ダイヤモンド）	4.01	5.17	5.50 [6]	6.0×10^-8	0.014
3C-SiC	1.30	2.18	2.36 [6]	3.7×10^-8	0.010
AlP	1.60	2.31	2.45 [6]	1.0×10^-9	0.008
AlAs	1.45	2.06	2.16 [6]	4.5×10^-9	0.003
MgO	4.49	6.28	7.83 [7]	1.3×10^-8	0.058

※バンドギャップは、L–Γ–X線上のバンド構造計算から求めた値です（擬ポテンシャルやk点などの計算条件が異なるため、解析事例 [5] の値とは最大0.09 eV異なります）。差はいずれも |ACE法 − 従来法| です。バンド端エネルギーの差は、Γ点・X点・L点それぞれの価電子帯上端と、最も低い伝導帯下端（バンドギャップを決める伝導帯下端）のエネルギーについて、ACE法のSCF計算に続くバンド構造計算の値と、従来法のSCF計算の値を比べた最大値です。

SiのバンドPBE/HSE06
図2. ACE法によるHSE06計算の結果を用いたSiのバンド構造（青の実線）と、PBEによるバンド構造（灰色の破線）。価電子帯上端をエネルギーの原点としています。HSE06では伝導帯が押し上げられ、間接バンドギャップが実験値（1.17 eV）に近づきます。

計算時間：5～10倍の高速化#

図3と表3に、HSE06のSCF計算に要した1コアあたりの時間を、従来法を1とした相対値で示します。ACE法を用いることで、すべての物質で計算時間が大幅に短縮され、全体で5.1～9.9倍、SCF反復部分だけでみると6.5～23.1倍の高速化が得られました。ACE法ではSCF反復の回数が33～69回と、従来法（10～22回）より増えますが、1反復あたりの計算時間が従来法の約35～55分の1に減るため、全体として大きく短縮されます。一方、初期化と収束後の処理（原子間力の計算と、ACE法での固有値・交換エネルギーの再計算）には、ACE法では短縮されない交換項の計算が含まれるため、全体の高速化率は、どの物質でもSCF反復部分の高速化率より小さくなります。特にウルトラソフト擬ポテンシャルを含む物質（C、3C-SiC、MgO）では、これらの時間がACE法の計算時間の約6～8割を占めるため、その差が大きくなっています。

計算時間の比較
図3. HSE06のSCF計算に要した1コアあたりの時間の比較（従来法を1とした相対値。橙：従来法、青：ACE法）。青い棒の上の数字はACE法の計算時間の相対値、各物質の上の太字の数字は高速化率（従来法の計算時間 ÷ ACE法の計算時間）です。

表3. HSE06のSCF反復回数と計算時間（1コアあたり、従来法を1とした相対値）

物質	SCF反復回数（従来法 → ACE法）	ACE法の計算時間（従来法 = 1）	高速化率（全体）	高速化率（SCF反復部分）
Si	10 → 69	0.198	5.1 倍	6.5 倍
C（ダイヤモンド）	20 → 65	0.177	5.7 倍	14.4 倍
3C-SiC	20 → 37	0.161	6.2 倍	23.1 倍
AlP	16 → 57	0.102	9.8 倍	13.1 倍
AlAs	13 → 33	0.101	9.9 倍	14.7 倍
MgO	22 → 69	0.135	7.4 倍	18.3 倍

高速化率は、従来法の計算時間 ÷ ACE法の計算時間です。SCF反復部分：全体の時間から、初期化（交換項の初期計算や擬ポテンシャルの準備など。ACE法では最初のACE演算子の作成を含む）と、収束後の処理（原子間力の計算。ACE法では固有値と交換エネルギーの再計算を含む）を除いた時間。

図4は、Siについて、SCF反復ごとの全エネルギーの変化を経過時間に対して示したものです。ACE法では、数反復ごとに交換演算子を作り直すたびにエネルギーの変化が一時的に大きくなるため、反復回数は従来法の10回から69回に増えました。それでも、1反復あたりの計算時間が平均で従来法の約45分の1と短いため、収束後の処理を含めても、従来法の計算時間を1とすると0.198で計算が終了しました（5.1倍の高速化）。

SCF収束の比較
図4. SiのHSE06計算における、SCF反復ごとの全エネルギーの変化と経過時間の関係（経過時間は1コアあたりの時間で、従来法の計算時間を1とした相対値）。点線は収束判定値（1原子あたり1×10^-9 Hartree、2原子の単位胞では2×10^-9 Hartree）で、全エネルギーの変化がこの値を3回連続で下回ると収束とみなします。ACE法の◆は、収束後の処理（固有値と交換エネルギーの再計算）を含めた計算の終了時刻です。

ACE法の効果が大きくなる計算#

ACE法が短縮するのは、SCF反復の中で交換演算子を波動関数に作用させる計算です。従来法では、この計算に必要なFFTの回数が（バンド数）×（占有軌道数）×（k点数）×（q点数）にほぼ比例し、FFT 1回の計算量も単位胞が大きいほど増えます。このため、k点メッシュを変えずに原子数（電子数）を増やした場合や、同じ系でk点・q点のメッシュを細かくした場合には、この計算の計算量がほかの部分よりも速く増え、全体の計算時間に占める割合が大きくなります（単位胞を大きくするときにk点の密度を一定に保つ場合は、k点数とq点数が減るため、この割合は必ずしも大きくなりません）。したがって、ほかの条件が同じであれば、このような計算ほどACE法による高速化率は大きくなります。

ただし、ACE演算子の作成にも同じ種類のFFT計算が必要なため、高速化率には、ACE演算子を作り直す頻度やSCF反復回数の増え方などで決まる上限があります。従来法で、交換演算子を作用させる計算が全体の計算時間の大部分を占める場合には、高速化率はその上限に近い値になります。

まとめ#

本解析では、第一原理計算ソフトウェアAdvance/PHASEに実装した圧縮交換演算子（ACE）法を用いて、ハイブリッド汎関数HSE06計算の高速化を検証しました。代表的な半導体・絶縁体6物質について、バンドギャップなどの計算結果が従来法と一致したまま、計算を5.1～9.9倍（SCF反復部分では6.5～23.1倍）高速化できることを示しました。入力ファイルに1行を追加するだけで利用できるため、これまで計算コストの面で適用が難しかった大きな系や、多数の材料を対象としたバンドギャップの高精度予測にも、ハイブリッド汎関数を活用しやすくなります。

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

お問い合わせ

参考文献#

J. Heyd, G. E. Scuseria, and M. Ernzerhof, "Hybrid functionals based on a screened Coulomb potential", J. Chem. Phys. 118, 8207 (2003); Erratum: J. Chem. Phys. 124, 219906 (2006).
A. V. Krukau, O. A. Vydrov, A. F. Izmaylov, and G. E. Scuseria, "Influence of the exchange screening parameter on the performance of screened hybrid functionals", J. Chem. Phys. 125, 224106 (2006).
L. Lin, "Adaptively Compressed Exchange Operator", J. Chem. Theory Comput. 12, 2242 (2016).
I. Carnimeo, S. Baroni, and P. Giannozzi, "Fast hybrid density-functional computations using plane-wave basis sets", Electron. Struct. 1, 015009 (2019).
ハイブリッド汎関数法を用いた半導体のバンドギャップ計算
M. Levinstein, S. Rumyantsev, and M. Shur, Handbook Series on Semiconductor Parameters, World Scientific (1999).
R. C. Whited, C. J. Flaten, and W. C. Walker, "Exciton thermoreflectance of MgO and CaO", Solid State Commun. 13, 1903 (1973).