特許 7467892 分子の化学反応の解析方法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-04-08

(45)【発行日】2024-04-16

(54)【発明の名称】分子の化学反応の解析方法

(51)【国際特許分類】

   G16C 20/10 20190101AFI20240409BHJP

   G16C 10/00 20190101ALI20240409BHJP

【ＦＩ】

G16C20/10

G16C10/00

【請求項の数】 6

(21)【出願番号】P 2019209777

(22)【出願日】2019-11-20

(65)【公開番号】P2021082066

(43)【公開日】2021-05-27

【審査請求日】2022-09-16

(73)【特許権者】

【識別番号】000183233

【氏名又は名称】住友ゴム工業株式会社

(74)【代理人】

【識別番号】100104134

【弁理士】

【氏名又は名称】住友 慎太郎

(74)【代理人】

【識別番号】100156225

【弁理士】

【氏名又は名称】浦 重剛

(74)【代理人】

【識別番号】100168549

【弁理士】

【氏名又は名称】苗村 潤

(74)【代理人】

【識別番号】100200403

【弁理士】

【氏名又は名称】石原 幸信

(74)【代理人】

【識別番号】100206586

【弁理士】

【氏名又は名称】市田 哲

(72)【発明者】

【氏名】図師 知文

(72)【発明者】

【氏名】角田 昌也

【審査官】渡邉 加寿磨

(56)【参考文献】

【文献】特開２０００－１０７５９３（ＪＰ，Ａ）

【文献】特開２０１８－８１４７３（ＪＰ，Ａ）

【文献】特開２０１１－２２１８６８（ＪＰ，Ａ）

【文献】中国特許出願公開第１０３７９３６２２（ＣＮ，Ａ）

【文献】特開２００６－１９０２３４（ＪＰ，Ａ）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ１６Ｃ １０／００－９９／００

Ｇ０６Ｑ １０／００－９９／００

Ｇ１６Ｚ ９９／００

Ｂ０１Ｊ １９／００

(57)【特許請求の範囲】

【請求項1】

分子間の化学反応を解析するための方法であって、
コンピュータが、密度汎関数強束縛法に基づいて、前記分子間の化学反応の過程を示す反応経路を特定する工程と、
前記コンピュータに、前記分子をモデル化した複数の分子モデルを入力する工程と、
前記コンピュータに、前記特定された反応経路が分子動力学計算中で反映されるように、前記分子動力学計算用のパラメータを設定する工程と、
前記コンピュータが、前記分子モデルを対象として、前記パラメータを用いて前記分子動力学計算を行う工程とを含み、
前記パラメータは、隣接する前記分子モデルの離間距離に基づいて、前記隣接する分子モデルで化学反応が生じるか否かを判断するための反応距離を含む、
分子の化学反応の解析方法。

【請求項2】

前記分子動力学計算を行う工程は、前記離間距離が前記反応距離よりも大きい前記分子モデルの化学結合を開裂する工程と、
前記離間距離が前記反応距離以下である前記分子モデルを化学結合する工程とを含む、請求項１に記載の分子の化学反応の解析方法。

【請求項3】

前記パラメータは、隣接する前記分子モデルで化学反応が生じる確率を示す反応確率を含む、請求項１又は２記載の分子の化学反応の解析方法。

【請求項4】

分子間の化学反応を解析するための方法であって、
コンピュータが、密度汎関数強束縛法に基づいて、前記分子間の化学反応の過程を示す反応経路を特定する工程と、
前記コンピュータに、前記分子をモデル化した複数の分子モデルを入力する工程と、
前記コンピュータに、前記特定された反応経路が分子動力学計算中で反映されるように、前記分子動力学計算用のパラメータを設定する工程と、
前記コンピュータが、前記分子モデルを対象として、前記パラメータを用いて前記分子動力学計算を行う工程とを含み、
前記パラメータは、隣接する前記分子モデルで化学反応が生じる確率を示す反応確率を含む、
分子の化学反応の解析方法。

【請求項5】

前記分子動力学計算を行う工程は、前記反応確率に基づいて、前記隣接する分子モデルで化学反応が生じるか否かを判断する工程を含む、請求項４記載の分子の化学反応の解析方法。

【請求項6】

前記反応経路を特定する工程は、複数の反応経路を取得する工程と、
前記複数の反応経路について、活性化エネルギーをそれぞれ計算する工程と、
前記複数の反応経路のうち、前記活性化エネルギーが相対的に小さい反応経路を、前記反応経路として選択する工程とを含む、請求項１ないし５のいずれかに記載の分子の化学反応の解析方法。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、分子の化学反応の解析するための方法に関する。

【背景技術】

【0002】

近年、分子間の化学反応を解析するための方法として、例えば、密度汎関数強束縛法が知られている。この方法では、分子間の化学反応の過程を示す反応経路が取得される。関連する技術としては、下記の非特許文献１及び非特許文献２がある。

【先行技術文献】

【特許文献】

【0003】

【文献】西村 好史、中井 浩巳 著、「分割統治型密度汎関数強束縛分子動力学（DC-DFTB-MD）法によるナノスケール系化学反応シミュレーション」、アンサンブル、分子シミュレーション学会、２０１６年４月、Ｖｏｌ．１８、Ｎｏ．２、ｐ９５－１０１

【文献】西村 好史、中井 浩巳 著、「分割統治型密度汎関数強束縛分子動力学 (DC-DFTB-MD)法の最近の展開」、Journal of Computer Chemistry, Japan、日本コンピュータ化学会、２０１５年、Ｖｏｌ．１４、Ｎｏ．３、ｐ４３－４６

【発明の概要】

【発明が解決しようとする課題】

【0004】

密度汎関数強束縛法では、量子化学に基づいて、数十分子レベルの反応経路を、比較的精度良く追跡することができる。しかしながら、密度汎関数強束縛法は、数百以上といった多数の分子を対象とする大規模な化学反応を解析するには、現実的ではなかった。

【0005】

本発明は、以上のような実状に鑑み案出されたもので、未知の反応経路を有する大規模な分子間の化学反応を解析することが可能な方法を提供することを主たる目的としている。

【課題を解決するための手段】

【0006】

本発明は、分子間の化学反応を解析するための方法であって、コンピュータが、密度汎関数強束縛法に基づいて、前記分子間の化学反応の過程を示す反応経路を特定する工程と、前記コンピュータに、前記分子をモデル化した複数の分子モデルを入力する工程と、前記コンピュータに、前記特定された反応経路が分子動力学計算中で反映されるように、前記分子動力学計算用のパラメータを設定する工程と、前記コンピュータが、前記分子モデルを対象として、前記パラメータを用いて前記分子動力学計算を行う工程とを含むことを特徴とする。

【0007】

本発明に係る前記分子の化学反応の解析方法において、前記反応経路を特定する工程は、複数の反応経路を取得する工程と、前記複数の反応経路について、活性化エネルギーをそれぞれ計算する工程と、前記複数の反応経路のうち、前記活性化エネルギーが相対的に小さい反応経路を、前記反応経路として選択する工程とを含んでもよい。

【0008】

本発明に係る前記分子の化学反応の解析方法において、前記パラメータは、隣接する前記分子モデルの離間距離に基づいて、前記隣接する分子モデルで化学反応が生じるか否かを判断するための反応距離を含んでもよい。

