JP2018194990A - 原子間の結合解離ポテンシャルの計算方法及び高分子材料のシミュレーション方法 - Google Patents

Abstract

【課題】 原子間の結合及び解離を精度よく表現できるポテンシャルを定義する。【解決手段】 コンピュータを用いて、分子を構成する原子間の結合解離ポテンシャルを計算するための方法である。この計算方法は、原子をモデル化した粒子モデルを互いに結合した分子モデルを、コンピュータに定義する工程Ｓ１と、コンピュータが、分子モデルを用いた量子力学計算によって波動関数を計算する計算工程Ｓ４と、コンピュータが、波動関数が安定しているか否かを判定する評価工程Ｓ５とを含む。波動関数が安定していると判定した場合、波動関数に基づいて、粒子モデル間の結合解離ポテンシャルを求める工程Ｓ６が実施される。波動関数が安定していないと判定した場合、波動関数のエネルギーが低くなるように、波動関数を再計算する再計算工程Ｓ７が実施される。【選択図】図２

Description

本発明は、原子間の結合解離ポテンシャルを計算する方法、及び、高分子材料のシミュレーション方法に関する。

近年、ゴム等の高分子材料の開発のために、高分子材料の性質を、コンピュータを用いて評価するためのシミュレーション方法が種々提案されている（例えば、下記特許文献１参照）。

この種のシミュレーション方法では、先ず、コンピュータ上に定義された空間（セル）に、複数の粒子モデルを含む分子鎖モデルが配置された高分子材料モデルが定義される。この粒子モデルは、高分子鎖の原子に対応している。また、高分子材料モデルの粒子モデル間には、粒子モデル間の距離の関数であるポテンシャルが定義される。このような高分子材料モデルは、上記空間内で分子動力学に基づく緩和計算が行われる。その後、例えば、高分子材料モデルを用いた変形計算が実施され、高分子材料の性能等が解析される。

特開２０１３−１９５２２０号公報

実際の高分子材料では、例えば大きな力が作用すると、分子鎖が切断することが知られている。しかしながら、上記シミュレーション方法では、分子鎖モデルを切断する方法が未だ十分に確立されていない。このような分子鎖モデルの切断を計算するには、粒子モデル間に、高分子鎖の原子間の結合及び解離を精度よく表現しうるポテンシャルを定義することが重要である。

本発明は、以上のような実状に鑑み案出されたもので、原子間の結合及び解離を精度よく表現できるポテンシャルを定義しうる結合解離ポテンシャルの計算方法及び高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、分子を構成する原子間の結合解離ポテンシャルを計算するための方法であって、前記原子をモデル化した粒子モデルを互いに結合した分子モデルを、前記コンピュータに定義する工程と、前記コンピュータが、前記分子モデルを用いた量子力学計算によって波動関数を計算する計算工程と、前記コンピュータが、前記波動関数が安定しているか否かを判定する評価工程とを含み、前記波動関数が安定していると判定した場合、前記波動関数に基づいて、前記粒子モデル間の結合解離ポテンシャルを求める工程が実施され、前記波動関数が安定していないと判定した場合、前記波動関数のエネルギーが低くなるように、前記波動関数を再計算する再計算工程が実施されることを特徴とする。

本発明に係る前記原子間の結合解離ポテンシャルの計算方法において、前記量子力学計算は、密度汎関数法に基づいて行われるのが望ましい。

本発明に係る前記原子間の結合解離ポテンシャルの計算方法において、前記計算工程は、前記波動関数のエネルギーを計算する工程を含み、前記評価工程は、前記エネルギーを前記波動関数で２次変分した固有値に、負の固有値が含まれない場合に、前記波動関数が安定していると判定するのが望ましい。

本発明に係る前記原子間の結合解離ポテンシャルの計算方法において、前記再計算工程は、前記負の固有値が無くなるように、前記波動関数を再計算するのが望ましい。

本発明は、請求項１乃至４のいずれかに記載の計算方法で求められた前記結合解離ポテンシャルを用いて、高分子材料の分子鎖の切断を解析するためのシミュレーション方法であって、前記分子鎖に基づいて、複数の前記粒子モデルと、前記粒子モデル間を結合するボンドモデルとを含む分子鎖モデルを、前記コンピュータに定義する工程、前記コンピュータに定義された空間に、少なくとも一つの前記分子鎖モデルを配置した高分子材料モデルを定義する工程、及び前記コンピュータが、予め定められた条件に基づいて、前記分子鎖モデルの切断を計算する切断計算工程を含み、前記切断計算工程は、前記高分子材料モデルの前記粒子モデル間に、前記結合解離ポテンシャルに近似する結合ポテンシャルを定義する工程を含むことを特徴とすることを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記結合ポテンシャルは、前記結合解離ポテンシャルのエネルギーの下限値から、前記エネルギーの上限値の５０％〜８０％までの範囲において、前記結合解離ポテンシャルに近似する関数で定義されるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記切断計算工程は、前記分子力学計算工程後の前記分子鎖モデルについて、前記ボンドモデルを介して隣り合う前記粒子モデル間の距離が、予め定められた第１距離よりも大きい離反粒子モデル部が存在するか否かを判断する工程をさらに含み、前記離反粒子モデル部が存在すると判断された場合、前記離反粒子モデル部を含む非全体領域である第１領域を対象に量子力学計算を行う量子力学計算工程と、前記量子力学計算工程後の前記第１領域について、前記粒子モデル間の距離が、予め定められた第２距離よりも大きいときに、それらの間を切断する工程とを含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１距離は、前記第２距離よりも小であるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２距離は、前記結合解離ポテンシャルについて、そのエネルギーの増分が小さくなる変曲点での前記粒子モデル間の距離に基づいて定義されるのが望ましい。

本願の第１の発明の原子間の結合解離ポテンシャルの計算方法は、コンピュータが、分子モデルを用いた量子力学計算によって波動関数を計算する計算工程を含んでいる。

発明者らは、鋭意研究を重ねた結果、計算工程において、安定しない波動関数が計算される場合があり、安定しない波動関数から結合解離ポテンシャルを求めると、原子間の結合及び解離を精度良く表現できないことを知見した。さらに、発明者らは、量子力学計算で用いられる条件を緩和することにより、安定した波動関数を求めることができることを知見した。

これらの知見に基づいて、本願の第１の発明では、結合解離ポテンシャルを求める工程に先立ち、波動関数が安定しているか否かを判定する評価工程を含んでいる。そして、波動関数が安定していると判定した場合、波動関数に基づいて、粒子モデル間の結合解離ポテンシャルを求める工程が実施され、波動関数が安定していないと判定した場合、前記波動関数のエネルギーが低くなるように、波動関数を再計算する再計算工程が実施される。

このように、本願の第１の発明では、安定した波動関数に基づいて、結合解離ポテンシャルが求められる。しかも、波動関数が安定していない場合には、前記波動関数のエネルギーが低くなるように波動関数が再計算されることで、安定した波動関数を求めることができる。このため、本願の第１の発明では、原子間の結合及び解離を精度良く表現しうる結合解離ポテンシャルを確実に求めることができる。

本願の第２の発明の高分子材料のシミュレーション方法は、コンピュータが、予め定められた条件に基づいて、分子鎖モデルの切断を計算する切断計算工程を含んでいる。切断計算工程は、高分子材料モデルの粒子モデル間に、結合解離ポテンシャルに近似する結合ポテンシャルを定義する工程を含んでいる。

結合解離ポテンシャルは、本願の第１の発明の計算方法で求められるため、原子間の結合及び解離を精度良く表現しうる。このため、本願の第２の発明では、粒子モデル間の切断を、実際の分子鎖の切断に近似させることができる。

本実施形態の原子間の結合解離ポテンシャルの計算方法を実行するコンピュータの一例を示す斜視図である。 原子間の結合解離ポテンシャルの計算方法の処理手順の一例を示すフローチャートである。 分子モデルの一例を示す概念図である。 分子モデル定義工程の処理手順の一例を示すフローチャートである。 （ａ）〜（ｃ）は、ポテンシャルを説明する図である。 図３の分子モデルの側面図である。 結合解離ポテンシャルの一例を示すグラフである。 ポリイソプレンの構造式である。 高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。 分子鎖モデルの一例を示す概念図である。 分子鎖モデル設定工程の処理手順の一例を示すフローチャートである。 高分子材料モデルの一例を示す概念図である。 切断計算工程の処理手順の一例を示すフローチャートである。 高分子材料モデルの伸長シミュレーションの一例を説明する側面図である。 結合ポテンシャル定義工程の処理手順の一例を示すフローチャートである。 結合解離ポテンシャルを求めるのに使用した分子モデルの一例を示す図である。 結合解離ポテンシャル及び結合ポテンシャルを示すグラフである。 量子力学計算工程の処理手順の一例を示すフローチャートである。 第１領域を示す図である。 実施例及び比較例の結合解離ポテンシャルのグラフである。 実施例の粒子モデル間の切断個数と、高分子材料モデルの伸長率との関係を示したグラフである。 実施例の高分子材料の応力と、高分子材料モデルの伸長率との関係を示したグラフである。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の原子間の結合解離ポテンシャルの計算方法（以下、単に「計算方法」ということがある）は、コンピュータを用いて、分子を構成する原子間の結合及び解離を表現するポテンシャルを求めるためのものである。

