JP2023110211A - 反応経路探索プログラム、反応経路探索システム及び反応経路探索方法 - Google Patents

Abstract

【課題】化学的な反応経路の探索において、原子数の増加に応じた計算コストの増加を抑制しやすい反応経路探索プログラム、反応経路探索システム及び反応経路探索方法を提供する。【解決手段】反応経路の探索対象となる構造中の原子からなる複数のフラグメントペアのそれぞれについて、当該構造のポテンシャルエネルギー曲面（ＰＥＳ）を示すＥ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応するＱ＝Ｑ０におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する。次に、複数のフラグメントペアのうち、ａ及びｂの少なくともいずれかの大きさに基づく優先順位の高いフラグメントペアに関して、それより優先順位の低いその他のフラグメントペアに先行して、第１の平衡状態から他の平衡状態である第２の平衡状態までの遷移における反応経路探索を実行する。【選択図】図１

Description

本発明は、化学的な反応経路を探索する反応経路探索プログラム、反応経路探索システム及び反応経路探索方法に関する。

与えられた構造に対する化学的な反応経路を探索する方法として非特許文献１に記載の人工力誘起反応法（ＡＦＩＲ法）がある。ＡＦＩＲ法は、構造中の複数の原子の組み合わせからなるフラグメントペアに力を加えることで構造の変形を誘起し、別の構造へと遷移させる手法である。このとき、遷移の過程において系が辿る座標変化から、対応する反応経路が得られる。この操作を、構造中の様々なフラグメントペアに対して系統的に適用することにより、様々な別の構造への反応経路が自動探索される。

前田 理，原渕 祐、「分子および周期系における化学変換の経路探索：力を利用したアプローチ（Exploring paths of chemical transformations in molecular and periodic systems: An approach utilizing force）」、Wiley Interdisciplinary Reviews: Computational Molecular Science, 11, e1538

ある構造からなる系において考慮すべきフラグメントペアの数は、系に含まれる原子数Ｎの二乗に比例して増加する。ＡＦＩＲ法における反応経路の探索はフラグメントペアごとに実施される。このため、探索対象となる構造における原子数が増加すると、計算すべきフラグメントペアの組み合わせが飛躍的に増加し、反応経路を探索するための計算コストが膨大なものとなるおそれがある。

本発明の目的は、原子数の増加に応じた計算コストの増加を抑制しやすい反応経路探索プログラム、反応経路探索システム及び反応経路探索方法を提供することにある。

本発明に係る反応経路探索プログラムは、複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出する反応経路探索システムとして１又は複数のコンピュータを機能させるプログラムであって、前記反応経路探索システムが、前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出手段と、互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出手段と、前記複数のフラグメントペアのうち、前記微分係数算出手段が算出したａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出手段に算出させる平衡状態間変化算出手段とを備えている。

（数式１）

（数式２）

また、本発明の別の観点に係る反応経路探索システムは、複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出するシステムであって、前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出手段と、互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出手段と、前記複数のフラグメントペアのうち、前記微分係数算出手段が算出したａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出手段に算出させる平衡状態間変化算出手段とを備えている。

また、本発明の別の観点に係る反応経路探索方法は、複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出する方法であって、前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出ステップと、互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出ステップと、前記複数のフラグメントペアのうち、前記微分係数算出ステップにおいて算出されたａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出ステップによって算出させる平衡状態間変化算出ステップとを備えている。

本発明の反応経路探索プログラム、反応経路探索システム及び反応経路探索方法によると、複数個の原子からなる構造体における反応経路を、原子の位置関係によって表される構造の変化として算出する。このとき、第１の平衡状態に対応する原子の位置におけるＥ（Ｑ）の二階の微係数ｂ及び三階の微係数ａの少なくともいずれかの大きさに基づく優先順位の高いフラグメントペアに関して、優先順位の低いフラグメントペアに先行して構造の変化が算出される。例えば、二階の微係数ｂが小さいことは、局所的なＥ（Ｑ）の上昇が小さい傾向を示す。このため、例えば、ｂが小さいフラグメントペアを優先した場合には第１の平衡状態付近におけるＥ（Ｑ）の上昇が小さい経路を得ることができる。また、三階の微係数ａが小さいことは、Ｅ（Ｑ）において下に凸から上に凸への変化が起こる傾向を示す。このため、別の例としてａが小さいフラグメントペアを優先した場合には第１の平衡状態付近において遷移状態が見つかりやすいと期待される。したがって、Ｅ（Ｑ）の二階の微係数ｂ及び三階の微係数ａの少なくともいずれかの大きさに基づいて一部のフラグメントペアに関してその他のフラグメントペアよりも優先的に計算を実行することにより、計算コストをそれほど掛けることなく目的となる反応経路を取得できる可能性が高くなる。このため、原子数の増加に応じた計算コストの増加が抑制されやすい。

また、本発明においては、前記反応経路探索システムが、前記平衡状態間変化算出手段による過去の算出結果に基づいて得られた全ての前記第２の平衡状態と初期状態とからなるＮ個（Ｎ：２以上の自然数）の前記平衡状態からいずれかを前記第１の平衡状態として選択する状態選択手段をさらに備えており、前記状態選択手段が選択した前記第１の平衡状態から前記第２の平衡状態までの遷移における前記フラグメントペアに関する前記構造の変化を前記平衡状態間変化算出手段が前記構造変化算出手段に算出させることを繰り返すことが好ましい。これによると、すでに取得された平衡状態のいずれかを第１の平衡状態とし、第１の平衡状態から第２の平衡状態までの構造の変化の算出、つまり、反応経路の取得を繰り返す。これにより、初期状態及び過去に取得された第２の平衡状態からなる取得済みの平衡状態並びにこれらの平衡状態間を結ぶ反応経路からなるネットワークを取得可能である。

また、本発明においては、前記反応経路探索システムが、前記Ｎ個の平衡状態間の遷移に係るＮ＊Ｎ個の速度定数からなるＮ行Ｎ列の速度定数行列から、前記Ｎ個の平衡状態の加重和で表されるｌ個（ｌ：ｌ＜Ｎを満たす自然数）の超状態間の遷移に係るｌ＊ｌ個の速度定数からなる速度定数行列を、ＲＣＭＣ法に基づいて速度定数行列を縮約することをｍ回（ｍ：ｍ＝Ｎ－ｌを満たす自然数）繰り返すことで取得する速度定数縮約手段を含んでおり、前記状態選択手段が、ｍ＝１，２，…，Ｍ（Ｍ：Ｍ＜Ｎを満たす自然数）のそれぞれについて、縮約された速度定数行列を前記速度定数縮約手段が取得した結果に基づいて、前記Ｎ個の平衡状態のそれぞれをＥＱ_ｉ（ｉ：Ｎ以下の自然数）とし、前記ｌ個の超状態のそれぞれをＳＳ_ｊ（ｊ：ｌ以下の自然数）とし、ｍ回の縮約後のＥＱ_ｉのポピュレーションをｐ_ｉ（ｍ）とし、ｍ回の縮約後のＳＳ_ｊのポピュレーションをＱ_ｊ（ｍ）とし、ｍ回の縮約後のＳＳ_ｊに対するＥＱ_ｉの寄与をχ_ｊｉ（ｍ）とし、ＥＱ_ｉの相対ギブズエネルギーをΔＧ_ｉとし、気体定数をＲとし、モデル温度パラメータをＴ_Ｒとするときに、数式３に基づいて算出されるｐ_ｉ（ｍ）及び数式４に基づいて算出されるΛ_ｉが大きいほど、ＥＱ_ｉが選択されやすくなるように、前記Ｎ個の平衡状態からいずれかを前記第１の平衡状態として選択することが好ましい。これによると、反応経路のトラフィック量を表すΛ_ｉが大きい平衡状態が反応経路の算出において選択されやすくなる。