【0009】

本発明に係る前記分子の化学反応の解析方法において、前記分子動力学計算を行う工程は、前記離間距離が前記反応距離よりも大きい前記分子モデルの化学結合を開裂する工程と、前記離間距離が前記反応距離以下である前記分子モデルを化学結合する工程とを含んでもよい。

【0010】

本発明に係る前記分子の化学反応の解析方法において、前記パラメータは、隣接する前記分子モデルで化学反応が生じる確率を示す反応確率を含んでもよい。

【0011】

本発明に係る前記分子の化学反応の解析方法において、前記分子動力学計算を行う工程は、前記反応確率に基づいて、前記隣接する分子モデルで化学反応が生じるか否かを判断する工程を含んでもよい。

【発明の効果】

【0012】

本発明の分子の化学反応の解析方法は、前記密度汎関数強束縛法に基づいて特定された前記反応経路を、前記分子動力学計算に反映させて、分子間の化学反応を解析することができる。前記密度汎関数強束縛法は、前記反応経路を高い精度で特定することができる。一方、前記分子動力学計算は、前記密度汎関数強束縛法に比べて、多数の分子モデルを対象とする大規模系での化学反応を、短時間で計算することができる。したがって、本発明の分子の化学反応の解析方法は、未知の反応経路を有する大規模な分子間の化学反応を、短時間かつ高い精度で計算することが可能となる。

【図面の簡単な説明】

【0013】

【図1】分子の化学反応の解析方法を実施するためのコンピュータの一例を示すブロック図である。

【図2】（ａ）は、ＭＢＴＳの構造式であり、（ｂ）は、八硫黄の構造式である。

【図3】分子の化学反応の解析方法の処理手順の一例を示すフローチャートである。

【図4】（ａ）は、ＭＢＴＳモデルを示す図であり、（ｂ）は、八硫黄モデルを示す図である。

【図5】（ａ）は、第１セルに配置されたＭＢＴＳモデル及び八硫黄モデルの一例を示す図であり、

【図6】（ａ）～（ｃ）は、反応経路の一例を説明する構造式である。

【図7】分子動力学計算での仮想空間である第２セル及び分子モデルの一例を示す概念図である。

【図8】（ａ）～（ｃ）は、化学結合の開裂、及び、化学結合の一例を説明するための概念図である。

【図9】分子動力学計算工程の処理手順の一例を示すフローチャートである。

【図10】本発明の他の実施形態の反応経路特定工程の処理手順の一例を示すフローチャートである。

【発明を実施するための形態】

【0014】

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の分子の化学反応の解析方法（以下、単に「解析方法」ということがある。）では、コンピュータを用いて、分子間の化学反応が解析される。本実施形態の解析方法では、未知の反応経路を有する大規模な分子間の化学反応が計算される。

【0015】

図１は、分子の化学反応の解析方法を実施するためのコンピュータの一例を示すブロック図である。本実施形態のコンピュータ１は、入力デバイスとしての入力部１１、出力デバイスとしての出力部１２、及び、種々の計算を実行する演算処理装置１３を有しており、分子の化学反応を解析するための装置（解析装置）として構成されている。

【0016】

入力部１１には、例えば、キーボード又はマウス等が用いられる。出力部１２には、例えば、ディスプレイ装置又はプリンタ等が用いられる。演算処理装置１３は、各種の演算を行う演算部（ＣＰＵ）１４、データやプログラム等が記憶される記憶部１５、及び、作業用メモリ１６を含んで構成されている。

【0017】

記憶部１５は、例えば、磁気ディスク、光ディスク又はＳＳＤ等からなる不揮発性の情報記憶装置である。記憶部１５には、データ部１７、及び、プログラム部１８が設けられている。

【0018】

データ部１７には、分子の化学反応を解析するのに必要な情報が記憶される。データ部１７は、分子に関する情報が記憶される初期データ記憶部１７ａ、分子モデルが記憶される分子モデル記憶部１７ｂ、計算条件が記憶される計算条件記憶部１７ｃ、分子間の反応経路が記憶される反応経路記憶部１７ｄ、分子動力学計算用のパラメータが記憶されるパラメータ記憶部１７ｅ、及び、解析結果が記憶される解析結果記憶部１７ｆを含んで構成されている。

【0019】

プログラム部１８は、演算部１４によって実行される分子の化学反応の解析プログラム（以下、単に「プログラム」ということがある。）である。このコンピュータ用のプログラム（プログラム部１８）は、コンピュータ１を、分子の化学反応の解析装置として機能させることができる。

【0020】

プログラム部１８は、分子間の反応経路を特定する反応経路特定部１８ａ、分子モデルを設定するための分子モデル設定部１８ｂ、分子モデル間のポテンシャルを定義するためのポテンシャル定義部１８ｃ、分子動力学計算用のパラメータを設定するためのパラメータ設定部１８ｄ、及び、分子動力学計算を行うための分子動力学計算部１８ｅを含んで構成される。このようなプログラム部１８によって実施される解析方法の処理手順については、後述する。

【0021】

本実施形態の解析方法において、化学反応が解析される分子については、特に限定されるわけではない。本実施形態の解析方法では、種々の分子の化学反応が解析されうる。本実施形態では、ジ－２－ベンゾチアゾリルジスルフィド（以下、単に「ＭＢＴＳ」ということがある。）、及び、環状の八硫黄（以下、単に「八硫黄」ということがある。）Ｓ₈の化学反応が解析される。なお、ＭＢＴＳは、チアゾール系加硫促進剤である。図２（ａ）は、ＭＢＴＳの構造式である、図２（ｂ）は、八硫黄Ｓ₈の構造式である。

【0022】

図３は、分子の化学反応の解析方法の処理手順の一例を示すフローチャートである。本実施形態の解析方法では、先ず、コンピュータ１が、分子間の反応経路を特定する（反応経路特定工程Ｓ１）。本実施形態の反応経路特定工程Ｓ１では、先ず、図１に示した初期データ記憶部１７ａに記憶されている分子に関する情報、及び、反応経路特定部１８ａが、作業用メモリ１６に読み込まれる。そして、反応経路特定部１８ａが、演算部１４によって実行される。この反応経路特定部１８ａの実行により、プログラム（プログラム部１８）は、コンピュータ１を、分子間の反応経路を特定するための手段として機能させることができる。

【0023】

反応経路とは、分子間の化学反応の過程を示すものである。本実施形態の反応経路特定工程Ｓ１では、密度汎関数強束縛法に基づいて、図２（ａ）に示したＭＢＴＳと、図２（ｂ）に示した八硫黄Ｓ₈との未知の反応経路が特定される。

【0024】

密度汎関数強束縛法は、密度汎関数法に基づく半経験的手法である。この密度汎関数強束縛法に基づく分子動力学計算（以下、「ＱＭ／ＭＤ計算」ということがある。）では、量子力学に基づく結合の生成・開裂を伴うダイナミクスの解析が可能となる。

【0025】