図１は、本実施形態の計算方法を実行するコンピュータの一例を示す斜視図である。コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成されている。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２が設けられている。また、記憶装置には、本実施形態の計算方法、及び、後述のシミュレーション方法を実行するための処理手順（プログラム）が予め記憶されている。

本実施形態の計算方法では、メタンチオール（ＣＨ₃ＳＨ）の炭素原子と硫黄原子との間の結合解離ポテンシャルが計算される態様が例示される。なお、メタンチオールに限定されるわけではなく、様々な分子を構成している原子間の結合解離ポテンシャルを計算することができる。図２は、本実施形態の計算方法の処理手順の一例を示すフローチャートである。

本実施形態の計算方法では、先ず、分子モデルが、コンピュータ１に定義される（分子モデル定義工程Ｓ１）。図３は、分子モデルの一例を示す概念図である。

分子モデル２は、全原子モデルとして構成されている。分子モデル２は、複数の粒子モデル３を含んでおり、これらの粒子モデル３が互いに結合されている。本実施形態の分子モデル２には、粒子モデル３間を結合するボンドモデル４が含まれている。図４は、分子モデル定義工程Ｓ１の処理手順の一例を示すフローチャートである。

本実施形態の分子モデル定義工程Ｓ１では、先ず、図３に示されるように、コンピュータ１に、粒子モデル３が設定される（工程Ｓ１１）。粒子モデル３は、後述する量子力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。

本実施形態の粒子モデル３は、分子（本実施形態では、メタンチオール）について、炭素原子をモデル化した炭素粒子モデル３Ｃ、水素原子をモデル化した水素粒子モデル３Ｈ、及び、硫黄原子をモデル化した硫黄粒子モデル３Ｓを含んでいる。

次に、本実施形態の分子モデル定義工程Ｓ１では、コンピュータ１に、ボンドモデル４が設定される（工程Ｓ１２）。ボンドモデル４は、粒子モデル３、３間を拘束するものである。本実施形態のボンドモデル４は、炭素粒子モデル３Ｃと硫黄粒子モデル３Ｓとの間、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとの間にそれぞれ設けられている。

次に、本実施形態の分子モデル定義工程Ｓ１では、隣り合う粒子モデル３、３間に、相互作用（斥力及び引力を含む）が生じさせるポテンシャルが定義される（工程Ｓ１３）。図５（ａ）〜（ｃ）は、ポテンシャルを説明する図である。

図５（ａ）〜（ｂ）に示されるように、分子モデル２には、各粒子モデル３、３間の結合長さである結合距離（結合長）ｒ、及び、ボンドモデル４を介して連続する３つの粒子モデル３がなす角度である結合角θが定義されている。さらに、図５（ｃ）に示されるように、分子モデル２には、ボンドモデル４を介して連続する４つの粒子モデル３において、隣り合う３つの粒子モデル３が作る一方の平面５Ａと他方の平面５Ｂとのなす角度ある二面角φが定義される。結合距離ｒ、結合角θ及び二面角φは、分子モデル２に作用する外力又は内力によって変化する。

結合距離ｒ、結合角θ及び二面角φは、下記式（１）で定義される結合ポテンシャルＵ_bond（ｒ）、下記式（２）で定義される結合角ポテンシャルＵ_angle（θ）、及び、下記式（３）で定義される結合二面角ポテンシャルＵ_torsion（φ）によって設定される。
ここで、各定数及び変数は、次のとおりである。
ｒ：結合距離
ｒ₀：平衡距離（平衡長）
ｋ₁、ｋ₂：ばね定数
θ：結合角
θ₀：平衡角度
ｋ₃：二面角ポテンシャルの強度
Ｎ−1：二面角ポテンシャル多項式の次数
φ：二面角
Ａ_n：二面角定数
なお、結合距離ｒ及び平衡距離ｒ₀は、各粒子モデル３の中心３ｃ（図６に示す）間の距離として定義される。

各結合ポテンシャルＵ_bond（ｒ）、結合角ポテンシャルＵ_angle（θ）、及び、結合二面角ポテンシャルＵ_torsion（φ）の各定数については、適宜設定することができる。なお、これらのポテンシャルは、論文（J. Phys. Chem. 94, 8897（1990））に基づいて、分子鎖Ｍｃの構造に応じて設定されるのが望ましい。

これらのポテンシャルが、隣り合う粒子モデル３、３間に定義されることにより、図３に示した分子モデル２が定義される。分子モデル２は、コンピュータに記憶される。

次に、本実施形態の計算方法では、粒子モデル３、３間に、後述の量子力学計算で使用する理論及び基底関数が定義される（工程Ｓ２）。本実施形態で定義される理論としては、例えば、密度汎関数法（汎関数：Ｂ３ＬＹＰ）が定義される。また、基底関数としては、例えば、基底関数（６−３１Ｇ（ｄ））が定義される。これらの量子力学計算に使用する理論及び基底関数は、Ｇａｕｓｓｉａｎ社製の量子化学計算プログラムＧａｕｓｓｉａｎ０３、又は、Ｇａｕｓｓｉａｎ０９を用いることによって定義することができる。これらの理論及び基底関数は、コンピュータ１に記憶される。

次に、本実施形態の計算方法では、量子力学計算で用いられる条件（制約）が定義される（工程Ｓ３）。条件は、分子軌道に関するものである。条件の一例としては、電荷、スピン多重度、空間軌道、及び、対称性に関する制限等が含まれる。これらの条件は、量子力学計算に用いられるパラメータであり、コンピュータ１に記憶される。

次に、本実施形態の計算方法では、コンピュータ１が、分子モデル２を用いた量子力学計算によって波動関数を計算する（計算工程Ｓ４）。計算工程Ｓ４では、上記工程Ｓ３で予め定められた条件（制約）に基づいて、量子力学計算が行われる。

量子力学計算では、分子モデル２が、原子核と電子とに基づく相互作用に従うものとして、シミュレーションの単位時間ステップ毎に、粒子モデル３の電子状態が追跡される。従って、量子力学計算は、分子力学計算に比べて計算精度が高い。このような量子力学計算は、上記した量子化学計算プログラムを用いて処理することができる。

図６は、図３の分子モデルの側面図である。図６に示されるように、分子モデル２は、上記理論及び基底関数により、原子核と電子とに基づく相互作用のバランスによって、粒子モデル３、３の平衡距離ｒ₀が決定されている。本実施形態の計算工程Ｓ４では、例えば、粒子モデル３、３が当接した状態から、粒子モデル３、３間の距離ｒが平衡距離ｒ₀よりも大（例えば、４．５Å）になるまでの間、粒子モデル３、３を単位時間ステップ毎に０．１Åずつ離間させる量子力学計算が実施され、波動関数のエネルギーが計算される。

本実施形態の計算工程Ｓ４では、スピン多重度が一重項（singlet）に設定されている。これにより、計算工程Ｓ４では、不対電子を有さない波動関数のエネルギーが計算される。波動関数のエネルギーは、コンピュータ１に入力される。

次に、図２に示されるように、本実施形態の計算方法では、コンピュータ１が、波動関数が安定しているか否かを判定する（評価工程Ｓ５）。波動関数は、そのエネルギーが低くなるほど安定している。一般に、現在の波動関数のエネルギーをその波動関数で２次変分した固有値に、負の固有値が含まれる場合、よりエネルギーの低い安定した波動関数が存在することを意味している。このような場合、現在の波動関数は、不安定であると判定することができる。このような波動関数の不安定性は、とりわけ単一の電子配置しか考慮しない密度汎関数法において、原子間距離が非常に長い分子を計算する場合に発生しやすい。本実施形態の評価工程Ｓ５は、波動関数のエネルギーをその波動関数で２次変分した固有値に、負の固有値が含まれない（即ち、全て正の固有値である）場合に、波動関数が安定していると判定している。これにより、本実施形態の評価工程Ｓ５では、安定した波動関数を確実に判定することができる。

評価工程Ｓ５において、コンピュータ１が、波動関数が安定していると判定した場合（評価工程Ｓ５で、「Ｙ」）、安定した波動関数に基づいて、粒子モデル３、３（図６に示す）間の結合解離ポテンシャルが求められる（工程Ｓ６）。他方、コンピュータ１が、波動関数が安定していないと判定した場合（評価工程Ｓ５で、「Ｎ」）、現在の波動関数のエネルギーが低くなるように、波動関数を再計算する再計算工程Ｓ７が実施される。