（数式３）

（数式４）

また、本発明においては、前記状態選択手段が、Λ_ｉが大きいＥＱ_ｉほど選択されやすく、且つ、ＥＱ_ｉを前記第１の平衡状態として前記構造変化算出手段が前記構造の変化を算出した回数ｎ_ｉが小さいＥＱ_ｉほど選択されやすく、且つ、前記Ｎ個の平衡状態について前記構造変化算出手段が前記構造の変化を算出した全体の回数ｎ_ａｌｌが大きいほどΛ_ｉの違いが前記平衡状態の選択されやすさの違いに表れにくくなるように、前記Ｎ個の平衡状態からいずれかを前記第１の平衡状態として選択することが好ましい。これによると、ｎ_ａｌｌが比較的小さい場合には、Λ_ｉが大きく且つｎ_ｉが小さい平衡状態ほど計算対象として選択されやすくなる。これによって、構造の変化の算出開始からある程度ｎ_ａｌｌが大きくなるまで、Λ_ｉが大きい（トラフィック量が大きい）平衡状態が適切に選択されやすい。一方、ｎ_ａｌｌが比較的大きくなっていくと、それに応じてΛ_ｉの違いが平衡状態の選択されやすさの違いに表れにくくなってくる。したがって、トラフィック量が小さい平衡状態も選択されるようになってくる。このように、上記構成によると、計算の進度によって変動する優先度に適切に応じた平衡状態の選択が実施される。

また、本発明においては、前記状態選択手段が、前記Ｎ個の平衡状態について前記構造変化算出手段が前記構造の変化を算出した全体の回数をｎ_ａｌｌとし、ＥＱ_ｉを前記第１の平衡状態として前記構造変化算出手段が前記構造の変化を算出した回数をｎ_ｉとし、ξ_ｉを０以上１以下の実数とし、α及びβをそれぞれ０以上１以下の実数とするときに、数式５のν_ｉが最も大きいＥＱ_ｉを前記第１の平衡状態として選択することが好ましい。数式５によって求められるν_ｉは、Λ_ｉが大きいほど大きく、且つ、ｎ_ａｌｌが比較的小さい場合にはｎ_ｉが小さいほど大きい値を有するが、ｎ_ａｌｌが比較的大きくなるとΛ_ｉの違いがν_ｉの違いに現れにくくなる。したがって、ν_ｉが最も大きいＥＱ_ｉを選択することで、計算の進度によって変動する優先度に適切に応じた平衡状態が選択される。

（数式５）

本発明の一実施形態に関して用いられるＡＦＩＲ法における複数の原子からなるフラグメントペアの概念図である。 本発明の一実施形態に関して用いられるＡＦＩＲ法に基づいて取得される反応経路のネットワークを示す概念図である。 本発明の一実施形態に係る反応経路探索システムの機能的構成を示すブロック図である。 図１のフラグメントペアにおける第０層原子同士の相対的な位置をこれらの原子を通る直線に沿って変更した場合におけるＰＥＳの変化を、平衡状態ＥＱ_０の周辺に関して概念的に示すグラフである。 図３の反応経路探索システムが実行する一連の処理の流れを示すフロー図である。

本発明の一実施形態に係る反応経路探索システム１について図面を参照しつつ説明する。反応経路探索システム１は、複数個の原子からなる分子等の構造体における反応経路を探索するものである。反応経路探索システム１は、人工力誘起反応法（ＡＦＩＲ法）、特に、単一成分ＡＦＩＲ法（ＳＣ－ＡＦＩＲ法）に基づいて反応経路を探索する。まず、ＳＣ－ＡＦＩＲ法に基づく反応経路の探索方法の概要について説明する。続いて、反応経路探索システム１の構成や機能等の詳細について説明する。

［反応経路探索方法の概要］
ＡＦＩＲ法は、次の数式１及び２に示すＡＦＩＲ関数Ｆ^ＡＦＩＲを使用して、図１に示すようなフラグメントＡ及びＢからなるフラグメントペア間に力を加えることによって構造変換を誘発する方法である。

（数式１）

（数式２）

ここで、フラグメントＡは、上記構造体中のＮ_１個（Ｎ_１：自然数）の原子からなる構造体の断片である。フラグメントＢは、上記構造体中の別のＮ_２個（Ｎ_２：自然数）の原子からなる構造体の断片である。Ｅ（Ｑ）は、上記構造体中の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギー曲面（ＰＥＳ）に対応する。r_ｓｔは、フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子とフラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離である。Ｒ_ｓ及びＲ_ｔは、それぞれ上記ｓ番目の原子及び上記ｔ番目の原子の共有結合半径である。ρは＋１及び－１のいずれかに設定される。ｐは、４以上且つ８以下の範囲から適宜（例えば、ｐ＝６に）設定される。フラグメント間に作用させる力の強さを示すαは次の数式１－１に基づく定数である。数式１－１のＲ_０とεとして、Ａｒ－Ａｒ Ｌｅｎｎａｒｄ－Ｊｏｎｅｓパラメータ、つまり、それぞれ３．８１６４＊１０^－１０ｍと１．００６１ｋＪ／ｍｏｌが採用される。数式１－１が表すαは、Ａｒ－Ａｒペアが互いに衝突エネルギーγで直接衝突する場合において、最小点から転換点まで移動するときに受ける平均力に対応する。したがって、γはモデル衝突エネルギーパラメータと呼ばれ、ＡＦＩＲ法における力が克服できるバリアのおおよその上限に対応する。γはユーザによって適宜設定されればよい。

（数式１－１）

反応経路は、標準の準ニュートン法によってＡＦＩＲ関数を最小化して得られる経路である（非特許文献１参照）。この経路を「ＡＦＩＲ経路」と呼ぶ。ＡＦＩＲ経路に沿って、近似平衡（ＥＱ）構造及び遷移状態（ＴＳ）構造が、それぞれＡＦＩＲ関数が極小値を取る構造及び極大値を取る構造として取得される。これらの構造は、実際のＥＱ及びＴＳの初期推定として使用可能である。ただし、これらのエネルギーの最大値が実際のＴＳから大きく外れている場合がある。この状況は、特にγが大きすぎる場合に発生する可能性がある。このような場合でも、ＡＦＩＲ経路の緩和計算を実行することで、これらのＴＳの見落としが回避される。この目的のために、パスポイント全体を最適化する任意のダブルエンド方式を使用できる。このための方法として、局所平面更新（ＬＵＰ）法が採用されている（Ｊ．Ｃｈｅｍ．Ｐｈｙｓ．１９９１，９４，７５１）。ＬＵＰ法を使用した一実装例によると、パスポイントは特定のパスに沿って均等に分散され、パスの接線に垂直な超平面内の低エネルギーポイントに移動する。これによって得られた経路を「ＬＵＰ経路」と呼ぶ。ここで、２つのエンドポイントとすべての極大値は、それぞれ最小値と最大値に直接最適化される。ＡＦＩＲ経路又はＬＵＰ経路上の近似ＥＱ及びＴＳは、標準の準ニュートン法によって実際のＥＱ及びＴＳに最適化される。次に、取得したすべてのＴＳから開始して、固有反応座標（ＩＲＣ）の計算が実行される。最後に、取得したすべてのＥＱ及びＴＳにおいて基準振動解析が実行される。