密度汎関数強束縛法の概要、計算条件、及び、計算方法等については、例えば、非特許文献１及び非特許文献２等に記載のように公知である。また、密度汎関数強束縛法で計算対象とする分子の個数については、コンピュータ１の性能等に基づいて適宜設定される。本実施形態の密度汎関数強束縛法に基づく計算では、２～５００個程度の分子を対象とした計算が行われる。なお、分子の個数が２個未満であると、反応経路を特定できないおそれがある。一方、分子の個数が５００個を超えると、計算時間が増大し、現実的な時間で計算を終了できないおそれがある。

【0026】

本実施形態の反応経路特定工程Ｓ１において、反応経路の特定には、例えば、密度汎関数強束縛法に基づくＱＭ／ＭＤ計算が可能なソフトウェアが用いられる。このソフトウェアの一例としては、例えば、ダッソー・システムズ株式会社製の「ＢＩＯＶＩＡ Ｍａｔｅｒｉａｌｓ Ｓｔｕｄｉｏ ＤＦＴＢ＋」、株式会社クロスアビリティ製の「ＤＣＤＦＴＢＭＤ」、及び、オープンソースプログラムの「ＤＦＴＢ＋」等である。このようなソフトウェアの少なくとも１つが、反応経路特定部１８ａ（図１に示す）に含まれている。

【0027】

本実施形態の反応経路特定工程Ｓ１では、先ず、図１に示したコンピュータ１（本例では、分子モデル記憶部１７ｂ）に、解析対象の分子（本例では、図２（ａ）及び（ｂ）に示したＭＢＴＳ及び八硫黄Ｓ₈）をモデル化した分子モデルが入力される。本実施形態では、図１に示した分子モデル設定部１８ｂが演算部１４によって実行されることにより、図１に示した初期データ記憶部１７ａに記憶されている分子に関数情報に基づいて、分子モデルが設定される。図４（ａ）及び（ｂ）は、分子モデル２の一例を示す図である。図４（ａ）は、ＭＢＴＳモデル３を示す図である。図４（ｂ）は、八硫黄モデル４を示す図である。

【0028】

本実施形態の分子モデル２は、ＭＢＴＳ（図２（ａ）に示す）をモデル化したＭＢＴＳモデル３（図４（ａ）に示す）と、八硫黄Ｓ₈（図２（ｂ）に示す）をモデル化した八硫黄モデル４（図４（ｂ）に示す）とが含まれる。これらの分子モデル２は、全原子モデルとして定義されており、複数の粒子モデル５と、粒子モデル５、５間を結合するボンドモデル６とを含んでいる。

【0029】

粒子モデル５は、前述のＱＭ／ＭＤ計算や、後述の分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粒子モデル５には、質量、直径、電荷、又は、初期座標などのパラメータが定義される。粒子モデル５は、炭素原子をモデル化した炭素粒子モデル５ｃ、水素原子をモデル化した水素粒子モデル５ｈ、窒素原子をモデル化した窒素粒子モデル５ｎ、及び、硫黄原子をモデル化した硫黄粒子モデル５ｓを含んでいる。

【0030】

ボンドモデル６は、粒子モデル５、５間を拘束するものである。本実施形態のボンドモデル６は、主鎖６ａと側鎖６ｂ（図４（ａ）に示す）とを含んでいる。主鎖６ａには、例えば、単結合や、二重結合が含まれてもよい。

【0031】

ボンドモデル６を介して隣り合う粒子モデル５、５間には、相互作用（斥力及び引力を含む）が生じさせるポテンシャル（図示省略）が定義される。これらのポテンシャルは、例えば、結合ポテンシャル、結合角ポテンシャル、及び、結合二面角ポテンシャルを含んでいる。このようなポテンシャルは、例えば、特許文献（特開２０１８－０３２０７７号公報）の記載に基づいて適宜定義することができる。これにより、分子モデル２（本例では、図４（ａ）に示したＭＢＴＳモデル３、及び、図４（ｂ）に示した八硫黄モデル４）が定義される。これらの分子モデル２は、図１に示したコンピュータ１（分子モデル記憶部１７ｂ）に記憶される。

【0032】

次に、本実施形態の反応経路特定工程Ｓ１では、図１に示したコンピュータ１（本例では、計算条件記憶部１７ｃ）に、密度汎関数強束縛法に基づくＱＭ／ＭＤ計算をするための計算条件が定義される。計算条件には、密度汎関数強束縛法に基づくエネルギー、温度、体積、圧力、原子の初速度、及び、ＭＤ計算の時間刻み幅等が含まれる。これらの計算条件は、前述のソフトウェアの仕様に応じて、適宜定義される。

【0033】

次に、本実施形態の反応経路特定工程Ｓ１では、ＱＭ／ＭＤ計算での仮想空間である第１セルの内部に、ＭＢＴＳモデル３及び八硫黄モデル４が配置される。図５（ａ）は、第１セル９に配置されたＭＢＴＳモデル３及び八硫黄モデル４の一例を示す図である。

【0034】

第１セル９は、少なくとも互いに向き合う一対の面２１、２１、本実施形態では、互いに向き合う三対の面２１、２１を有している。第１セル９は、直方体又は立方体（本例では、直方体）として定義されている。各面２１、２１には、周期境界条件が定義されている。このような第１セル９を用いた分子動力学計算では、例えば、一方側の面２１から出て行ったＭＢＴＳモデル３又は八硫黄モデル４の一部が、他方側の面２１から入ってくるように計算することができる。第１セル９の大きさは、例えば、ＭＢＴＳモデル３及び八硫黄モデル４の個数等に基づいて、適宜設定することができる。

【0035】

第１セル９の内部には、ＭＢＴＳモデル３及び八硫黄モデル４がランダムに配置される。ＭＢＴＳモデル３及び八硫黄モデル４の各個数の合計は、前述した密度汎関数強束縛法での分子の個数の範囲内となるように設定される。

【0036】

そして、本実施形態の反応経路特定工程Ｓ１では、図１に示したコンピュータ１（本例では、反応経路特定部１８ａ）が、図５（ａ）に示した分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）を対象に、前述した計算条件に基づいて、密度汎関数強束縛法に基づくＱＭ／ＭＤ計算が行われる。

【0037】

密度汎関数強束縛法に基づくＱＭ／ＭＤ計算は、例えば、一つの反応経路が取得されるまで行われてもよいし、任意の計算時間が経過するまで行われてもよい。図５（ｂ）は、ＱＭ／ＭＤ計算が終了した後の分子モデル２の一例を示す図である。

【0038】

本実施形態の密度汎関数強束縛法に基づくＱＭ／ＭＤ計算では、図４（ａ）及び図５（ａ）に示したＭＢＴＳモデル３での硫黄粒子モデル５ｓ－５ｓの化学結合が開裂して、図５（ｂ）に示した一対のＢｔ－Ｓモデル７が生成される。一方、図４（ｂ）及び図５（ｂ）に示した八硫黄モデル４は、硫黄粒子モデル５ｓ－５ｓの化学結合が開裂して、図５（ｂ）に示した二硫黄モデル８が生成される。そして、図５（ｂ）に示されるように、Ｂｔ－Ｓモデル７の硫黄粒子モデル５ｓと、二硫黄モデル８の硫黄粒子モデル５ｓとが化学結合し、Ｂｔ－Ｓ4－Ｂｔモデル２０が生成される。

