Implementation of the hybrid exchange-correlation functionals in the SIESTA code

本論文は、数値原子軌道のガウス関数近似表現を活用して、大幅に改善されたバンドギャップ予測を可能にする拡張系の大規模かつスケーラブルなシミュレーションを実現する、SIESTA コードにおけるハイブリッド交換相関汎関数の効率的かつ正確な実装を提示する。

要約

原子で構成された都市の振る舞いをシミュレーションしようとしている状況を想像してください。あなたは、原子を結びつけている電子（微小な粒子）がどのように移動し、相互作用するかを知りたいのです。何十年もの間、科学者たちはこれを達成するために**密度汎関数理論（DFT）**と呼ばれるツールを用いてきました。DFTを非常に高速で効率的な地図だと考えてみてください。それは都市の全体的な構造を把握するには優れていますが、盲点があります。それは「エネルギーギャップ」（建物の1 階と 2 階の間の距離）を誤って評価しがちだということです。実際よりもギャップが小さいと判断する傾向があり、その結果、実際には絶縁体である材料が導体のように見えることがあります。

これを修正するために、科学者たちはハイブリッド汎関数を開発しました。これらは、地図を高解像度の衛星画像にアップグレードするようなものです。これらは盲点を補正する特定の「正確な交換」計算を追加し、正しいエネルギーギャップをもたらします。しかし、欠点があります。この高解像度の視点は計算が非常に遅いのです。巨大な都市のすべての車について同時に交通の流れを計算しようとするようなものです。コンピュータは圧倒され、シミュレーションには永遠に時間がかかります。

問題：「4 中心」のボトルネック
ハイブリッド計算がこれほど遅い主な理由は、「4 中心積分」と呼ばれる数学的な問題にあります。部屋にいる 4 人の異なる人々の間の相互作用を計算しようとしている状況を想像してください。1,000 人がいれば、可能な 4 人グループの数は天文学的なものになります。原子の世界では、すべての可能なグループについてこれらの相互作用を計算することが計算上のボトルネックとなります。

解決策：「ガウス」翻訳者
この論文の著者たちは、材料シミュレーション用の人気ソフトウェアであるSIESTAコードと協力して、これを高速化するための巧妙な方法を見つけました。

1. ネイティブ言語（NAO）： SIESTA は通常、「数値原子軌道（NAO）」という言語で話します。これらは、特定の距離で突然終わる厳格で局所化された地図のようなものです。これらは標準的な計算には効率的ですが、ハイブリッド汎関数に必要な複雑な「4 中心」の数学には非常に使いにくいものです。
2. 翻訳（GTO）： チームは翻訳者を作成しました。彼らは、その厳格で局所化された地図（NAO）を取り出し、「ガウス型軌道（GTO）」を用いて近似しました。GTO を数学的に扱いやすい滑らかなベル曲線の形状だと考えてください。
3. ライブラリ（Libint）： GTO は数学的に滑らかであるため、それらの間の相互作用を瞬時に計算できる、事前に最適化された「ライブラリ」（libintと呼ばれます）が既に存在します。これは、4 人間のあらゆる可能な会話のための事前計算された辞書を持っているようなものです。

どのように機能させたか
チームは単に言語を交換しただけではなく、橋を架けました。

• フィッティング： 彼らは数学的に、厳格な SIESTA の地図を滑らかなガウス形状に「適合」させました。これは、元の画像の詳細を失うことなく、高品質なプリンターが処理できるように、ギザギザしたピクセル化された画像を滑らかにするようなものです。
• スクリーニング： 彼らは入り口に「用心棒」を追加しました。ほとんどの原子は十分に離れており、有意に相互作用しないため、コードは遠くのペアを無視します。これにより、計算回数は数十億から管理可能な数百万に削減されます。
• 並列処理の力： 彼らは、数千のコンピュータプロセッサが互いの邪魔をすることなく都市の異なる部分に同時に作業できるシステムを構築しました。

結果：高速かつ高精度
この論文は、この新しい方法を、シリコンチップからグラフェンなどの 2 次元材料まで、さまざまな材料でテストしました。

• 精度： 新しい方法は「盲点」を修正しました。例えば、黒リンが金属ではなく（ギャップを持つ）半導体であると正しく予測し、シリコンとダイヤモンドのエネルギーギャップを実験的な現実とほぼ同じに計算しました。
• 速度： ガウス変換とスクリーニングによる「用心棒」を使用することで、以前は実行に時間がかかりすぎた大規模なシステム（数百、あるいは数千の原子）に対して、これらの高精度計算を実行可能にしました。

