(57)【要約】【課題】  分子鎖モデルの切断を安定して計算する。【解決手段】  コンピュータを用いて、高分子材料の分子鎖の切断を解析するためのシミュレーション方法である。このシミュレーション方法では、高分子材料の一部に対応する仮想空間であるセルに、分子鎖モデルを配置した高分子材料モデルを定義する工程、及び、分子鎖モデルの切断を計算するための切断計算工程を含む。切断計算工程では、高分子材料モデルを用いて、分子力学計算を行う第１分子力学計算工程Ｓ４６、第１分子力学計算工程Ｓ４６後の分子鎖モデルについて、粒子モデル間の距離が第１距離よりも大きい粒子モデル対を含む第１領域を対象に量子力学計算を行う工程Ｓ４８と、第１領域について、第１距離よりも大きく設定され、かつ、結合解離ポテンシャルの変曲点での粒子モデル間の距離に１．０〜２．５Åを加えた第２距離よりも、粒子モデル間の距離が大きいときに、粒子モデル対に対応する分子鎖モデルの粒子モデル間の結合を切断する工程Ｓ５４を含む。【選択図】図８
2018032077.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータを用いて、高分子材料の分子鎖の切断を解析するためのシミュレーション方法であって、
前記分子鎖に基づいて、複数の粒子モデルと、前記粒子モデル間を結合するボンドモデルとを含む分子鎖モデルを、前記コンピュータに定義する工程、
前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに定義する工程、
前記セルに、少なくとも一つの前記分子鎖モデルを配置した高分子材料モデルを定義する工程、及び
前記コンピュータが、予め定められた条件に基づいて、前記分子鎖モデルの切断を計算するための切断計算工程を含み、
前記切断計算工程は、前記分子鎖モデルから選択された一対の粒子モデルと、前記粒子モデル間を連結するボンドモデルとを含む分子モデルを設定する工程、
前記分子モデルを用いた量子力学計算に基づいて、エネルギーの増分が小さくなる変曲点を有する結合解離ポテンシャルを求める工程、
前記高分子材料モデルを用いて、分子力学計算を行う第１分子力学計算工程、及び
前記第１分子力学計算工程後の前記分子鎖モデルについて、前記粒子モデル間の距離が、予め定められた第１距離よりも大きい粒子モデル対が存在するか否かを判断する工程を含み、
前記粒子モデル対が存在すると判断された場合、前記粒子モデル対を含む第１領域を対象に量子力学計算を行う量子力学計算工程と、
前記量子力学計算工程後の前記第１領域について、前記第１距離よりも大きく設定されかつ前記変曲点での前記粒子モデル間の距離に１．０〜２．５Åを加えた第２距離よりも、前記粒子モデル間の距離が大きいときに、前記粒子モデル対に対応する前記分子鎖モデルの前記粒子モデル間の結合を切断する切断工程とを含むことを特徴とする高分子材料のシミュレーション方法。
【請求項２】
前記第１距離は、前記変曲点での前記粒子モデル間の距離よりも小さい請求項１記載の高分子材料のシミュレーション方法。
【請求項３】
前記切断計算工程は、前記第１領域を対象に分子力学計算を行う第２分子力学計算工程をさらに含む請求項１又は２記載の高分子材料のシミュレーション方法。
【請求項４】
前記切断計算工程は、前記第１分子力学計算工程、前記第２分子力学計算工程、及び、前記量子力学計算工程の計算結果に基づいて、前記高分子材料モデルのエネルギーＥｔ、及び、前記高分子材料モデルに作用する力Ｆｔを含む物理量を計算する工程、並びに
前記物理量に基づいて、前記粒子モデルの座標値を更新する工程をさらに含み、
前記切断工程は、前記座標値から計算された前記粒子モデル間の距離に基づいて、前記粒子モデル間の結合を切断する請求項３記載の高分子材料のシミュレーション方法。
【請求項５】
前記エネルギーＥｔは、下記式に基づいて計算される請求項４記載の高分子材料のシミュレーション方法。
Ｅｔ＝Ｅａ−Ｅｂ＋Ｅｃ
ここで、
Ｅａ：第１分子力学計算工程で求められた高分子材料モデルのエネルギー
Ｅｂ：第２分子力学計算工程で求められた第１領域のエネルギー
Ｅｃ：量子力学計算工程で求められた第１領域のエネルギー
【請求項６】
前記力Ｆｔは、下記式に基づいて計算される請求項４又は５記載の高分子材料のシミュレーション方法。
Ｆｔ＝Ｆａ−Ｆｂ＋Ｆｃ
ここで、
Ｆａ：第１分子力学計算工程で求められた高分子材料モデルに作用する力
Ｆｂ：第２分子力学計算工程で求められた第１領域に作用する力
Ｆｃ：量子力学計算工程で求められた第１領域に作用する力
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、高分子材料の解析に役立つ高分子材料のシミュレーション方法に関する。
【背景技術】
【０００２】
近年、ゴム等の高分子材料の開発のために、高分子材料の性質を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている（例えば、下記特許文献１参照）。
【０００３】
この種のシミュレーション方法では、先ず、高分子材料の一部に対応する仮想空間であるセルに、複数の粒子モデルを含む分子鎖モデルを配置した高分子材料モデルが、コンピュータに定義される。高分子材料モデルの粒子モデル間には、粒子モデル間の距離の関数であるポテンシャルが定義されている。このような高分子材料モデルは、上記空間内で分子動力学に基づく緩和計算が行われる。その後、例えば、高分子材料モデルを用いた変形計算が実施され、高分子材料の性能等が解析される。
【先行技術文献】
【特許文献】
【０００４】
【特許文献１】特開２０１３−１９５２２０号公報
【発明の概要】
【発明が解決しようとする課題】
【０００５】
実際の高分子材料では、例えば大きな力が作用すると、分子鎖が切断することが知られている。しかしながら、上記シミュレーション方法では、分子鎖モデルを切断する方法が未だ十分に確立されていない。このような分子鎖モデルの切断を、安定して計算できる方法が求められていた。
【０００６】
本発明は、以上のような実状に鑑み案出されたもので、分子鎖モデルの切断を安定して計算することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【０００７】
本発明は、コンピュータを用いて、高分子材料の分子鎖の切断を解析するためのシミュレーション方法であって、前記分子鎖に基づいて、複数の粒子モデルと、前記粒子モデル間を結合するボンドモデルとを含む分子鎖モデルを、前記コンピュータに定義する工程、前記高分子材料の一部に対応する仮想空間であるセルを、前記コンピュータに定義する工程、前記セルに、少なくとも一つの前記分子鎖モデルを配置した高分子材料モデルを定義する工程、及び前記コンピュータが、予め定められた条件に基づいて、前記分子鎖モデルの切断を計算するための切断計算工程を含み、前記切断計算工程は、前記分子鎖モデルから選択された一対の粒子モデルと、前記粒子モデル間を連結するボンドモデルとを含む分子モデルを設定する工程、前記分子モデルを用いた量子力学計算に基づいて、エネルギーの増分が小さくなる変曲点を有する結合解離ポテンシャルを求める工程、前記高分子材料モデルを用いて、分子力学計算を行う第１分子力学計算工程、及び前記第１分子力学計算工程後の前記分子鎖モデルについて、前記粒子モデル間の距離が、予め定められた第１距離よりも大きい粒子モデル対が存在するか否かを判断する工程を含み、前記粒子モデル対が存在すると判断された場合、前記粒子モデル対を含む第１領域を対象に量子力学計算を行う量子力学計算工程と、前記量子力学計算工程後の前記第１領域について、前記第１距離よりも大きく設定されかつ前記変曲点での前記粒子モデル間の距離に１．０〜２．５Åを加えた第２距離よりも、前記粒子モデル間の距離が大きいときに、前記粒子モデル対に対応する前記分子鎖モデルの前記粒子モデル間の結合を切断する切断工程とを含むことを特徴とする。
【０００８】
本発明に係る前記高分子材料のシミュレーション方法において、前記第１距離は、前記変曲点での前記粒子モデル間の距離よりも小さいのが望ましい。
【０００９】
本発明に係る前記高分子材料のシミュレーション方法において、前記切断計算工程は、前記第１領域を対象に分子力学計算を行う第２分子力学計算工程をさらに含むのが望ましい。
【００１０】
本発明に係る前記高分子材料のシミュレーション方法において、前記切断計算工程は、前記第１分子力学計算工程、前記第２分子力学計算工程、及び、前記量子力学計算工程の計算結果に基づいて、前記高分子材料モデルのエネルギーＥｔ、及び、前記高分子材料モデルに作用する力Ｆｔを含む物理量を計算する工程、並びに前記物理量に基づいて、前記粒子モデルの座標値を更新する工程をさらに含み、前記切断工程は、前記座標値から計算された前記粒子モデル間の距離に基づいて、前記粒子モデル間の結合を切断するのが望ましい。
【００１１】
本発明に係る前記高分子材料のシミュレーション方法において、前記エネルギーＥｔは、下記式に基づいて計算されるのが望ましい。
Ｅｔ＝Ｅａ−Ｅｂ＋Ｅｃ
ここで、
Ｅａ：第１分子力学計算工程で求められた高分子材料モデルのエネルギー
Ｅｂ：第２分子力学計算工程で求められた第１領域のエネルギー
Ｅｃ：量子力学計算工程で求められた第１領域のエネルギー
【００１２】
本発明に係る前記高分子材料のシミュレーション方法において、前記力Ｆｔは、下記式に基づいて計算されるのが望ましい。
Ｆｔ＝Ｆａ−Ｆｂ＋Ｆｃ
ここで、
Ｆａ：第１分子力学計算工程で求められた高分子材料モデルに作用する力
Ｆｂ：第２分子力学計算工程で求められた第１領域に作用する力
Ｆｃ：量子力学計算工程で求められた第１領域に作用する力
【発明の効果】
【００１３】
本発明の高分子材料のシミュレーション方法は、コンピュータが、予め定められた条件に基づいて、分子鎖モデルの切断を計算するための切断計算工程を含んでいる。
【００１４】
切断計算工程は、前記分子鎖モデルから選択された一対の粒子モデルと、粒子モデル間を連結するボンドモデルとを含む分子モデルを設定する工程、及び、分子モデルを用いた量子力学計算に基づいて、エネルギーの増分が小さくなる変曲点を有する結合解離ポテンシャルを求める工程を含んでいる。
【００１５】
変曲点は、粒子モデル間の結合と解離との境界と考えることができる。即ち、粒子モデル間の距離が変曲点よりも大きくなると、粒子モデル間で切断、又は、結合次数の減少が生じていると考えることができる。従って、結合解離ポテンシャルは、分子鎖モデルの切断を計算するのに役立つ。
【００１６】
切断計算工程は、高分子材料モデルを用いて、分子力学計算を行う第１分子力学計算工程、及び、第１分子力学計算工程後の分子鎖モデルについて、粒子モデル間の距離が、予め定められた第１距離よりも大きい粒子モデル対が存在するか否かを判断する工程を含んでいる。
【００１７】
粒子モデル対が存在すると判断された場合、粒子モデル対を含む第１領域を対象に量子力学計算を行う量子力学計算工程と、量子力学計算工程後の第１領域について、粒子モデル間の距離が、第１距離よりも大きく設定されかつ変曲点での粒子モデル間に距離に１．０〜２．５Åを加えた第２距離よりも大きいときに、粒子モデル対に対応する分子鎖モデルの粒子モデル間の結合を切断する切断工程とを含んでいる。
【００１８】
粒子モデル対は、分子鎖モデルにおいて、粒子モデル間で切断する可能性が高い部分である。この粒子モデル対を含んだ第１領域を対象に、分子力学計算よりも計算精度の高い量子力学計算が行われている。このため、本発明のシミュレーション方法によれば、粒子モデル間の切断を、実際の分子鎖の切断に近似させることができ、高分子材料を精度良く解析することができる。
【００１９】
また、粒子モデル対の有無は、量子力学計算よりも計算時間が短い分子力学計算に基づいて判断される。このため、本発明のシミュレーション方法によれば、高分子材料モデルの全体領域を対象とした量子力学計算を実施する必要がないため、計算時間を短縮できる。
【００２０】
一般に、量子力学計算で求められる第１領域のエネルギーと、分子力学計算で求められる第１領域のエネルギーとは互いに異なる。また、量子力学計算と分子力学計算とのエネルギーの差異は、粒子モデル間の距離が近いほど大きくなる一方、粒子モデル間の距離が大きいほど小さくなる。発明者らは、鋭意研究を重ねた結果、粒子モデル間の距離を十分に確保せずに、粒子モデル間の結合を切断すると、切断された粒子モデル間において、前記差異が大きくなり、計算が異常終了しやすくなることを知見した。
【００２１】
本発明では、変曲点での粒子モデル間の距離に１．０〜２．５Åを加えた第２距離よりも、粒子モデル間の距離が大きいときに、分子鎖モデルの粒子モデル間の結合を切断しているため、前記差異が大きくなるのを防ぐことができる。従って、本発明のシミュレーションによれば、分子鎖モデルの切断を安定して計算することができる。
【図面の簡単な説明】
【００２２】
【図１】高分子材料のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。
【図２】ポリイソプレンの構造式である。
【図３】シミュレーション方法の処理手順の一例を示すフローチャートである。
【図４】分子鎖モデルの一例を示す概念図である。
【図５】分子鎖モデル設定工程の処理手順の一例を示すフローチャートである。
【図６】（ａ）〜（ｃ）は、第１ポテンシャルを説明する分子鎖モデルの部分図である。
【図７】高分子材料モデルの一例を示す概念図である。
【図８】切断計算工程の処理手順の一例を示すフローチャートである。
【図９】高分子材料モデルの伸長シミュレーションを説明する側面図である。
【図１０】分子モデルの一例を示す図である。
【図１１】結合解離ポテンシャルの一例を示すグラフである。
【図１２】結合解離ポテンシャル及び結合ポテンシャルの一例を示すグラフである。
【図１３】量子力学計算工程の処理手順の一例を示すフローチャートである。
【図１４】第１領域の一例を示す図である。
【図１５】（ａ）は、粒子モデル間で切断する前の分子鎖モデル５一部を示す図、（ｂ）は、粒子モデル間で切断した分子モデルの一部を示す図である。
【発明を実施するための形態】
【００２３】
以下、本発明の実施の一形態が図面に基づき説明される。
本発明の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、コンピュータを用いて、高分子材料の分子鎖の切断を解析するための方法である。
【００２４】
図１は、本発明のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成されている。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２が設けられている。また、記憶装置には、本実施形態のシミュレーション方法を実行するための処理手順（プログラム）が予め記憶されている。
【００２５】
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態の高分子材料は、cis-１，４ポリイソプレン（以下、単に「ポリイソプレン」ということがある。）である。図２は、ポリイソプレンの構造式である。
【００２６】
ポリイソプレンを構成する分子鎖２Ａは、メチン基等（例えば、−ＣＨ＝、＞Ｃ＝）、メチレン基（−ＣＨ２−）、及び、メチル基（−ＣＨ３）によって構成されるイソプレンのモノマー（イソプレン分子）３が、重合度ｎで連結されて構成されている。なお、高分子材料には、ポリイソプレン以外の高分子材料が用いられてもよい。
【００２７】
図３は、シミュレーション方法の処理手順の一例を示すフローチャートである。本実施形態のシミュレーション方法では、先ず、図２に示した分子鎖２Ａに基づいて、分子鎖モデルがコンピュータ１に定義される（分子鎖モデル設定工程Ｓ１）。図４は、分子鎖モデル５の一例を示す概念図である。
【００２８】
分子鎖モデル５は、全原子モデルとして構成されている。本実施形態の分子鎖モデル５は、一つの硫黄原子が架橋されたポリイソプレンとして設定されている。分子鎖モデル５は、複数の粒子モデル６と、粒子モデル６、６間を結合するボンドモデル７とを含んで構成されている。
【００２９】
粒子モデル６及びボンドモデル７は、図２に示した分子鎖２Ａのモノマー３をなす単位構造に基づいて連結される。これにより、モノマーモデル９が設定される。このモノマーモデル９は、分子量（重合度）Ｍnに基づいて連結される。これにより、分子鎖モデル５が設定される。図５は、分子鎖モデル設定工程Ｓ１の処理手順の一例を示すフローチャートである。
【００３０】
本実施形態の分子鎖モデル設定工程Ｓ１では、先ず、図４に示されるように、粒子モデル６が設定される（工程Ｓ１１）。粒子モデル６は、後述する分子動力学計算、分子力学計算、及び、量子力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル６は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。
【００３１】
本実施形態の粒子モデル６は、図２に示した分子鎖２Ａの炭素原子をモデル化した炭素粒子モデル６Ｃ、水素原子をモデル化した水素粒子モデル６Ｈ、及び、硫黄原子をモデル化した硫黄粒子モデル６Ｓを含んでいる。これらの粒子モデル６は、コンピュータ１に記憶される。
【００３２】
次に、分子鎖モデル設定工程Ｓ１では、ボンドモデル７が設定される（工程Ｓ１２）。ボンドモデル７は、粒子モデル６、６間を拘束するものである。本実施形態のボンドモデル７は、主鎖７ａと側鎖７ｂとを含んでいる。主鎖７ａは、炭素粒子モデル６Ｃ、６Ｃを連結するものである。側鎖７ｂは、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとの間、及び、炭素粒子モデル６Ｃと硫黄粒子モデル６Ｓとの間を連結するものである。これらのボンドモデル７は、コンピュータ１に記憶される。
【００３３】
次に、分子鎖モデル設定工程Ｓ１では、ボンドモデル７を介して隣り合う粒子モデル６、６間には、相互作用（斥力及び引力を含む）が生じさせるポテンシャルが定義される（工程Ｓ１３）。図４に示されるように、ポテンシャルには、ボンドモデル７を介して隣り合う粒子モデル６、６間に定義される第１ポテンシャルＰ１と、ボンドモデル７を介さずに隣り合う粒子モデル６、６間に定義される第２ポテンシャルＰ２とが定義される。
【００３４】
図６（ａ）〜（ｃ）は、第１ポテンシャルＰ１を説明する分子鎖モデル５の部分図である。図６（ａ）〜（ｂ）に示されるように、分子鎖モデル５には、各粒子モデル６、６間の結合長さである結合距離（結合長）ｒ、及び、ボンドモデル７を介して連続する３つの粒子モデル６がなす角度である結合角θが定義されている。さらに、図６（ｃ）に示されるように、分子鎖モデル５には、ボンドモデル７を介して連続する４つの粒子モデル６において、隣り合う３つの粒子モデル６が作る一方の平面８Ａと他方の平面８Ｂとのなす角度ある二面角φが定義される。結合距離ｒ、結合角θ及び二面角φは、分子鎖モデル５に作用する外力又は内力によって変化する。
【００３５】
結合距離ｒ、結合角θ及び二面角φは、下記式（１）で定義される結合ポテンシャルＵbond（ｒ）、下記式（２）で定義される結合角ポテンシャルＵangle（θ）、及び、下記式（３）で定義される結合二面角ポテンシャルＵtorsion（φ）によって設定される。
【数１】
2018032077.tif 000003
【数２】
2018032077.tif 000004
【数３】
2018032077.tif 000005
ここで、各定数及び変数は、次のとおりである。
ｒ：結合距離
ｒ0：平衡距離（平衡長）
ｋ1、ｋ2：ばね定数
θ：結合角
θ0：平衡角度
ｋ3：二面角ポテンシャルの強度
Ｎ−1：二面角ポテンシャル多項式の次数
φ：二面角
Ａn：二面角定数
なお、結合距離ｒ及び平衡距離ｒ0は、各粒子モデル６の中心（図示省略）間の距離として定義される。
【００３６】
各結合ポテンシャルＵbond（ｒ）、結合角ポテンシャルＵangle（θ）、及び、結合二面角ポテンシャルＵtorsion（φ）の各定数については、適宜設定することができる。なお、これらのポテンシャルは、論文（J. Phys. Chem. 94, 8897（1990））に基づいて、分子鎖Ｍｃの構造に応じて設定されるのが望ましい。
【００３７】
第２ポテンシャルＰ２は、上述したように、ボンドモデル７を介さずに隣り合う粒子モデル６、６間に定義されるものである。本実施形態の第２ポテンシャルＰ２は、下記式（４）で定義されるＬＪポテンシャルＵLJ（ｒij）である。
【００３８】
【数４】
2018032077.tif 000006
ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒij：粒子モデル間の距離
ε：粒子モデル間に定義されるＬＪポテンシャルの強度
σ：粒子モデルの直径に相当
なお、距離ｒijは、各粒子モデル６、６の中心（図示省略）間の距離として定義される。
【００３９】
第２ポテンシャルＰ２は、粒子モデル間の距離ｒijがσよりも小さくなるほど、粒子モデル６、６間に作用する斥力が大きくなる。また、第２ポテンシャルＰ２は、粒子モデル間の距離ｒijがσになるときに最小となり、粒子モデル６、６間に斥力や引力は働かない。さらに、第２ポテンシャルＰ２は、粒子モデル間の距離ｒijがσよりも大になると、粒子モデル６、６間に作用する引力が働く。このように、第２ポテンシャルＰ２は、粒子モデル間の距離ｒijに応じて、斥力及び引力を定義することができる。
【００４０】
本実施形態では、硫黄粒子モデル６Ｓ、６Ｓ間、炭素粒子モデル６Ｃ、６Ｃ間、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとの間、硫黄粒子モデル６Ｓと炭素粒子モデル６Ｃとの間、及び、硫黄粒子モデル６Ｓと水素粒子モデル６Ｈとの間に、それぞれ異なる第２ポテンシャルＰ２が設定されている。各第２ポテンシャルＰ２は、上記式（４）の定数がそれぞれ異なっている。なお、各定数は、例えば、論文（J. Phys. Chem. 94, 8897（1990））に基づいて、適宜設定することができる。
【００４１】
次に、本実施形態のシミュレーション方法では、高分子材料の一部に対応する仮想空間であるセルが、コンピュータ１に定義される（工程Ｓ２）。図７は、高分子材料モデルの一例を示す概念図である。
【００４２】
本実施形態のセル１７は、分子鎖モデル５を配置して、高分子材料モデル１６を定義するためのものである。セル１７は、互いに向き合う三対の平面１７ａ、１７ｂを有する直方体として定義されている。各平面１７ａ、１７ｂには、周期境界条件が定義されている。このようなセル１７は、例えば、一方の平面１７ａから出て行った分子鎖モデル５の一部が、反対側の平面１７ｂから入ってくるように計算することができる。従って、一方の平面１７ａと、反対側の平面１７ｂとが連続している（繋がっている）ものとして取り扱うことができる。
【００４３】
セル１７の一辺の長さＬ１、Ｌ２、Ｌ３は、適宜設定することができる。本実施形態の長さＬ１は、分子鎖モデル５の慣性半径（図示省略）の２倍以上が望ましい。慣性半径は、分子鎖モデル５の拡がりを示す量である。このようなセル１７では、分子動力学計算において、周期境界条件による自己のイメージとの衝突が起こりにくいため、分子鎖モデル５の空間的拡がりを適切に計算することができる。さらに、セル１７の大きさは、例えば１気圧で安定な体積に設定される。このようなセル１７は、高分子材料の少なくとも一部の体積を定義することができる。セル１７は、コンピュータ１に記憶される。
【００４４】
次に、本実施形態のシミュレーション方法では、コンピュータ１に、高分子材料モデルが定義される（工程Ｓ３）。高分子材料モデル１６は、コンピュータ１上に定義されたセル１７に、少なくとも一つの分子鎖モデル５を配置することによって定義される。セル１７に配置される分子鎖モデル５の個数については、適宜設定することができる。本実施形態の分子鎖モデル５の個数としては、例えば、１個以上、好ましくは３０個以上に設定されるのが望ましい。
【００４５】
本実施形態では、分子鎖モデル５が配置されたセル１７を用いて、分子動力学（ Molecular Dynamics ： ＭＤ ）計算による緩和計算が実施される。これにより、高分子材料モデル１６が設定される。
【００４６】
本実施形態の分子動力学計算では、例えば、セル１７について所定の時間、分子鎖モデル５が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粒子モデル６の動きが、単位時間ステップ毎に追跡される。このような構造緩和の計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ−ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。
【００４７】
分子動力学計算では、セル１７において、圧力及び温度が一定、又は体積及び温度が一定に保たれる。また、分子動力学計算では、分子鎖モデル５の人為的な初期配置が十分に排除されたとみなせるまで行われる。これにより、工程Ｓ３では、実際の高分子材料の分子運動に近似させて、分子鎖モデル５の初期配置を精度よく緩和することができる。このような緩和計算を経て、高分子材料モデル１６が定義される。高分子材料モデル１６は、コンピュータ１に記憶される。
【００４８】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、予め定められた条件に基づいて、分子鎖モデル５の切断を計算する（切断計算工程Ｓ４）。本実施形態の切断計算工程Ｓ４では、図７に示した高分子材料モデル１６を、一軸方向（例えば、Ｚ軸方向）に伸長させる伸長シミュレーションが実施される。
【００４９】
本実施形態の伸長シミュレーションでは、分子力学計算（Molecular Mechanics ： ＭＭ）、及び、量子力学計算（Quantum Mechanics ： ＱＭ）を用いて、単位時間ステップ毎に、高分子材料モデル１６の伸長が計算される。そして、単位時間ステップを進展させることによって、予め定められた長さ（上限値）まで伸長した高分子材料モデル１６が計算される。この高分子材料モデル１６の伸長によって、分子鎖モデル５は、粒子モデル６、６間が離間する。この粒子モデル６、６間の離間に伴って、分子鎖モデル５が切断される状態が計算される。分子力学計算は、分子動力学計算と同様の上記ソフトウェアによって計算されうる。
【００５０】
量子力学計算では、分子鎖モデル５が、原子核と電子とに基づく相互作用に従うものとして、シミュレーションの単位時間ステップ毎に、粒子モデル６の電子状態と動きとが追跡される。従って、量子力学計算は、分子力学計算に比べて計算精度が高い。このような量子力学計算は、Ｇａｕｓｓｉａｎ社製の量子化学計算プログラムＧａｕｓｓｉａｎ０３、又は、Ｇａｕｓｓｉａｎ０９を用いて処理することができる。
【００５１】
本実施形態では、高分子材料モデル１６の全体領域を対象に、分子力学計算を用いて、伸長シミュレーションが実施される。そして、分子鎖モデル５のうち粒子モデル６、６間で切断しそうな部分のみを対象に、量子力学計算を用いて、分子鎖モデル５の切断が詳細に計算される。従って、本実施形態の切断計算工程Ｓ４では、例えば、高分子材料モデル１６の全体領域を対象に量子力学計算を実施する場合に比べて、計算時間を大幅に短縮することができる。図８は、切断計算工程Ｓ４の処理手順の一例を示すフローチャートである。
【００５２】
本実施形態の切断計算工程Ｓ４では、先ず、高分子材料モデル１６の伸長シミュレーションを実施するための初期条件が設定される（工程Ｓ４１）。図９は、高分子材料モデル１６の伸長シミュレーションを説明する側面図である。初期条件としては、高分子材料モデル１６のポアソン比に基づいて、各単位時間ステップでの高分子材料モデル１６の長さＬ１、Ｌ２、Ｌ３が定義される。これにより、切断計算工程Ｓ４では、高分子材料モデル１６のポアソン比に基づいて、高分子材料モデルを一軸方向（例えば、Ｚ軸方向）に伸長させるシミュレーションを実施することができる。
【００５３】
高分子材料モデル１６の一軸方向の伸長速度Ｖａについては、解析対象の高分子材料や、シミュレーション条件に応じて適宜設定することができる。本実施形態の伸長速度Ｖａは、例えば、２００ｍ／秒〜３００ｍ／秒に設定されている。これにより、本実施形態の切断計算工程Ｓ４では、高分子材料が高速に引き伸ばされた状態を計算することができる。これらの初期条件は、コンピュータ１に記憶される。
【００５４】
次に、本実施形態の切断計算工程Ｓ４では、分子鎖モデル５から選択された一対の粒子モデル６と、粒子モデル６、６間を連結するボンドモデル７とを含む分子モデルが設定される（工程Ｓ４２）。図１０は、分子モデル１１の一例を示す図である。分子モデルは、後述する結合解離ポテンシャルを求める際に用いられる。
【００５５】
本実施形態の分子モデル１１は、分子鎖モデル５において、一対の粒子モデル６、６の全ての組み合わせが設定される。本実施形態では、硫黄粒子モデル６Ｓ、６Ｓを結合した分子モデル１１Ａ、炭素粒子モデル６Ｃ、６Ｃを結合した分子モデル（図示省略）、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとを結合した分子モデル（図示省略）、硫黄粒子モデル６Ｓと炭素粒子モデル６Ｃとを連結した分子モデル（図示省略）、及び、硫黄粒子モデル６Ｓと水素粒子モデル６Ｈとを連結した分子モデル（図示省略）が設定される。図１０は、硫黄粒子モデル６Ｓ、６Ｓを結合した分子モデル１１Ａを代表して示している。
【００５６】
本実施形態の分子モデル１１は、一対の粒子モデル６、６（本実施形態では、硫黄粒子モデル６Ｓ、６Ｓ）と、粒子モデル６、６間を連結するボンドモデル７とを含んでいる。なお、各粒子モデル６には、粒子モデル６でモデル化された原子の原子価（例えば、硫黄原子の場合：２）から、粒子モデル６の結合数（例えば、１）を減じた数の水素粒子モデル６Ｈ（例えば、１個）が、ボンドモデル７を介して連結される。これにより、分子モデル１１は、量子力学計算において、安定した構造を維持することができる。このような方法に基づいて、上記した各分子モデル１１がそれぞれ設定される。これらの分子モデル１１は、コンピュータ１に記憶される。
【００５７】
次に、本実施形態の切断計算工程Ｓ４では、分子モデル１１を用いた量子力学計算に基づいて、エネルギーの増分が小さくなる変曲点を有する結合解離ポテンシャルが求められる（工程Ｓ４３）。結合解離ポテンシャルＰｕ（図示省略）は、結合する粒子モデル６の種類に応じてそれぞれ異なる。本実施形態では、工程Ｓ４３で求められた分子モデル１１を用いて、硫黄粒子モデル６Ｓ、６Ｓ間、炭素粒子モデル６Ｃ、６Ｃ間、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとの間、硫黄粒子モデル６Ｓと炭素粒子モデル６Ｃとの間、及び、硫黄粒子モデル６Ｓと水素粒子モデル６Ｈとの間について、それぞれ異なる結合解離ポテンシャルＰｕが求められる。
【００５８】
工程Ｓ４３では、先ず、各分子モデル１１の粒子モデル６、６間に、量子力学計算に使用する理論及び基底関数が定義される。本実施形態で定義される理論としては、例えば、密度汎関数法（汎関数：Ｂ３ＬＹＰ）が定義される。また、基底関数としては、例えば、基底関数（６−３１Ｇ（ｄ））が定義される。これらの量子力学計算に使用する理論及び基底関数は、上記した量子力学計算用のプログラムを用いることによって定義することができる。これらの理論及び基底関数は、コンピュータ１に記憶される。
【００５９】
次に、工程Ｓ４３では、分子モデル１１を用いた量子力学計算を行うシミュレーションによって、各粒子モデル６、６間の結合解離ポテンシャルＰｕが求められる。図１０に示されるように、分子モデル１１では、上記理論及び基底関数により、原子核と電子とに基づく相互作用のバランスによって、粒子モデル６、６の平衡距離ｒ0が決定されている。本実施形態の工程Ｓ４３では、例えば、粒子モデル６、６が当接した状態から、粒子モデル６、６間の距離ｒが平衡距離ｒ0よりも大（例えば、４．５Å）になるまでの間、粒子モデル６、６を単位ステップ毎に０．１Åずつ離間させる量子力学計算が実施される。また、本実施形態の量子力学計算では、非制限密度汎関数法に基づいて、分子モデル１１の全エネルギーが、単位ステップ毎に計算される。そして、粒子モデル６、６間の距離ｒと、分子モデル１１の全エネルギーとの関係に基づいて、結合解離ポテンシャルＰｕを求めることができる。図１１は、結合解離ポテンシャルＰｕの一例を示すグラフである。平衡距離ｒ0及び距離ｒは、粒子モデル６、６の中心６ｃ、６ｃ（図１０に示す）間の距離として定義される。
【００６０】
図１１のグラフにおいて、結合解離ポテンシャルＰｕは、粒子モデル６、６間の距離ｒが平衡距離ｒ0（本実施形態では、２．１）よりも小さくなるほど、分子モデル１１の全エネルギー（粒子モデル６、６間の斥力）が大きくなる。また、結合解離ポテンシャルＰｕは、距離ｒが２．１になるときにエネルギーが最小となり、粒子モデル６、６間に斥力や引力は作用しない。さらに、結合解離ポテンシャルＰｕは、距離ｒが平衡距離ｒ0よりも大きくなるほど、分子モデル１１の全エネルギー（粒子モデル６、６間の引力）が大きくなる。このように、結合解離ポテンシャルＰｕは、距離ｒに応じて、斥力及び引力を定義することができる。
【００６１】
結合解離ポテンシャルＰｕは、粒子モデル６、６間の距離ｒが平衡距離ｒ0よりも大きい領域において、分子モデル１１の全エネルギーの増分が小さくなる変曲点１３を有している。変曲点１３は、変曲点１３よりも距離ｒが大きいと予想される領域、及び、変曲点１３よりも距離ｒが小さいと予想される領域において、結合解離ポテンシャルＰｕに近似する一対の直線１５ａ、１５ｂの交点で特定することができる。
【００６２】
変曲点１３よりも距離ｒが小さい領域Ｔｂでは、変曲点１３よりも粒子モデル６、６間の距離ｒが大きい領域Ｔａに比べて、分子モデル１１の全エネルギーの増分が大きい。これは、変曲点１３よりも小さい領域Ｔｂにおいて、粒子モデル６、６間の結合が強固に維持されることを示している。また、変曲点１３よりも距離ｒが大きい領域Ｔａでは、粒子モデル６、６間で切断、又は、結合次数の減少が生じているため、分子モデル１１の全エネルギーの増分が小さくなっている。従って、変曲点１３は、粒子モデル６、６間の結合と解離との境界と考えることができる。
【００６３】
工程Ｓ４３では、上記した各分子モデル１１の全てについて、量子力学計算がそれぞれ実施される。これにより、硫黄粒子モデル６Ｓ、６Ｓ間の結合解離ポテンシャルＰｕ、炭素粒子モデル６Ｃ、６Ｃ間の結合解離ポテンシャルＰｕ、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとの間の結合解離ポテンシャルＰｕ、硫黄粒子モデル６Ｓと炭素粒子モデル６Ｃとの間の結合解離ポテンシャルＰｕ、及び、硫黄粒子モデル６Ｓと水素粒子モデル６Ｈとの間の結合解離ポテンシャルＰｕがそれぞれ求められる。各結合解離ポテンシャルＰｕは、コンピュータ１に記憶される。
【００６４】
次に、本実施形態の切断計算工程Ｓ４では、分子力学計算で用いられる結合ポテンシャルＰｎが定義される（工程Ｓ４４）。図４に示されるように、結合ポテンシャルＰｎは、分子鎖モデル５において、ボンドモデル７を介して隣り合う粒子モデル６、６間に定義されるものである。
【００６５】
本実施形態の結合ポテンシャルＰｎとしては、結合解離ポテンシャルＰｕ（図示省略）に近似するポテンシャルが定義される。このような結合ポテンシャルＰｎは、量子力学計算でのポテンシャルに近似するため、分子力学計算の結果と量子力学計算の結果との差異を、最小限にすることができる。図１２は、結合解離ポテンシャルＰｕ及び結合ポテンシャルＰｎの一例を示すグラフである。
【００６６】
本実施形態の工程Ｓ４４は、フックの法則に基づいた調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させている。これにより、結合ポテンシャルＰｎが求められる。調和振動子型の関数は、下記式（５）で定義される。
【００６７】
【数５】
2018032077.tif 000007
ここで、
Ｖ：分子モデルのエネルギー
ｋ：バネ定数
ｒ：粒子モデル間の距離
ｒ0：粒子モデル間の平衡距離
【００６８】
上記式（５）の調和振動子型の関数は、粒子モデル６、６間の距離ｒと、分子モデルのエネルギーとの関係を示す二次関数である。本実施形態では、変曲点１３よりも距離ｒが小さい領域Ｔｂ（図１１に示す）において、調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させることによって、結合ポテンシャルＰｎを求めている。従って、本実施形態の結合ポテンシャルＰｎは、粒子モデル６、６間で切断、又は、結合次数の減少が生じる前の結合解離ポテンシャルＰｕに近似するものである。
【００６９】
調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させる方法としては、例えば、最小二乗法を用いることができる。このような結合ポテンシャルＰｎでは、複雑な結合解離ポテンシャルＰｕを、簡単な調和振動子型の関数で表すことができる。
【００７０】
工程Ｓ４４では、結合解離ポテンシャルＰｕのエネルギー（分子モデル１１のエネルギー）の下限値（即ち、平衡距離ｒ0）から、エネルギーの上限値の５０％〜８０％までの範囲において、調和振動子型の関数を、結合解離ポテンシャルＰｕに近似させるのが望ましい。これにより、結合ポテンシャルＰｎは、変曲点１３よりも距離ｒ0が小さい領域Ｔｂ（図１１に示す）において、結合解離ポテンシャルＰｕに精度良く近似することができる。
【００７１】
なお、調和振動子型の関数に近似させる範囲において、結合解離ポテンシャルＰｕのエネルギーの上限値の８０％を超えると、変曲点１３よりも距離が大きい領域Ｔａ（図１１に示す）を含めて、調和振動子型の関数を近似させるおそれがある。逆に、結合解離ポテンシャルＰｕのエネルギーの上限値の５０％未満であると、結合解離ポテンシャルＰｕに近似させる範囲が小さくなり、結合解離ポテンシャルＰｕに精度良く近似できないおそれがある。このような観点より、調和振動子型の関数に近似させる範囲において、より好ましくは、結合解離ポテンシャルＰｕのエネルギーの上限値の７０％以下であり、より好ましくは６０％以上である。
【００７２】
工程Ｓ４４では、硫黄粒子モデル６Ｓ、６Ｓ間の結合解離ポテンシャルＰｕ、炭素粒子モデル６Ｃ、６Ｃ間の結合解離ポテンシャルＰｕ、炭素粒子モデル６Ｃと水素粒子モデル６Ｈとの間の結合解離ポテンシャルＰｕ、硫黄粒子モデル６Ｓと炭素粒子モデル６Ｃとの間の結合解離ポテンシャルＰｕ、及び、硫黄粒子モデル６Ｓと水素粒子モデル６Ｈとの間の結合解離ポテンシャルＰｕについて、結合ポテンシャルＰｎがそれぞれ求められる。
【００７３】
工程Ｓ４４では、図４に示した分子鎖モデル５において、ボンドモデル７を介して隣り合う粒子モデル６、６間に、第１ポテンシャルＰ１に代えて、結合ポテンシャルＰｎがそれぞれ定義される。なお、ボンドモデル７を介さずに隣り合う粒子モデル６、６間に定義されている第２ポテンシャルＰ２については、維持される。
【００７４】
次に、本実施形態の切断計算工程Ｓ４では、高分子材料モデル１６（図７に示す）を伸長させる変形計算が実施される（伸長工程Ｓ４５）。本実施形態では、初期条件として設定されたポアソン比に基づいて、高分子材料モデル１６の変形が計算される。
【００７５】
伸長工程Ｓ４５では、現在の単位時間ステップにおいて、一軸方向に伸長した高分子材料モデル１６が計算される。この高分子材料モデル１６の伸長により、図７に示した分子鎖モデル５も一軸方向に引き伸ばされ、ボンドモデル７を介して隣り合う粒子モデル６、６間の距離ｒ（図４に示す）の少なくとも一部が大きくなる。
【００７６】
次に、本実施形態の切断計算工程Ｓ４では、伸長工程Ｓ４５で変形した高分子材料モデル１６に対して、分子力学計算（Molecular Mechanics ： ＭＭ）が実施される（第１分子力学計算工程Ｓ４６）。第１分子力学計算工程Ｓ４６では、高分子材料モデル１６のエネルギーＥａと、高分子材料モデル１６に作用する力Ｆａとを含む物理量が計算される。このようなエネルギーＥａと、力Ｆａとを含む物理量は、コンピュータ１に記憶される。
【００７７】
次に、本実施形態の切断計算工程Ｓ４では、第１分子力学計算工程Ｓ４６後の分子鎖モデル５について、ボンドモデル７を介して隣り合う粒子モデル６、６間の距離ｒ（図４に示す）が、予め定められた第１距離Ｌａよりも大きい粒子モデル対（離反粒子モデル部）２１が存在するか否かが判断される（工程Ｓ４７）。
【００７８】
工程Ｓ４７では、セル１７に配置されている全ての分子鎖モデル５について、ボンドモデル７を介して隣り合う粒子モデル６、６間の距離ｒが計算される。なお、工程Ｓ４７では、ボンドモデル７で拘束されていない粒子モデル６、６間の距離は計算されない。そして、粒子モデル６、６間の距離ｒが、第１距離Ｌａよりも大きい一対の粒子モデル６、６及びボンドモデル７を、粒子モデル対２１として判断される。
【００７９】
本実施形態の第１距離Ｌａは、図１２に示した結合解離ポテンシャルＰｕの変曲点１３での粒子モデル６、６間の距離ｒ（Ｌｎ）よりも小に設定されている。また、第１距離Ｌａは、粒子モデル６、６の平衡距離ｒ0よりも大に設定されている。上述したように、変曲点１３は、粒子モデル６、６間の結合と解離との境界点と考えることができる。このため、第１距離Ｌａよりも大きく離間した粒子モデル対２１では、粒子モデル６、６間で切断する可能性が高いと考えることができる。従って、工程Ｓ４７では、粒子モデル６、６間で切断する可能性が高い粒子モデル対２１（図１４に示す）の有無を判断することができる。
【００８０】
また、第１距離Ｌａについては、変曲点１３での距離ｒよりも小さい値であれば、適宜設定することができる。本実施形態では、図１２に示されるように、結合解離ポテンシャルＰｕと、結合ポテンシャルＰｎとの交点１８での距離ｒを、第１距離Ｌａとして求めている。これにより、量子力学計算よりも計算時間を短縮しうる分子力学計算に基づいて、粒子モデル対２１の存在の有無を、精度良く判断することができる。
【００８１】
工程Ｓ４７において、粒子モデル対２１が存在すると判断された場合（工程Ｓ４７で、「Ｙ」）、粒子モデル対２１のみを対象に、次の量子力学計算工程Ｓ４８が実施される。他方、粒子モデル対２１が存在しないと判断された場合（工程Ｓ４７で、「Ｎ」）、後述する工程Ｓ５２が実施される。
【００８２】
次に、本実施形態の切断計算工程Ｓ４では、粒子モデル対２１を含む第１領域２２を対象に、量子力学計算が行われる（量子力学計算工程Ｓ４８）。図１３は、量子力学計算工程Ｓ４８の処理手順の一例を示すフローチャートである。
【００８３】
本実施形態の量子力学計算工程Ｓ４８では、先ず、粒子モデル対２１を含む第１領域２２が設定される（工程Ｓ４８１）。図１４は、第１領域２２の一例を示す図である。第１領域２２は、高分子材料モデル１６の分子鎖モデル５から、粒子モデル対２１を抜き出した非全体領域である。なお、高分子材料モデル１６に、複数の粒子モデル対２１が存在する場合、全ての粒子モデル対２１が、第１領域２２に含められる。
【００８４】
また、粒子モデル対２１の各粒子モデル６には、粒子モデル６でモデル化された原子の原子価（例えば、炭素原子の場合：４）から、粒子モデル６の結合数（例えば、１）を減じた数の水素粒子モデル６Ｈ（例えば、３個）が、ボンドモデル７を介して連結される。これにより、量子力学計算において、安定した構造を有する粒子モデル対２１を設定することができる。
【００８５】
次に、本実施形態の量子力学計算工程Ｓ４８では、粒子モデル対２１の粒子モデル６に、量子力学計算に使用する理論及び基底関数が定義される（工程Ｓ４８２）。工程Ｓ４８２では、粒子モデル対２１の粒子モデル６、６間に設定されている結合ポテンシャルＰｎを無効にして、量子力学計算に使用する理論及び基底関数が定義される。なお、量子力学計算に使用する理論及び基底関数については、上述のとおりである。
【００８６】
次に、本実施形態の量子力学計算工程Ｓ４８では、初期条件として設定されたポアソン比に基づいて、第１領域２２のみを対象に、量子力学計算が実施される（工程Ｓ４８３）。上述したように、量子力学計算では、原子核と電子とに基づく相互作用に従うものとして、シミュレーションの単位時間ステップ毎に、粒子モデル６の電子状態と動きとが追跡される。従って、工程Ｓ４８３では、分子力学計算に比べて、第１領域２２及び粒子モデル対２１の伸長を、高い精度で計算することができる。また、工程Ｓ４８３では、高分子材料モデル１６から抜き出された第１領域２２のみを対象に、量子力学計算が実施されている。このため、図７に示した高分子材料モデル１６の全体領域を対象に量子力学計算が実施される場合に比べて、計算時間を大幅に短縮することができる。
【００８７】
工程Ｓ４８３では、第１領域２２のエネルギーＥｃと、第１領域２２に作用する力Ｆｃとを含む物理量が計算される。このようなエネルギーＥｃ及び力Ｆｃを含む物理量は、コンピュータ１に記憶される。
【００８８】
次に、本実施形態の切断計算工程Ｓ４では、第１領域２２を対象に分子力学計算が行われる（第２分子力学計算工程Ｓ４９）。図１４に示されるように、第１領域２２は、伸長工程Ｓ４５で変形した高分子材料モデル１６の分子鎖モデル５から、粒子モデル対２１を抜き出した非全体領域である。従って、第２分子力学計算工程Ｓ４９は、第１領域２２のエネルギーＥｂと、第１領域２２に作用する力Ｆｂとを含む物理量が計算される。このようなエネルギーＥｂと力Ｆｂとを含む物理量は、後述するＯＮＩＯＭ法に従って、高分子材料モデル１６のエネルギーＥｔ、及び、高分子材料モデル１６に作用する力Ｆｔを求めるのに利用される。エネルギーＥｂと力Ｆｂとを含む物理量は、コンピュータ１に記憶される。また、第２分子力学計算工程Ｓ４９では、粒子モデル対２１の粒子モデル６、６間に、結合ポテンシャルＰｎが設定される。
【００８９】
次に、本実施形態の切断計算工程Ｓ４では、高分子材料モデル１６のエネルギーＥｔが計算される（工程Ｓ５０）。本実施形態のエネルギーＥｔは、ＯＮＩＯＭ（ Our own N-layered Integrated molecular Orbital and molecular Mechanics ） 法に従って、下記式（６）に基づいて計算される。
Ｅｔ＝Ｅａ−Ｅｂ＋Ｅｃ  …  （６）
ここで、
Ｅａ：第１分子力学計算工程で求められた高分子材料モデルのエネルギー
Ｅｂ：第２分子力学計算工程で求められた第１領域のエネルギー
Ｅｃ：量子力学計算工程で求められた第１領域のエネルギー
【００９０】
上記式（６）において、第１分子力学計算工程Ｓ４６で求められたエネルギーＥａは、高分子材料モデル１６の全体領域のエネルギーである。また、第２分子力学計算工程Ｓ４９で求められたエネルギーＥｂは、第１領域２２でのエネルギーである。これらのエネルギーＥａ、Ｅｂは、分子力学計算で求められているため、量子力学計算に比べて、計算精度が低い。量子力学計算工程Ｓ４８で求められたエネルギーＥｃは、第１領域２２のエネルギーである。このエネルギーＥｃは、量子力学計算で求められているため、分子力学計算に比べて、計算精度が高い。
【００９１】
ＯＮＩＯＭ法は、計算精度の異なる階層構造を組み合わせることによって最適解を得る方法である。本実施形態では、計算精度が比較的低い高分子材料モデルのエネルギーＥａのうち、第１領域のエネルギーＥｂを、計算精度が比較的高い第１領域のエネルギーＥｃに置き換えることによって、高分子材料モデル１６のエネルギーＥｔを求めている。このようなエネルギーＥｔは、分子鎖モデル５の切断による影響を詳細に考慮することができるため、計算精度を高めることができる。また、本実施形態では、例えば、量子力学計算に基づいた高分子材料モデル１６の全体領域のエネルギーを求めることなく、高分子材料モデル１６のエネルギーＥｔを求めることができるため、計算時間を短縮することができる。エネルギーＥｔは、コンピュータ１に記憶される。
【００９２】
次に、本実施形態の切断計算工程Ｓ４では、高分子材料モデル１６に作用する力Ｆｔが計算される（工程Ｓ５１）。本実施形態の力Ｆｔは、高分子材料モデルのエネルギーＥｔと同様に、ＯＮＩＯＭ法に従って、下記式（７）に基づいて計算される。
Ｆｔ＝Ｆａ−Ｆｂ＋Ｆｃ  …  （７）
ここで、
Ｆａ：第１分子力学計算工程で求められた高分子材料モデルに作用する力
Ｆｂ：第２分子力学計算工程で求められた第１領域に作用する力
Ｆｃ：量子力学計算工程で求められた第１領域に作用する力
【００９３】
本実施形態では、計算精度が比較的低い高分子材料モデル１６の力Ｆａのうち、第１領域２２の力Ｆｂを、計算精度が比較的高い第１領域の力Ｆｃに置き換えることによって、高分子材料モデル１６の力Ｆｔを求めている。このような力Ｆｔは、分子鎖モデル５の切断による影響を詳細に考慮することができるため、計算精度を高めることができる。また、本実施形態では、量子力学計算に基づいた高分子材料モデル１６の全体領域に作用する力を計算する必要がないため、計算時間を短縮することができる。この力Ｆｔは、コンピュータ１に記憶される。
【００９４】
次に、本実施形態の切断計算工程Ｓ４では、現在の単位時間ステップにおいて、分子鎖モデル５の各粒子モデル６の座標値及び速度が更新される（工程Ｓ５２）。本実施形態では、ＯＮＩＯＭ法で計算されたエネルギーＥｔ及び力Ｆｔを含む物理量から加速度を求める。そして、エネルギーＥｔ及び力Ｆｔを含む物理量から求められた加速度に基づいて、第１分子力学計算工程Ｓ４６が実施される前の分子鎖モデル５の各粒子モデル６の座標値から、各粒子モデル６の新たな座標値（即ち、次の単位時間ステップで用いられる座標値）が計算される。
【００９５】
このような各粒子モデル６の座標値は、切断される可能性の高い第１領域２２を量子力学計算した座標値が含められているため、分子鎖モデル５の切断を精度よく解析することができる。また、切断される可能性の低い部分については、比較的精度の低い分子力学計算が実施されるため、計算時間を短縮することができる。
【００９６】
本実施形態では、第１分子力学計算工程Ｓ４６、及び、第２分子力学計算工程Ｓ４９において、結合解離ポテンシャルＰｕに近似する結合ポテンシャルＰｎが定義されている。このため、第２分子力学計算工程Ｓ４９後の各粒子モデル６に作用する力と、量子力学計算工程Ｓ４８後の各粒子モデル６に作用する力との間に生じる差異を最小限に抑えることができる。従って、ＯＮＩＯＭ法に基づく、座標値の計算精度を高めることができる。また、各粒子モデル６の速度は、分子鎖モデル５の各粒子モデル６の座標値、高分子材料モデル１６のエネルギーＥｔ、及び、高分子材料モデル１６に作用する力Ｆｔに基づいて計算される。
【００９７】
なお、上述した工程Ｓ４７において、粒子モデル対２１が存在しないと判断された場合（工程Ｓ４７で、「Ｎ」）は、第１領域２２を対象とした量子力学計算（工程Ｓ４８）及び分子力学計算（工程Ｓ４９）が実施されていない。この場合、工程Ｓ５２では、高分子材料モデル１６を対象とした分子力学計算によって求められた粒子モデル６の座標値、第１分子力学計算工程で求められたエネルギーＥａ、及び、第１分子力学計算工程で求められた力Ｆａから、座標値及び速度が更新される。
【００９８】
次に、切断計算工程Ｓ４では、量子力学計算工程Ｓ４８（本実施形態では量子力学計算工程Ｓ４８及び第２分子力学計算工程Ｓ４９）が行われた後の第１領域２２について、粒子モデル６、６間の距離ｒが、予め定められた第２距離Ｌｂよりも大きい粒子モデル対２１が存在するか否かが判断される（工程Ｓ５３）。工程Ｓ５３では、工程Ｓ５２で更新された分子鎖モデル５の各粒子モデル６の座標値に基づいて、第１領域２２に配置されている全ての粒子モデル対２１について、粒子モデル６、６間の距離ｒが計算される。
【００９９】
本実施形態の第２距離Ｌｂは、粒子モデル６、６間で切断される粒子モデル対２１の有無を判断するためのパラメータである。図１２に示されるように、本実施形態の第２距離Ｌｂは、第１距離Ｌａよりも大きく設定されている。さらに、第２距離Ｌｂは、変曲点１３での粒子モデル間の距離Ｌｎ以上に設定されている。上述したように、変曲点１３は、粒子モデル６、６（図４に示す）間の結合と解離との境界点と考えることができる。このため、粒子モデル６、６間の距離ｒが、第２距離Ｌｂよりも大きい場合、粒子モデル６、６間で切断すると判断できる。
【０１００】
図１５（ａ）は、粒子モデル６、６間で切断する前の分子鎖モデル５の一部を示す図である。図１５（ｂ）は、粒子モデル６、６間で切断した分子鎖モデル５の一部を示す図である。後述の切断工程Ｓ５４では、切断対象の粒子モデル対２１に対応する分子鎖モデル５の粒子モデル６、６間の結合が切断される。切断の定義は、先ず、粒子モデル６、６間を拘束しているボンドモデル７（図１５（ａ）に示す）を無効に設定する。次に、切断された粒子モデル６、６間に定義されていた結合ポテンシャルＰｎ（図１５（ａ）に示す）が無効とされ、図１５（ｂ）に示されるように、ボンドモデル７を介さずに隣り合う粒子モデル６、６間に、第２ポテンシャルＰ２が定義される。
【０１０１】
さらに、切断された各粒子モデル６、６には、例えば、粒子モデル６でモデル化された原子の原子価（例えば、炭素原子の場合：４）から、粒子モデル６の結合数（例えば、３）を減じた数の水素粒子モデル６Ｈ（例えば、１個）が、ボンドモデル７を介して連結される。
【０１０２】
このような切断の定義により、分子鎖モデル５の粒子モデル６、６間の切断を計算することができるが、シミュレーションが異常終了しやすいという問題がある。発明者らは、鋭意研究を重ねた結果、粒子モデル６、６（例えば、炭素粒子モデル６Ｃ、６Ｃ）間の距離を十分に確保せずに切断すると、切断された粒子モデル６、６間において、量子力学計算で求められる第１領域２２のエネルギーと、分子力学計算で求められる第１領域２２のエネルギーとの差異が大きくなり、異常終了が発生しやすくなることを知見した。
【０１０３】
量子力学計算でのエネルギーと、分子力学計算でのエネルギーとは互いに異なる。また、量子力学計算でのエネルギーと、分子力学計算でのエネルギーとの差異は、粒子モデル６、６間の距離が近いほど大きくなる一方、粒子モデル６、６間の距離が大きいほど小さくなる。発明者らは、種々検証した結果、切断された各粒子モデル６、６に新たに連結された水素粒子モデル６Ｈ、６Ｈ間の距離が平衡距離ｒ0よりも小さくなり、水素粒子モデル６Ｈ、６Ｈ間の第２ポテンシャルＰ２の斥力（即ち、上記式（４）の１２乗の項）が過度に大きくなることで異常終了が発生すると解明した。
【０１０４】
このような検証結果に基づいて、図１２に示されるように、本実施形態の第２距離Ｌｂは、変曲点１３での粒子モデル６、６間の距離Ｌｎに１．０〜２．５Åを加えた値に設定されている。これにより、後述の切断工程Ｓ５４において、切断対象の粒子モデル６、６間の距離を十分に確保したのちに、粒子モデル６、６間で切断されるため、新たに連結された水素粒子モデル６Ｈ、６Ｈ（図１５（ｂ）に示す）も大きく離間させることができる。従って、本実施形態では、前記エネルギーの差異の増大に起因する計算の異常終了を確実に防ぐことができるため、分子鎖モデル５の切断を安定して計算することができる。
【０１０５】
なお、第２距離Ｌｂにおいて、変曲点１３での前記距離Ｌｎに加える値が１．０Å未満であると、量子力学計算と分子力学計算とのエネルギーの差異の増大を防げないおそれがある。逆に、変曲点１３での前記距離Ｌｎに加える値が２．５Åを超えると、量子力学計算が不安定になって異常終了しやすくなる。さらに、第１領域２２に含まれる粒子モデル対が切断されるまでの時間が長くなるため、量子力学計算での計算負荷が増大しやすくなる。このような観点より、変曲点１３での前記距離Ｌｎに加える値は、好ましくは１．５以上であり、また、好ましくは、２．０以下である。
【０１０６】
工程Ｓ５３において、粒子モデル６、６間の距離ｒが第２距離Ｌｂ以上の粒子モデル対２１が存在すると判断された場合（工程Ｓ５３で、「Ｙ」）、次の切断工程Ｓ５４が実施される。他方、全ての粒子モデル対２１において、粒子モデル６、６間の距離が第２距離Ｌｂ未満であると判断された場合（工程Ｓ５３で、「Ｎ」）、高分子材料モデル１６が予め定められた条件まで伸長したか否かを判断する後述の工程Ｓ５５が実施される。
【０１０７】
次に、本実施形態の切断計算工程Ｓ４では、切断対象の粒子モデル対２１に対応する分子鎖モデル５の粒子モデル６、６間の結合が切断される（切断工程Ｓ５４）。粒子モデル６、６間の切断の定義については、上述の手順で行われる。切断工程Ｓ５４では、工程Ｓ５２で更新された座標値に基から計算された粒子モデル６、６間の距離に基づいて、粒子モデル６、６間の結合が切断される。更新された座標値は、量子力学計算で計算された座標値が反映されているため、分子鎖モデル５の切断を精度よく解析することができる。
【０１０８】
このように、本実施形態のシミュレーション方法では、粒子モデル対２１を含んだ第１領域２２を対象として、分子力学計算よりも計算精度の高い量子力学計算の結果に基づき、粒子モデル６、６間を切断している。従って、本実施形態では、分子鎖モデル５の切断を、実際の分子鎖の切断に近似させることができ、高分子材料を精度良く解析するのに役立つ。
【０１０９】
また、量子力学計算が実施される粒子モデル対２１の有無は、工程Ｓ４７において、量子力学計算よりも計算時間が短い分子力学計算の結果に基づいて判断される。このため、高分子材料モデル１６の全体領域を対象とした量子力学計算を実施する必要がない。従って、本実施形態のシミュレーション方法では、計算時間を短縮することができる。さらに、本実施形態では、第２距離Ｌｂが上記範囲に設定されるため、計算の異常終了を防ぐことができ、分子鎖モデル５の切断を安定して計算することができる。
【０１１０】
次に、本実施形態の切断計算工程Ｓ４では、高分子材料モデル１６が、予め定められた条件まで伸長したか否かが判断される（工程Ｓ５５）。工程Ｓ５５では、高分子材料モデル１６が上限値まで伸長したと判断された場合（工程Ｓ５５で、「Ｙ」）、次の工程Ｓ５が実施される。他方、高分子材料モデル１６が上限値まで伸長していないと判断された場合は（工程Ｓ５５で、「Ｎ」）、単位時間ステップを一つ進展させて（工程Ｓ５６）、工程Ｓ４５〜工程Ｓ５５が再度実施される。これにより、切断計算工程Ｓ４では、高分子材料モデルを上限値まで伸長させて、分子鎖モデル５の切断を解析することができる。
【０１１１】
次に、図３に示されるように、本実施形態のシミュレーション方法では、高分子材料モデル１６の物理量が、許容範囲内であるか否かが判断される（工程Ｓ５）。この工程Ｓ５では、高分子材料モデル１６の物理量が許容範囲内であると判断された場合（工程Ｓ５で、「Ｙ」）、高分子材料モデル１６に基づいて、高分子材料が製造される（工程Ｓ６）。他方、高分子材料モデル１６の物理量が、許容範囲外と判断された場合は（工程Ｓ５で、「Ｎ」）、分子鎖モデル５の諸条件が変更されて（工程Ｓ７）、工程Ｓ１〜工程Ｓ５が再度実施される。このように、本実施形態では、高分子材料モデル１６の物理量が許容範囲内になるまで、分子鎖モデル５の諸条件が変更されるため、所望の性能を有する高分子材料を、効率よく設計することができる。
【０１１２】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【０１１３】
図３及び図５に示した処理手順に従って、分子鎖モデル及び高分子モデルが定義され、分子鎖モデルを切断する切断計算工程が実施された（実施例、比較例）。
【０１１４】
実施例及び比較例の切断計算工程では、図８及び図１３に示した処理手順に従って、高分子材料モデルに基づいて分子力学計算を行う工程、粒子モデル間の距離が第１距離Ｌａよりも  大きい離反粒子モデル部を含む第１領域を対象に、量子力学計算を行う工程、及び、第１領域を対象に分子力学計算を行う工程が実施され、上記式（６）、（７）に基づいて、高分子材料モデルのエネルギーＥｔ及び高分子材料モデルに作用する力Ｆｔが計算された。
【０１１５】
実施例及び比較例の分子力学計算を行う工程では、分子鎖モデルの粒子モデル間に、結合解離ポテンシャルに近似する結合ポテンシャルが定義された。さらに、実施例及び比較例の切断工程では、粒子モデル間の距離が第２距離Ｌｂよりも大きいときに、粒子モデル対に対応する前記分子鎖モデルの粒子モデル間の結合が切断された。そして、高分子材料モデルが下記のセル変形率の上限値になるまでの計算時間が測定された。
【０１１６】
実施例１〜４の第２距離Ｌｂは、結合解離ポテンシャルの変曲点での粒子モデル間の距離に１．０〜２．５Åを加えた値に設定された。他方、比較例１の第２距離は、結合解離ポテンシャルの変曲点での粒子モデル間の距離に０．５Åを加えた値に設定された。また、比較例２の第２距離は、結合解離ポテンシャルの変曲点での粒子モデル間の距離に３．０Åを加えた値に設定された。
【０１１７】
なお、粒子モデル間に設定されるポテンシャルやソフトウェアは、明細書に記載のとおりである。共通仕様は、次のとおりである。
コンピュータ：ＳＧＩ社のワークステーション
ＣＰＵのコア数：１２コア
搭載メモリ：６４ＧＢ
高分子モデル：
セルの１辺の長さＬ１：２．３nm
分子鎖モデルの個数：１個
分子鎖モデル：
ポリイソプレン（イソプレン分子：１００重合）に硫黄（１原子）を架橋
粒子モデル（原子）の個数：１３０１個
変曲点での距離Ｌｎ：３．０９Å
第１距離Ｌａ：２．６９Å
伸長シミュレーション：
伸長方向：一軸方向（Ｚ軸方向）
伸長速度Ｖａ：２５０ｍ／sec
セル変形率の上限値：１１００％
テスト結果を、表１に示す。
【０１１８】
【表１】
2018032077.tif 000008
【０１１９】
テストの結果、実施例では、セル変形率が１１００％になるまで伸長することができた。他方比較例では、セル変形率が７００％の時点で異常終了した。また、実施例のうち、結合解離ポテンシャルの変曲点での粒子モデル間の距離に１．０〜２．０Åを加えた実施例２及び実施例３が、実施例１及び実施例２に比べて、計算精度を維持しつつ、計算時間の増大を防ぐことができた。
【符号の説明】
【０１２０】
Ｓ４６      第１分子力学計算工程
Ｓ４８      量子力学計算を行う工程
Ｓ５４      分子鎖モデルの粒子モデル間の結合を切断する工程
【図１】
2018032077.tif 000009
【図２】
2018032077.tif 000010
【図３】
2018032077.tif 000011
【図４】
2018032077.tif 000012
【図５】
2018032077.tif 000013
【図６】
2018032077.tif 000014
【図７】
2018032077.tif 000015
【図８】
2018032077.tif 000016
【図９】
2018032077.tif 000017
【図１０】
2018032077.tif 000018
【図１１】
2018032077.tif 000019
【図１２】
2018032077.tif 000020
【図１３】
2018032077.tif 000021
【図１４】
2018032077.tif 000022
【図１５】
2018032077.tif 000023