ＳＣ－ＡＦＩＲ法では、ＡＦＩＲ経路はＥＱから計算される（Ｊ．Ｃｏｍｐｕｔ．Ｃｈｅｍ．２０１４，３５，１６６）。ＥＱに対応する構造は、分子、分子複合体、有機金属複合体、クラスター、表面の吸着構造、結晶構造等である。単一のＥＱから、様々なフラグメントＡ及びＢからなるフラグメントペアのＡＦＩＲ関数を最小化することにより、様々なＡＦＩＲ経路が取得される。各ＥＱに対応する構造では、フラグメントペアは、原子ペアに焦点を当てることによって体系的に定義される。つまり、図１に示す通り、フラグメントＡとフラグメントＢでは、２つの原子が第０層原子として割り当てられている。次に、第０層の原子に接続された第１層の原子と、第１層の原子に接続された第２層の原子が、フラグメントＡとＢに追加される。各フラグメント中の第１層及び第２層の原子において、他方のフラグメント内の原子との距離が第０層原子間の距離未満の原子は、各フラグメントから削除される。一連のフラグメントペアは、距離が非常に長いものを除く全ての原子ペアにこの手順を適用することによって生成される。次に、各フラグメントペアのＡＦＩＲ経路が計算される。仮に、全てのフラグメントペアについてＡＦＩＲ経路が計算される場合、Ｎ_{ａｔｏｍｓ}を原子の数とすると、単一のＥＱの周りの計算を完了するために必要なコストはＮ_{ａｔｏｍｓ} ^２とスケーリングされる。既に取得済みのＥＱを出発点としてＡＦＩＲ経路を算出し、もって、新たなＥＱやＥＱ間の経路の取得を繰り返すことにより、Ｎ個（Ｎ：２以上の自然数）のＥＱで構成される図２に示すような反応経路ネットワークが取得される。以下、Ｎ個のＥＱからなる反応経路ネットワークにおいて、ｉ番目（ｉ：Ｎ以下の自然数）のＥＱをＥＱ_ｉと示す場合がある。このうち、ＥＱ_１は初期構造としての反応物を示す。図２は、Ｎ＝６の場合の反応経路ネットワークの一例を示す。

次に、反応経路の計算コストを低下させるために採用されるキネティックベースのナビゲーションについて説明する。ＥＱの数は、原子の数Ｎ_{ａｔｏｍｓ}に関して指数関数的に増加することが示されている。したがって、Ｎ_{ａｔｏｍｓ}が小さくない限り、上記のＳＣ－ＡＦＩＲ法に基づく上記の計算手順を全てのＥＱに適用することは困難である。一方、穏やかな実験条件下で動的にアクセスできるＥＱの数は限られている。そこで、速度論的方法である速度定数行列縮約（ＲＣＭＣ）法を使用し、特定の実験条件下でＥＱが速度論的にアクセス可能かどうかを判断することにより、速度論的にアクセスできないＥＱの探索コストを排除する（Ｃｈｅｍ．Ｌｅｔｔ．２０１９，４８，４７）。

ＲＣＭＣ法では、縮約と呼ばれる操作を繰り返すことにより、Ｎ個のＥＱで構成される反応経路ネットワーク全体を考慮した速度論的分析を可能にする。Ｎ個のＥＱ同士の遷移に関するＮ＊Ｎ個の速度定数から、Ｎ行Ｎ列の速度定数行列が得られる。ＲＣＭＣ法では、Ｍ回の縮約によって、より小さなｎ行ｎ列の速度定数行列（ｎ＝Ｎ-Ｍ）が得られる。Ｍ回の縮約後のｎ行ｎ列の速度定数行列は、ｎ個の超状態（ＳＳ）間の遷移における速度定数を表す。ＳＳは、全てのＥＱの加重和として表される。１回の縮約により速度定数行列のサイズが１つ減少し、準定常状態近似とボルツマン分布保存条件に基づいて、速度定数行列のすべての非対角要素が更新される。縮約によって得た速度定数行列の非対角要素は、あるＳＳから別のＳＳへの遷移の速度定数に対応する。縮約は、ユーザが適宜設定可能なしきい値ｔ_ＭＡＸよりもライフタイムが短いすべての状態に適用される。Ｍ回の縮約により、ｔ_ＭＡＸよりも長いタイムスケールで相互に遷移するｎ個のＳＳが取得される。ＲＣＭＣ法のより詳細な説明は後述の通りである。

ｍ回縮約後のＥＱ_ｉのポピュレーションをｐ_ｉ（ｍ）とする。ポピュレーションとは、ある時刻におけるそのＥＱ_ｉの存在確率を示す。キネティックベースのナビゲーションでは、ユーザは実験条件を考慮して初期ポピュレーションを提供する必要がある。たとえば、ＥＱ_１が反応物である場合、ＥＱ_１の初期ポピュレーションは１（ｐ_１（０）＝１）に設定され、他のすべてのＥＱのポピュレーションは０（ｐ_ｉ≠１（０）＝０）に設定される。次に、全てのＥＱ_ｉについて、数式３及び４に基づき、インデックスΛ_i、いわゆるトラフィック量が計算される。

（数式３）

（数式４）

数式３、４において、Ｑ^（ｍ） _ｊは、ｎ個のＳＳのそれぞれをＳＳ_ｊ（ｊ：ｊ≦ｎを満たす自然数）としたときに、ｍ回の縮約後のＳＳ_ｊのポピュレーションを示す。χ^（ｍ） _ｊｉは、ｍ回の縮約後のＳＳ_ｊへのＥＱ_ｉの寄与である。ΔＧ_ｉは、量子化学計算に基づく基準振動解析によって得られるＥＱ_ｉの相対ギブズエネルギーである。Ｒは気体定数、Ｔ_Ｒはモデル温度パラメータである。数式４は、ｔ_ＭＡＸ内で発生するＥＱ_ｉへのポピュレーション流入とＥＱ_ｉからのポピュレーション流出の合計を表す。例えば、反応物ＥＱ、反応物ＥＱより十分安定な中間体ＥＱ、および中間体ＥＱより十分安定な生成物ＥＱで構成される３状態システムでは、実験条件として反応物ＥＱから生成物ＥＱへの経路が中間体ＥＱを経由する場合にのみ開いている場合には、反応物ＥＱと生成物ＥＱはΛ_ｉ＝１．０であり、中間体ＥＱはΛ_ｉ＝２．０である。一方、実験条件からＥＱ_ｉにアクセスできない場合、Λ_ｉはゼロになる。本実施形態においては、後述の通り、かかるΛ_ｉ値を使用してＮ個のＥＱから１つが選択される。そして、選択されたＥＱにおいてＡＦＩＲ経路が算出される。ＡＦＩＲ経路は、まだ処理されていないフラグメントペアから適宜選択されたフラグメントペアに対して計算される。Ｔ_Ｒは、ゼロ以外のΛ_ｉを持つＥＱの中から高エネルギーＥＱが選択される頻度を決定する。通常は、例えばＴ_Ｒ＝３０００Ｋのような高い値に設定される。