結合解離ポテンシャルを求める工程Ｓ６では、波動関数に基づいて、単位時間ステップの結合解離ポテンシャルが求められている。即ち、結合解離ポテンシャルは、粒子モデル３、３（図６に示す）間の距離ｒと、波動関数のエネルギーとの関係に基づいて求められる。結合解離ポテンシャルは、コンピュータ１に入力される。

本実施形態の再計算工程Ｓ７では、現在の波動関数のエネルギーをその波動関数で２次変分した固有値に、負の固有値が無くなる（即ち、全て正の固有値になる）ように、電子の配置又は軌道の順序の入れ替えが行われる。このような電子の配置及び軌道の順序の入れ替えは、上記した量子化学計算プログラムにおいて、オプション「Stable=Opt」を追加することにより計算することができる。これにより、再計算工程Ｓ７では、よりエネルギーの低い波動関数を計算することができる。

次に、図２に示されるように、本実施形態の計算方法では、コンピュータ１が、粒子モデル３、３（図６に示す）が、予め定められた距離（例えば、４．５Å）まで離間したか否かが判断する（工程Ｓ８）。工程Ｓ８において、粒子モデル３、３が予め定められた距離まで離間したと判断された場合（工程Ｓ８で、「Ｙ」）、本実施形態の計算方法の一連の処理が終了する。

他方、工程Ｓ８において、粒子モデル３、３（図６に示す）が予め定められた距離まで離間していないと判断された場合（工程Ｓ６で、「Ｎ」）、単位時間ステップを一つ進展させて（工程Ｓ９）、計算工程Ｓ４〜工程Ｓ８が再度実施される。これにより、粒子モデル３、３が当接した状態から、粒子モデル３、３間の距離ｒが平衡距離ｒ₀よりも大（例えば、４．５Å）になるまでの結合解離ポテンシャルを求めることができる。なお、計算工程Ｓ４では、再計算工程Ｓ７で緩和された条件を、元の条件に再定義した後に、波動関数が計算される。

図７は、結合解離ポテンシャルＰｕの一例を示すグラフである。図７のグラフにおいて、結合解離ポテンシャルＰｕは、図６に示した粒子モデル３、３間の距離ｒが平衡距離ｒ₀（本実施形態では、２．１）よりも小さくなるほど、波動関数のポテンシャルエネルギー（粒子モデル３、３間の斥力）が大きくなる。また、結合解離ポテンシャルＰｕは、距離ｒが２．１になるときにポテンシャルエネルギーが最小となり、粒子モデル３、３間に斥力や引力は作用しない。さらに、結合解離ポテンシャルＰｕは、距離ｒが平衡距離ｒ₀よりも大きくなるほど、波動関数のポテンシャルエネルギー（粒子モデル３、３間の引力）が大きくなる。このように、結合解離ポテンシャルＰｕは、距離ｒに応じて、斥力及び引力を定義することができる。

結合解離ポテンシャルＰｕは、粒子モデル３、３（図６に示す）間の距離ｒが平衡距離ｒ₀よりも大きい領域において、波動関数のエネルギーの増分が小さくなる変曲点８を有している。変曲点８は、変曲点８よりも距離ｒが大きいと予想される領域、及び、変曲点８よりも距離ｒが小さいと予想される領域において、結合解離ポテンシャルＰｕに近似する一対の直線９ａ、９ｂの交点で特定することができる。

変曲点８よりも距離ｒが小さい領域Ｔｂでは、変曲点８よりも粒子モデル３、３（図６に示す）間の距離ｒが大きい領域Ｔａに比べて、波動関数のエネルギーの増分が大きい。これは、変曲点８よりも小さい領域Ｔｂにおいて、粒子モデル３、３間の結合が強固に維持されることを示している。また、変曲点８よりも距離ｒが大きい領域Ｔａでは、粒子モデル３、３間で切断、又は、結合次数の減少が生じており、波動関数のエネルギーの増分が小さくなっている。従って、変曲点８は、粒子モデル３、３間の結合と解離との境界と考えることができる。従って、本実施形態の結合解離ポテンシャルＰｕは、原子間の結合及び解離を精度良く表現することができる。

このように、本実施形態の計算方法は、安定した波動関数のみを用いて、結合解離ポテンシャルＰｕを求めることができる。しかも、波動関数が安定していない場合には、波動関数が安定するまで再計算される。このため、本実施形態の計算方法は、原子間の結合及び解離を精度良く表現しうる結合解離ポテンシャルＰｕを確実に求めることができる。

また、本実施形態では、密度汎関数法（ＤＦＴ）に基づいて量子力学計算が行われているため、他の理論（例えば、ＭＰ２、ＣＣＳＤＴ、ＣＡＳＳＣＦ及びＭＲＣＩ等）に比べて、計算時間を短縮することができる。従って、本実施形態の計算方法では、原子間の結合及び解離を精度良く表現しうる結合解離ポテンシャルＰｕを、短時間で求めることができる。

次に、本実施形態の計算方法で求められた結合解離ポテンシャルを用いて、高分子材料の分子鎖の切断を解析するための高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）について説明する。

本実施形態のシミュレーション方法は、コンピュータ１（図１に示す）を用いて、高分子材料の分子鎖の切断を解析するための方法である。分子鎖の切断を計算するには、高分子鎖の原子間の結合及び解離を精度よく表現しうるポテンシャルを定義することが重要である。本実施形態では、上記した計算方法で求められた結合解離ポテンシャルＰｕが用いられるため、分子鎖の切断を、精度よく解析するのに役立つ。

高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態の高分子材料は、cis-１，４ポリイソプレン（以下、単に「ポリイソプレン」ということがある。）である。図８は、ポリイソプレンの構造式である。

ポリイソプレンを構成する分子鎖１１Ａは、メチン基等（例えば、−ＣＨ＝、＞Ｃ＝）、メチレン基（−ＣＨ_２−）、及び、メチル基（−ＣＨ_３）によって構成されるイソプレンのモノマー（イソプレン分子）１２が、重合度ｎで連結されて構成されている。なお、高分子材料１１には、ポリイソプレン以外の高分子材料が用いられてもよい。

図９は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、図８に示した分子鎖１１Ａに基づいて、分子鎖モデルがコンピュータ１に定義される（分子鎖モデル設定工程Ｓ２１）。図１０は、分子鎖モデル１３の一例を示す概念図である。

分子鎖モデル１３は、全原子モデルとして構成されている。本実施形態の分子鎖モデル１３は、一つの硫黄原子が架橋されたポリイソプレンとして設定されている。分子鎖モデル１３は、複数の粒子モデル３と、粒子モデル３、３間を結合するボンドモデル４とを含んで構成されている。

粒子モデル３及びボンドモデル４は、図８に示した分子鎖１１Ａのモノマー１２をなす単位構造に基づいて連結される。これにより、モノマーモデル１４が設定される。このモノマーモデル１４は、分子量（重合度）Ｍnに基づいて連結される。これにより、分子鎖モデル１３が設定される。図１１は、分子鎖モデル設定工程Ｓ２１の処理手順の一例を示すフローチャートである。

本実施形態の分子鎖モデル設定工程Ｓ２１では、先ず、図１０に示されるように、粒子モデル３が設定される（工程Ｓ２１１）。粒子モデル３は、後述する分子動力学計算、分子力学計算、及び、量子力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル３は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。

本実施形態の粒子モデル３は、炭素粒子モデル３Ｃ、水素粒子モデル３Ｈ及び硫黄粒子モデル３Ｓを含んでいる。これらの粒子モデル３は、コンピュータ１に記憶される。

次に、分子鎖モデル設定工程Ｓ２１では、ボンドモデル４が設定される（工程Ｓ２１２）。ボンドモデル４は、粒子モデル３、３間を拘束するものである。本実施形態のボンドモデル４は、主鎖４ａと側鎖４ｂとを含んでいる。主鎖４ａは、炭素粒子モデル３Ｃ、３Ｃを連結するものである。側鎖４ｂは、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間、及び、炭素粒子モデル３Ｃと硫黄粒子モデル３Ｓとの間を連結するものである。これらのボンドモデル４は、コンピュータ１に記憶される。

次に、分子鎖モデル設定工程Ｓ２１では、ボンドモデル４を介して隣り合う粒子モデル３、３間には、相互作用（斥力及び引力を含む）が生じさせるポテンシャルが定義される（ポテンシャル定義工程Ｓ２１３）。図１０に示されるように、ポテンシャルには、ボンドモデル４を介して隣り合う粒子モデル３、３間に定義される第１ポテンシャルＰ１と、ボンドモデル４を介さずに隣り合う粒子モデル３、３間に定義される第２ポテンシャルＰ２とが定義される。第１ポテンシャルＰ１は、図３及び図５に示した分子モデル２と同様に、上記式（１）〜（３）に基づいて設定される。