トレードオフ
著者たちはまた、速度と精度の間の最適なバランスをどのように得るか分析しました。彼らは以下を発見しました。

• 各原子を表すために中程度の数の「ガウス形状」（約 4 から 6）を使用するのが通常で十分です。
• 相互作用のための特定の「カットオフ」距離を設定することは、すべての遠くの原子を計算する必要なくうまく機能します。
• このバランスにより、科学者は最も高価な方法とほぼ同じ精度の結果を、その数分の一の時間で得ることができます。

まとめ
この論文は、SIESTA ソフトウェアのための新しいエンジンを紹介しています。これにより、科学者たちは、ソフトウェアのネイティブ言語を数学的に滑らかなものに変換し、瞬時に処理できるようにすることで、大規模な材料に対して高精度な「ハイブリッド」シミュレーションを実行できます。これにより、コンピュータが作業を完了するのを数週間待つことなく、半導体や 2 次元シートなどの複雑な材料の電子特性を正確に予測することが可能になります。

技術概要

以下は、論文「siesta コードにおけるハイブリッド交換相関汎関数の実装」の詳細な技術的概要です。

1. 問題提起

密度汎関数理論（DFT）は材料科学の基盤ですが、標準的な半局所汎関数（LDA や GGA/PBE など）は自己相互作用誤差に悩まされています。これにより、バンドギャップの系統的な過小評価、電子の非局在化の過大評価、欠陥準位や反応障壁の記述における不正確さが生じます。
半局所交換と正確なハートリー - フック（HF）交換の一定割合を混合するハイブリッド汎関数（HSE06、PBE0 など）は、これらの問題を解決します。しかし、SIESTAコード内での大規模な拡張系（固体、表面、欠陥）への適用は、歴史的に以下の理由により制限されてきました：

• 計算コスト: 正確な交換に必要な 4 中心電子反発積分（ERI）の評価は、$O(N^4)$ と悪くスケーリングし、計算的に実行不可能です。
• 基底関数の不一致: SIESTA は厳密に局在した**数値原子軌道（NAO）**を使用しており、優れた空間的局在性と疎行列構造を提供します。しかし、解析的な ERI 評価のための標準的な効率的なライブラリ（libint など）は、**ガウス型軌道（GTO）**向けに設計されています。
• 解析的力の欠如: 以前の実装では、完全な解析的力計算が欠けており、幾何構造最適化や分子動力学を妨げていました。

2. 手法

著者らは、NAO と GTO の両方の長所を活用し、SIESTA 内で効率的なハイブリッド汎関数計算を可能にするハイブリッド手法を提示します。

A. NAO から GTO へのフィッティング戦略

数値的に ERI を評価する（遅い）か、NAO 基底を放棄するのではなく、著者らはネイティブ NAO の半径部分をガウスの線形結合として近似します：
$$\chi_\mu(r) \approx \sum_{k=1}^{N_G} C_k \exp(-\alpha_k r^2)$$

• 内部フィッティング: SIESTA 内に新しい堅牢なフィッティングアルゴリズムを実装しました。これは、係数を最小二乗法で解くために指数を固定し、その後指数を最適化する交互最適化戦略を用いて、NAO とガウス展開の差を最小化します。
• 厳密な局在化: ガウス関数は無限遠でゼロにならないため、著者らは切断戦略を導入します。ガウスの和に滑らかなカットオフ関数（$1 - e^{-\alpha(r-r_c)^2}$）を乗算し、軌道が厳密に局在し、切断半径で連続であることを保証することで、SIESTA の疎行列の利点を維持します。
• 角変換: SIESTA で使用される実球面調和関数を、libint ライブラリとインターフェースするために直交ガウス調和関数に変換します。

B. 積分評価とスクリーニング

• 解析的評価: NAO が縮約ガウスとして表現されると、4 中心 ERI は Obara-Saika アルゴリズムを採用する**libint**ライブラリ（バージョン 1.1.4）を用いて解析的に計算されます。
• 多段階スクリーニング: 計算複雑度を $O(N^4)$ から準線形スケーリングに削減するため、著者らは階層的なスクリーニング手法を採用します：
  1. シュワルツ不等式: 積分の大きさを評価し、無視できる対を除外します。
  2. シェル四重項フィルタリング: 軌道をシェルにグループ化し、空間的重なりと距離に基づいてフィルタリングします。
  3. 密度行列スクリーニング: シュワルツの上限に、シェル間を接続する最大密度行列要素を乗算します。積が閾値を下回る場合、積分はスキップされます。
  4. 範囲分離: HSE06 汎関数の短距離性（遮蔽されたクーロンポテンシャル）を利用し、長距離相互作用を自然に切断することで、エワルド総和の必要性を排除します。