【0039】

図６（ａ）～（ｃ）は、反応経路の一例を説明する構造式である。本実施形態の反応経路特定工程Ｓ１では、上述の密度汎関数強束縛法に基づくＱＭ／ＭＤ計算により、図６（ａ）～（ｃ）に示されるように、ＭＢＴＳ及び八硫黄Ｓ₈が化学反応して、Ｂｔ－Ｓ4－Ｂｔが生成されるまでの反応経路を特定することができる。このように、本実施形態の解析方法では、密度汎関数強束縛法に基づいて、分子間の未知の反応経路を特定することができる。特定された反応経路は、図１に示したコンピュータ１（本例では、反応経路記憶部１７ｄ）に記憶される。

【0040】

次に、本実施形態の解析方法では、図１に示したコンピュータ１（本例では、分子モデル記憶部１７ｂ）に、分子（図２（ａ）及び（ｂ）に示す）をモデル化した複数の分子モデルが入力される（工程Ｓ２）。工程Ｓ２では、先ず、初期データ記憶部１７ａに記憶されている分子に関する情報、及び、分子モデル設定部１８ｂが、作業用メモリ１６に読み込まれる。そして、工程Ｓ２では、分子モデル設定部１８ｂが、演算部１４によって実行される。この分子モデル設定部１８ｂの実行により、プログラム（プログラム部１８）は、コンピュータ１を、分子モデルを設定するための手段として機能させることができる。

【0041】

工程Ｓ２で入力される分子モデルは、後述の分子動力学計算に用いるためのものである。本実施形態の工程Ｓ２では、図２（ａ）及び（ｂ）に示した解析対象のＭＢＴＳ、及び、八硫黄Ｓ₈がモデル化される。

【0042】

分子モデルは、複数の原子の集団を一つのビーズに置き換えた粗視化分子モデル（図示省略）ではなく、後述の分子動力学計算において、化学反応（即ち、分子を構成する原子の結合の組換えによる変化）の詳細な解析が可能な全原子モデルとして定義される。本実施形態では、図４（ａ）及び（ｂ）に示されるように、反応経路特定工程Ｓ１で定義された分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）が用いられる。分子モデル２の詳細については、前述のとおりである。

【0043】

図７は、分子動力学計算での仮想空間である第２セル１０及び分子モデル２の一例を示す概念図である。本実施形態の工程Ｓ２では、第２セル１０の内部に、複数の分子モデル２が配置される。図７では、分子モデル２の一部を代表して示している。

【0044】

第２セル１０は、図５（ａ）に示した第１セル９と同様に設定することができる。第２セル１０の大きさは、例えば、分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）の個数等に基づいて、適宜設定することができる。

【0045】

分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）は、第２セル１０の内部に、ランダムに配置される。後述の分子動力学計算では、前述の密度汎関数強束縛法に比べて、多数の分子モデル２を対象とする大規模系での化学反応を、短時間で計算することができる。このため、工程Ｓ２では、前述の密度汎関数強束縛法で計算対象とする分子の個数に比べて、多くの分子モデル２が入力される。

【0046】

分子モデル２の個数については、後述の分子動力学計算を行うコンピュータ１の性能等に基づいて、適宜設定することができる。なお、分子モデル２の個数が多いと、計算時間が増大し、現実的な時間で計算を終了できないおそれがある。一方、分子モデル２の個数が少ないと、大規模系での分子間の化学反応を解析できないおそれがある。このような観点より、本実施形態の分子モデル２の個数（ＭＢＴＳモデル３の個数と八硫黄モデル４の個数との合計）は、４～１０００個に設定される。

【0047】

ＭＢＴＳモデル３の個数と、八硫黄モデル４の個数との比については、例えば、実際に予定されているＭＢＴＳ及び八硫黄Ｓ₈の配合割合等を考慮して、適宜設定することができる。複数の分子モデル２及び第２セル１０は、図１に示したコンピュータ１（本例では、分子モデル記憶部１７ｂ）に入力される。

【0048】

本実施形態の解析方法では、図３に示されるように、反応経路特定工程Ｓ１の後に、複数の分子モデルを入力する工程Ｓ２が実施されているが、特に限定されない。例えば、工程Ｓ２の後に、反応経路特定工程Ｓ１が実施されてもよい。この場合、反応経路特定工程Ｓ１では、工程Ｓ２で定義された分子モデル２を用いることができる。

【0049】

次に、本実施形態の解析方法では、隣接する分子モデル２間に、ポテンシャルＰが定義される（工程Ｓ３）。工程Ｓ３では、先ず、図１に示されるように、分子モデル記憶部１７ｂに記憶されている分子モデル２（図７に示す）、及び、ポテンシャル定義部１８ｃが、作業用メモリ１６に読み込まれる。そして、工程Ｓ３では、ポテンシャル定義部１８ｃが、演算部１４によって実行される。このポテンシャル定義部１８ｃの実行により、プログラム（プログラム部１８）は、コンピュータ１を、ポテンシャルＰ（図７に示す）を定義するための手段として機能させることができる。

【0050】

図７に示されるように、ポテンシャルＰは、後述の分子動力学計算において、隣接する分子モデル２の粒子モデル５、５間に、斥力及び引力を含む相互作用を定義するためのものである。ポテンシャルＰについては、適宜定義することができる。本実施形態のポテンシャルＰとしては、例えば、特許文献（特開２０１８－０３２０７７号公報）の記載に基づいて、ＬＪポテンシャルが定義される。ＬＪポテンシャルのパラメータについては、前述の特許文献の記載に基づいて、適宜定義することができる。ポテンシャルＰは、図１に示したコンピュータ１（計算条件記憶部１７ｃ）に記憶される。

【0051】

次に、本実施形態の解析方法では、図１に示したコンピュータ１に、特定された反応経路が、後述の分子動力学計算中で反映されるように、分子動力学計算用のパラメータが設定される（工程Ｓ４）。工程Ｓ４では、先ず、図１に示されるように、反応経路記憶部１７ｄに記憶されている分子間の反応経路（図６に示す）、分子モデル記憶部１７ｂに記憶されている分子モデル２（図７に示す）、及び、パラメータ設定部１８ｄが、作業用メモリ１６に読み込まれる。そして、工程Ｓ４では、パラメータ設定部１８ｄが演算部１４によって実行される。このパラメータ設定部１８ｄの実行により、プログラム（プログラム部１８）は、コンピュータ１を、分子動力学計算用のパラメータを設定する手段として機能させることができる。

【0052】

パラメータについては、反応経路特定工程Ｓ１で特定された反応経路（図６に示す）が、分子動力学計算中で反映されるものであれば、適宜設定することができる。本実施形態のパラメータには、反応距離及び反応確率が含まれる。なお、パラメータは、反応距離及び反応確率のいずれか一方のみでもよいし、その他のパラメータがさらに含まれてもよい。図８（ａ）～（ｃ）は、化学結合の開裂、及び、化学結合の一例を説明するための概念図である。

【0053】