第２ポテンシャルＰ２は、下記式（４）で定義される ＬＪポテンシャルＵ_LJ（ｒ_ij）である。

ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒ_ij：粒子モデル間の距離
ε：粒子モデル間に定義されるＬＪポテンシャルの強度
σ：粒子モデルの直径に相当
なお、距離ｒ_ijは、各粒子モデル３、３の中心（図示省略）間の距離として定義される。

第２ポテンシャルＰ２は、粒子モデル３、３間の距離ｒ_ijがｒ_ｃよりも小さくなるほど、粒子モデル３、３間に作用する斥力が大きくなる。また、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijがｒ_ｃになるときに最小となり、粒子モデル３、３間に斥力や引力は働かない。さらに、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijがｒ_ｃよりも大になると、粒子モデル３、３間に作用する引力が働く。このように、第２ポテンシャルＰ２は、粒子モデル間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

本実施形態では、硫黄粒子モデル３Ｓ、３Ｓ間、炭素粒子モデル３Ｃ、３Ｃ間、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間、硫黄粒子モデル３Ｓと炭素粒子モデル３Ｃとの間、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとの間に、それぞれ異なる第２ポテンシャルＰ２が設定されている。各第２ポテンシャルＰ２は、上記式（４）の定数がそれぞれ異なっている。なお、各定数は、例えば、論文（J. Phys. Chem. 94, 8897（1990））に基づいて、適宜設定することができる。

次に、本実施形態のシミュレーション方法では、コンピュータ１に、高分子材料モデルが定義される（工程Ｓ２２）。図１２は、高分子材料モデルの一例を示す概念図である。高分子材料モデル１５は、コンピュータ１上に定義された空間１６に、少なくとも一つの分子鎖モデル１３が配置されることによって定義される。

本実施形態の空間１６は、互いに向き合う三対の平面１６ａ、１６ｂを有する直方体として定義されている。各平面１６ａ、１６ｂには、周期境界条件が定義されている。これにより、空間１６では、例えば、一方の平面１６ａから出て行った分子鎖モデル１３の一部が、反対側の平面１６ｂから入ってくるように計算することができる。従って、一方の平面１６ａと、反対側の平面１６ｂとが連続している（繋がっている）ものとして取り扱うことができる。

空間１６の一辺の長さＬ１、Ｌ２、Ｌ３は、適宜設定することができる。本実施形態の長さＬ１は、分子鎖モデル１３の慣性半径（図示省略）の２倍以上が望ましい。慣性半径は、分子鎖モデル１３の拡がりを示す量である。このような空間１６では、分子動力学計算において、周期境界条件による自己のイメージとの衝突が起こり難くなり、分子鎖モデル１３の空間的拡がりを適切に計算することができる。さらに、空間１６の大きさは、例えば１気圧で安定な体積に設定される。このような空間１６は、高分子材料の少なくとも一部の体積を定義することができる。また、空間１６に配置される分子鎖モデル１３の個数については、適宜設定することができる。本実施形態の分子鎖モデル１３の個数としては、例えば、１個以上、好ましくは３０個以上に設定されるのが望ましい。

本実施形態では、分子鎖モデル１３が配置された空間１６を用いて、分子動力学（ Molecular Dynamics ： ＭＤ ）計算による緩和計算が実施される。これにより、高分子材料モデル１５が設定される。

本実施形態の分子動力学計算では、例えば、空間１６について所定の時間、分子鎖モデル１３が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル３の動きが、単位時間ステップ毎に追跡される。このような構造緩和の計算は、例えばソフトマテリアル総合シミュレーター（ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。

分子動力学計算では、空間１６において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。また、分子動力学計算では、分子鎖モデル１３の人為的な初期配置が十分に排除されたとみなせるまで行われる。これにより、工程Ｓ２２では、実際の高分子材料の分子運動に近似させて、分子鎖モデル１３の初期配置を精度よく緩和することができる。このような緩和計算を経て、高分子材料モデル１５が定義される。高分子材料モデル１５は、コンピュータ１に記憶される。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、予め定められた条件に基づいて、分子鎖モデル１３の切断を計算する（切断計算工程Ｓ２３）。本実施形態の切断計算工程Ｓ２３では、図１２に示した高分子材料モデル１５を、一軸方向（例えば、Ｚ軸方向）に伸長させる伸長シミュレーションが実施される。

本実施形態の伸長シミュレーションでは、分子力学計算（Molecular Mechanics ： ＭＭ）及び量子力学計算（Quantum Mechanics ： ＱＭ）を用いて、単位時間ステップ毎に、高分子材料モデル１５の伸長が計算される。そして、単位時間ステップを進展させることによって、予め定められた長さ（上限値）まで伸長した高分子材料モデル１５が計算される。この高分子材料モデル１５の伸長によって、分子鎖モデル１３は、粒子モデル３、３間が離間する。この粒子モデル３、３間の離間に伴って、分子鎖モデル１３が切断される状態が計算される。分子力学計算は、分子動力学計算と同様の上記ソフトウェアによって計算されうる。

量子力学計算では、分子鎖モデル１３が、原子核と電子とに基づく相互作用に従うものとして、シミュレーションの単位時間ステップ毎に、粒子モデル３の電子状態と動きとが追跡される。従って、量子力学計算は、分子力学計算に比べて計算精度が高い。このような量子力学計算は、上記した量子化学計算プログラムを用いて処理することができる。

本実施形態では、高分子材料モデル１５の全体領域を対象に、分子力学計算を用いて、伸長シミュレーションが実施される。そして、分子鎖モデル１３のうち粒子モデル３、３間で切断しそうな部分のみを対象に、量子力学計算を用いて、分子鎖モデル１３の切断が詳細に計算される。従って、本実施形態の切断計算工程Ｓ２３では、例えば、高分子材料モデル１５の全体領域を対象に量子力学計算を実施する場合に比べて、計算時間を大幅に短縮することができる。図１３は、本実施形態の切断計算工程Ｓ２３の処理手順の一例を示すフローチャートである。

本実施形態の切断計算工程Ｓ２３では、先ず、高分子材料モデル１５の伸長シミュレーションを実施するための初期条件が設定される（工程Ｓ３１）。図１４は、高分子材料モデル１５の伸長シミュレーションの一例を説明する側面図である。初期条件としては、高分子材料モデル１５のポアソン比に基づいて、各単位時間ステップでの高分子材料モデル１５の長さＬ１、Ｌ２、Ｌ３が定義される。これにより、切断計算工程Ｓ２３では、高分子材料モデル１５のポアソン比に基づいて、高分子材料モデルを一軸方向（例えば、Ｚ軸方向）に伸長させるシミュレーションを実施することができる。

高分子材料モデル１５の一軸方向の伸長速度Ｖａについては、解析対象の高分子材料や、シミュレーション条件に応じて適宜設定することができる。本実施形態の伸長速度Ｖａは、例えば、５０〜２００ｍ／秒に設定されている。これにより、本実施形態の切断計算工程Ｓ２３では、高分子材料が高速に引き伸ばされた状態を計算することができる。これらの初期条件は、コンピュータ１に記憶される。

次に、本実施形態の切断計算工程Ｓ２３では、分子力学計算で用いられる結合ポテンシャルＰｎが定義される（結合ポテンシャル定義工程Ｓ３２）。結合ポテンシャルＰｎは、分子鎖モデル１３において、ボンドモデル４を介して隣り合う粒子モデル３、３間に定義されるものである。本実施形態の結合ポテンシャルＰｎとしては、上記した計算方法で求められた結合解離ポテンシャルＰｕ（図７に示す）に近似するポテンシャルが定義される。このような結合ポテンシャルＰｎは、量子力学計算でのポテンシャルに近似するため、分子力学計算の結果と量子力学計算の結果との差異を、最小限にすることができる。

結合解離ポテンシャルＰｕ（図示省略）は、結合する粒子モデル３の種類に応じてそれぞれ異なる。本実施形態では、硫黄粒子モデル３Ｓ、３Ｓ間、炭素粒子モデル３Ｃ、３Ｃ間、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間、硫黄粒子モデル３Ｓと炭素粒子モデル３Ｃとの間、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとの間に、それぞれ異なる結合解離ポテンシャルＰｕが求められる。図１５は、本実施形態の結合ポテンシャル定義工程の処理手順の一例を示すフローチャートである。図１６は、結合解離ポテンシャルＰｕを求めるのに使用した分子モデル２の一例を示す図である。

