GPR_calculator: An On-the-Fly Surrogate Model to Accelerate Massive Nudged Elastic Band Calculations

技術要約：GPR Calculator

問題提起
原子論的シミュレーション、特に幾何構造最適化、障壁計算、および分子動力学を伴うものには、コンパクトな位相空間内での原子構成の広範なサンプリングが必要となる。密度汎関数理論（DFT）は高い精度を提供するが、その計算コストは、表面拡散や反応のための能動的（NEB）計算のような大規模なシミュレーションにおいては、極めて高くなる。典型的なNEBシミュレーションでは、複数のイメージ（構成）に対してエネルギーと力の評価を必要とする数百の最適化ステップが含まれる。各DFT評価に数十から数百のCPU分（minutes）を要する場合があるため、単一のNEBシミュレーションの完了に数日を要することもある。さらに、表面化学反応は多くの場合、複数の素反応ステップを伴い、それがコストをさらに増大させる。ニューラルネットワークや「ワンショット」のガウス過程回帰（GPR）モデルといった既存の機械学習力場（MLFF）アプローチは、インタイム（実行時）の検証メカニズムを欠いていることが多く、あるいは、探索的なタスクにおいて高い精度を要求する動的なデータ収集手順と統合されていないことが多く、そのため、高度な精度を必要とする作業への適合性が制限されている。

手法
著者らは、高価な電子状態計算を近似するために、ガウス過程回帰（GPR）を用いた「オンザフライ（on-the-fly）」サロゲートモデリング戦略を採用した、Atomic Simulation Environment (ASE) のアドオンモジュールとして設計された、PythonとC++のハイブリッドパッケージであるGPR calculatorを提案している。

• コアアルゴリズム: 本手法は、ベースとなる計算機（グラウンドトゥルースとなるDFTデータを提供する）とサロゲートGPRモデルからなる、ハイブリッド計算機アーキテクチャを採用している。NEB最適化の間、GPRモデルは原子構成のエネルギーと力を予測する。決定的なのは、モデルが不確実性定量化を行う点である。特定のイメージに対する予測不確実性が、ユーザー定義の閾値（例：エネルギーについては1構造あたり0.05 eV、力については0.075 eV/Å）を超えた場合、正確なエネルギーと力を計算するためにベースとなる計算機が呼び出される。この新しいデータポイントは、動的にGPRモデルを更新するために使用される。
• 記述子とカーネル: システムは、並進、回転、および置換対称性を保証するために、球関数展開係数のパワースペクトルから導出されたSO(3)記述子を利用して局所的な原子環境を特徴付ける。GPRモデルは、径方向基底関数（RBF）およびドット積カーネルを採用している。共分散行列は、異なる構造間でのエネルギーと力の相関を扱うために、特定のカーネル関数（$k_{ee}$, $k_{ef}$, $k_{ff}$）を利用して構築される。
• 実装と最適化: GPRの$O(N^3)$の計算複雑性とPythonにおける記述子計算のボトルネックに対処するため、著者らはcffiを介してPythonとインターフェースされたC++によるカーネル計算を実装した。本パッケージは、mpi4pyを利用して並列化を行い、共分散行列の計算と記述子の計算を複数のプロセッサに分散させている。スーパーコンピュータ環境では、物理コアはGPR計算機（共分散と予測を担当）とベース計算機（電子状態計算を担当）の間で分割される。

主な貢献

1. 動的なオンザフライ学習: 静的な「ワンショット」MLFFモデルとは異なり、GPR calculatorは学習プロセスをNEB最適化ループに直接統合している。不確実性が高い場合にのみ高価なベース計算機を選択的に呼び出すことで、モデルを適応的に洗練させる。
2. 不確実性定量化: フレームワークは予測不確実性を明示的に定量化し、計算効率と精度の厳密なバランスを可能にする。
3. ASEとの統合: 本パッケージはASEとのシームレスな統合を目的として設計されており、DFT以外の様々な電子状態手法をサポートしている。
4. 性能最適化: ハイブリッドC++/Python実装とMPI並列化により、大規模なデータセットの処理を可能にし、計算負荷の高いカーネルおよび記述子の計算を加速させている。

結果
著者らは、表面拡散および解離反応のNEBシミュレーションを用いてパッケージを検証し、GPR加速アプローチを純粋なDFT（VASP）計算と比較した。

• MgO(100)上のPd$_4$クラスター拡散: GPR calculatorは、純粋なVASPと比較して5倍の高速化（6時間 vs 36.5時間）を達成した。この手法は、収束までに337回のVASP呼び出しを必要とした純粋なDFTアプローチに対し、わずか80回のVASP呼び出しと340回のGPR呼び出しを必要とした。得られた活性化障壁（0.48 eV）は、純粋なVASPの結果（0.44 eV）とほぼ一致した。
• Pd(100)上でのH$_2$S解離: 結合切断を伴うより複雑なこの反応において、GPR calculatorは3.5倍の高速化（4.5時間 vs 推定約16時間）を達成した。これには395回のVASP呼び出しと1350回のGPR呼び出しを要した。エネルギープロファイルは、純粋なVASPの結果とイメージあたり0.015 eV以内で一致した。
• 遷移金属（Au, Ag, Cu, Pt）上でのH$_2$S解離: 本パッケージは4種類の(100)遷移金属表面に適用された。結果は、これらのシステム全体で3〜5倍の高速化を示した。GPR calculatorは、Ptが解離に最も有利なエネルギー特性を示し、Agが最も反応性が低いといった既知の傾向を正確に再現した。
• 偶発的な発見: Pd(100)上でのH$_2$S解離において、GPRの不確実性の確率的な性質により、モデルは曲がった反応経路を探索することができ、それが純粋なDFTで見出された線形経路（0.46 eV）よりも低い活性化障壁（0.31 eV）をもたらした。これは、最適化中の浅い局所解を回避するモデルの能力に起因すると著者らは考えている。

意義
本論文は、GPR calculatorが、鞍点や反応障壁を特定するために必要な精度を損なうことなく、NEB計算の計算コストを大幅に削減（3〜10倍の加速）できると主張している。シミュレーション中にサロゲートモデルを動的に更新することで、このアプローチは静的なMLFFの限界を克服し、高精度が極めて重要となる探索的なシミュレーションに適している。著者らは、このツールを計算材料科学および触媒分野における貴重な資産として位置づけており、複雑な化学プロセスのより効率的な探索を可能にするものであるとしている。なお、サロゲート駆動のシミュレーションは、固有のノイズにより収束までにより多くのイテレーションを必要とする場合があるが、これは力の許容誤差を調整することで緩和可能であり、最終的な結果は、GPRによって最適化された軌跡を最終的な純粋DFT計算の初期推測値として用いることで、迅速に精緻化できる。