図６（ａ）及び（ｂ）に示されるように、本実施形態の反応経路では、ＭＢＴＳの硫黄原子Ｓ－Ｓ間の化学結合が開裂している。さらに、本実施形態の反応経路では、八硫黄Ｓ₈の硫黄原子Ｓ－Ｓ間の化学結合も開裂している。このため、本実施形態の工程Ｓ４では、図８（ａ）に示されるように、ＭＢＴＳモデル３の一対の硫黄粒子モデル５ｓ－５ｓ間、八硫黄モデル４の硫黄粒子モデル５ｓ－５ｓ間に、分子動力学計算用のパラメータ（本例では、反応距離及び反応確率）が設定される。このパラメータは、図８（ｂ）に示されるように、ＭＢＴＳモデル３が開裂したＢｔ－Ｓモデル７、及び、八硫黄モデル４が開裂した二硫黄モデル８において、無効となる（即ち、後述の分子動力学計算において、反映されなくなる）ように定義される。

【0054】

さらに、図６（ｂ）及び（ｃ）に示されるように、本実施形態の反応経路では、Ｂｔ－Ｓの硫黄原子Ｓと、二硫黄Ｓ₂の硫黄原子Ｓとが化学結合している。このため、図８（ｂ）に示されるように、本実施形態の工程Ｓ４では、Ｂｔ－Ｓモデル７の硫黄粒子モデル５ｓと、二硫黄モデル８の硫黄粒子モデル５ｓとの間に、分子動力学計算用のパラメータ（本例では、反応距離及び反応確率）が設定される。このパラメータは、後述の分子動力学計算において、ＭＢＴＳモデル３（図８（ａ）に示す）が開裂したＢｔ－Ｓモデル７、及び、八硫黄モデル４（図８（ａ）に示す）が開裂した二硫黄モデル８が生成されてから有効となるように定義される。さらに、このパラメータは、Ｂｔ－Ｓモデル７と二硫黄モデル８とが化学結合したＢｔ－Ｓ4－Ｂｔモデル２０（図８（ｃ）に示す）において、無効となる（即ち、後述の分子動力学計算において、反映されなくなる）ように定義される。

【0055】

反応距離は、隣接する分子モデル２の離間距離に基づいて、その隣接する分子モデル２で化学反応（本例では、化学結合、及び、化学結合の開裂）が生じるか否かを判断するためのものである。本実施形態の離間距離は、互いに化学結合された分子モデル２の離間距離と、互いに化学結合されていない分子モデル２の離間距離とが含まれる。

【0056】

互いに化学結合された分子モデル２の離間距離には、図８（ａ）に示されるように、ＭＢＴＳモデル３内でボンドモデル６を介して隣接するＢｔ－Ｓモデル７、７の第１離間距離Ｌ１、及び、八硫黄モデル４内でボンドモデル６を介して隣接する二硫黄モデル８、８の第２離間距離Ｌ２が含まれる。一方、互いに化学結合されていない分子モデル２の離間距離には、図８（ｂ）に示されるように、Ｂｔ－Ｓモデル７と、二硫黄モデル８との間の第３離間距離Ｌ３が含まれる。

【0057】

反応距離は、前述の離間距離（本例では、第１離間距離Ｌ１～第３離間距離Ｌ３）に対応して定義される。本実施形態の反応距離は、第１離間距離Ｌ１に対応する第１反応距離Ｒ１（図示省略）、第２離間距離Ｌ２に対応する第２反応距離Ｒ２（図示省略）、及び、第３離間距離Ｌ３に対応する第３反応距離Ｒ３（図示省略）を含んでいる。本実施形態の離間距離Ｌ１～Ｌ３及び反応距離Ｒ１～Ｒ３は、化学結合の開裂、及び、化学結合される粒子モデル５（本実施形態では、硫黄粒子モデル５ｓ）の中心間の距離として定義される。

【0058】

図８（ａ）に示されるように、本実施形態では、後述の分子動力学計算において、ＭＢＴＳモデル３内で互いに隣接する一方のＢｔ－Ｓモデル７と、他方のＢｔ－Ｓモデル７との間の離間距離（硫黄粒子モデル５ｓ、Ｓの第１離間距離）Ｌ１が、それらの間に定義された第１反応距離Ｒ１（図示省略）よりも大きくなった場合に、図８（ｂ）に示されるように、Ｂｔ－Ｓモデル７、７の化学結合（ＭＢＴＳモデル）を開裂させる。この化学結合の開裂は、例えば、硫黄粒子モデル５ｓ、５ｓ間のボンドモデル６（図８（ａ）に示す）を無効（又は削除）することによって定義することができる。

【0059】

同様に、図８（ａ）に示されるように、八硫黄モデル４を構成する一方の二硫黄モデル８と、他方の二硫黄モデル８との間の離間距離（硫黄粒子モデル５ｓの第２離間距離）Ｌ２が、それらの間に定義された第２反応距離Ｒ２（図示省略）よりも大きくなった場合に、図８（ｂ）に示されるように、二硫黄モデル８、８の化学結合（八硫黄モデル４）を開裂させる。この化学結合の開裂の定義については、前述のとおりである。

【0060】

図８（ｂ）に示されるように、本実施形態では、後述の分子動力学計算において、隣接するＢｔ－Ｓモデル７と、二硫黄モデル８との間の離間距離（Ｂｔ－Ｓモデル７の硫黄粒子モデル５ｓと、二硫黄モデル８の硫黄粒子モデル５ｓとの離間距離）Ｌ３が、それらの間に定義された第３反応距離Ｒ３（図示省略）以下である場合に、図８（ｃ）に示されるように、Ｂｔ－Ｓモデル７と、二硫黄モデル８とを化学結合させる。この化学結合は、例えば、Ｂｔ－Ｓモデル７の硫黄粒子モデル５ｓと、二硫黄モデル８の硫黄粒子モデル５ｓとの間に、ボンドモデル６（図８（ｃ）に示す）が設定されることによって定義することができる。

【0061】

各反応距離Ｒ１～Ｒ３（図示省略）については、特定された反応経路（図６（ａ）～（ｃ）に示す）等を考慮して、適宜設定することができる。本実施形態の反応距離Ｒ１～Ｒ３は、化学結合、及び、化学結合の開裂が生じる原子（本例では、硫黄原子）の結合距離の２～６倍（本例では、４倍）にそれぞれ設定されている。第１反応距離Ｒ１～第３反応距離Ｒ３は、図１に示したコンピュータ１（パラメータ記憶部１７ｅ）に入力される。

【0062】

反応確率は、隣接する分子モデル２、２で化学反応が生じる確率を示している。したがって、後述の分子動力学計算では、反応確率が相対的に大きい分子モデル２、２において、化学反応を生じやすくなる。反応確率については、特定された反応経路（図６（ａ）～（ｃ）に示す）等を考慮して、適宜設定することができる。本実施形態では、下記式（１）で定義されるアレニウスの反応確率ｐが用いられる。

【0063】

【数1】

ここで、
Ｅ_a：活性化エネルギー
Ｒ：気体定数
Ｔ：温度
Ａ：頻度因子

【0064】