［システムの詳細］
反応経路探索システム１は、１又は複数台のコンピュータを含むハードウェアと、反応経路探索システム１の各種機能を実現するようにハードウェアを機能させるプログラムを含むソフトウェアとを備えている。各コンピュータは、ＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ）、ＲＯＭ（Ｒｅａｄ－Ｏｎｌｙ Ｍｅｍｏｒｙ）及びＲＡＭ（Ｒａｎｄｏｍ Ａｃｃｅｓｓ Ｍｅｍｏｒｙ）等のメモリ装置、並びに、入出力インタフェース等の各種インタフェース等からなる。また、コンピュータには、ユーザ入力を受け付ける入力装置やディスプレイ等の出力装置、ハードディスク等の記録装置等が外部装置として接続されるか、内部装置として搭載されている場合がある。ソフトウェアは、メモリ装置やハードディスク等の記録装置に記録されたプログラムデータ等からなる。そして、コンピュータ等のハードウェアがソフトウェアに従って演算処理、入出力処理等の各種の情報処理を実行する。ソフトウェアに従ったハードウェアの機能により、以下において説明する反応経路探索システム１の各種の機能が実現される。

なお、ハードウェアが複数台のコンピュータを含んでいる場合、複数台のコンピュータが協働して上記各種の機能に係る処理を実行してもよい。かかる協働には様々な態様があり得る。例えば、コンピュータ同士がインターネット等の通信ネットワークを通じて接続され、各コンピュータが必要な情報を互いに通信しつつ必要な演算処理を分担して実行してもよい。また、一のコンピュータはユーザインタフェース機能を担当し、他のコンピュータは、一のコンピュータにおいてユーザ入力され、通信ネットワークを通じて送信された入力情報に基づいて演算処理を実行してもよい。

反応経路探索システム１は、図３に示すように、記憶部１０、フラグメントペア導出部２０、フラグメントペアランキング取得部３０、状態選択部４０及び反応経路算出部５０を備えている。記憶部１０は、取得済みＥＱのリストを示すデータ、原子の座標を示す構造データ、フラグメントＡ及びＢにおける原子の組み合わせを示すフラグメントデータ、探索回数データ等を記憶している。取得済みＥＱのリストは、反応物に対応するＥＱ_１と反応経路探索によって取得されたその他のＥＱ（ＥＱ_ｉ≠１）とからなる。構造データは、取得済みＥＱのリストに含まれる各ＥＱについて原子の座標を示すデータを含んでいる。原子の座標は、構造体を構成する各原子の種類及び位置を表している。原子の位置は、例えば、原子の３次元の位置座標（例えば、デカルト座標系上の座標）によって示される。フラグメントデータは、取得済みＥＱのリストについて、各ＥＱにおいて設定されたフラグメントペア（フラグメントＡ及びＢ）における原子の組み合わせを示すデータである。記憶部１０は、ＥＱと関連付けてフラグメントデータを記憶している。フラグメントデータは、フラグメントペアに関してランキング（本発明における優先順位）が取得されている場合には、そのフラグメントペアのランキングを示すデータを含んでいる。ランキングの高さは、反応経路探索が行われる優先順位の高さを示す。探索回数データは、ＥＱ_ｉを出発点とする反応経路探索が実行された回数ｎ_ｉを取得済みのＥＱリストに含まれる各ＥＱ_ｉ（ｉ＝１，２，…）について示す。さらに、フラグメントデータは、各フラグメントペアに関して反応経路探索が実行済みであるか否かを示すデータを含んでいる。このデータは、デフォルトでは実行済みでないことを示す。また、記憶部１０は、反応経路探索システム１の各部が以下の演算を行う際に使用する各種の数値条件を示すデータを記憶している。

フラグメントペア導出部２０は、ＳＣ－ＡＦＩＲ法に従った上述の方法に基づいて、構造体を構成する全ての原子からフラグメントＡ及びＢのそれぞれを構成する原子の組み合わせを抽出する。この抽出結果を示すフラグメントデータは記憶部１０に記憶される。フラグメントペア導出部２０は、各ＥＱ_ｉについてフラグメントペアを構成する原子の組み合わせを生成する。

フラグメントペアランキング取得部３０は、記憶部１０のフラグメントデータを参照し、取得済みＥＱのリストから、フラグメントペアのランキングが取得されていないＥＱ（以下、対象ＥＱとする。）を把握する。そして、対象ＥＱのそれぞれについて以下の通りランキングを取得する。

まず、フラグメントペアランキング取得部３０は、対象ＥＱに関してフラグメントペア導出部２０が導出したフラグメントペアのそれぞれについて、Ｑ＝Ｑ_０におけるＥ（Ｑ）の二階微分係数ｂ及び三階微分係数ａを導出する。ここで、Ｅ（Ｑ）は上述の通りＰＥＳを示す。Ｑ_０は対象ＥＱに対応する構造における原子の位置を示す。二階微分係数ｂ及び三階微分係数ａは、例えば、図４に示すＱ＝Ｑ_０、Ｑ＝Ｑ’における曲線Ｃの接線の傾きとＥ（Ｑ）の値とに基づいて算出される。Ｑ’（＝Ｑ_０＋ΔＱ）は、Ｑ＝Ｑ_０におけるフラグメントペア中の一対の第０層原子間の距離をＤ（図１参照）とするとき、かかる距離がＤからＤ＋ΔＤ（例えば、ΔＤ＝０．０５オングストローム）に微小に変化するように、一対の第０層原子間を結ぶ直線に沿って、一対の第０層原子における少なくともいずれか一方の位置を他方に対して微小に相対移動させた後の構造における原子の位置を示す。また、曲線Ｃは、上記直線に沿って上記一対の第０層原子同士の相対位置を変化させた際のＥ（Ｑ）の変化を示すＰＥＳに沿った曲線である。フラグメントペアランキング取得部３０は、Ｑ’又はΔＱを次のような条件を満たすように最適化する。その条件は、一対の第０層原子間の距離をＤ＋ΔＤに固定しつつ、全原子間の距離を要素とする距離行列に関し、Ｑ＝Ｑ’における距離行列の第０層原子間に対応する要素以外の行列要素が、Ｑ＝Ｑ_０における距離行列の第０層原子間に対応する要素以外の行列要素をできる限りよく再現することである。Ｅ（Ｑ）を三次関数ｆ（ｘ）＝ａ_３＊ｘ^３＋ａ_２＊ｘ^２＋ａ_１＊ｘ＋ａ_０で表されるとみなし、４つの数式：ｆ（ｘ（Ｑ_０））＝Ｅ（Ｑ_０）、ｆ（ｘ（Ｑ’））＝Ｅ（Ｑ’）、ｆ’（ｘ（Ｑ_０））＝０、ｆ’（ｘ（Ｑ’））＝（Ｑ＝Ｑ’における曲線Ｃの接線ｌの傾き）からなる４連立方程式をａ_０、ａ_１、ａ_２及びａ_３について解くことでｆ（ｘ）の各係数を求める。なお、ｘは一対の第０層原子間を結ぶ直線に沿った座標軸を示し、ｆ’（ｘ）はｆ（ｘ）の一次導関数（ｄｆ／ｄｘ）であり、ｆ’’（ｘ）はｆ（ｘ）の二次導関数（ｄ^２ｆ／ｄｘ^２）であるとする。また、ｆ（ｘ（Ｑ））及びｆ’（ｘ（Ｑ））は、位置Ｑに対応するｘにおけるｆ（ｘ）及びｆ’（ｘ）の値を示すものとする。そして、求めたａ_０、ａ_１、ａ_２及びａ_３を代入したｆ（ｘ）を用いてｂ＝ｆ’’（Ｑ_０）及びａ＝ｆ’’’（Ｑ_０）としてａ及びｂを求める。なお、ｆ’’’（ｘ）はｆ（ｘ）の三次導関数（ｄ^３ｆ／ｄｘ^３）であるとする。なお、フラグメントペアランキング取得部３０による二階の微分係数ｂ及び三階の微分係数ａを算出する機能及びステップが、本発明における微分係数算出手段の機能及び微分係数算出ステップに対応する。

