ポテンシャルエネルギー曲面の新記述子:反応経路地図のパーシステント・ホモロジー

Abstract

A reaction route map (RRM), which is a collection of elementary reaction pathways, contracts the potential energy surface (PES) with 3N − 6 variables (N: the number of atoms) into a weighted graph representation. Although the automated construction of RRMs has greatly contributed to the accurate understanding of chemical reaction mechanisms, only a small fraction of networks with low activation energies are relevant to actual chemical reactions, and thus studies focusing on the entire RRM have not been conducted. In this letter, we summarize our recent approach to applying the persistent homology (PH) analysis to the graph structure of an RRM.

Translated Abstract

A reaction route map (RRM), which is a collection of elementary reaction pathways, contracts the potential energy surface (PES) with 3N − 6 variables (N: the number of atoms) into a weighted graph representation. Although the automated construction of RRMs has greatly contributed to the accurate understanding of chemical reaction mechanisms, only a small fraction of networks with low activation energies are relevant to actual chemical reactions, and thus studies focusing on the entire RRM have not been conducted. In this letter, we summarize our recent approach to applying the persistent homology (PH) analysis to the graph structure of an RRM.

1 背景

化学反応の静的・動的な描像を思考する上で，ポテンシャルエネルギー曲面(PES)は欠かせない概念である [1]．しかし，3N-6 変数(N: 原子数)の関数であるPESは，計算機で扱うことを考えても複雑すぎる．それだけでなく, 異なる原子数の系のPES同士を同じ土俵で直接比較することもできない． Mirthらは数学的に定義された「穴」に注目した解析手法であるパーシステント・ホモロジー(PH)を用いたPESのトポロジー解析法を提案している [2]．特に低次のパーシステンス・バーコード(PB)がPESの特徴を表す記述子として共通に利用できることを示唆している．しかし，Mirthらの手法では3N− 6変数関数であるPESのすべての定常点を解析的に求める必要があり，実在系に対してそのまま用いることは難しい．一方，前田らが開発するGRRMプログラム [3]等により, 系の反応経路地図(RRM)を自動構築する手法がめざましく発展している RRMは，PESの中で化学的に特に重要な極小点(EQ)，一次鞍点(TS)と固有反応座標(IRC)を集めたものである．RRMの自動構築は，化学反応メカニズムの正確な理解に大きく貢献したが，実際の化学反応に関係するのは活性化エネルギーの低いほんの一部のネットワークのみであり，RRM全体に着目した研究は行われていない そこで我々は，RRMに対してPH解析を適用することを着想し，報告した [4]．ある種PESの縮図ともいえるRRMに内在するグローバルな情報から，PESのトポロジー情報を再構成する現実的な手法を構築することにより，異なる系や化学反応の横断的比較・分類，さらには予測を行うことを目指している．

2 手法

Figure 3.

Persistence barcodes for RRMs of (a) Au₃Ag₂, (b) Au₂Ag₃, (c) Au₃Cu₂, (d) Au₂Cu₃, (e) Ag₃Cu₂, and (f) Ag₂Cu₃. (Reprinted with permission from Ref [4]. Fig.16. Copyright 2023 American Chemical Society.)

RRMは，数学的には頂点(EQ)と辺(TS)の双方に重み(エネルギー値)を持つ重み付きグラフで表現される．例として，Figure 1(a)に4つのEQ(v_i)と5つのTS($\overline{v_i{v}_j}$)から成る模式的なRRMを示したので，これのPBがどのように得られるか概説しておこう．EQとTSに付随する数字は重み(エネルギー値)を示しており，TS $\overline{v_i{v}_j}$の重みは常にEQ v_iとv_jの重みよりも大きいと仮定する．PBを求めるには，Figure 1(a)で示したRRMグラフ全体に対し，閾値以下の重みを持つ要素のみから成る部分グラフを考え，閾値を徐々に大きくしていく操作を行う．すると，Figure 1(b)のような増大列が得られる． 0次のホモロジーは独立要素の数を表しており，ある閾値で独立成分が増えたら「生成」(birth, b)を，独立成分が減ったら「消滅」(death, d)を記録する．消滅する成分は，直近で生成した成分とする，というlast-in first-out方式(elder rule)を適用するのがPHの流儀である． 区間[b,d]で表されるバーでこれを表現すると，Figure 1(b)の増大列に対しては，Figure 1(c)のβ₀で示したような4つのバーで0次のホモロジーが表現される． これがPBである． 一方，1次のホモロジーは周回する経路(輪っか)の数を表しており，周回経路ができたら「生成」を記録する．グラフ構造で輪っかが消えてなくなることはないが，今回は便宜的に，3辺から成る輪っか(三角形)ができた時点で「消滅」を記録することとした．Figure 1(b)の増大列に対しては，閾値が4で四角形ができ，5でそれが三角形2つに分かれて，すぐに「消滅」が記録されるので，Figure 1(c)のβ₁で示したような2つのバーで1次のホモロジーが表現される．このようなルールを作ると，任意の重み付き単純グラフに対して0次と1次のPBを得ることができる．