上記式（１）において、活性化エネルギーＥ_aは、化学反応させる分子モデル２ごとに取得される。本実施形態では、図８（ａ）に示されるように、ＭＢＴＳモデル３を構成する一方のＢｔ－Ｓモデル７と、他方のＢｔ－Ｓモデル７との間の活性化エネルギーＥ_a、八硫黄モデル４を構成する一方の二硫黄モデル８と、他方の二硫黄モデル８との間の活性化エネルギーＥ_a、及び、図８（ｂ）に示した隣接するＢｔ－Ｓモデル７と、二硫黄モデル８との間の活性化エネルギーＥ_aがそれぞれ取得される。

【0065】

この活性化エネルギーＥ_aは、例えば、実験やシミュレーション等に基づいて、適宜取得することができる。本実施形態では、ＤＦＴ法に基づく量子化学計算に基づいて、活性化エネルギーＥ_aを取得している。

【0066】

ＤＦＴ法に基づく量子化学計算で用いられる基底関数としては、例えば、基底関数（6-31G(d)）が定義される。ＤＦＴ法に基づく基底関数の定義、及び、量子化学計算には、量子化学計算プログラム（例えば、Gaussian社製のGaussian03、Gaussian09又はGaussian16等）が用いられる。なお、ＤＦＴ法の概要や計算方法については、例えば、論文（田中章夫、前川健典、鈴木机倫 著「反応機構解析における理論計算の役割」、住友化学技術誌２０１３、住友化学株式会社、２０１３年７月３１日、ｐ４３－５２）に記載されているとおりである。

【0067】

上記式（１）において、頻度因子Ａは、化学反応が生じるか否かを判断する時間間隔に基づいて、適宜設定することができる。本実施形態の頻度因子Ａは、例えば、１×１０^５～１×１０^２０に設定される。上記式（１）において、温度Ｔは、後述の分子動力学計算での分子モデル２の温度である。

【0068】

本実施形態の反応確率ｐは、図８（ａ）に示したＭＢＴＳモデル３を構成する一方のＢｔ－Ｓモデル７と、他方のＢｔ－Ｓモデル７との間の第１反応確率、八硫黄モデル４を構成する一方の二硫黄モデル８と、他方の二硫黄モデル８との間の第２反応確率、及び、図８（ｂ）に示した隣接するＢｔ－Ｓモデル７と、二硫黄モデル８との間の第３反応確率が含まれる。これらの反応確率は、図１に示したコンピュータ１（パラメータ記憶部１７ｅ）に記憶される。

【0069】

次に、本実施形態の解析方法では、コンピュータ１（図１に示す）が、図７に示した分子モデル２（本例では、ＭＢＴＳモデル３及び八硫黄モデル４）を対象として、パラメータ（本例では、反応距離及び反応確率）を用いて分子動力学計算を行う（分子動力学計算工程Ｓ５）。

【0070】

分子動力学計算工程Ｓ５では、先ず、図１に示されるように、分子モデル記憶部１７ｂに記憶されている分子モデル２及び第２セル１０（図７に示す）、計算条件記憶部１７ｃに記憶されている計算条件（図７に示したポテンシャルＰ等）、パラメータ記憶部１７ｅに記憶されているパラメータ（反応距離及び反応確率）、並びに、分子動力学計算部１８ｅが、作業用メモリ１６に読み込まれる。そして、分子動力学計算工程Ｓ５では、分子動力学計算部１８ｅが演算部１４によって実行される。この分子動力学計算部１８ｅの実行により、プログラム（プログラム部１８）は、コンピュータ１を、分子動力学計算を行う手段として機能させることができる。

【0071】

分子動力学計算では、例えば、図７に示した第２セル１０について所定の時間、分子モデル２が古典力学に従うものとして、ニュートンの運動方程式が適用される。第２セル１０において、圧力及び温度が一定、又は、体積及び温度が一定に保たれる。そして、各時刻での粒子モデル５の動きが、単位時間ステップ毎に追跡される。これにより、分子動力学計算は、分子モデル２を、実際の分子の動きに近似させることができる。

【0072】

分子動力学計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ－ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。図９は、分子動力学計算工程Ｓ５の処理手順の一例を示すフローチャートである。

【0073】

本実施形態の分子動力学計算工程Ｓ５では、先ず、図７に示した分子モデル２（本例では、ＭＢＴＳモデル３及び八硫黄モデル４）を対象とする分子動力学計算が開始され（工程Ｓ５１）、隣接する分子モデル２、２の離間距離（図８（ａ）～（ｂ）に示した離間距離Ｌ１～Ｌ３）と、反応距離（図示しない反応距離Ｒ１～Ｒ３）とが比較される（工程Ｓ５２）。

【0074】

工程Ｓ５２では、図７に示した第２セル１０に配置されている全ての分子モデル２から選択された隣接する分子モデル２、２の離間距離（図８（ａ）～（ｂ）に示した離間距離Ｌ１～Ｌ３のいずれか）と、反応距離（図示しない反応距離Ｒ１～Ｒ３のいずれか）とが比較される。前述したように、本実施形態の離間距離は、互いに化学結合された分子モデル２の第１離間距離Ｌ１及び第２離間距離Ｌ２（図８（ａ）に示す）と、互いに化学結合されていない分子モデル２の第３離間距離Ｌ３（図８（ｂ）に示す）とが含まれている。

【0075】

工程Ｓ５２では、図８（ａ）に示されるように、選択された隣接する分子モデル２が互いに化学結合されている場合に、隣接する分子モデル２の離間距離（図８（ａ）に示した第１離間距離Ｌ１及び第２離間距離Ｌ２）が、反応距離よりも大きいか否かが判断される。

【0076】

工程Ｓ５２において、ＭＢＴＳモデル３内で互いに隣接する（化学結合された）一方のＢｔ－Ｓモデル７（図８（ａ）において二点鎖線で示す）と、他方のＢｔ－Ｓモデル７（図８（ａ）において二点鎖線で示す）との間の第１離間距離Ｌ１が、第１反応距離Ｒ１（図示省略）よりも大きい場合（工程Ｓ５２で、「＞反応距離」）に、次の工程Ｓ５３が実施される。

【0077】

本実施形態の工程Ｓ５３では、前述の反応確率（本例では、第１反応確率）に基づいて、隣接する分子モデル２、２（本例では、ＭＢＴＳモデル３内で互いに隣接する一対のＢｔ－Ｓモデル７、７）で化学反応（本例では、図８（ｂ）に示した化学結合の開裂）が生じるか否かが判断される。

【0078】

本例において、化学反応が生じるか否かの判断は、第１反応確率（０～１の範囲）と、乱数（０～１の範囲）とが比較される。乱数は、例えば、モンテカルロ法に基づいて取得される。そして、工程Ｓ５３において、第１反応確率が乱数よりも大きい場合に、図８（ａ）に示した一対のＢｔ－Ｓモデル７、７間で化学反応が生じると判断され（工程Ｓ５３で、「Ｙ」）、図８（ｂ）に示されるように、分子モデル２（一対のＢｔ－Ｓモデル７、７）の化学結合を開裂させる（工程Ｓ５４）。

【0079】

一方、工程Ｓ５３において、第１反応確率が乱数以下である場合には、分子モデル２（図８（ａ）に示した一対のＢｔ－Ｓモデル７、７）で化学反応が生じないと判断され（工程Ｓ５３で、「Ｎ」）、次の工程Ｓ５５が実施される。

【0080】