フラグメントペアランキング取得部３０は、上記のように取得した二階微分係数ｂ及び三階微分係数ａの少なくともいずれかによって表される数式４―１～４－４のζ^Second、ζ^Third、ζ^Scaled及びζ^Averageのいずれかに基づき、フラグメントペアのランキングを決定する。ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageのいずれを使用するかは、入力装置を通じたユーザ入力に従って選択される。ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageのいずれかのうち、選択されたものをζとする。このとき、フラグメントペアランキング取得部３０は、ζが小さいフラグメントペアほどランキングを高く設定する。取得されたランキングは記憶部１０に記憶される。

（数式４－１）

（数式４－２）

（数式４－３）

（数式４－４）

ζ^Secondが小さい方向は、Ｅ（Ｑ）の局所的な上昇が小さい方向と言える。このため、ζ^Secondが小さいフラグメントペアについて反応経路探索を実行することで、ＥＱ付近におけるＥ（Ｑ）の上昇が小さい経路を得ることができると期待される。ζ^Thirdが小さい方向は、Ｅ（Ｑ）において下に凸から上に凸への変化が起こり、遷移状態が見つかると期待される。このとき、三次関数上での極大点の高さはａとｂのバランスで決まる。ζ^Thirdは、ｂが大きな方向においてＥ（Ｑ）が急激に下に凸から上に凸に変化する際に、特に小さくなる。そのため、ζ^Thirdが小さなフラグメントペアについて反応経路探索を実行することにより、化学結合のような硬い結合が低障壁で組み変わる経路を得ることができると期待される。一方、ζ^Scaledが小さいほど、三次関数上の極大点は低くなる。そのため、ζ^Scaledが小さなフラグメントペアについて反応経路探索を実行することにより、内部回転や水素結合組み換えなどの障壁の低い経路を得ることができると期待される。多くの化学反応では、コンフォメーション変化などの低障壁のイベントがナノ秒やマイクロ秒といった短い時間スケールで繰り返され、局所的な熱平衡が達成される。化学結合の組み換えのような高い障壁を経るいわゆるレア・イベントは、数秒や数時間などの長い時間スケールで進行する。多段階反応では、局所的な熱平衡化とレア・イベントが繰り返されることによって反応全体が進行する。このような状況を再現するには、ζ^Third値が小さなフラグメントペアとζ^Scaled値が小さなフラグメントペアの両方を計算する必要がある。そこで、ζ^Thirdとζ^Scaledの相乗平均を取ったζ^Averageが導入されている。ζ^Averageは、ζ^Third又はζ^Scaledが小さいフラグメントペアにおいて小さくなる。そのため、ζ^Averageが小さいフラグメントペアについて反応経路探索を実行することにより、内部回転、水素結合組み換え、化学結合組み換えなど、多様な経路を得ることができると期待される。

状態選択部４０（本発明における状態選択手段）は、反応経路探索の出発点となるＥＱを選択する。具体的には、状態選択部４０は、上述の数式４に示すΛ_ｉを算出すると共に、以下の数式５のν_ｉを算出する。そして、状態選択部４０は、ν_ｉが最も大きいＥＱ_ｉを、次に実行される反応経路探索の出発点として選択する。なお、状態選択部４０は、Λ_ｉを算出する際に、上記の通りＲＣＭＣ法に基づく速度定数行列の縮約演算を実行する。かかる状態選択部４０の機能が、本発明における速度定数縮約手段の機能に対応する。

（数式５）

ｎ_ｉは、記憶部１０が記憶している探索回数データによって示される。ｎ_ａｌｌは、その時点までの全てのＥＱ_ｉにおけるｎ_ｉの合計である。α及びβは、それぞれ、０以上１以下の実数である。例えば、α＝０．５、β＝０．５と設定されてもよい。αは、Λ_ｉが小さいＥＱ_ｉであってもそれらのＥＱ_ｉが何回かは選択されるようにする役割を担う。βは、αを含んだ項によって選択されようとするＥＱ_ｉを選択されにくくする役割を担う。ξ_ｉは、０以上１以下のランダムな実数である。なお、ξ_ｉは、数式５に基づくν_ｉの算出ごとに設定される。

反応経路算出部５０は、状態選択部４０が選択したＥＱを出発点のＥＱ（本発明における第１の平衡状態）とする反応経路探索を実行する。反応経路算出部５０は、このＥＱに関連付けて記憶部１０が記憶しているフラグメントデータが示すフラグメントペアに関して反応経路探索を実行する。この際、反応経路算出部５０は、かかるフラグメントペアのうち、ランキングが高いものをランキングが低いものに先行して反応経路探索を実行する。つまり、ランキングが高いフラグメントペアはランキングが低いフラグメントペアより先に反応経路探索が実行される。かかる反応経路探索が実行されたフラグメントペアに関し、反応経路探索が実行済みであることを示すように、記憶部１０に記憶されたフラグメントデータが更新される。

反応経路探索においては、上述の通り、ＡＦＩＲ経路及びＬＵＰ経路の算出、遷移状態の最適化、並びに、ＩＲＣの算出が実行される。これにより、状態選択部４０が選択したＥＱを出発点とする新たな反応経路が取得される。反応経路算出部５０は、記憶部１０の記憶内容を参照して、取得された反応経路が通る、又は到達するＥＱ（本発明における第２の平衡状態）において、取得済みＥＱのリストに含まれない新たなＥＱがある場合に、記憶部１０の記憶内容が示す取得済みＥＱのリストにその新たなＥＱを追加する。また、新たなＥＱに対応する原子の座標を示す構造データが記憶部１０の記憶内容に追加される。この新たなＥＱについては、フラグメントペア導出部２０がフラグメントペアを構成する原子の組み合わせを抽出すると共に、これによって導出されたフラグメントペアについてフラグメントペアランキング取得部３０がランキングを取得する。なお、反応経路算出部５０によるＡＦＩＲ経路を算出する機能及びステップが、本発明における構造変化算出手段の機能及び構造変化算出ステップに対応する。また、状態選択部４０が選択したＥＱを出発点とするＡＦＩＲ経路の算出により別のＥＱまでの遷移に対応する反応経路を算出する反応経路算出部５０の機能及びステップが、本発明における平衡状態間変化算出手段の機能及び平衡状態間変化算出ステップに対応する。

以下、反応経路探索システム１が実行する反応経路探索方法の一連の流れを図５に基づいて説明する。まず、ユーザによって与えられた初期状態ＥＱ_１について、フラグメントペア導出部２０及びフラグメントペアランキング取得部３０がフラグメントペア及びそのランキングを取得する（ステップＳ１）。初期状態ＥＱ_１における構造データは、例えば、ユーザが生成したデータファイルをシステムに取り込むことによって与えられる。

次に、反応経路算出部５０が、ステップＳ１において取得されたランキングのうち、最も高いランキングに対応するフラグメントペアについて、ＥＱ_１を出発点とする反応経路探索を実行する（ステップＳ２）。これによって取得された反応経路に基づき、新たなＥＱであるＥＱ_２が取得済みＥＱのリストに追加される。