C. 解析的力と並列化

• 力: 基底関数が原子位置に依存することから生じるプルアイ項を含む、完全な解析的力を導出します。これにより、幾何構造最適化と第一原理分子動力学が可能になります。
• 並列化: ERI の計算を MPI プロセス間で分配する動的な並列分配方式を使用し、優れた負荷分散とスケーラビリティを確保します。

3. 主な貢献

1. SIESTA における最初の効率的なハイブリッド実装: NAO ベースの SIESTA フレームワークに正確な交換を統合し、コード固有の線形スケーリング能力を犠牲にすることなく成功させました。
2. 堅牢な内部フィッティングアルゴリズム: NAO をガウスにフィッティングするための自己完結型で体系的な手法を開発し、外部フィッティングツールや複雑なワークフローの必要性を排除しました。
3. 解析的力: SIESTA NAO フレームワーク内のハイブリッド汎関数に対する、完全な解析的力定式化を初めて提供し、構造緩和と動力学を可能にしました。
4. 包括的なベンチマーク: 半導体、絶縁体、2 次元材料、分子など、幅広い材料にわたり、PBE、$G_0W_0$、および実験データと比較して実装を検証しました。

4. 結果

論文は、広範なベンチマークを通じて実装を検証しています：

• バンドギャップ: HSE06 汎関数は、PBE と比較してバンドギャップ予測を大幅に改善します。
  • 例: シリコンのバンドギャップは、PBE の 0.59 eV から1.20 eV（HSE06）に改善され、実験値 1.12 eV と一致します。
  • 例: PBE が金属として誤って予測する黒リンは、0.29 eV のギャップを持つ半導体として正しく予測され（$G_0W_0$および実験に近い）、正しい性質を示します。
  • 結果は、Si、ダイヤモンド、LiF、c-BN、TiO$_2$、ZnO、h-BN において、$G_0W_0$計算および実験データと優れた一致を示しています。
• 精度と効率のトレードオフ:
  • ガウスの数（$N_G$）: $N_G$を増やすと精度は向上しますが、コストは劇的に増加します（$O(N_G^4)$）。著者らは、$N_G = 4$が最適なバランスを提供し、計算時間の一部で完全に収束した値から数十 meV 以内のギャップを提供することを見出しました。
  • カットオフ半径（$\epsilon_{Gauss}$）: $5 \times 10^{-3}$の閾値は、疎性を維持しつつ、ほとんどの系で収束したバンドギャップを確保する十分にコンパクトな軌道を提供します。
  • スクリーニング閾値: 積分スクリーニングのための$10^{-5}$から$10^{-6}$の閾値は、フィルタリングされていない計算と比較して大幅な高速化（桁違い）をもたらすほぼ収束した結果を提供します。
• 並列スケーラビリティ: 5$\times$5$\times$5 の SrTiO$_3$スーパーセル（数百原子）でのベンチマークでは、4 コアから 64 コアへのスケーリング時に11.8 倍の高速化を示し、並列効率は 73% 以上を維持しました。
• 分子検証: この手法は O$_2$分子における一重項 - 三重項分裂を正しく再現し、スピン分解交換実装の精度を確認しました。

5. 意義

この研究は、大規模な拡張系の第一原理シミュレーションにおける主要なボトルネックを除去します。SIESTA 内で数百から数千原子を含む系に対するハイブリッド汎関数計算を可能にすることで、著者らは、以前は半局所汎関数やより小規模な系に限定されていた電子構造（バンドギャップ、欠陥準位）および構造特性の正確な予測を、より大規模なスケールで利用可能にしました。

ハイブリッド汎関数を用いた幾何構造最適化と分子動力学の実行能力は、以下の分野の研究に対する新たな道を開きます：

• ポラロン挙動と電荷の局在化。
• 半導体および絶縁体における欠陥物理学。
• ヘテロ構造における表面および界面特性。
• 複雑な材料に対する正確な熱化学的性質。

この実装は、局在基底関数の効率性とハイブリッド汎関数の精度の間のギャップを埋めることで、SIESTA を高精度な材料モデリングのための競争力があり多用途なツールとして位置づけています。