同様に、工程Ｓ５２において、八硫黄モデル４内の一方の二硫黄モデル８（図８（ａ）において二点鎖線で示す）と、他方の二硫黄モデル８（図８（ａ）において二点鎖線で示す）との間の離間距離Ｌ２が、第２反応距離Ｒ２（図示省略）よりも大きい場合（工程Ｓ５２で、「＞反応距離」）に、次の化学反応（本例では、化学結合の開裂）が生じるか否かを判断する工程Ｓ５３が実施される。

【0081】

工程Ｓ５３において、第２反応確率が、乱数よりも大きい場合に、化学反応が生じると判断され（工程Ｓ５３で、「Ｙ」）、図８（ｂ）に示されるように、分子モデル２（一対の二硫黄モデル８、８）の化学結合を開裂させる（工程Ｓ５４）。一方、工程Ｓ５３において、第２反応確率が乱数以下である場合には、分子モデル２（図８（ａ）に示した一対の二硫黄モデル８、８）で化学反応が生じないと判断され（工程Ｓ５３で、「Ｎ」）、次の工程Ｓ５５が実施される。

【0082】

また、工程Ｓ５２では、図８（ｂ）に示されるように、選択された隣接する分子モデル２、２との間で化学結合されていない場合に、隣接する分子モデル２の離間距離（図８（ｂ）に示した第３離間距離Ｌ３）が、反応距離以下であるか否かが判断される。

【0083】

工程Ｓ５２において、図８（ｂ）に示した隣接するＢｔ－Ｓモデル７と、二硫黄モデル８との間の第３離間距離Ｌ３が、第３反応距離Ｒ３（図示省略）以下である場合（工程Ｓ５２で、「≦反応距離」）に、次の工程Ｓ５６が実施される。

【0084】

本実施形態の工程Ｓ５６では、前述の反応確率（本例では、第３反応確率）に基づいて、隣接する分子モデル２、２（本例では、隣接するＢｔ－Ｓモデル７、及び、二硫黄モデル８）で化学反応（本例では、化学結合）が生じるか否かが判断される。化学反応が生じるか否かについては、工程Ｓ５３と同様に判断することができる。

【0085】

工程Ｓ５６において、第３反応確率が乱数よりも大きい場合に、化学反応（化学結合）が生じると判断され（工程Ｓ５６で、「Ｙ」）、図８（ｃ）に示されるように、分子モデル２（Ｂｔ－Ｓモデル７及び二硫黄モデル８）を化学結合させる（工程Ｓ５７）。一方、工程Ｓ５６において、第３反応確率が乱数以下である場合に、化学反応（化学結合）が生じないと判断されて（工程Ｓ５６で、「Ｎ」）、次の工程Ｓ５５が実施される。

【0086】

また、工程Ｓ５２では、前述の条件を満たさない場合（工程Ｓ５２で、「その他」）、次の工程Ｓ５５が実施される。

【0087】

工程Ｓ５５では、第２セル１０に配置されている全ての分子モデル２について、隣接する分子モデル２、２の離間距離（本例では、図８（ａ）～（ｂ）に示した離間距離Ｌ１～Ｌ３）が比較されたか否かが判断される。工程Ｓ５５において、全ての分子モデル２の離間距離と、反応距離とが比較されたと判断された場合（工程Ｓ５５で、「Ｙ」）、次の工程Ｓ５８が実施される。

【0088】

一方、工程Ｓ５５において、全ての分子モデル２の離間距離が、反応距離と比較されていないと判断された場合（工程Ｓ５５で、「Ｎ」）、他の分子モデル２が選択され（工程Ｓ５９）、工程Ｓ５２～工程Ｓ５７が再度実施される。これにより、分子動力学計算工程Ｓ５では、第２セル１０に配置されている全ての分子モデル２について、隣接する分子モデル２、２の離間距離と、反応距離とを比較して、化学結合の開裂、及び、化学結合を計算することができるため、図６（ａ）～（ｃ）に示した反応経路を、分子動力学計算中に反映させることができる。

【0089】

次に、工程Ｓ５８では、予め定められた終了時間が経過したか否かが判断される。計算終了時間については、例えば、特定された反応経路（図６（ａ）～（ｃ）に示す）等に基づいて、適宜設定することができる。

【0090】

工程Ｓ５８において、終了時間が経過したと判断された場合（工程Ｓ５８で、「Ｙ」）、分子動力学計算工程Ｓ５の一連の処理が終了し、次の工程Ｓ６（図３に示す）が実施される。一方、工程Ｓ５８において、終了時間が経過していないと判断された場合（工程Ｓ５８で、「Ｎ」）、単位時間ステップを一つ進めて（工程Ｓ６０）、工程Ｓ５２～工程Ｓ５８が再度実施される。したがって、本実施形態の分子動力学計算では、終了時間が経過するまで、分子間の化学反応を計算することができる。分子動力学計算での解析結果は、図１に示したコンピュータ１（解析結果記憶部１７ｆ）に入力される。

【0091】

このように、本実施形態の解析方法では、密度汎関数強束縛法に基づいて特定された反応経路（図６（ａ）～（ｃ）に示す）を、分子動力学計算に反映させて、分子間の化学反応を解析することができる。したがって、本実施形態の解析方法は、未知の反応経路を有する大規模な分子間の化学反応を、短時間かつ高い精度で計算することが可能となる。さらに、本実施形態の解析方法では、反応距離（図示しない反応距離Ｒ１～Ｒ３）と、反応確率との双方を考慮することができるため、実際の分子間の化学反応を、精度良く解析することができる。

【0092】

本実施形態の分子動力学計算工程Ｓ５では、図７に示されるように、化学反応前の分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）を対象として、分子動力学計算が行われたが、このような態様に限定されない。例えば、中間生成物の分子を対象として、分子動力学計算が行われてもよい。このような分子動力学計算の一例としては、全ての八硫黄Ｓ₈（図６（ａ）に示す）は、中間生成物である二硫黄Ｓ₂（図６（ｂ）に示す）に開裂すると仮定して、図７に示した八硫黄モデル４に代えて、図５（ｂ）に示した二硫黄モデル８を第２セル１０に配置して、ＭＢＴＳモデル３及び二硫黄モデル８を対象とする分子動力学計算が行われても良い。これにより、図８（ａ）及び（ｂ）に示したように、八硫黄モデル４（図８（ａ）に示す）が開裂した二硫黄モデル８を計算する必要がないため、計算時間を短縮することができる。

【0093】

次に、本実施形態の解析方法では、分子間の化学反応が、良好か否かが判断される（工程Ｓ６）。化学反応の良否については、図１に示した解析結果記憶部１７ｆに記憶されている分子動力学計算工程Ｓ５での解析結果に基づいて、適宜判断することができる。本実施形態では、生成されたＢｔ－Ｓ4－Ｂｔモデル２０（図８（ｃ）に示す）の割合が多いほど、分子間の化学反応が良好であると判断している。化学反応の良否については、コンピュータ１が判断してもよいし、オペレータが判断してもよい。

【0094】