本実施形態の結合ポテンシャル定義工程Ｓ３２では、先ず、分子モデル２が設定される（工程Ｓ３２１）。分子モデル２は、分子鎖モデル１３において、一対の粒子モデル３、３の全ての組み合わせが設定される。本実施形態では、硫黄粒子モデル３Ｓ、３Ｓを結合した分子モデル２Ａ、炭素粒子モデル３Ｃ、３Ｃを結合した分子モデル（図示省略）、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとを結合した分子モデル（図示省略）、硫黄粒子モデル３Ｓと炭素粒子モデル３Ｃとを連結した分子モデル（図示省略）、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとを連結した分子モデル（図示省略）が設定される。図１６では、硫黄粒子モデル３Ｓ、３Ｓを結合した分子モデル２Ａが代表して示されている。

本実施形態の分子モデル２は、一対の粒子モデル３、３（本実施形態では、硫黄粒子モデル３Ｓ、３Ｓ）と、ボンドモデル４とを含んでいる。なお、各粒子モデル３には、粒子モデル３でモデル化された原子の原子価（例えば、硫黄原子の場合：２）から、粒子モデル３の結合数（例えば、１）を減じた数の水素粒子モデル３Ｈ（例えば、１個）が、ボンドモデル４を介して連結される。これにより、分子モデル２は、量子力学計算において、安定した構造を維持することができる。このような方法に基づいて、上記した各分子モデル２がそれぞれ設定される。これらの分子モデル２は、コンピュータ１に記憶される。

次に、本実施形態の結合ポテンシャル定義工程Ｓ３２では、各分子モデル２の粒子モデル３、３間に、量子力学計算に使用する理論及び基底関数が定義される（工程Ｓ３２２）。理論及び基底関数は、上記した計算方法と同様のものが定義され、コンピュータ１に記憶される。

次に、本実施形態の結合ポテンシャル定義工程Ｓ３２では、上記した計算方法に基づいて、分子モデル２の各粒子モデル３、３間の結合解離ポテンシャルＰｕ（図７に示す）が求められる（工程Ｓ３２３）。工程Ｓ３２３では、上記した各分子モデル２の全てについて、量子力学計算がそれぞれ実施される。これにより、硫黄粒子モデル３Ｓ、３Ｓ間の結合解離ポテンシャルＰｕ、炭素粒子モデル３Ｃ、３Ｃ間の結合解離ポテンシャルＰｕ、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間の結合解離ポテンシャルＰｕ、硫黄粒子モデル３Ｓと炭素粒子モデル３Ｃとの間の結合解離ポテンシャルＰｕ、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとの間の結合解離ポテンシャルＰｕがそれぞれ求められる。各結合解離ポテンシャルＰｕは、コンピュータ１に記憶される。

次に、本実施形態の結合ポテンシャル定義工程Ｓ３２では、結合解離ポテンシャルＰｕに近似するポテンシャル（結合ポテンシャル）Ｐｎが求められる（工程Ｓ３２４）。図１７は、結合解離ポテンシャルＰｕ及び結合ポテンシャルＰｎを示すグラフである。本実施形態の工程Ｓ３２４は、フックの法則に基づいた調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させている。これにより、結合ポテンシャルＰｎが求められる。調和振動子型の関数は、下記式（５）で定義される。

ここで、
Ｖ：分子モデルのエネルギー
ｋ：バネ定数
ｒ：粒子モデル間の距離
ｒ₀：粒子モデル間の平衡距離

上記式（５）の調和振動子型の関数は、粒子モデル３、３間の距離ｒと、分子モデルのエネルギー（波動関数のエネルギー）との関係を示す二次関数である。本実施形態では、図７に示した変曲点８よりも距離ｒが小さい領域Ｔｂにおいて、調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させることによって、結合ポテンシャルＰｎを求めている。従って、本実施形態の結合ポテンシャルＰｎは、粒子モデル３、３間で切断、又は、結合次数の減少が生じる前の結合解離ポテンシャルＰｕに近似するものである。

調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させる方法としては、例えば、最小二乗法を用いることができる。このような結合ポテンシャルＰｎでは、複雑な結合解離ポテンシャルＰｕを、簡単な調和振動子型の関数で表すことができる。

工程Ｓ３２４では、結合解離ポテンシャルＰｕのエネルギー（分子モデル２のエネルギー）の下限値（即ち、平衡距離ｒ₀）から、エネルギーの上限値の５０％〜８０％までの範囲において、調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させるのが望ましい。これにより、結合ポテンシャルＰｎは、図７に示した変曲点８よりも距離ｒ₀が小さい領域Ｔｂにおいて、結合解離ポテンシャルＰｕに精度良く近似することができる。

なお、調和振動子型の関数に近似させる範囲において、結合解離ポテンシャルＰｕのエネルギーの上限値の８０％を超えると、図７に示した変曲点８よりも距離が大きい領域Ｔａを含めて、調和振動子型の関数を近似させるおそれがある。逆に、結合解離ポテンシャルＰｕのエネルギーの上限値の５０％未満であると、結合解離ポテンシャルＰｕに近似させる範囲が小さくなり、結合解離ポテンシャルＰｕに精度良く近似できないおそれがある。このような観点より、調和振動子型の関数に近似させる範囲において、より好ましくは、結合解離ポテンシャルＰｕのエネルギーの上限値の７０％以下であり、より好ましくは６０％以上である。

工程Ｓ３２４では、図１６に示した硫黄粒子モデル３Ｓ、３Ｓ間の結合解離ポテンシャルＰｕ、炭素粒子モデル３Ｃ、３Ｃ間の結合解離ポテンシャルＰｕ、炭素粒子モデル３Ｃと水素粒子モデル３Ｈとの間の結合解離ポテンシャルＰｕ、硫黄粒子モデル３Ｓと炭素粒子モデル３Ｃとの間の結合解離ポテンシャルＰｕ、及び、硫黄粒子モデル３Ｓと水素粒子モデル３Ｈとの間の結合解離ポテンシャルＰｕについて、結合ポテンシャルＰｎがそれぞれ求められる。

工程Ｓ３２４では、図１０に示した分子鎖モデル１３において、ボンドモデル４を介して隣り合う粒子モデル３、３間に、第１ポテンシャルＰ１に代えて、結合ポテンシャルＰｎがそれぞれ定義される。なお、ボンドモデル４を介さずに隣り合う粒子モデル３、３間に定義されている第２ポテンシャルＰ２については、維持される。

次に、本実施形態の切断計算工程Ｓ２３では、高分子材料モデル１５（図１２に示す）を伸長させる変形計算が実施される（伸長工程Ｓ３３）。本実施形態では、初期条件として設定されたポアソン比に基づいて、高分子材料モデル１５の変形が計算される。

伸長工程Ｓ３３では、現在の単位時間ステップにおいて、一軸方向に伸長した高分子材料モデル１５（図１４に示す）が計算される。この高分子材料モデル１５の伸長により、図１２に示した分子鎖モデル１３も一軸方向に引き伸ばされ、ボンドモデル４を介して隣り合う粒子モデル３、３間の距離ｒ（図１０に示す）の少なくとも一部が大きくなる。

次に、本実施形態の切断計算工程Ｓ２３では、伸長工程Ｓ３３で変形した高分子材料モデル１５に対して、分子力学計算（Molecular Mechanics ： ＭＭ）が実施される（第１分子力学計算工程Ｓ３４）。第１分子力学計算工程Ｓ３４では、高分子材料モデル１５のエネルギーＥａと、高分子材料モデル１５に作用する力Ｆａとを含む物理量が計算される。このようなエネルギーＥａと、力Ｆａとを含む物理量は、コンピュータ１に記憶される。

次に、本実施形態の切断計算工程Ｓ２３では、第１分子力学計算工程Ｓ３４後の分子鎖モデル１３について、ボンドモデル４を介して隣り合う粒子モデル３、３間の距離ｒ（図１０に示す）が、予め定められた第１距離Ｌａよりも大きい離反粒子モデル部２１が存在するか否かが判断される（工程Ｓ３５）。

工程Ｓ３５では、空間１６に配置されている全ての分子鎖モデル１３について、ボンドモデル４を介して隣り合う粒子モデル３、３間の距離ｒが計算される。なお、工程Ｓ３５では、ボンドモデル４で拘束されていない粒子モデル３、３間の距離は計算されない。そして、粒子モデル３、３間の距離ｒが、第１距離Ｌａよりも大きい一対の粒子モデル３、３及びボンドモデル４を、離反粒子モデル部２１として判断される。

本実施形態の第１距離Ｌａは、図１７に示した結合解離ポテンシャルＰｕの変曲点８での粒子モデル３、３間の距離ｒ（Ｌｎ）よりも小に設定されている。上述したように、変曲点８は、粒子モデル３、３間の結合と解離との境界点と考えることができる。このため、距離ｒが第１距離Ｌａよりも大きい離反粒子モデル部２１では、粒子モデル３、３間で切断する可能性が高いと考えることができる。従って、工程Ｓ３５では、粒子モデル３、３間で切断する可能性が高い離反粒子モデル部２１（図１９に示す）の有無を判断することができる。

