機械学習で加速する遷移状態探索：Ru(0001)表面での窒素固定シミュレーション#

PHASE logo

アンモニア合成（ハーバー・ボッシュ法）において、ルテニウム(Ru)は非常に優れた触媒活性を示すことが知られています。この反応全体の律速段階は、強固な三重結合を持つ窒素分子（N₂）が金属表面上で解離し、窒素原子（N）になる過程です。本事例では、第一原理計算ソフトウェアAdvance/PHASEを用いて、Ru(0001)面上の窒素分子解離反応経路および反応障壁を算出しました。遷移状態探索には機械学習NEB（ML-NEB）法を適用し、計算コストの大幅な削減と高精度な解析を両立しています。

Keywords: 第一原理計算 (DFT), 窒素固定, 触媒反応, 解離吸着, Ru(0001), 遷移状態探索, ML-NEB, RPBE, 構造敏感反応

計算モデルと計算条件#

Ruは六方最密充填構造（hcp）を持つ金属であり、(0001)面はその最密充填面として最も安定で、表面の大部分を占める面です。本解析ではこのRu(0001)テラス面をp(2x4)スラブモデルで表現し、窒素分子が表面に垂直吸着した「始状態」から、解離して2つのN原子がhcp hollowサイトに吸着する「終状態」までの反応経路を計算しました。

交換相関汎関数には、金属表面上の吸着エネルギー評価において実績のあるRPBEを使用しました。これは、標準的なPBE汎関数が吸着エネルギーを過大評価（過剰結合）する傾向を改善したものです。

Ru(0001)上のN2解離反応モデル
図1.　Ru(0001)テラス面上でのN₂解離反応モデル（表面スーパーセルp(2x4)）。(左) 始状態: Topサイトへの垂直な分子吸着, (右) 終状態: Hollowサイトへの解離吸着。

表1.　計算条件の概要

項目	設定
計算手法	平面波基底・擬ポテンシャル法 (DFT) ※ウルトラソフト擬ポテンシャル使用
交換相関汎関数	GGA-RPBE
遷移状態探索法	ML-NEB (Machine Learning NEB)
波動関数のカットオフ	25 Ry (約340 eV)
k点サンプリング	6x3x1 【表面スーパーセル： p(2x4)】
スラブ層数	4層、下部2層固定

ML-NEB機能による計算コストの大幅削減#

従来のNEB（Nudged Elastic Band）法は、反応経路上の多数のイメージ（構造）を同時に最適化するため、多大な計算コストを要しました。本解析で用いたML-NEB（機械学習NEB）[1, 2] は、ポテンシャル曲面の予測に機械学習モデル（ガウス過程回帰）を導入することで、DFT計算の回数を大幅に削減します。最終的な遷移状態に対してはDFT計算による力の収束判定を行うため、遷移状態の計算精度は通常のCI-NEBと同等に保たれます。

解析結果#

計算された反応経路におけるエネルギー変化を図2に、各状態のエネルギー値を表2に示します。

Reaction Energy Profile
図2.　Ru(0001)上における窒素解離反応のエネルギープロファイル（基準：N₂分子が吸着した表面）。横軸のPath distanceは、反応経路に沿った全原子の変位の累積距離（反応座標）を表します。ML-NEBの計算には、反応経路全体の不確かさの最大値（閾値）max_uncertainty = 0.05 eVとしています。遷移状態はDFT計算による力の収束判定で確認されています。

表2.　各反応状態のエネルギー（分子あたり）

状態	構造の特徴	エネルギー [eV]	備考
始状態 (IS)	N₂分子がRu Topサイトに垂直吸着	0	ガス分子基準の吸着エネルギーは -0.64 eV（ISはガス状態より安定）
遷移状態 (TS)	N-N結合解離の鞍点	+1.59	ISからの障壁: 1.59 eV
終状態 (FS)	N原子がhcp hollowサイトに解離吸着	-1.17	発熱反応（安定化）