工程Ｓ６において、分子間の化学反応が良好であると判断された場合（工程Ｓ６で、「Ｙ」）、図６（ａ）～（ｃ）に示した分子（本例では、ＭＢＴＳ及び八硫黄Ｓ₈）間の化学反応による生成物（Ｂｔ－Ｓ4－Ｂｔ）を含有する製品（例えば、ゴム製品）が製造される（工程Ｓ７）。一方、工程Ｓ６において、分子間の化学反応が良好でないと判断された場合（工程Ｓ６で、「Ｎ」）、化学反応させる分子を変更して（工程Ｓ８）、反応経路特定工程Ｓ１～工程Ｓ６が再度実施される。これにより、本実施形態では、分子間の化学反応が良好な生成物を含有する製品を確実に製造することができる。

【0095】

これまでの実施形態の解析方法では、図９に示した分子動力学計算工程Ｓ５において、反応距離（図示しない反応距離Ｒ１～Ｒ３）と、反応確率との双方を考慮して、分子間の化学反応が解析されたが、このような態様に限定されない。例えば、反応確率を考慮せずに、反応距離のみを考慮して、分子間の化学反応が解析されてもよい。このような解析方法では、分子間の化学反応を、短時間で計算することができる。

【0096】

これまでの実施形態の解析方法では、反応経路特定工程Ｓ１において、図６（ａ）～（ｃ）に示した一つの反応経路のみが取得されたが、このような態様に限定されない。反応経路特定工程Ｓ１では、例えば、複数の反応経路を取得して、それらの反応経路から相対的に反応しやすい反応経路（即ち、活性化エネルギーが相対的に小さい反応経路）が選択されてもよい。図１０は、本発明の他の実施形態の反応経路特定工程Ｓ１の処理手順の一例を示すフローチャートである。この実施形態において、これまでの実施形態と同一の構成については、同一の符号を付し、説明を省略することがある。

【0097】

この実施形態の反応経路特定工程Ｓ１は、複数の反応経路が取得される（工程Ｓ１１）。この実施形態の工程Ｓ１１では、これまでの実施形態の反応経路特定工程Ｓ１と同様に、図５（ａ）に示した分子モデル２（ＭＢＴＳモデル３及び八硫黄モデル４）、及び、計算条件に基づいて、密度汎関数強束縛法に基づくＱＭ／ＭＤ計算が行われる。この工程Ｓ１１では、複数の反応経路が取得されるまで、ＱＭ／ＭＤ計算が行われる。取得される反応経路の個数については、解析対象の分子の構造に基づいて、適宜設定される。取得された複数の反応経路は、コンピュータ１に記憶される。

【0098】

次に、この実施形態の反応経路特定工程Ｓ１は、取得された複数の反応経路（一つの反応経路を図６（ａ）～（ｃ）に示す）について、活性化エネルギーがそれぞれ計算される（工程Ｓ１２）。工程Ｓ１２では、工程Ｓ１１で取得された反応経路について、化学結合が開裂した分子間の活性化エネルギー、及び、新たな化学結合が生じた分子間の活性化エネルギーがそれぞれ計算される。図６（ａ）～（ｃ）に示した反応経路では、ＭＢＴＳを構成する一対のＢｔ－Ｓ、Ｂｔ－Ｓ（図６（ａ）において二点鎖線で示す）間の活性化エネルギー、八硫黄Ｓ₈を構成する一対の二硫黄Ｓ₂－Ｓ₂（図６（ａ）において二点鎖線で示す）間の活性化エネルギー、及び、図６（ｂ）に示した隣接するＢｔ－Ｓと二硫黄Ｓ₂との間の活性化エネルギーがそれぞれ計算される。そして、これらの活性化エネルギーの合計値を計算することで、図６（ａ）～（ｃ）に示した反応経路に必要な活性化エネルギーを求めることができる。

【0099】

活性化エネルギーは、ＤＦＴ法に基づく量子化学計算に基づいて計算することができる。活性化エネルギーの計算には、例えば、図５（ａ）及び（ｂ）に示したＭＢＴＳモデル３、八硫黄モデル４、Ｂｔ－Ｓモデル７、及び、二硫黄モデル８を用いることができる。活性化エネルギーの計算方法については、前述のとおりである。活性化エネルギーは、図１に示したコンピュータ１に記憶される。

【0100】

次に、この実施形態の反応経路特定工程Ｓ１は、複数の反応経路のうち、活性化エネルギーが相対的に小さい反応経路を選択する（工程Ｓ１３）。一般に、活性化エネルギーが相対的に小さいほど、反応経路に基づく化学反応が生じやすい傾向がある。したがって、この実施形態の反応経路特定工程Ｓ１では、工程Ｓ１３において、活性化エネルギーが相対的に小さい反応経路が、前述の分子動力学計算で反映させる反応経路として選択される。これにより、この実施形態の解析方法では、実際に生じやすい化学反応を、分子動力学計算で確実に再現することができる。したがって、この実施形態の解析方法は、未知の反応経路を有する大規模な分子間の化学反応を、より高い精度で計算することが可能となる。

【0101】

本実施形態では、１つの反応経路のみが選択されているが、２つ以上の反応経路が選択されてもよい。また、取得された全ての反応経路について、活性化エネルギーが高く計算された場合には、化学反応が生じやすい反応経路が存在しないため、化学反応させる分子を変更して、反応経路特定工程Ｓ１が再度実施されてもよい。これにより、この実施形態の解析方法では、実際に生じやすい化学反応を、分子動力学計算で確実に再現することができる。

【0102】

以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。

【実施例】

【0103】

分子間の化学反応が解析された（実施例及び比較例）。実施例では、先ず、図３に示した処理手順に基づいて、密度汎関数強束縛法に基づいて、分子間の反応経路が特定された。反応経路の特定には、合計１５個の分子モデルが用いられた。また、実施例では、１つの反応経路のみが特定された。

【0104】

次に、実施例では、特定された反応経路が分子動力学計算中で反映されるように、分子動力学計算用のパラメータが設定された。そして、実施例では、図９に示した処理手順に基づいて、分子モデルを対象として、パラメータを用いて分子動力学計算が行われた。この分子動力学計算では、合計２００個の分子モデルがそれぞれ用いられた。

【0105】

一方、比較例では、密度汎関数強束縛法に基づいて、分子間の化学反応が解析された。この比較例では、実施例の分子動力学計算で用いられた個数分の分子モデルが用いられた。共通仕様は、次のとおりである。

【0106】

化学反応が解析される分子：
ＭＢＴＳ（ジ－２－ベンゾチアゾリルジスルフィド）
八硫黄Ｓ₈
ＴＭＥ（テトラメチルエチレン）

【0107】

テストの結果、実施例の計算時間は、５．８時間であった。一方、比較例の計算時間を見積もると０．４年であり、大規模な化学反応を現実的な時間で解析できなかった。また、実施例の分子動力学計算では、密度汎関数強束縛法に基づいて特定された反応経路が反映されているため、比較例と同様の精度で解析することができた。したがって、実施例は、未知の反応経路を有する大規模な分子間の化学反応を解析することができた。

【符号の説明】

【0108】

Ｓ１ 分子間の反応経路を特定する工程
Ｓ２ 複数の分子モデルを入力する工程
Ｓ４ 分子動力学計算用のパラメータを設定する工程
Ｓ５ パラメータを用いて分子動力学計算を行う工程

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】