また、第１距離Ｌａについては、変曲点８での距離ｒよりも小さい値であれば、適宜設定することができる。本実施形態では、図１７に示されるように、結合解離ポテンシャルＰｕと、結合ポテンシャルＰｎとの交点２２での距離ｒを、第１距離Ｌａとして求めている。これにより、量子力学計算よりも計算時間を短縮しうる分子力学計算に基づいて、離反粒子モデル部２１の存在の有無を、精度良く判断することができる。

工程Ｓ３５において、離反粒子モデル部２１が存在すると判断された場合（工程Ｓ３５で、「Ｙ」）、離反粒子モデル部２１を対象に、次の量子力学計算工程Ｓ３６が実施される。他方、離反粒子モデル部２１が存在しないと判断された場合（工程Ｓ３５で、「Ｎ」）、後述する工程Ｓ４０が実施される。

次に、本実施形態の切断計算工程Ｓ２３では、離反粒子モデル部２１を含む第１領域２３を対象に、量子力学計算が行われる（量子力学計算工程Ｓ３６）。図１８は、量子力学計算工程Ｓ３６の処理手順の一例を示すフローチャートである。

本実施形態の量子力学計算工程Ｓ３６では、先ず、離反粒子モデル部２１を含む第１領域２３が設定される（工程Ｓ３６１）。図１９は、第１領域２３を示す図である。第１領域２３は、高分子材料モデル１５の分子鎖モデル１３から、離反粒子モデル部２１のみを抜き出した非全体領域である。なお、高分子材料モデル１５に、複数の離反粒子モデル部２１が存在する場合、全ての離反粒子モデル部２１が、第１領域２３に含められる。

また、離反粒子モデル部２１の各粒子モデル３には、粒子モデル３でモデル化された原子の原子価（例えば、炭素原子の場合：４）から、粒子モデル３の結合数（例えば、１）を減じた数の水素粒子モデル３Ｈ（例えば、３個）が、ボンドモデル４を介して連結される。これにより、量子力学計算において、安定した構造を有する離反粒子モデル部２１を設定することができる。

次に、本実施形態の量子力学計算工程Ｓ３６では、離反粒子モデル部２１の粒子モデル３に、量子力学計算に使用する理論及び基底関数が定義される（工程Ｓ３６２）。工程Ｓ３６２では、離反粒子モデル部２１の粒子モデル３、３間に設定されている結合ポテンシャルＰｎを無効にして、量子力学計算に使用する理論及び基底関数が定義される。なお、量子力学計算に使用する理論及び基底関数については、上述のとおりである。

次に、本実施形態の量子力学計算工程Ｓ３６では、初期条件として設定されたポアソン比に基づいて、第１領域２３を対象に、量子力学計算が実施される（工程Ｓ３６３）。上述したように、量子力学計算では、原子核と電子とに基づく相互作用に従うものとして、シミュレーションの単位時間ステップ毎に、粒子モデル３の電子状態と動きとが追跡される。従って、工程Ｓ３６３では、分子力学計算に比べて、第１領域２３及び離反粒子モデル部２１の伸長を、高い精度で計算することができる。また、工程Ｓ３６３では、高分子材料モデル１５から抜き出された第１領域２３のみを対象に、量子力学計算が実施されている。このため、図１２に示した高分子材料モデル１５の全体領域を対象に量子力学計算が実施される場合に比べて、計算時間を大幅に短縮することができる。

工程Ｓ３６３では、第１領域２３のエネルギーＥｃと、第１領域２３に作用する力Ｆｃとを含む物理量が計算される。このようなエネルギーＥｃ及び力Ｆｃを含む物理量は、コンピュータ１に記憶される。

次に、本実施形態の切断計算工程Ｓ２３では、第１領域２３を対象に分子力学計算が行われる（第２分子力学計算工程Ｓ３７）。図１９に示されるように、第１領域２３は、伸長工程Ｓ３３で変形した高分子材料モデル１５の分子鎖モデル１３から、離反粒子モデル部２１のみを抜き出した非全体領域である。従って、第２分子力学計算工程Ｓ３７は、第１領域２３のエネルギーＥｂと、第１領域２３に作用する力Ｆｂとを含む物理量が計算される。このようなエネルギーＥｂと力Ｆｂとを含む物理量は、後述するＯＮＩＯＭ法に従って、高分子材料モデル１５のエネルギーＥｔ、及び、高分子材料モデル１５に作用する力Ｆｔを求めるのに利用される。エネルギーＥｂと力Ｆｂとを含む物理量は、コンピュータ１に記憶される。また、第２分子力学計算工程Ｓ３７では、離反粒子モデル部２１の粒子モデル３、３間に、結合ポテンシャルＰｎが設定される。

次に、本実施形態の切断計算工程Ｓ２３では、高分子材料モデル１５のエネルギーＥｔが計算される（工程Ｓ３８）。本実施形態のエネルギーＥｔは、ＯＮＩＯＭ（ Our own N-layered Integrated molecular Orbital and molecular Mechanics ） 法に従って、下記式（６）に基づいて計算される。
Ｅｔ＝Ｅａ−Ｅｂ＋Ｅｃ … （６）
ここで、
Ｅａ：第１分子力学計算工程で求められた高分子材料モデルのエネルギー
Ｅｂ：第２分子力学計算工程で求められた第１領域のエネルギー
Ｅｃ：量子力学計算工程で求められた第１領域のエネルギー

上記式（６）において、第１分子力学計算工程Ｓ３４で求められたエネルギーＥａは、高分子材料モデル１５の全体領域のエネルギーである。また、第２分子力学計算工程Ｓ３７で求められたエネルギーＥｂは、第１領域２３でのエネルギーである。これらのエネルギーＥａ、Ｅｂは、分子力学計算で求められているため、量子力学計算に比べて、計算精度が低い。量子力学計算工程Ｓ３６で求められたエネルギーＥｃは、第１領域２３のエネルギーである。このエネルギーＥｃは、量子力学計算で求められているため、分子力学計算に比べて、計算精度が高い。

ＯＮＩＯＭ法は、計算精度の異なる階層構造を組み合わせることによって最適解を得る方法である。本実施形態では、計算精度が比較的低い高分子材料モデルのエネルギーＥａのうち、第１領域のエネルギーＥｂを、計算精度が比較的高い第１領域のエネルギーＥｃに置き換えることによって、高分子材料モデル１５のエネルギーＥｔを求めている。このようなエネルギーＥｔは、分子鎖モデル１３の切断による影響を詳細に考慮することができるため、計算精度を高めることができる。また、本実施形態では、例えば、量子力学計算に基づいた高分子材料モデル１５の全体領域のエネルギーを求めることなく、高分子材料モデル１５のエネルギーＥｔを求めることができるため、計算時間を短縮することができる。エネルギーＥｔは、コンピュータ１に記憶される。

次に、本実施形態の切断計算工程Ｓ２３では、高分子材料モデル１５に作用する力Ｆｔが計算される（工程Ｓ３９）。本実施形態の力Ｆｔは、高分子材料モデルのエネルギーＥｔと同様に、ＯＮＩＯＭ法に従って、下記式（７）に基づいて計算される。
Ｆｔ＝Ｆａ−Ｆｂ＋Ｆｃ … （７）
ここで、
Ｆａ：第１分子力学計算工程で求められた高分子材料モデルに作用する力
Ｆｂ：第２分子力学計算工程で求められた第１領域に作用する力
Ｆｃ：量子力学計算工程で求められた第１領域に作用する力

本実施形態では、計算精度が比較的低い高分子材料モデル１５の力Ｆａのうち、第１領域２３の力Ｆｂを、計算精度が比較的高い第１領域２３の力Ｆｃに置き換えることによって、高分子材料モデル１５の力Ｆｔを求めている。このような力Ｆｔは、分子鎖モデル１３の切断による影響を詳細に考慮することができるため、計算精度を高めることができる。また、本実施形態では、量子力学計算に基づいた高分子材料モデル１５の全体領域に作用する力を計算する必要がないため、計算時間を短縮することができる。この力Ｆｔは、コンピュータ１に記憶される。

次に、本実施形態の切断計算工程Ｓ２３では、現在の単位時間ステップにおいて、分子鎖モデル１３の各粒子モデル３の座標値及び速度が更新される（工程Ｓ４０）。本実施形態では、ＯＮＩＯＭ法で計算されたエネルギーＥｔ及び力Ｆｔから加速度を求める。そして、エネルギーＥｔ及び力Ｆｔから加速度に基づいて、第１分子力学計算工程Ｓ３４が実施される前の分子鎖モデル１３の各粒子モデル３の座標値から、各粒子モデル３の新たな座標値（即ち、次の単位時間ステップで用いられる座標値）が計算される。