次に、状態選択部４０が、取得済みＥＱのリストに含まれる各ＥＱ_ｉ（ｉ＝１，２，…）について数式５に示すν_ｉを算出する（ステップＳ３）。次に、フラグメントペアランキング取得部３０が、ステップＳ３で算出したν_ｉのうち、最も大きいものに対応するＥＱ_ｉについて、フラグメントペアのランキングが取得済みであるか否かを判定する（ステップＳ４）。この判定は、記憶部１０に記憶されたフラグメントデータを参照することで実施される。ランキングが取得済みであるとフラグメントペアランキング取得部３０が判定した場合には（ステップＳ４、Ｙｅｓ）、ステップＳ６の処理が実行される。ランキングが取得済みでない（つまり、当該ＥＱ_ｉが上述の対象ＥＱである）と判定した場合には（ステップＳ４、Ｎｏ）、ステップＳ３で算出したν_ｉのうち、最も大きいものに対応するＥＱ_ｉに係るフラグメントペアについて、フラグメントペアランキング取得部３０がランキングを取得する（ステップＳ５）。

次に、反応経路算出部５０が、ステップＳ３で算出したν_ｉのうち、最も大きいものに対応するＥＱ_ｉについて記憶部１０の記憶内容を参照し、反応経路の探索が未だ実行されていないフラグメントペアのうち、最もランキングが高いものに関して反応経路探索を実行する（ステップＳ６）。そして、かかる反応経路探索によって、取得済みＥＱのリストにない新たなＥＱが取得された場合には、その新たなＥＱに係るデータが記憶部１０の記憶内容に追加される。

次に、反応経路探索システム１は、反応経路の探索条件が満たされたか否かを判定する（ステップＳ７）。探索条件は、例えば、「過去ｎ回のＡＦＩＲ経路の計算において、トラフィック量Λ_ｉが大きい方からｍ個のＥＱが更新されていないこと」である。この場合、自然数ｎ及びｍはそれぞれユーザが設定可能としてもよい。かかる探索条件が満たされていないと判定した場合（ステップＳ７、Ｎｏ）、ステップＳ３からの処理が実行される。探索条件が満たされていると判定した場合（ステップＳ７、Ｙｅｓ）、反応経路探索システム１は一連の処理を終了する。

以上説明した本発明の反応経路探索システム１によると、ある平衡状態ＥＱに対応する原子の位置におけるＥ（Ｑ）の二階の微係数ｂ及び三階の微係数ａの少なくともいずれかの大きさに基づく優先順位の高いフラグメントペアに関して、優先順位の低いフラグメントペアに先行して反応経路探索が実行される。具体的には、上記数式４－１～４－４に示すζ^Second、ζ^Third、ζ^Scaled及びζ^Averageのいずれかからユーザが選択したζに基づいてフラグメントペアにランキングが設定される。そして、ランキングが高いフラグメントペアから順に反応経路探索が実行される。このようにζ^Second、ζ^Third、ζ^Scaled及びζ^Averageに基づくランキングが高いフラグメントペアについて優先的に反応経路探索を実行する場合、上記の通り、遷移状態が見つかりやすい、化学結合のような硬い結合が低障壁で組み変わる経路を得やすい、内部回転や水素結合組み換えなどの障壁の低い経路を得やすい、内部回転、水素結合組み換え、化学結合組み換えなど、多様な経路を得やすい等、目的となる反応経路を取得しやすくなる。これにより、計算コストをそれほど掛けることなく目的となる反応経路を取得できる可能性が高くなる。よって、原子数の増加に応じた計算コストの増加が抑制されやすい。

また、本実施形態においては、ν_ｉが最も大きいＥＱ_ｉが選択され、そのＥＱ_ｉについて反応経路探索が実行される。ν_ｉは、上記数式５に示されるように、ｎ_ａｌｌが比較的小さい場合に、数式５におけるα＾｛ｎ_ｉ／ｌｏｇ（ｎ_ａｌｌ）｝の値がΛ_ｉに対して比較的小さくなり、Λ_ｉの違いが表れやすい。したがって、トラフィック量が大きい平衡状態ほど計算対象として適切に選択されやすくなる。一方、ｎ_ａｌｌが比較的大きくなっていくと、数式５におけるα＾｛ｎ_ｉ／ｌｏｇ（ｎ_ａｌｌ）｝の値がΛ_ｉに対して比較的大きくなってくる。つまり、Λ_ｉの違いがν_ｉの違いに表れにくくなっていく。したがって、ｎ_ａｌｌが比較的大きくなると平衡状態の選択されやすさにΛ_ｉの違いが表れにくくなってくる。したがって、トラフィック量が小さい平衡状態も選択されるようになってくる。このように、上記構成によると、計算の進度に応じて適切に優先度の高い平衡状態が選択される。

［実施例］
以下、本発明の一実施例について説明する。本実施例は、最も基本的な有機合成反応であるストレッカー反応について、上述の実施形態に係る図５に示す反応経路探索方法を、一部工程を変更しつつ実施した。ストレッカー反応は、ケトン（またはアルデヒド）、アンモニア、および、シアン化水素が反応し、アミノニトリルが生成する反応である。本実施例では、ケトンとしてアセトンを用いた。また、水溶液中でのプロトン移動効率を考慮して、余計な水１分子を系に追加した。ＳＣ－ＡＦＩＲ法によって加える力の衝突エネルギーパラメータγは５００ｋＪ/ｍｏｌとした。分子同士が離れすぎないよう、全原子ペアに対してγ＝１００／｛Ｎ（Ｎ-１）／２｝ｋＪ／ｍｏｌとの弱い力を付加した（ここでのＮは系に含まれる原子数）。計算量を抑えるために、上述の工程変更として、ＬＵＰ経路の計算後の遷移状態の最適化とＩＲＣの算出は実施せず、ＬＵＰ経路のエネルギー極大点を近似遷移状態、ＬＵＰ経路をＥＱ同士のコネクションを定義する経路として使用した。ギブスエネルギーを算出する際の調和振動解析では、５０ｃｍ^-1以下の振動数を持つモードの振動数は５０ｃｍ^-1とした。電子状態計算は、Ｇａｕｓｓｉａｎ１６を用い、ＲＨＦ／ＳＶレベルで行った。

ν_ｉ値の計算では、２００Ｋ、３００Ｋ及び４００Ｋの３つの温度でＲＣＭＣ法を実行することで得たΛ_ｉ値の内、最も大きなものを用いた。なお、３つの温度全てでΛ_ｉ値がゼロとなった場合には、対応するν_ｉ値も強制的にゼロとした。また、反応時間は１日、気圧は１ａｔｍ、Ｔ_Ｒは３０００Ｋとした。ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageを用いて付与されたランキングに基づく優先順位で反応経路探索の対象となるフラグメントペアを選択する場合に加え、対象となるフラグメントペアをランダムに選択した場合についても反応経路探索を実行した。探索は、生成物が得られた時点で打ち切り、それまでに要した計算コストを、ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageを用いた場合の間で比較した。これらの計算では同一の初期構造（反応物の構造）と生成物の構造を用いた。表１はその比較結果である。表２は初期構造における各原子のｘ座標、ｙ座標及びｚ座標、表３は生成物の構造における各原子のｘ座標、ｙ座標及びｚ座標である。

（表１）

（表２）

（表３）