Figure 1.

(a) RRM toy model, (b) its increasing sequence, and (c) its persistence barcode. (Reprinted with permission from Ref [4]. Fig.9. Copyright 2023 American Chemical Society.)

Figure 2.

(a) CNPI-contracted RRM of Au₅ and (b) its reconstructed RRM in (3N − 6)-dimensional shape space. Each blue dot corresponds to a distinct permutation isomer in the shape space, and each edge corresponds to a reaction path. (Reprinted with permission from Ref [6]. cover art. Copyright 2023 American Chemical Society.)

ところで，一般にRRMを構築するプログラムでは，探索の効率化のために，同種核の置換と反転操作からなる群(CNPI群)の作用により一致する配座は同一視される [5]． すなわち3N− 6次元よりも縮約されたPESの部分空間に対してRRMを構築する． 異なるCNPI群を持つ分子間でRRMを比較したければ，元の3N− 6次元のPES全体のRRMを復元してから行うべきである． そこで，縮約されたRRMから3N− 6次元空間におけるRRMを再現するアルゴリズムを開発した [6]． このアルゴリズムを用いると，例えば金5量体について, Figure 2(a)に示す縮約されたRRMから，Figure 2(b)に示す3N− 6次元空間のRRMが得られる．

3 結果

ここからは，貨幣金属ナノクラスターに適用した結果をいくつか紹介する． まず，貨幣金属ナノクラスターX₃Y₂ (X_,Y = Cu, Ag, Au)のRRMに対するPBをFigure 3に示す． これらの系のCNPI群は同一なため，RRMの展開は行っていない． 横に並んだ同じ元素の組合せのPBは似た傾向を示すことが確認できる． 具体的には，Au/Cuの組合せではバーが50 kcal/molまでの広い範囲に見られるが，他の組合せでは，30 kcal/mol程度までにしか見られない． また，Ag/Cuの組合せでは，他よりもバーの数自体が少ない． さらに，Au/Agの組合せでは，赤の半直線で表される1次のバーの数が他よりも多い． これらの結果は，開発した手法により2つの金属元素の親和性や合金反応の特性の抽出に利用できる可能性を示唆している．

次に，金5量体について，CNPI群で縮約されたGRRMプログラムが出力するRRM(Figure 2(a))と，3N-6次元空間で再構成したRRM(Figure 2(b))に対してPBを構築して比較した(Figure 4(a)と(b))．再構成されたRRMに対するPBに表れるバーは，縮約されたRRMに対するものよりも多い．また，縮約RRMでは自己ループとなる辺や二重辺の一部が再構成RRMでは単純グラフの辺で表現されるため，再構成RRMのPBには縮約RRMのPBにはないバーが現れる． さらに，となる1次のバーが多く現れた．これは，RRMの展開により多数の周回経路が生じることに起因する．

Figure 4.

Persistence barcodes of Au₅ for (a) CNPI-contracted RRM and (b) reconstructed RRM in 3N− 6-dimensional shape space.

4 結論

本稿では, RRMに対するPH解析法を提案し, PESの横断的比較に向けた予備検討を行った. 特定の反応にのみ興味を持つのであれば, その系のRRMを詳細に分析すれば十分だが, 本手法により抽出されたPBという特徴量は, 複数のRRMを比較・分類するのに有用だと期待される.

謝辞

本研究は文部科学省科研費(JP21H05544, JP23H04093),日本学術振興会科研費(JP23H01915)およびIQCE研究助成の支援を受けて実施された．計算の一部は岡崎の計算科学研究センターの計算機を利用して行った(21-IMS-C018, 22-IMS-C019, 23-IMS-C016)．村山は北海道大学博士人材フェローシップ(情報・AI)による支援を受けた．この場を借りて謝意を表したい．

参考文献

• [1]   D. Wales, Energy Landscapes: Applications to Clusters, Biomolecules and Glasses, Cambridge University Press (2003).
• [2]   J. Mirth, Y. Zhai, J. Bush, E. G. Alvarado, H. Jordan, M. Heim, et al., J. Chem. Phys., 154, 114114 (2021). doi:10.1063/5.0036747 PMID:33752361
• [3]   S. Maeda, K. Ohno, K. Morokuma, Phys. Chem. Chem. Phys., 15, 3683 (2013). doi:10.1039/c3cp44063j PMID:23389653
• [4]   B. Murayama, M. Kobayashi, M. Aoki, S. Ishibashi, T. Saito, T. Nakamura, et al., J. Chem. Theory Comput., 19, 5007 (2023). doi:10.1021/acs.jctc.2c01204 PMID:37395411
• [5]   K. Ohno, H. Satoh, Exploration on Quantum Chemical Potential Energy Surface, Royal Society of Chemistry (2023).
• [6]   H. Teramoto, T. Saito, M. Aoki, B. Murayama, M. Kobayashi, T. Nakamura, et al., J. Chem. Theory Comput., 19, 5886 (2023). doi:10.1021/acs.jctc.3c00500 PMID:37642714