このような各粒子モデル３の座標値は、切断される可能性の高い第１領域２３を量子力学計算した座標値が含められているため、分子鎖モデル１３の切断を精度よく解析することができる。また、切断される可能性の低い部分については、比較的精度の低い分子力学計算が実施されるため、計算時間を短縮することができる。

本実施形態では、第１分子力学計算工程Ｓ３４、及び、第２分子力学計算工程Ｓ３７において、結合解離ポテンシャルＰｕに近似する結合ポテンシャルＰｎが定義されている。このため、第２分子力学計算工程Ｓ３７後の各粒子モデル３に作用する力と、量子力学計算工程Ｓ３６後の各粒子モデル３に作用する力との間に生じる差異を最小限に抑えることができる。従って、ＯＮＩＯＭ法に基づく、座標値の計算精度を高めることができる。また、各粒子モデル３の速度は、分子鎖モデル１３の各粒子モデル３の座標値、高分子材料モデル１５のエネルギーＥｔ、及び、高分子材料モデル１５に作用する力Ｆｔに基づいて計算される。

なお、上述した工程Ｓ３５において、離反粒子モデル部２１が存在しないと判断された場合（工程Ｓ３５で、「Ｎ」）は、第１領域２３を対象とした量子力学計算（工程Ｓ３６）及び分子力学計算（工程Ｓ３７）が実施されていない。この場合、工程Ｓ４０では、高分子材料モデル１５を対象とした分子力学計算によって求められた粒子モデル３の座標値、第１分子力学計算工程で求められたエネルギーＥａ、及び、第１分子力学計算工程で求められた力Ｆａから、座標値及び速度が更新される。

次に、本実施形態の切断計算工程Ｓ２３では、量子力学計算工程Ｓ３６後、本実施形態では量子力学計算工程Ｓ３６及び第２分子力学計算工程Ｓ３７後の第１領域２３について、粒子モデル３、３間の距離ｒが、予め定められた第２距離Ｌｂよりも大きい離反粒子モデル部２１が存在するか否かが判断される（工程Ｓ５１）。工程Ｓ５１では、工程Ｓ４０で更新された分子鎖モデル１３の各粒子モデル３の座標値に基づいて、第１領域２３に配置されている全ての離反粒子モデル部２１について、粒子モデル３、３間の距離ｒが計算される。

本実施形態の第２距離Ｌｂは、図１７に示した結合解離ポテンシャルＰｕの変曲点８での粒子モデル３、３間の距離Ｌｎに基づいて定義されている。本実施形態の第２距離Ｌｂは、第１距離Ｌａよりも大に設定されている。さらに、第２距離Ｌｂは、変曲点８での粒子モデル間の距離Ｌｎ以上に設定されている。

上述したように、変曲点８は、粒子モデル３、３間の結合と解離との境界点と考えることができる。このため、粒子モデル３、３間の距離ｒが第２距離Ｌｂよりも大きい場合、粒子モデル３、３間で切断すると判断できる。従って、工程Ｓ５１では、粒子モデル３、３間の距離ｒが第２距離Ｌｂよりも大きいか否かが判断されることにより、粒子モデル３、３間で切断される離反粒子モデル部２１の有無を判断できる。なお、粒子モデル３、３間で切断される離反粒子モデル部２１を精度良く判断するために、第２距離Ｌｂは、変曲点８での粒子モデル３、３間の距離Ｌｎよりも大が望ましく、変曲点８での距離Ｌｎの１００％より大が望ましく、また、１３０％以下が望ましい。

工程Ｓ５１において、粒子モデル３、３間の距離ｒが第２距離Ｌｂ以上の離反粒子モデル部２１が存在すると判断された場合（工程Ｓ５１で、「Ｙ」）、次の切断工程Ｓ５２が実施される。他方、全ての離反粒子モデル部２１において、粒子モデル３、３間の距離が第２距離Ｌｂ未満であると判断された場合（工程Ｓ５１で、「Ｎ」）、高分子材料モデル１５が予め定められた条件まで伸長したか否かを判断する後述の工程Ｓ５３が実施される。

次に、本実施形態の切断計算工程Ｓ２３では、粒子モデル３、３間が切断される（切断工程Ｓ５２）。切断工程Ｓ５２では、距離ｒが第２距離Ｌｂ以上であると判断された粒子モデル３、３について、粒子モデル３、３間のボンドモデル４の定義を無効にする。これにより、粒子モデル３、３間の切断が定義される。なお、粒子モデル３、３間の切断が定義されるのは、図１９に示した第１領域２３の粒子モデル３ではなく、第１領域２３の粒子モデル３に対応する分子鎖モデル１３（図１０に示す）の粒子モデル３、３間において、切断が定義される。

このように、本実施形態のシミュレーション方法では、離反粒子モデル部２１を含んだ第１領域２３を対象として、分子力学計算よりも計算精度の高い量子力学計算の結果に基づき、粒子モデル３、３間を切断している。従って、本実施形態では、分子鎖モデル１３の切断を、実際の分子鎖の切断に近似させることができ、高分子材料を精度良く解析するのに役立つ。

また、結合解離ポテンシャルＰｕは、本実施形態の計算方法に基づいて求められるため、原子間の結合及び解離を精度良く表現しうる。このため、本実施形態のシミュレーション方法では、粒子モデル３、３間の切断を、実際の分子鎖の切断により近似させることができる。

さらに、量子力学計算が実施される離反粒子モデル部２１の有無は、工程Ｓ３５において、量子力学計算よりも計算時間が短い分子力学計算の結果に基づいて判断される。このため、高分子材料モデル１５の全体領域を対象とした量子力学計算を実施する必要がない。従って、本実施形態のシミュレーション方法では、計算時間を短縮することができる。

なお、切断の定義は、粒子モデル３、３間を拘束しているボンドモデル４を無効にすることによって設定することができる。また、切断された粒子モデル３、３間に定義されていた結合ポテンシャルＰｎが無効とされ、ボンドモデル４を介さずに隣り合う粒子モデル３、３間に定義される第２ポテンシャルＰ２が定義される。

切断された粒子モデル３、３には、例えば、粒子モデル３でモデル化された原子の原子価（例えば、硫黄原子の場合：２）から、粒子モデル３の結合数（例えば、１）を減じた数の水素粒子モデル３Ｈ（例えば、１個）が、ボンドモデル４を介して連結される。これにより、分子鎖モデル１３は、安定した構造を維持することができる。

次に、本実施形態の切断計算工程Ｓ２３では、高分子材料モデル１５が、予め定められた条件まで伸長したか否かが判断される（工程Ｓ５３）。工程Ｓ５３では、高分子材料モデル１５が上限値まで伸長したと判断された場合（工程Ｓ５３で、「Ｙ」）、次の工程Ｓ２４が実施される。他方、高分子材料モデル１５が上限値まで伸長していないと判断された場合は（工程Ｓ５３で、「Ｎ」）、単位時間ステップを一つ進展させて（工程Ｓ５４）、工程Ｓ３３〜工程Ｓ５３が再度実施される。これにより、切断計算工程Ｓ２３では、高分子材料モデルを上限値まで伸長させて、分子鎖モデル１３の切断を解析することができる。

次に、本実施形態のシミュレーション方法では、高分子材料モデル１５の物理量が、許容範囲内であるか否かが判断される（工程Ｓ２４）。物理量の許容範囲としては、例えば、高分子材料に求められる性能に基づいて、適宜設定されうる。工程Ｓ２４において、高分子材料モデル１５の物理量が許容範囲内であると判断された場合（工程Ｓ２４で、「Ｙ」）、高分子材料モデル１５に基づいて、高分子材料が製造される（工程Ｓ２５）。他方、高分子材料モデル１５の物理量が、許容範囲外と判断された場合は（工程Ｓ２４で、「Ｎ」）、分子鎖モデル１３の諸条件が変更されて（工程Ｓ２６）、工程Ｓ２１〜工程Ｓ２４が再度実施される。このように、本実施形態では、高分子材料モデル１５の物理量が許容範囲内になるまで、分子鎖モデル１３の諸条件が変更されるため、所望の性能を有する高分子材料を、効率よく設計することができる。

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

［実施例Ａ］
図２及び図４に示した処理手順に従って、分子を構成する原子間の結合解離ポテンシャルが計算された（実施例１）。実施例１では、密度汎関数法（ＤＦＴ）の理論に基づいて量子力学計算が行われ、安定した波動関数のみから結合解離ポテンシャルが計算された。