表１において、Ｎ^gradient、N^Hessian、Ｎ^EQ、Ｎ^pathはそれぞれ、生成物の取得までに行った勾配計算の数、Ｈｅｓｓｉａｎ計算（２階微分に係る計算）の数、生成物の取得までに得たＥＱの数、及び、経路の数に対応する。ここで、すべての計算において、イミニウムイオンへのシアノアニオンの求核付加によってアミノニトリルが生成する反応機構が得られていることを確認した。

表１から分かるように、ランダムにフラグメントペアを選択した場合に最も計算コストが高かった。つまり、本発明が有効であることが確認できた。ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageの中では、ζ^Averageが最も良いパフォーマンスを示した。これは、局所的な熱平衡化の経路とレア・イベントに対応する遅い過程の経路の両方が効率よくサンプリングされたためであると考えられる。反応物分子が弱く会合した状態では、１０ｋＪ／ｍｏｌ程度以下の非常に低い障壁で分子同士の配向が変化する経路が多数存在する。そういった経路はζ^Scaledを指標とすることで効率よく探索できる。一方、ヘミアミナール中間体から水分子が脱離してイミニウムイオンが生成する経路は１００ｋＪ／ｍｏｌ程度と比較的高い障壁を有し、ζ^Scaledを指標とした場合には、その計算は後回しにされてしまう。一方、ζ^Scaledとζ^Thirdの相乗平均であるζ^Averageを指標とした場合、ζ^Scaledが小さいフラグメントペアとζ^Thirdが小さいフラグメントペアの両方が探索される。そのため、ヘミアミナール中間体からイミニウムイオンへの変換のような、比較的障壁の高い化学結合の組み換え経路についても、優先探索が行われたものと考えられる。

本実施例で反応経路探索の対象としたのは、高々２０原子（Ｎ＝２０）からなる系である。一方、フラグメントの数はＮの二乗に比例して増加するため、ランダムとζ^Averageとのコスト比もＮと共に急激に増加すると予想される。本実施例では有機合成反応を例に数値検証を行ったが、化学結合の組み換えを伴わない過程、例えば、ポリマーのコンフォメーション変化や凝集構造の相転移などでは、ζ^Scaledを指標として低障壁の経路のみを優先探索する方が良いと予想される。そのため、ユーザが目的に応じてζ^Averageとζ^Scaledを使い分けることが推奨される。

［ＲＣＭＣ法の詳細］
ＲＣＭＣ法の詳細について説明する。ｎ番目の縮約における状態の集合をＮ^（ｎ）、状態ｉからｊへの速度定数をｋ_ｉ→ｊ ^（ｎ）、状態ｉのボルツマン分布をＰ_ｉ ^（ｎ）とする。つまり、Ｎ^（０）は反応経路ネットワークに含まれる全安定構造（ＥＱ）の集合、ｋ_ｉ→ｊ ^（０）は安定構造ｉからｊへの遷移の速度定数、Ｐ_ｉ ^（ｎ）は安定構造ｉのボルツマン分布に相当する。これらの量を用いた縮約操作は以下の通りである。

１．速度定数の上限ｋ_ＭＡＸ（＝１／ｔ_ＭＡＸ)を与え、ｎ＝０にセットする。
２．全ての素過程に対してｋ_ｉ→ｊ ^（０）を計算し、速度定数行列を得る。このとき、安定構造ｉとｊを直接つなぐ素過程が二つ以上ある場合、ｋ_ｉ→ｊ ^（０）はそれらに対する速度定数の和とする。安定構造ｉとｊを直接つなぐ素過程が存在しない場合、ｋ_ｉ→ｊ ^（０）はゼロとする。安定構造ｉとｉをつなぐ素過程に対するｋ_ｉ→ｉ ^（０）もゼロとする。
３．全安定構造のボルツマン分布Ｐ_ｉ ^（０）を求める。
４．全安定構造の（初期）分布Ｑ_ｉ ^（０）を与える。
５．χ_ｉｉ ^（０）＝１、χ_ｊｉ ^（０）＝０（ｊ≠ｉ）とする。以下の６～１４の処理がループとして繰り返し実行される。
６．全てのｉとｊのペアの中から、ｋ_ｉ→ｊ ^（ｎ）が最大のものを見つける。
７．もし、６で得た最大のｋ_ｉ→ｊ ^（ｎ）についてｋ_ｉ→ｊ ^（ｎ）＜ｋ_ＭＡＸであればループを抜けて終了する。
８．定常状態近似を状態ｉに適用する。その際、速度定数行列内の全ての速度定数を数式６と数式７を用いて補正する。
（数式６）

（数式７）

９．状態ｉが関連する速度定数を全てゼロにする（ｋ´_ｉ→ｍ ^{（ｎ＋１）}＝ｋ´_ｍ→ｉ ^{（ｎ＋１）}＝０）。
１０．数式８の通り、状態ｉのボルツマン分布をゼロとし（Ｐ_ｉ ^{（ｎ＋１）}＝０）、別の状態ｋ（ｋ≠ｉ）へ速度定数の重みで分配する。
（数式８）

１１．数式９の通り、状態ｉの分布をゼロとし（Ｑ_ｉ ^{（ｎ＋１）}＝０）、別の状態ｋ（ｋ≠ｉ）へ速度定数の重みで分配する。
（数式９）

１２．状態ｋへの安定構造ｊの寄与を数式１０を用いて更新する。
（数式１０）

１３．更新されたボルツマン分布を再現するために、ｋ´_ｋ→ｌ ^{（ｎ＋１）}を数式１１の通り補正する。
（数式１１）

１４．状態ｉを消去（重み付きで分配）した残りＮ^{（ｎ＋１）}を状態の集合とし、ステップ６へ戻る。

上記４及び１１におけるＱ_ｉ ^（０）及びＱ_ｉ ^（ｎ）が数式３のＱ^（ｍ） _ｊとして用いられる。また、上記５及び１２におけるχ_ｉｉ ^（０）、χ_ｊｉ ^（０）（ｊ≠ｉ）及びχ_ｊｉ ^（ｎ）が数式３のχ_ｊｉ ^（ｍ）として用いられる。

＜変形例＞
以上は、本発明の好適な実施形態についての説明であるが、本発明は上述の実施形態に限られるものではなく、課題を解決するための手段に記載された範囲の限りにおいて様々な変更が可能なものである。

例えば、本実施形態では、ζ^Second、ζ^Third、ζ^Scaled及びζ^Averageのいずれかからユーザが選択したζが小さいほどランキングが高くなるようにフラグメントペアごとにランキングが設定される。そして、ランキングが高いフラグメントペアから順に反応経路探索が実行される。しかし、ランキングの設定や反応経路探索を実行する順序は、その他の方法が採用されてもよい。例えば、ζが所定の値以下か否かに基づき、各フラグメントペアに２段階でランキングが設定されてもよい。そして、ランキングが高い方のフラグメントペアについてランダムな順に反応経路探索が実行された後に、ランキングが低い方のフラグメントペアについてランダムな順に探索が実行されてもよい。この場合にも、ランキングが高いフラグメントペアについてランキングが低いものに先行して探索が実行されることで、全てのフラグメントペアについてランダムな順で探索が実行される場合と比べて目的の反応経路が先行して取得されやすい。