transition state
図3.　Ru(0001)上における窒素解離反応の遷移状態（左：上面図、右：斜視図）。

先行研究との比較・検証#

本計算の信頼性を確認するため、得られた解離障壁（1.59 eV）を過去の文献値と比較検証しました。

Mortensenらの計算 [3]: Ru(0001)テラス面上の解離障壁として約 1.36 ～ 1.48 eV を報告しています。汎関数やモデルの大きさの違いを考慮すれば、本計算の 1.59 eV はこれに近い値であり、計算手法が妥当であることを示しています。
Logadottirらの解析 [4]: BEP則 (Brønsted–Evans–Polanyi relation)に基づく解析において、Ru(0001)テラス面の障壁が約 1.8 eV と示されており、本計算結果はこのトレンドとも整合します。

考察：平坦面の不活性と真の活性中心#

本解析で得られた 1.59 eV という高いエネルギー障壁は、Ru(0001)テラス面が窒素解離に対して実質的に「不活性」であることを示しています。

1. 反応確率の低さ（実験事実との一致）#

アレニウス式 $k \propto \exp(-E_a/k_BT)$ に基づき、この障壁高さにおける室温付近での反応確率を見積もると、約 $10^{-12}$ 以下となります。これは、Mortensenらが指摘している「Ru(0001)平坦面でのN₂付着確率は $10^{-12}$ 以下である」という実験事実 [3] と良く一致します。すなわち、本計算結果はRuテラス面が反応にほとんど寄与しないことを正しく再現しています。

2. 構造敏感性とステップサイトの役割#

テラス面が不活性であるにもかかわらず、なぜ実際のRu触媒は高活性を示すのでしょうか。この理由は、Ru上のアンモニア合成が、表面の局所構造に強く依存する「構造敏感反応」であることで説明されます。

Logadottirら [4]: 平坦面では高い障壁を持つ一方、ステップ（段差）サイトでは障壁が大幅に（約0.5 eV付近まで）に低下することを示しています。
Medfordら [5]: 同様に、ステップサイトにおける解離障壁は計算および実験による推定値の範囲として 0.4 ～ 0.8 eV 程度と報告しており、ここが活性中心であると結論づけています。

以上の比較から、本解析は「テラス面は不活性であり、実際の反応はステップや欠陥サイトで進行する」という触媒設計上の重要な指針を、第一原理計算によって理論的に裏付ける結果となりました。

まとめ#

本解析では、第一原理計算ソフトウェアAdvance/PHASEを用いたML-NEB法により、Ru(0001)テラス面の窒素解離障壁を 1.59 eV と算出しました。この値は先行研究と良く一致し、平坦面での反応確率が極めて低い実験事実を定量的に再現しました。高活性な触媒設計には、テラス面ではなくステップサイトのような特殊な反応場が不可欠であることが、理論計算からも強く示唆されました。

参考文献#

J. A. Garrido Torres, P. C. Jennings, M. H. Hansen, J. R. Boes, and T. Bligaard, "Low-scaling algorithm for nudged elastic band calculations using a surrogate machine learning model", Phys. Rev. Lett. 122, 156001 (2019).
グラフェン単層膜の水素原子透過性
J. J. Mortensen, Y. Morikawa, B. Hammer, and J. K. Nørskov, "Density Functional Calculations of N2 Adsorption and Dissociation on a Ru(0001) Surface", J. Catal. 169, 85 (1997).
A. Logadottir, T. H. Rod, J. K. Nørskov, B. Hammer, S. Dahl, and C. J. H. Jacobsen, "The Brønsted-Evans-Polanyi Relation and the Volcano Plot for Ammonia Synthesis over Transition Metal Catalysts", J. Catal. 197, 229 (2001).
A. J. Medford, J. Wellendorff, A. Vojvodic, F. Studt, F. Abild-Pedersen, K. W. Jacobsen, T. Bligaard, and J. K. Nørskov, "Assessing the reliability of calculated catalytic ammonia synthesis rates", Science 345, 197 (2014).