比較のために、従来の方法と同様の手順に基づいて、波動関数が安定しているか否かを判定する評価工程を行わずに、結合解離ポテンシャルが計算された（比較例１〜５）。即ち、比較例１〜５の結合解離ポテンシャルは、安定した波動関数、及び、不安定な波動関数の双方から計算される。比較例１〜５では、表１に示した理論に基づいて、量子力学計算が実施された。共通仕様は次のとおりである。
コンピュータ：ＳＧＩ社のワークステーション
ＣＰＵのコア数：１２コア
搭載メモリ：６４ＧＢ
分子：メタンチオール
図２０は、実施例１及び比較例１〜５の結合解離ポテンシャルのグラフである。計算時間を表１に示す。

テストの結果、実施例１の結合解離ポテンシャルは、実施例１と同一の理論（ＤＦＴ）に基づいて計算された比較例１に比べて、最も計算精度が高い理論（ＭＲＣＩ）に基づいて計算された比較例２の結合解離ポテンシャルに近似させることができた。従って、実施例は、原子間の結合及び解離を精度よく表現できるポテンシャルを定義できた。表１に示されるように、実施例１は、比較例２〜５に比べて、計算時間を大幅に短縮できた。従って、実施例１は、計算時間を短縮しつつ、結合解離ポテンシャルを精度よく求めることができた。

［実施例Ｂ］
図９及び図１１に示した処理手順に従って、分子鎖モデル及び高分子モデルが定義され、分子鎖モデルを切断する切断工程が実施された（実施例２、比較例６）。

実施例２の切断工程では、図１３、図１５及び図１８に示した処理手順に従って、高分子材料モデルに基づいて分子力学計算を行う工程、粒子モデル間の距離が第１距離Ｌａよりも大きい離反粒子モデル部を含む第１領域を対象に、量子力学計算を行う工程、及び、第１領域を対象に分子力学計算を行う工程が実施され、上記式（６）、（７）に基づいて、高分子材料モデルのエネルギーＥｔ及び高分子材料モデルに作用する力Ｆｔが計算された。

分子力学計算を行う工程では、分子鎖モデルの粒子モデル間に、結合解離ポテンシャルに近似する近似ポテンシャルが定義された。なお、結合解離ポテンシャルは、図２及び図４に示した処理手順に従って、安定した波動関数のみから求められた。実施例２の切断工程では、粒子モデル間の距離が第２距離Ｌｂよりも大きいときに、それらの間を切断する工程が実施された。

比較例６の切断工程では、高分子材料モデルの全体領域について、量子力学計算が実施され、粒子モデル間の距離に応じて、分子鎖モデルの切断が計算された。なお、粒子モデル間に設定されるポテンシャルは、明細書に記載のとおりである。共通仕様は、次のとおりであり、使用されたコンピュータは、実施例Ａと同一である。
高分子モデル：
１辺の長さＬ１：２．３nm
分子鎖モデルの個数：１個
分子鎖モデル：
ポリイソプレン（イソプレン分子：１００重合）に硫黄（１原子）を架橋
粒子モデル（原子）の個数：１３０１個
変曲点での距離Ｌｎ：３．０９Å
第１距離Ｌａ：２．６９Å
第２距離Ｌｂ：３．３５Å
伸長シミュレーション：
伸長方向：一軸方向（Ｚ軸方向）
伸長速度Ｖａ：２５０ｍ／sec

図２１は、実施例２の粒子モデル間の切断個数と、高分子材料モデルの伸長率との関係を示したグラフである。図２１に示されるように、実施例２では、高分子材料モデルが伸長するに従って、粒子モデル間で切断する個数が増加する状態を示すことができた。また、実施例２では、伸長率が３００％程度でＣ−Ｓ結合の切断が計算され、伸長率が５６０％程度からＣ−Ｃ結合の切断が計算された。このように、実施例２では、粒子モデル毎に、切断が開始する伸長率を計算することができた。

図２２は、実施例２の高分子材料の応力と、高分子材料モデルの伸長率との関係を示したグラフである。図１５に示されるように、伸長率が５００％程度から応力が急激に大きくなっているが、伸長率が６００％程度から応力が急激に減少している。これは、Ｃ−Ｃ結合の切断によって、応力が減少したためである。

このように、実施例２では、高分子材料モデルを用いた伸長シミュレーションにより、分子鎖モデルの切断を精度良く計算できたため、高分子材料を精度良く解析することができた。また、実施例２では、本発明の計算方法に基づいて、安定した波動関数のみから結合解離ポテンシャルが計算されるため、実際の分子鎖の切断に近似させることができた。さらに、実施例２では、粒子モデル間で切断する可能性が高い離反粒子モデル部を含む第１領域のみを対象に、量子力学計算が実施されるため、計算時間を短縮することができた。

他方、比較例６では、高分子材料モデルの全体領域について、量子力学計算が実施されるため、実施例２と同様の変形シミュレーションを実施するには、上記コンピュータの処理能力において、事実上不可能であると考えられる。従って、比較例６では、分子鎖モデルの切断を計算することができなかった。

Ｓ１ 分子モデルを定義する工程
Ｓ４ 計算工程
Ｓ５ 評価工程
Ｓ６ 結合解離ポテンシャルを求める工程
Ｓ７ 再計算工程

Claims (9)

1. コンピュータを用いて、分子を構成する原子間の結合解離ポテンシャルを計算するための方法であって、
   前記原子をモデル化した粒子モデルを互いに結合した分子モデルを、前記コンピュータに定義する工程と、
   前記コンピュータが、前記分子モデルを用いた量子力学計算によって波動関数を計算する計算工程と、
   前記コンピュータが、前記波動関数が安定しているか否かを判定する評価工程とを含み、
   前記波動関数が安定していると判定した場合、前記波動関数に基づいて、前記粒子モデル間の結合解離ポテンシャルを求める工程が実施され、
   前記波動関数が安定していないと判定した場合、前記波動関数のエネルギーが低くなるように、前記波動関数を再計算する再計算工程が実施されることを特徴とする原子間の結合解離ポテンシャルの計算方法。
2. 前記量子力学計算は、密度汎関数法に基づいて行われる請求項１記載の原子間の結合解離ポテンシャルの計算方法。
3. 前記計算工程は、前記波動関数のエネルギーを計算する工程を含み、
   前記評価工程は、前記エネルギーを前記波動関数で２次変分した固有値に、負の固有値が含まれない場合に、前記波動関数が安定していると判定する請求項１又は２記載の原子間の結合解離ポテンシャルの計算方法。
4. 前記再計算工程は、前記負の固有値が無くなるように、前記波動関数を再計算する請求項３記載の原子間の結合解離ポテンシャルの計算方法。
5. 請求項１乃至４のいずれかに記載の計算方法で求められた前記結合解離ポテンシャルを用いて、高分子材料の分子鎖の切断を解析するためのシミュレーション方法であって、
   前記分子鎖に基づいて、複数の前記粒子モデルと、前記粒子モデル間を結合するボンドモデルとを含む分子鎖モデルを、前記コンピュータに定義する工程、
   前記コンピュータに定義された空間に、少なくとも一つの前記分子鎖モデルを配置した高分子材料モデルを定義する工程、及び
   前記コンピュータが、予め定められた条件に基づいて、前記分子鎖モデルの切断を計算する切断計算工程を含み、
   前記切断計算工程は、前記高分子材料モデルの前記粒子モデル間に、前記結合解離ポテンシャルに近似する結合ポテンシャルを定義する工程を含むことを特徴とする高分子材料のシミュレーション方法。
6. 前記結合ポテンシャルは、前記結合解離ポテンシャルのエネルギーの下限値から、前記エネルギーの上限値の５０％〜８０％までの範囲において、前記結合解離ポテンシャルに近似する関数で定義される請求項５記載の高分子材料のシミュレーション方法。
7. 前記切断計算工程は、前記分子力学計算工程後の前記分子鎖モデルについて、前記ボンドモデルを介して隣り合う前記粒子モデル間の距離が、予め定められた第１距離よりも大きい離反粒子モデル部が存在するか否かを判断する工程をさらに含み、
   前記離反粒子モデル部が存在すると判断された場合、前記離反粒子モデル部を含む非全体領域である第１領域を対象に量子力学計算を行う量子力学計算工程と、
   前記量子力学計算工程後の前記第１領域について、前記粒子モデル間の距離が、予め定められた第２距離よりも大きいときに、それらの間を切断する工程とを含む請求項５又は６記載の高分子材料のシミュレーション方法。
8. 前記第１距離は、前記第２距離よりも小である請求項７記載の高分子材料のシミュレーション方法。
9. 前記第２距離は、前記結合解離ポテンシャルについて、そのエネルギーの増分が小さくなる変曲点での前記粒子モデル間の距離に基づいて定義される請求項７又は８記載の高分子材料のシミュレーション方法。