また、上述の実施形態において、状態選択部４０が、上記数式５で表されるν_ｉが最も大きいＥＱ_ｉを反応経路探索の出発点として選択する。しかし、かかるν_ｉ以外の基準値に基づいて反応経路探索の出発点が選択されてもよい。この場合、Λ_ｉが大きい平衡状態ほど選択されやすいと共に、ｎ_ａｌｌが大きいほどΛ_ｉの違いが平衡状態の選択されやすさの違いに表れにくくなるように反応経路探索の出発点が選択されればよい。

１ 反応経路探索システム
１０ 記憶部
２０ フラグメントペア導出部
３０ フラグメントペアランキング取得部
４０ 状態選択部
５０ 反応経路算出部

Claims (7)

1. 複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出する反応経路探索システムとして１又は複数のコンピュータを機能させるプログラムであって、
   前記反応経路探索システムが、
   前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出手段と、
   （数式１）
   

   

   （数式２）
   

   

   互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出手段と、
   前記複数のフラグメントペアのうち、前記微分係数算出手段が算出したａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出手段に算出させる平衡状態間変化算出手段とを備えていることを特徴とする反応経路探索プログラム。
2. 前記反応経路探索システムが、
   前記平衡状態間変化算出手段による過去の算出結果に基づいて得られた全ての前記第２の平衡状態と初期状態とからなるＮ個（Ｎ：２以上の自然数）の前記平衡状態からいずれかを前記第１の平衡状態として選択する状態選択手段をさらに備えており、
   前記状態選択手段が選択した前記第１の平衡状態から前記第２の平衡状態までの遷移における前記フラグメントペアに関する前記構造の変化を前記平衡状態間変化算出手段が前記構造変化算出手段に算出させることを繰り返すことを特徴とする請求項１に記載の反応経路探索プログラム。
3. 前記反応経路探索システムが、
   前記Ｎ個の平衡状態間の遷移に係るＮ＊Ｎ個の速度定数からなるＮ行Ｎ列の速度定数行列から、前記Ｎ個の平衡状態の加重和で表されるｌ個（ｌ：ｌ＜Ｎを満たす自然数）の超状態間の遷移に係るｌ＊ｌ個の速度定数からなる速度定数行列を、ＲＣＭＣ法に基づいて速度定数行列を縮約することをｍ回（ｍ：ｍ＝Ｎ－ｌを満たす自然数）繰り返すことで取得する速度定数縮約手段を含んでおり、
   前記状態選択手段が、
   ｍ＝１，２，…，Ｍ（Ｍ：Ｍ＜Ｎを満たす自然数）のそれぞれについて、縮約された速度定数行列を前記速度定数縮約手段が取得した結果に基づいて、前記Ｎ個の平衡状態のそれぞれをＥＱ_ｉ（ｉ：Ｎ以下の自然数）とし、前記ｌ個の超状態のそれぞれをＳＳ_ｊ（ｊ：ｌ以下の自然数）とし、ｍ回の縮約後のＥＱ_ｉのポピュレーションをｐ_ｉ（ｍ）とし、ｍ回の縮約後のＳＳ_ｊのポピュレーションをＱ_ｊ（ｍ）とし、ｍ回の縮約後のＳＳ_ｊに対するＥＱ_ｉの寄与をχ_ｊｉ（ｍ）とし、ＥＱ_ｉの相対ギブズエネルギーをΔＧ_ｉとし、気体定数をＲとし、モデル温度パラメータをＴ_Ｒとするときに、数式３に基づいて算出されるｐ_ｉ（ｍ）及び数式４に基づいて算出されるΛ_ｉが大きいほど、ＥＱ_ｉが選択されやすくなるように、前記Ｎ個の平衡状態からいずれかを前記第１の平衡状態として選択することを特徴とする請求項２に記載の反応経路探索プログラム。
   （数式３）
   

   

   （数式４）
   

   
4. 前記状態選択手段が、
   Λ_ｉが大きいＥＱ_ｉほど選択されやすく、且つ、ＥＱ_ｉを前記第１の平衡状態として前記構造変化算出手段が前記構造の変化を算出した回数ｎ_ｉが小さいＥＱ_ｉほど選択されやすく、且つ、前記Ｎ個の平衡状態について前記構造変化算出手段が前記構造の変化を算出した全体の回数ｎ_ａｌｌが大きいほどΛ_ｉの違いが前記平衡状態の選択されやすさの違いに表れにくくなるように、前記Ｎ個の平衡状態からいずれかを前記第１の平衡状態として選択することを特徴とする請求項３に記載の反応経路探索プログラム。
5. 前記状態選択手段が、
   前記Ｎ個の平衡状態について前記構造変化算出手段が前記構造の変化を算出した全体の回数をｎ_ａｌｌとし、ＥＱ_ｉを前記第１の平衡状態として前記構造変化算出手段が前記構造の変化を算出した回数をｎ_ｉとし、ξ_ｉを０以上１以下の実数とし、α及びβをそれぞれ０以上１以下の実数とするときに、数式５のν_ｉが最も大きいＥＱ_ｉを前記第１の平衡状態として選択することを特徴とする請求項３又は４に記載の反応経路探索プログラム。
   （数式５）
   

   
6. 複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出するシステムであって、
   前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出手段と、
   （数式１）
   

   

   （数式２）
   

   

   互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出手段と、
   前記複数のフラグメントペアのうち、前記微分係数算出手段が算出したａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出手段に算出させる平衡状態間変化算出手段とを備えていることを特徴とする反応経路探索システム。
7. 複数個の原子からなる構造体における反応経路を、前記複数個の原子間の位置関係によって表される構造の変化として算出する方法であって、
   前記複数個の原子から抽出されたＮ_１個（Ｎ_１：自然数）の原子からなるフラグメントＡと前記複数個の原子から抽出された、前記フラグメントＡと別のＮ_２個（Ｎ_２：自然数）の原子からなるフラグメントＢとによって構成されるフラグメントペアに関し、前記複数個の原子の位置を示す幾何学的パラメータＱにおけるポテンシャルエネルギーをＥ（Ｑ）とし、前記フラグメントＡに属するｓ番目（ｓ：ｓ≦Ｎ_１を満たす自然数）の原子と前記フラグメントＢに属するｔ番目（ｔ：ｔ≦Ｎ_２を満たす自然数）の原子との距離をr_ｓｔとし、Ｒ_ｓ及びＲ_ｔをそれぞれ前記ｓ番目の原子及び前記ｔ番目の原子の共有結合半径とし、ρを＋１及び－１のいずれかとし、ｐ及びαを定数値としたときに、数式２に基づいて算出される数式１の関数Ｆ^ＡＦＩＲ（Ｑ）が最小値を取る前記構造の変化を算出する構造変化算出ステップと、
   （数式１）
   

   

   （数式２）
   

   

   互いに異なる原子の組み合わせを有する複数の前記フラグメントペアのそれぞれについて、Ｅ（Ｑ）が極小値を取る平衡状態の１つである第１の平衡状態に対応する前記複数個の原子の位置におけるＥ（Ｑ）の二階の微分係数ｂ及び三階の微分係数ａの少なくともいずれかを算出する微分係数算出ステップと、
   前記複数のフラグメントペアのうち、前記微分係数算出ステップにおいて算出されたａ及びｂの少なくともいずれかの大きさに基づく優先順位の高い前記フラグメントペアに関して、それより優先順位の低いその他の前記フラグメントペアに先行して、前記第１の平衡状態から他の前記平衡状態である第２の平衡状態までの遷移における前記構造の変化を前記構造変化算出ステップによって算出させる平衡状態間変化算出ステップとを備えていることを特徴とする反応経路探索方法。

Images