(57)【特許請求の範囲】
【請求項１】
コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、
前記コンピュータが、前記高分子鎖の粗視化モデルを用いて運動計算を行い、前記粗視化モデルの拡散の大きさを示す平均二乗変位を計算する第１計算工程、
前記コンピュータが、前記高分子鎖の全原子モデルを用いて運動計算を行い、前記全原子モデルの拡散の大きさを示す平均二乗変位を計算する第２計算工程、及び
前記コンピュータが、前記粗視化モデルの平均二乗変位と、前記全原子モデルの平均二乗変位とに基づいて、前記粗視化モデルの時間単位を、前記高分子鎖の時間単位に換算する工程を含むことを特徴とする高分子材料のシミュレーション方法。
【請求項２】
前記第１計算工程は、前記粗視化モデルの初期配置を決定する工程と、初期配置された前記粗視化モデルの緩和状態を計算する工程とを含み、
前記第２計算工程は、前記全原子モデルの構造緩和を計算することなく、前記第１計算工程で得られた緩和状態の前記粗視化モデルの配置に基づいて、緩和状態の前記全原子モデルを作成する工程を含む請求項１記載の高分子材料のシミュレーション方法。
【請求項３】
前記粗視化モデルは、前記高分子鎖のモノマーに対応する複数のビーズを含み、
前記緩和状態の全原子モデルを作成する工程は、前記コンピュータが、緩和状態の前記粗視化モデルの前記各ビーズに、前記モノマーを表すモノマーモデルを割り当てる請求項２記載の高分子材料のシミュレーション方法。
【請求項４】
前記粗視化モデルは、前記高分子鎖のモノマーに対応する複数のビーズを含み、
前記コンピュータが、前記全原子モデルのRouseのビーズ数を計算する第３計算工程、
前記コンピュータが、前記粗視化モデルのRouseのビーズ数を計算する第４計算工程、
及び
前記コンピュータが、前記全原子モデルのRouseのビーズ数と、前記粗視化モデルのRouseのビーズ数とに基づいて、前記粗視化モデルのビーズ数と、前記高分子鎖のモノマー数との比を計算する工程をさらに含み、
前記第３計算工程は、前記コンピュータに、モノマー数が異なる複数の全原子モデルを設定する工程と、
前記コンピュータが、前記モノマー数が異なる全原子モデル毎に、前記Rouseのビーズ数を計算する工程とを含み、
前記第４計算工程は、前記コンピュータに、ビーズ数が異なる複数の粗視化モデルを設定する工程と、
前記コンピュータが、前記ビーズ数が異なる粗視化モデル毎に、前記Rouseのビーズ数を計算する工程とを含む請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
【請求項５】
前記Rouseのビーズ数ＮRouseは、下記式（１）で計算される請求項４に記載の高分子材料のシミュレーション方法。
【数１】
0006050903.tif 000005
ここで、
Ｌ：粗視化モデル又は全原子モデルの全長
Ｒg：粗視化モデル又は全原子モデルの慣性半径
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粗視化モデルの時間単位を、高分子鎖の時間単位に換算することができる高分子材料のシミュレーション方法に関する。
【背景技術】
【０００２】
近年、ゴム配合の開発のために、高分子材料の反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。この種のシミュレーション方法では、例えば、高分子鎖を、複数のビーズでモデル化した粗視化モデルが、コンピュータに入力される。そして、粗視化モデルが予め定められた空間に配置され、コンピュータによる分子動力学計算が行われる。
【０００３】
粗視化モデルを用いた分子動力学計算では、現実の高分子鎖とは異なる時間単位が用いられている。このため、下記特許文献１では、粗視化モデルの時間単位を、高分子鎖の時間単位に換算する方法が提案されている。
【０００４】
下記特許文献１の換算方法では、先ず、粗視化モデル及び全原子モデルの末端間ベクトルの長さの時間相関関数がそれぞれ求められる。次に、粗視化モデル及び全原子モデルについて、末端間ベクトルの相関がなくなる時間として定義される時間のパラメータがそれぞれ求められる。そして、粗視化モデルの時間のパラメータと、全原子モデルの時間のパラメータとをフィッティグさせることにより、粗視化モデルの時間単位を、高分子鎖の時間単位に換算している。
【先行技術文献】
【特許文献】
【０００５】
【特許文献１】特開２０１４−２０６９１３号公報
【発明の概要】
【発明が解決しようとする課題】
【０００６】
上記特許文献１の換算方法では、全原子モデルの末端間ベクトルの長さの時間相関関数を求める際に、全原子モデルの配置が十分に緩和するまで平衡計算が行われる。この平衡計算では、時間単位を精度よく換算するために、全原子モデルの鎖長として、高分子鎖の物性を示す現実的な長さを用いると、非常に多くの計算時間を要する。従って、上記特許文献１の換算方法に基づいて、粗視化モデルの時間単位を高分子鎖の時間単位に換算することは、現実的ではない。
【０００７】
本発明は、以上のような実状に鑑み案出されたもので、粗視化モデルの時間単位を高分子鎖の時間単位に、短時間で換算することができる高分子材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【０００８】
本発明は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、前記コンピュータが、前記高分子鎖の粗視化モデルを用いて運動計算を行い、前記粗視化モデルの拡散の大きさを示す平均二乗変位を計算する第１計算工程、前記コンピュータが、前記高分子鎖の全原子モデルを用いて運動計算を行い、前記全原子モデルの拡散の大きさを示す平均二乗変位を計算する第２計算工程、及び前記コンピュータが、前記粗視化モデルの平均二乗変位と、前記全原子モデルの平均二乗変位とに基づいて、前記粗視化モデルの時間単位を、前記高分子鎖の時間単位に換算する工程を含むことを特徴とする。
【０００９】
本発明に係る前記高分子材料のシミュレーション方法において、前記第１計算工程は、前記粗視化モデルの初期配置を決定する工程と、初期配置された前記粗視化モデルの緩和状態を計算する工程とを含み、前記第２計算工程は、前記全原子モデルの構造緩和を計算することなく、前記第１計算工程で得られた緩和状態の前記粗視化モデルの配置に基づいて、緩和状態の前記全原子モデルを作成する工程を含んでもよい。
【００１０】
本発明に係る前記高分子材料のシミュレーション方法において、前記粗視化モデルは、前記高分子鎖のモノマーに対応する複数のビーズを含み、前記緩和状態の全原子モデルを作成する工程は、前記コンピュータが、緩和状態の前記粗視化モデルの前記各ビーズに、前記モノマーを表すモノマーモデルを割り当ててもよい。
【００１１】
本発明に係る前記高分子材料のシミュレーション方法において、前記粗視化モデルは、前記高分子鎖のモノマーに対応する複数のビーズを含み、前記コンピュータが、前記全原子モデルのRouseのビーズ数を計算する第３計算工程、前記コンピュータが、前記粗視化モデルのRouseのビーズ数を計算する第４計算工程、及び前記コンピュータが、前記全原子モデルのRouseのビーズ数と、前記粗視化モデルのRouseのビーズ数とに基づいて、前記粗視化モデルのビーズ数と、前記高分子鎖のモノマー数との比を計算する工程をさらに含み、前記第３計算工程は、前記コンピュータに、モノマー数が異なる複数の全原子モデルを設定する工程と、前記コンピュータが、前記モノマー数が異なる全原子モデル毎に、前記Rouseのビーズ数を計算する工程とを含み、前記第４計算工程は、前記コンピュータに、ビーズ数が異なる複数の粗視化モデルを設定する工程と、前記コンピュータが、前記ビーズ数が異なる粗視化モデル毎に、前記Rouseのビーズ数を計算する工程とを含でもよい。
【００１４】
本発明に係る前記高分子材料のシミュレーション方法において、前記Rouseのビーズ数ＮRouseは、下記式（１）で計算されてもよい。
【数１】
0006050903.tif 000002
ここで、
Ｌ：粗視化モデル又は全原子モデルの全長
Ｒg：粗視化モデル又は全原子モデルの慣性半径
【発明の効果】
【００１５】
本発明の高分子材料のシミュレーション方法は、コンピュータが、高分子鎖の粗視化モデルを用いて運動計算を行い、高分子鎖の拡散の大きさを示す平均二乗変位を計算する第１計算工程、高分子鎖の全原子モデルを用いて運動計算を行い、全原子モデルの拡散の大きさを示す平均二乗変位を求める第２計算工程、及び、粗視化モデルの平均二乗変位と、全原子モデルの平均二乗変位とに基づいて、粗視化モデルの時間単位を、高分子鎖の時間単位に換算する工程を含んでいる。
【００１６】
平均二乗変位は、任意の時間幅において、その時間幅での粒子の平均的な拡散の大きさ（即ち、移動距離の二乗平均）を示すものである。このような平均二乗変位を求める運動計算に基づいて、時間単位を換算するためには、全原子モデルの拡散挙動と、粗視化モデルの拡散挙動とが類似する転移域（即ち、分子鎖同士の絡み合いが起こるまで）の時間幅の計算で十分であるため、例えば、全原子モデルの構造が十分に緩和するまで（即ち、全原子モデルの分子鎖同士の絡み合いがほどけるまで）緩和計算を行う必要がない。従って、本発明の高分子材料のシミュレーション方法では、粗視化モデルの時間単位を高分子鎖の時間単位に、短時間で換算することができる。
【図面の簡単な説明】
【００１７】
【図１】本発明のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。
【図２】ポリブタジエンの構造式である。
【図３】全原子モデルの一例を示す概念図である。
【図４】粗視化モデルの一例を示す概念図である。
【図５】本実施形態の高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。
【図６】本実施形態のモノマー数換算工程の処理手順の一例を示すフローチャートである。
【図７】本実施形態の第３計算工程の処理手順の一例を示すフローチャートである。
【図８】全原子モデルが配置された高分子材料モデルの一例を示す概念図である。
【図９】全原子モデルのポテンシャルを説明する概念図である。
【図１０】全原子モデルのRouseのビーズ数とモノマー数との関係を示すグラフである。
【図１１】本実施形態の第４計算工程の処理手順の一例を示すフローチャートである。
【図１２】粗視化モデルが配置された高分子材料モデルの概念図である。
【図１３】粗視化モデルのポテンシャルを説明する概念図である。
【図１４】粗視化モデルのRouseのビーズ数とビーズ数との関係を示すグラフである。
【図１５】全原子モデルのRouseのビーズ数と粗視化モデルのRouseのビーズ数とをフィッティングしたグラフである。
【図１６】本実施形態の時間換算工程の処理手順の一例を示すフローチャートである。
【図１７】本実施形態の第１計算工程の処理手順の一例を示すフローチャートである。
【図１８】粗視化モデルの平均二乗変位と時間幅との関係を示すグラフである。
【図１９】本実施形態の第２計算工程の処理手順の一例を示すフローチャートである。
【図２０】全原子モデルの平均二乗変位と時間幅との関係を示すグラフである。
【図２１】全原子モデルと粗視化モデルとをフィッティングした平均二乗変位と時間幅との関係を示すグラフである。
【発明を実施するための形態】
【００１８】
以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、高分子鎖を有する高分子材料を解析するための方法である。高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。
【００１９】
図１は、本発明のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含んで構成されている。この本体１ａには、例えば、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置１ａ１、１ａ２が設けられている。また、記憶装置には、本実施形態のシミュレーション方法を実行するためのソフトウェア等が予め記憶されている。
【００２０】
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、高分子材料として、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある。）が例示される。図２は、ポリブタジエンの構造式である。このポリブタジエンを構成する高分子鎖２は、メチレン基（−ＣＨ２−）とメチン基（−ＣＨ−）とからなるモノマー５｛−［ＣＨ２−ＣＨ＝ＣＨ−ＣＨ２］−｝が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基（−ＣＨ２）に替えて、メチル基（−ＣＨ３）が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。
【００２１】
本実施形態のシミュレーション方法では、高分子材料の高分子鎖２をモデル化した全原子モデル、及び、粗視化モデルが用いられる。図３は、全原子モデルの一例を示す概念図である。図４は、粗視化モデルの一例を示す概念図である。
【００２２】
図３に示されるように、全原子モデル３は、図２に示した高分子鎖２の実際の構造に基づいて、原子をモデル化した粒子モデル７で表現したものである。図４に示されるように、粗視化モデル４は、図２に示した高分子鎖２のモノマー５に対応する複数のビーズ（粒子）１６を含んでいる。各ビーズ１６は、直径を持った球として表現される。
【００２３】
全原子モデル３（図３に示す）を用いた分子動力学計算では、現実の高分子鎖２（図２に示す）に基づいた時間単位が用いられている。他方、粗視化モデル４（図４に示す）を用いた分子動力学計算では、現実の高分子鎖２とは異なる時間単位が用いられている。このため、粗視化モデル４の分子動力学計算の計算結果を、実際の高分子鎖２の運動として取り扱うためには、粗視化モデル４の時間単位を、高分子鎖２の時間単位に換算する必要がある。
【００２４】
また、粗視化モデル４（図４に示す）の時間単位を、高分子鎖２（図２に示す）の時間単位に換算するには、粗視化モデル４の時間のパラメータと、全原子モデル３（図３に示す）の時間のパラメータとを精度よくフィッティングすることが重要である。このため、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比を予め計算し、このビーズ数とモノマー数との比に基づいて全原子モデル３及び粗視化モデル４が設定して、全原子モデル３を用いた計算結果と、粗視化モデル４を用いた計算結果とを精度よく対応させるのが望ましい。
【００２５】
このような観点より、本実施形態のシミュレーション方法では、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算するのに先立ち、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比（即ち、粗視化モデル４の１個のビーズ１６あたりの高分子鎖２のモノマー数）を計算している。図５は、本実施形態の高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。
【００２６】
図５に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を計算する（モノマー数換算工程Ｓ１）。
【００２７】
本実施形態のモノマー数換算工程Ｓ１では、全原子モデル３（図３に示す）のRouseのビーズ数と、粗視化モデル４（図４に示す）のRouseのビーズ数とに基づいて、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比を計算している。ここで、Rouseのビーズ数ＮRouseとは、高分子物理で知られているRouseモデルのパラメータの一つであり、下記式（１）で計算される。
【数１】
0006050903.tif 000003
ここで、
Ｌ：粗視化モデル又は全原子モデルの全長
Ｒg：粗視化モデル又は全原子モデルの慣性半径
【００２８】
上記式（１）では、例えば、全原子モデル３の末端間ベクトルＶ１（図３に示す）の長さのアンサンブル平均を用いる上記特許文献１とは異なり、慣性半径Ｒgを用いて、Rouseのビーズ数が求められる。これにより、本実施形態では、全原子モデル３の末端側の粒子モデル７の不安定な運動性の影響を小さくできるため、Rouseのビーズ数を精度よく求めることができる。
【００２９】
上記式（１）において、Rouseのビーズ数の値が大きいほど、全原子モデル３（図３に示す）又は粗視化モデル４（図４に示す）は、コンパクトに折り畳んだ形状となる。他方、Rouseのビーズ数の値が小さいほど、全原子モデル３又は粗視化モデル４の折れ曲がりが小さくなり、直線状に連続してのびる形状となる。従って、Rouseのビーズ数は、全原子モデル３又は粗視化モデル４の形状を示すパラメータとして用いられる。
【００３０】
全原子モデル３（図３に示す）のRouseのビーズ数と、粗視化モデル４（図４に示す）のRouseのビーズ数とをフィッティングさせることにより、互いに同一形状の全原子モデル３と粗視化モデル４とを対応づけて、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比を精度よく計算しうることが知られている（例えば、上記特許文献１参照）。このような知見に基づいて、本実施形態のシミュレーション方法では、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比（即ち、粗視化モデル４の１個のビーズ１６あたりの高分子鎖２のモノマー数）を求めている。図６は、本実施形態のモノマー数換算工程Ｓ１の処理手順の一例を示すフローチャートである。
【００３１】
本実施形態のモノマー数換算工程Ｓ１では、先ず、コンピュータ１が、全原子モデル３（図３に示す）のRouseのビーズ数を計算する（第３計算工程Ｓ１１）。図７は、本実施形態の第３計算工程Ｓ１１の処理手順の一例を示すフローチャートである。
【００３２】
本実施形態の第３計算工程Ｓ１１では、先ず、コンピュータ１に、モノマー数が異なる複数の全原子モデル３が設定される（工程Ｓ１１１）。図３に示されるように、全原子モデル３は、複数の粒子モデル７と、粒子モデル７、７間を結合するボンド８とを含んで構成されている。
【００３３】
粒子モデル７及びボンド８は、図２に示した高分子鎖２のモノマー５を表す単位構造に基づいて連結されることにより、モノマーモデル１１が設定される。このモノマーモデル１１が、予め定められたモノマー数（即ち、分子量（重合度））に基づいて連結されることにより、全原子モデル３が設定される。本実施形態では、予め定められた複数のモノマー数に基づいて、モノマーモデル１１がそれぞれ連結される。これにより、モノマー数が異なる複数の全原子モデル３が設定される。
【００３４】
モノマー数については、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ１の性能等に基づいて、緩和計算が現実的な計算時間で完了しうる範囲内で設定されるのが望ましい。なお、モノマー数は、計算精度を維持するために、極端に小さい値を除外するのが望ましい。モノマー数の一例としては、１０〜６０から選択されうる。なお、このモノマー数に限定されるわけではない。
【００３５】
粒子モデル７は、後述の分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル７には、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の粒子モデル７は、高分子鎖２の炭素原子をモデル化した炭素粒子モデル７Ｃ、及び、高分子鎖２の水素原子をモデル化した水素粒子モデル７Ｈを含んでいる。
【００３６】
ボンド８は、粒子モデル７、７間を拘束するものである。本実施形態のボンド８は、炭素粒子モデル７Ｃ、７Ｃを連結する主鎖８Ａ、及び、炭素粒子モデル７Ｃと水素粒子モデル７Ｈとの間を連結する側鎖８Ｂとを含んでいる。これらの主鎖８Ａ及び側鎖８Ｂは、例えば、平衡長とバネ定数とが定義されたバネとして取り扱われる。
【００３７】
全原子モデル３は、各粒子モデル７、７間の結合長さである結合長、ボンド８を介して連続する３つの粒子モデル７がなす角度である結合角、及び、ボンド８を介して連続する４つの粒子モデル７において、隣り合う３つの粒子モデル７が作る二面角などが定義される。これにより、全原子モデル３は、三次元構造を有する。全原子モデル３は、慣例に従い、外力又は内力を受けることによって、結合長、結合角及び二面角が変化する。これにより、全原子モデル３は、その三次元構造を変化させることができる。
【００３８】
結合長、結合角及び二面角は、例えば、論文１（J. Comput. Chem. 25, 1157-1174 (2004)）に基づいて設定されるポテンシャル（ＧＡＦＦ）によって定義されうる。ポテンシャルは、高分子鎖２の構造に応じて設定されるのが望ましい。このような全原子モデル３は、材料物性シミュレーションソフトウェア（例えば、（株）ＪＳＯＬ社製のＪ−ＯＣＴＡ）を用いて作成することができる。各全原子モデル３は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。
【００３９】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、モノマー数が異なる各全原子モデル３の初期配置を決定する（工程Ｓ１１２）。工程Ｓ１１２では、予め定められた空間に、全原子モデル３が配置された高分子材料モデルが設定される。図８は、全原子モデル３が配置された高分子材料モデル１０の一例を示す概念図である。
【００４０】
本実施形態の工程Ｓ１１２では、モノマー数が異なる各全原子モデル３が、独立して設けられた空間９にそれぞれ配置される。これにより、工程Ｓ１１２では、各全原子モデル３の初期配置が決定された高分子材料モデル１０がそれぞれ定義される。各高分子材料モデル１０には、同一のモノマー数の全原子モデル３が複数本配置されている。全原子モデル３は、モンテカルロ法に基づいて、空間９内に配置されるのが望ましい。
【００４１】
空間９は、解析対象の高分子材料の微小構造部分に相当する。本実施形態の空間９は、互いに向き合う三対の平面９Ｓ、９Ｓを有する立方体として定義されている。各平面９Ｓには、周期境界条件が定義されている。これにより、空間９では、例えば、一方の平面９Ｓａから出て行った全原子モデル３の一部が、反対側の平面９Ｓｂから入ってくるように計算することができる。従って、一方の平面９Ｓａと、反対側の平面９Ｓｂとが連続している（繋がっている）ものとして取り扱うことができる。
【００４２】
各空間９に配置される全原子モデル３の本数については、モノマー数や、空間９の大きさに基づいて適宜設定される。なお、全原子モデル３の本数が少ないと、後述の分子動力学計算において、一方の平面９Ｓａから出て行きかつ反対側の平面９Ｓｂから入ってきた全原子モデル３の一端と、平面９Ｓｂ側に配置されている全原子モデル３の他端とが絡まりやすくなり、計算落ちするおそれがある。さらに、全原子モデル３の慣性半径を精度良く計算するためには、平衡時において全原子モデル同士の絡まりを可能な限り防ぐ必要がある。このため、全原子モデル３の本数としては、平衡時の周期境界長が平衡時の慣性半径の３倍以上に長くなりうる本数が選択されるのが望ましい。逆に、全原子モデル３の本数が多くても、運動方程式の質点として取り扱われる粒子モデル７の粒子数が増大し、多くの計算時間を要するおそれがある。このような観点より、空間９に配置される全原子モデル３の本数は、好ましくは２０本以上であり、また、好ましくは２００本以下である。
【００４３】
空間９の一辺の長さＬａは、適宜設定することができる。本実施形態の長さＬａは、系内の粒子モデル７の初期密度が、例えば０．００１g／cm3となるように設定されるのが望ましい。また、空間９の大きさは、例えば１atmで安定な体積に設定されるのが望ましい。このような空間９では、後述の分子動力学計算において、全原子モデル３の回転運動をスムーズに計算することができる。
【００４４】
図９は、全原子モデル３のポテンシャルを説明する概念図である。工程Ｓ１１２では、全原子モデル３の初期配置が決定された各高分子材料モデル１０において、隣接する全原子モデル３、３の粒子モデル７、７間に、相互作用ポテンシャルＰ１が定義される。本実施形態の相互作用ポテンシャルＰ１は、ＬＪポテンシャルＵLJ（ｒij）であり、下記式（２）で定義される。
【００４５】
【数２】
0006050903.tif 000004
ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒij：粒子モデル間の距離
ｒc：カットオフ距離
ε：粒子モデル間に定義されるＬＪポテンシャルの強度
σ：粒子モデルの直径に相当
なお、距離ｒij及びカットオフ距離ｒcは、各粒子モデル７、７の中心間の距離として定義される。
【００４６】
相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒijがσよりも小さくなるほど、粒子モデル７、７間に作用する斥力が大きくなる。また、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒijがσになるときに最小となり、粒子モデル７、７間に斥力や引力は働かない。さらに、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒijがσよりも大になるほど、粒子モデル７、７間に作用する引力が働く。このように、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒijに応じて、斥力及び引力を定義することができる。
【００４７】
相互作用ポテンシャルＰ１は、炭素粒子モデル７Ｃ、７Ｃ間に設定される第１ポテンシャルＰ１ａ、水素粒子モデル７Ｈ、７Ｈ間に設定される第２ポテンシャルＰ１ｂ、及び、炭素粒子モデル７Ｃと水素粒子モデル７Ｈとの間に設定される第３ポテンシャルＰ１ｃを含んでいる。なお、上記式（２）中の各定数は、上記論文１に基づいて、適宜設定することができる。これらの全原子モデル３の初期配置が決定された高分子材料モデル１０は、コンピュータ１に入力される。
【００４８】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、初期配置された全原子モデル３の緩和状態を計算する（工程Ｓ１１３）。工程Ｓ１１３では、全原子モデル３の初期配置が決定された各高分子材料モデル１０（図８に示す）について、分子動力学計算が実施される。分子動力学計算は、例えば、空間９について所定の時間、配置した全ての全原子モデル３が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全ての粒子モデル７の動きが追跡され、コンピュータ１に記憶される。また、分子動力学計算の条件は、例えば、系内の粒子モデル７の個数、体積及び温度は一定で行われる。このような分子動力学計算は、例えば、分子動力学計算プログラムLAMMPSを用いて行うことができる。
【００４９】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、全原子モデル３の初期配置を十分に緩和できたか否かを判断する（工程Ｓ１１４）。なお、初期配置の緩和の判断基準については、全原子モデル３の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。工程Ｓ１１４では、全原子モデル３の初期配置を十分に緩和できたと判断された場合（工程Ｓ１１４で、「Ｙ」）、次の工程Ｓ１１５が実施される。他方、全原子モデル３の初期配置を十分に緩和できていないと判断された場合（工程Ｓ１１４で、「Ｎ」）は、単位時間を（例えば、工程Ｓ１１３で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ１１６）、工程Ｓ１１３（分子動力学計算）及び工程Ｓ１１４が再度実施される。これにより、全原子モデル３の初期配置を、確実に緩和することができる。
【００５０】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、全ての全原子モデル（即ち、モノマー数が異なる全ての全原子モデル）３の緩和状態が計算されたか否かを判断する（工程Ｓ１１５）。工程Ｓ１１５では、全ての全原子モデル３の緩和計算が計算されたと判断された場合（工程Ｓ１１５で、「Ｙ」）、次の工程Ｓ１１７が実施される。他方、全ての全原子モデル３の緩和計算が終了していないと判断された場合（工程Ｓ１１５で、「Ｎ」）、緩和状態が計算されていない全原子モデル３（即ち、図８に示した高分子材料モデル１０）が選択され（工程Ｓ１１８）、工程Ｓ１１３〜工程Ｓ１１６が実施される。これにより、全ての全原子モデル３の平衡状態（緩和状態）が計算される。
【００５１】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、全原子モデル３の全長Ｌを計算する（工程Ｓ１１７）。工程Ｓ１１７では、上記特許文献１の方法に基づいて、全原子モデル３を主鎖方向に強制的に引き伸ばす変形計算が実施される。これにより、全原子モデル３の全長Ｌを求めることができる。全原子モデル３の全長Ｌは、コンピュータ１に入力される。
【００５２】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、全原子モデル３の慣性半径Ｒgを計算する（工程Ｓ１１９）。全原子モデル３の慣性半径は、例えば、論文２（ Kurt Kremer & Gary S. Grest 著 「Dynamics of entangled linear polymer melts: A molecular-dynamics simulation」、J. Chem Phys. vol.92, No.8, 15 April 1990 ）に記載されている式２．５に基づいて、工程Ｓ１１３の分子動力学計算で求められた全原子モデル３の粒子モデル７の座標値から計算される。全原子モデル３の慣性半径は、コンピュータ１に入力される。
【００５３】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、全原子モデル３のRouseのビーズ数を計算する（工程Ｓ１２０）。上述したように、Rouseのビーズ数は、全原子モデル３の全長Ｌと、全原子モデル３の慣性半径Ｒgとに基づいて、上記式（１）で計算される。全原子モデル３のRouseのビーズ数は、コンピュータ１に入力される。
【００５４】
次に、本実施形態の第３計算工程Ｓ１１では、コンピュータ１が、モノマー数が異なる全ての全原子モデル３について、Rouseのビーズ数が計算されたか否かを判断する（工程Ｓ１２１）。工程Ｓ１２１において、全ての全原子モデル３のRouseのビーズ数が計算されたと判断された場合（工程Ｓ１２１において、「Ｙ」）、次の第４計算工程Ｓ１２が実施される。他方、全ての全原子モデル３のRouseのビーズ数が計算されていないと判断された場合（工程Ｓ１２１において、「Ｎ」）、Rouseのビーズ数が計算されていない他の全原子モデル３が選択され（工程Ｓ１２２）、工程Ｓ１１７、及び、工程Ｓ１１９〜工程Ｓ１２１が実施される。これにより、第３計算工程Ｓ１１では、モノマー数が異なる全原子モデル３毎に、Rouseのビーズ数が計算される。図１０は、全原子モデル３のRouseのビーズ数とモノマー数との関係を示すグラフである。図１０の例では、モノマー数が２０、３０及び４０の全原子モデル３のRouseのビーズ数が示されている。全原子モデル３のRouseのビーズ数は、コンピュータ１に入力される。
【００５５】
次に、本実施形態のモノマー数換算工程Ｓ１では、コンピュータ１が、粗視化モデル４のRouseのビーズ数を計算する（第４計算工程Ｓ１２）。図１１は、本実施形態の第４計算工程Ｓ１２の処理手順の一例を示すフローチャートである。
【００５６】
本実施形態の第４計算工程Ｓ１２は、先ず、コンピュータ１に、ビーズ数（鎖長）が異なる複数の粗視化モデル４が設定される（工程Ｓ１３１）。図４に示されるように、粗視化モデル４は、複数のビーズ１６と、ビーズ１６、１６間を結合する結合鎖１７とを含んで構成されている。
【００５７】
本実施形態では、予め定められたビーズ数に基づいて、ビーズ１６が連結されることにより、粗視化モデル４が設定される。本実施形態では、複数の異なるビーズ数が設定されている。これらのビーズ数に基づいて、ビーズ１６がそれぞれ連結されることにより、ビーズ数が異なる複数の粗視化モデル４が設定される。ビーズ数については、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ１の性能等に基づいて、適宜設定することができる。ビーズ数は、例えば、５〜１０００から選択されるのが望ましい。なお、このようなビーズ数に限定されるわけではない。
【００５８】
ビーズ１６は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ビーズ１６には、質量、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ１に記憶される。
【００５９】
結合鎖１７は、ビーズ１６、１６間に平衡長を定義した結合ポテンシャルとして構成される。ここで、「平衡長」とは、ビーズ１６、１６間の結合距離である。この結合距離が変化した場合は、結合鎖１７によって、元の平衡長が保持される。これにより、粗視化モデル４は、直鎖状の三次元構造を維持することができる。なお、平衡長は、隣り合うビーズ１６、１６の中心１６ａ、１６ａ間の距離として定義される。また、結合鎖１７の結合ポテンシャルには、例えば、上記論文２に基づいて設定されるのが望ましい。このような粗視化モデル４は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ１に入力される。
【００６０】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、ビーズ数が異なる各粗視化モデル４の初期配置を決定する（工程Ｓ１３２）。工程Ｓ１３２では、予め定められた空間に、粗視化モデル４が配置された高分子材料モデルが設定される。図１２は、粗視化モデル４が配置された高分子材料モデル１８の一例を示す概念図である。
【００６１】
本実施形態の工程Ｓ１３２では、ビーズ数が異なる粗視化モデル４が、独立して設けられた空間９にそれぞれ配置される。これにより、工程Ｓ１３２では、各粗視化モデル４の初期配置が決定された高分子材料モデル１８がそれぞれ定義される。各高分子材料モデル１８には、同一のビーズ数の全原子モデル３が複数本配置されている。粗視化モデル４は、モンテカルロ法に基づいて、空間９内に配置されるのが望ましい。
【００６２】
空間９は、解析対象の高分子材料の微小構造部分に相当し、図８に示した空間９と同様に設定される。各空間９に配置される粗視化モデル４の本数については、ビーズ数や、空間９の大きさに基づいて適宜設定される。なお、粗視化モデル４の本数が少ないと、全原子モデル３（図８に示す）と同様に、粗視化モデル４の一端と他端とが絡まりやすくなり、計算落ちするおそれがある。さらに、粗視化モデル４の慣性半径を精度良く計算するためには、平衡時において粗視化モデル４同士の絡まりを可能な限り防ぐ必要がある。このため、粗視化モデル４の本数としては、平衡時の周期境界長が平衡時の慣性半径の３倍以上に長くなりうる本数が選択されるのが望ましい。逆に、粗視化モデル４の本数が多くても、運動方程式の質点として取り扱われるビーズ１６の粒子数が増大し、多くの計算時間を要するおそれがある。このような観点より、空間９に配置される粗視化モデル４の本数は、好ましくは５０本以上であり、また、好ましくは２０００本以下である。
【００６３】
図１３は、粗視化モデル４のポテンシャルを説明する概念図である。工程Ｓ１３２では、各高分子材料モデル１８において、隣接する粗視化モデル４、４のビーズ１６、１６間に、相互作用ポテンシャルＰ２がそれぞれ定義される。本実施形態の相互作用ポテンシャルＰ２は、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（２）で定義される。なお、上記式（２）において、「粒子モデル間」は、「ビーズ間」に置き換えて適用される。
【００６４】
相互作用ポテンシャルＰ２の強度ε、相互作用ポテンシャルＰ２が作用する距離σ、及び、カットオフ距離ｒcは、例えば、上記論文１に基づいて設定されるのが望ましい。これらの粗視化モデル４の初期配置が決定された高分子材料モデル１８は、コンピュータ１に入力される。
【００６５】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、初期配置された粗視化モデル４の緩和状態を計算する（工程Ｓ１３３）。工程Ｓ１３３では、図１２に示した空間９に初期配置された粗視化モデル４について、分子動力学計算が実施される。分子動力学計算では、例えば、空間９について所定の時間、配置した全ての粗視化モデル４が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全てのビーズ１６の動きが追跡され、コンピュータ１に記憶される。また、分子動力学計算の条件は、例えば、系内のビーズ１６の個数、体積及び温度は一定で行われる。
【００６６】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、粗視化モデル４の初期配置を十分に緩和できたか否かを判断する（工程Ｓ１３４）。なお、初期配置の緩和の判断基準については、粗視化モデル４の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。工程Ｓ１３４では、粗視化モデル４の初期配置を十分に緩和できたと判断された場合（工程Ｓ１３４で、「Ｙ」）、次の工程Ｓ１３５が実施される。他方、粗視化モデル４の初期配置を十分に緩和できていないと判断された場合（工程Ｓ１３４で、「Ｎ」）、単位時間を（例えば、工程Ｓ１３３で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ１３６）、工程Ｓ１３３（分子動力学計算）及び工程Ｓ１３４が再度実施される。これにより、粗視化モデル４の初期配置を、確実に緩和することができる。
【００６７】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、全ての粗視化モデル（即ち、ビーズ数が異なる全ての粗視化モデル）４の緩和状態が計算されたか否かを判断する（工程Ｓ１３５）。工程Ｓ１３５では、全ての粗視化モデル４の緩和状態が計算されたと判断された場合（工程Ｓ１３５で、「Ｙ」）、次の工程Ｓ１３７が実施される。他方、全ての粗視化モデル４の緩和状態が計算されていないと判断された場合（工程Ｓ１３５で、「Ｎ」）、緩和状態が計算されていない粗視化モデル４（即ち、図１２に示した高分子材料モデル１８）が選択され（工程Ｓ１３８）、工程Ｓ１３３〜工程Ｓ１３６が実施される。これにより、全ての粗視化モデル４の平衡状態（緩和状態）が計算される。
【００６８】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、粗視化モデル４の全長Ｌを計算する（工程Ｓ１３７）。粗視化モデル４の全長Ｌは、粗視化モデル４内の隣接ビーズの平衡核間距離と、ビーズ数とを乗じることによって求めることができる。なお、平衡核間距離は、上記の分子動力学計算の結果に基づいて求められる。粗視化モデル４の全長Ｌは、コンピュータ１に入力される。
【００６９】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、粗視化モデル４の慣性半径Ｒgを計算する（工程Ｓ１３９）。粗視化モデル４の慣性半径は、例えば、上記論文２に記載されている式２．５に基づいて、工程Ｓ１３３の分子動力学計算で求められた粗視化モデル４のビーズ１６の座標値から計算される。粗視化モデル４の慣性半径は、コンピュータ１に入力される。
【００７０】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、粗視化モデルのRouseのビーズ数を計算する（工程Ｓ１４０）。上述したように、Rouseのビーズ数は、粗視化モデル４の全長Ｌと、粗視化モデル４の慣性半径Ｒgとに基づいて、上記式（１）で計算される。粗視化モデル４のRouseのビーズ数は、コンピュータ１に入力される。
【００７１】
次に、本実施形態の第４計算工程Ｓ１２では、コンピュータ１が、ビーズ数が異なる全ての粗視化モデル４について、Rouseのビーズ数が計算されたか否かを判断する（工程Ｓ１４１）。工程Ｓ１４１において、全ての粗視化モデル４のRouseのビーズ数が計算されたと判断された場合（工程Ｓ１４１において、「Ｙ」）、次の工程Ｓ１３（図６に示す）が実施される。他方、全ての粗視化モデル４のRouseのビーズ数が計算されていないと判断された場合（工程Ｓ１４１において、「Ｎ」）、Rouseのビーズ数が計算されていない粗視化モデル４が選択され（工程Ｓ１４２）、工程Ｓ１３７、工程Ｓ１３９〜工程Ｓ１４１が実施される。これにより、第４計算工程Ｓ１２では、ビーズ数が異なる粗視化モデル４毎に、Rouseのビーズ数が計算される。図１４は、粗視化モデル４のRouseのビーズ数とビーズ数との関係を示すグラフである。図１４の例では、ビーズ数（鎖長）が１０、１３、２０、２７、５０、８０、１００、２００及び５００の粗視化モデル４のRouseのビーズ数が示されている。粗視化モデル４のRouseのビーズ数は、コンピュータ１に入力される。
【００７２】
次に、図６に示されるように、本実施形態のモノマー数換算工程Ｓ１では、コンピュータ１が、粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を計算する（工程Ｓ１３）。工程Ｓ１３では、全原子モデル３（図３に示す）のRouseのビーズ数と、粗視化モデル４のRouseのビーズ数とに基づいて、粗視化モデル４の１個のビーズ１６あたりの高分子鎖２のモノマー数を求めている。
【００７３】
工程Ｓ１３では、先ず、統計誤差を考慮した非線形最小二乗法に基づいて、図１０に示した全原子モデル３のRouseのビーズ数とモノマー数との関係と、図１４に示した粗視化モデル４のRouseのビーズ数とビーズ数との関係とがフィッティングされる。図１５は、全原子モデル３のRouseのビーズ数と、粗視化モデル４のRouseのビーズ数とをフィッティングしたグラフである。
【００７４】
上記フィッティングは、先ず、粗視化モデル４の結果（図１４に示す）が、図１５に示した滑らかな曲線（破線で示している）でフィッティングされる。次に、この曲線と全原子モデル３の結果（図１０に示す）とが重なるように、粗視化モデル４のビーズ数と全原子モデル３のモノマー数との比についてフィッティングされる。これにより、全原子モデル３のRouseのビーズ数と、粗視化モデル４のRouseのビーズ数とがフィッティングされる。フィッティングの手法としては、例えば、論文３（K. Levenberg, Quarterly of Applied Mathematics 2, 164-168 (1944), D. Marquardt, SIAM Journal on Applied Mathematics 11(2), 431-441 (1963)）に記載のレーベンバーグ・マーカート法（Levenberg-Marquardt Method）を用いることができる。このレーベンバーグ・マーカート法を用いたフィッティングは、例えば、グラフ作成ソフトウェア（例えば、gnuplot）により、容易に行うことができる。
【００７５】
フィッティングに用いる統計誤差としては、全原子モデル３の慣性半径の統計誤差、及び、粗視化モデル４の慣性半径の統計誤差が計算される。全原子モデル３の慣性半径の統計誤差は、モノマー数が異なる全原子モデル３毎に、分子鎖毎の慣性半径の時間平均についての分散を、空間９に配置された全原子モデルの本数で除した値を求め、その値の平方根から求めることができる。また、粗視化モデル４の慣性半径の統計誤差は、ビーズ数が異なる粗視化モデル４毎に、分子鎖毎の慣性半径の時間平均についての分散を、空間９に配置された粗視化モデルの本数で除した値を求め、その値の平方根から求めることができる。
【００７６】
工程Ｓ１３では、図１５に示したグラフに基づいて、粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を計算している。なお、図１５に示したグラフの例において、粗視化モデルのビーズ数と、高分子鎖２のモノマー数との比（即ち、１個のビーズ１６あたりのモノマー数）は、１．５である。このような粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比は、コンピュータ１に入力される。
【００７７】
このように、本実施形態のモノマー数換算工程Ｓ１では、全原子モデル３（図３に示す）のRouseのビーズ数と、粗視化モデル４（図４に示す）のRouseのビーズ数とをフィッティングさせることにより、互いに略同一形状となる全原子モデル３と粗視化モデル４とを対応づけることができるため、粗視化モデル４のビーズ数と、高分子鎖２のモノマー数との比を精度よく計算することができる。さらに、本実施形態では、全原子モデル３のRouseのビーズ数、及び、粗視化モデル４のRouseのビーズ数がそれぞれ複数求められているため、精度良くフィッティングすることができる。
【００７８】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、粗視化モデル４（図４に示す）の時間単位を、高分子鎖２（図２に示す）の時間単位に換算する（時間換算工程Ｓ２）。本実施形態の時間換算工程Ｓ２では、粗視化モデル４の平均二乗変位と、全原子モデル３（図３に示す）の平均二乗変位とに基づいて、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算している。平均二乗変位は、任意の時間幅において、その時間幅での粗視化モデル４のビーズ１６又は全原子モデル３の粒子モデル７の拡散の大きさ（即ち、移動距離の二乗平均）を示すものである。従って、平均二乗変位は、時間のパラメータを含んでいる。また、本実施形態では、高分子材料のガラス域とゴム域との間の転移域について、粗視化モデル４の平均二乗変位、及び、全原子モデル３の平均二乗変位が計算される。図１６は、本実施形態の時間換算工程Ｓ２の処理手順の一例を示すフローチャートである。
【００７９】
本実施形態の時間換算工程Ｓ２では、先ず、コンピュータ１が、粗視化モデル４（図４に示す）の拡散の大きさを示す平均二乗変位を計算する（第１計算工程Ｓ２１）。図１７は、本実施形態の第１計算工程Ｓ２１の処理手順の一例を示すフローチャートである。
【００８０】
本実施形態の第１計算工程Ｓ２１では、先ず、コンピュータ１が、粗視化モデル４のビーズ数（鎖長）を決定する（工程Ｓ２１１）。工程Ｓ２１１では、現実の高分子鎖の物性を示すビーズ数が指定される。なお、ビーズ数（鎖長）が小さいと、粗視化モデル４同士の絡み合いが起きる前に流動してしまい、高分子鎖の物性を十分に示せないおそれがある。逆に、ビーズ数（鎖長）が大きいと、後述する第２計算工程Ｓ２２（図１６に示す）において、全原子モデル３の慣性半径の３倍程度の平均周期境界長とするために、全原子モデル３の本数を多くする必要がある。これにより、運動方程式の質点として取り扱われる粒子モデル７（図９に示す）が多くなり、計算時間が増加するおそれがある。例えば、論文２に基づいて設定された粗視化モデル４では、ビーズ数（鎖長）が８０〜２００程度に設定されるのが望ましい。
【００８１】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、決定されたビーズ数に基づいて、粗視化モデル４（図４に示す）を設定する（工程Ｓ２１２）。工程Ｓ２１２では、図４に示されるように、決定されたビーズ数分のビーズ１６が、結合鎖１７を介して連結される。これにより、決定されたビーズ数の粗視化モデル４が設定される。粗視化モデル４は、コンピュータ１に入力される。
【００８２】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、粗視化モデル４の初期配置を決定する（工程Ｓ２１３）。工程Ｓ２１３では、図１２に示されるように、予め定められた空間９に粗視化モデル４が配置された高分子材料モデル２１が設定される。空間９には、複数本の粗視化モデル４が配置される。粗視化モデル４は、モンテカルロ法に基づいて、空間９内に配置されるのが望ましい。なお、各空間９に配置される粗視化モデル４の本数や、相互ポテンシャル等については、上述したモノマー数換算工程Ｓ１と同一手順で設定される。粗視化モデル４の初期配置が決定された高分子材料モデル２１は、コンピュータ１に入力される。
【００８３】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、初期配置された粗視化モデル４の緩和状態を計算する（工程Ｓ２１４）。工程Ｓ２１４では、空間９に初期配置された粗視化モデル４について、分子動力学計算が実施される。分子動力学計算は、上述したモノマー数換算工程Ｓ１と同一手順で実施される。
【００８４】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、粗視化モデル４の初期配置を十分に緩和できたか否かを判断する（工程Ｓ２１５）。工程Ｓ２１５では、粗視化モデル４の初期配置を十分に緩和できたと判断された場合（工程Ｓ２１５で、「Ｙ」）、次の工程Ｓ２１６が実施される。他方、粗視化モデル４の初期配置を十分に緩和できていないと判断された場合（工程Ｓ２１５で、「Ｎ」）、単位時間を（例えば、工程Ｓ２１４で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ２１７）、工程Ｓ２１４（分子動力学計算）及び工程Ｓ２１５が再度実施される。これにより、粗視化モデル４の初期配置を、確実に緩和することができる。初期配置が緩和された粗視化モデル４は、コンピュータ１に入力される。
【００８５】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、初期配置が緩和された粗視化モデル４について、運動計算を行う（工程Ｓ２１６）。本実施形態の運動計算としては、上述した分子動力学計算が実施される。これにより、緩和状態の粗視化モデル４をさらに拡散させて、転移域の時間スケール（絡み合いが起こる時間程度、例えば、論文２に従う粗視化モデルの場合、１８００τ程度）の計算を行う。工程Ｓ２１６では、粗視化モデル４のビーズ１６の座標データが出力される。座標データは、コンピュータ１に入力される。
【００８６】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、粗視化モデル４の平均二乗変位を計算する（工程Ｓ２１８）。本実施形態の平均二乗変位は、工程Ｓ２１６の計算開始から計算終了までの単位時間毎のビーズ１６の座標値に基づいて計算される。平均二乗変位（ Mean-square Displacement：ＭＳＤ）は、粗視化モデル４のビーズ１６（図１２に示す）の拡散の大きさを示す指標である。平均二乗変位は、下記式（３）で定義される。
ＭＳＤ＝＜｜ｒ（ｔ'）−ｒ（０）｜2＞  …（３）
なお、符号は次のとおりである。
ｒ（０）：任意の時刻におけるビーズの座標値
ｒ（ｔ'）：任意の時刻から時間（時間幅）ｔ'後におけるビーズの座標値
＜＞：アンサンブル平均
【００８７】
上記式（３）に示されるように、平均二乗変位は、空間９に配置されている全ての粗視化モデル４について、各ビーズ１６（図１２に示す）の任意の時刻から時間後（任意の時間幅ｔ'）に移動した距離の二乗を、ビーズ１６でアンサンブル平均した値である。このような平均二乗変位により、ビーズ１６の平均的な拡散の大きさを把握することができる。
【００８８】
平均二乗変位は、粗視化モデル４を構成する全てのビーズ１６（図１２に示す）のうち、予め定められたビーズ１６について計算されるのが望ましい。これにより、平均二乗変位が計算されるビーズ１６が限定されるため、工程Ｓ２１６で出力される粗視化モデル４のビーズ１６の座標データの総量を削減し、平均二乗変位の計算時間を短縮することができる。また、平均二乗変位は、図４に示されるように、粗視化モデル４を構成する全てのビーズ１６のうち、粗視化モデル４の長手方向の中央側４Ｃに配置される一部のビーズ１６について計算されるのが望ましい。これにより、粗視化モデル４のうち、ビーズ１６の絡み合いが起こりにくい末端側４Ｅ、４Ｅの拡散の影響を除去できるため、転移域での平均二乗変位を精度よく計算することができる。
【００８９】
平均二乗変位が計算されるビーズ１６の個数については、適宜設定することができる。なお、平均二乗変位が計算されるビーズ１６の個数が少ないと、粗視化モデル４のビーズ１６のサンプル数が不足して、平均二乗変位の統計誤差が大きくなり、粗視化モデル４の時間単位を、高分子鎖の時間単位に精度よく換算できないおそれがある。逆に、平均二乗変位が計算されるビーズ１６の個数が多いと、ビーズ１６の絡み合いが起こりにくい末端側４Ｅ、４Ｅの拡散の影響を十分に除去できないおそれがある。このような観点より、平均二乗変位が計算されるビーズ１６の個数は、粗視化モデル４を構成する全てのビーズ１６の個数の３０〜７０％に設定されるのが望ましい。図１８は、粗視化モデル４の平均二乗変位と、時間幅との関係を示すグラフである。このような平均二乗変位の平均値は、コンピュータ１に入力される。
【００９０】
次に、本実施形態の第１計算工程Ｓ２１では、コンピュータ１が、粗視化モデル４同士の絡み合いが起き始めたか否かが判断される（工程Ｓ２１９）。粗視化モデル４同士の絡み合いが起き始めたか否かの判断は、例えば、工程Ｓ２１８で求められた平均二乗変位が、時間幅の約１／２乗に比例する領域（即ち、転移域の拡散に対応する領域）を含み、かつ、時間幅の約１／４乗に比例する領域（即ち、ゴム域の拡散に対応する領域）に向かって、拡散が約１／２乗より遅くなる領域を含んでいるか否かに基づいて判断されうる。なお、このような判断は、粗視化モデルの構造に応じて、適宜設定されうる。従って、第１計算工程Ｓ２１では、粗視化モデル４の配置が大きく変化する前に運動計算が終了される。
【００９１】
工程Ｓ２１９において、粗視化モデル４同士の絡み合いが起き始めたと判断された場合（工程Ｓ２１９で、「Ｙ」）、次の第２計算工程Ｓ２２（図１６及び図１９に示す）が実施される。他方、粗視化モデル４同士の絡み合いが起きていないと判断された場合（工程Ｓ２１９で、「Ｎ」）、単位時間を（例えば、工程Ｓ２１６で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ２２０）、工程Ｓ２１６及び工程Ｓ２１８が実施される。これにより、第１計算工程Ｓ２１では、粗視化モデル４同士の絡み合いが起き始めるまでの（即ち、転移域の拡散に対応する）粗視化モデル４の平均二乗変位が計算される。
【００９２】
次に、時間換算工程Ｓ２では、コンピュータ１が、全原子モデル３の拡散の大きさを示す平均二乗変位を計算する（第２計算工程Ｓ２２）。図１９は、本実施形態の第２計算工程Ｓ２２の処理手順の一例を示すフローチャートである。
【００９３】
本実施形態の第２計算工程Ｓ２２では、先ず、コンピュータ１が、緩和状態の全原子モデル３を作成する（工程Ｓ２２１）。工程Ｓ２２１では、図１７に示した第１計算工程Ｓ２１の工程Ｓ２１５で得られた緩和状態の粗視化モデル４の配置に基づいて、緩和状態の全原子モデル３が作成される。
【００９４】
工程Ｓ２２１では、リバースマッピングに基づいて、緩和状態の粗視化モデル４の各ビーズ１６（図１２に示す）に、モノマーモデル１１（図３に示す）を割り当てている。上述したように、モノマーモデル１１は、高分子鎖２のモノマー５（図２に示す）をなす単位構造に基づいて、粒子モデル７及びボンド８を連結したものである。このようなモノマーモデル１１が、粗視化モデル４の各ビーズ１６に割り当てられることにより、緩和状態の全原子モデル３が設定される。なお、リバースマッピングは、上記材料物性シミュレーションソフトウェアを用いて、容易に行うことができる。
【００９５】
なお、リバースマッピング後の粒子モデル７、７間の衝突（距離が極端に近づきすぎること）を避けるため、衝突が起こらないような十分に低い密度でリバースマッピングを行い、かつ、分子動力学計算プログラムLAMMPSなどを用いてエネルギー最小化を行った後に、全原子モデル３が平衡密度となるまで段階的にアフィン変形とエネルギー最小化とを繰り返し行うのが望ましい。なお、モノマーの立体異性体を保持する必要がある場合は、例えば、特開２０１４−２２５２２６に記載される方法と同様に、立体異性体に関する二面角や結合角のポテンシャル関数を大きくして、エネルギー最小化が行われるのが望ましい。また、全原子モデル３の平衡密度は、鎖長の短い複数の全原子モデル３の緩和計算後の平衡密度から、滑らかな曲線（図示省略）でフィッティングして外挿することによって推定することが望ましい。
【００９６】
本実施形態では、モノマー数換算工程Ｓ１で計算された粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比に基づいて、各ビーズ１６にモノマーモデル１１を割り当てる。これにより、緩和状態の全原子モデル３を精度よく作成することができる。そして、工程Ｓ２２１では、空間９に配置されている緩和状態の全ての粗視化モデル４に基づいて、緩和状態の全原子モデル３が作成される。
【００９７】
工程Ｓ２２１では、先ず、粗視化モデル４（図４に示す）のビーズ数と高分子鎖２（図２に示す）のモノマー数との比が１でない場合、モノマー数と同一のビーズ数になる（即ち、ビーズ数とモノマー数との比が１になる）ように、粗視化モデル４の主鎖上の両末端の位置を保持して、主鎖上を等間隔にリサンプリング（即ち、ビーズ１６の総数を変更して再配置）する。例えば、粗視化モデル４の１個のビーズ１６あたりの高分子鎖２のモノマー数が１．５であり、かつ、粗視化モデル４のビーズ数が８０である場合、対応する高分子鎖２のモノマー数は１２０（即ち、モノマー数１．５×ビーズ数８０）となるため、粗視化モデル４のビーズ数は、８０から１２０となるようにリサンプリングされる。
【００９８】
次に、工程Ｓ２２１では、上記材料物性シミュレーションソフトウェア等を用いてリバースマッピングが実行されることにより、リサンプリング後の粗視化モデル４の各ビーズ１６にモノマーモデル１１が１個割り当てられ、空間９から粗視化モデル４が取り除かれる。これにより、緩和状態の全原子モデル３が配置された高分子材料モデル（図示省略）を得ることができる。
【００９９】
このように、工程Ｓ２２１では、粗視化モデル４に比べて大幅に時間を要する全原子モデル３の構造緩和の計算（即ち、分子動力学計算）をすることなく、緩和状態の全原子モデル３を作成することができるため、計算時間を大幅に短縮することができる。緩和状態の全原子モデル３は、コンピュータ１に入力される。
【０１００】
次に、本実施形態の第２計算工程Ｓ２２では、コンピュータ１が、緩和状態の全原子モデル３について、運動計算を行う（工程Ｓ２２２）。本実施形態の運動計算としては、上述した分子動力学計算が実施される。分子動力学計算は、上述したモノマー数換算工程Ｓ１で実施された分子動力学計算と同一手順で実施される。これにより、緩和状態の全原子モデル３について、高分子材料の転移域の拡散の大きさを計算することができる。工程Ｓ２２２では、全原子モデル３の粒子モデル７の座標データが出力される。座標データは、コンピュータ１に入力される。
【０１０１】
次に、本実施形態の第２計算工程Ｓ２２では、コンピュータ１が、全原子モデル３の平均二乗変位を計算する（工程Ｓ２２３）。本実施形態の平均二乗変位は、工程Ｓ２２２の計算開始から計算終了までの本実施形態の平均二乗変位は、工程Ｓ２２２の計算開始から計算終了までの単位時間毎の粒子モデル７の座標値に基づいて計算される。平均二乗変位（ Mean-square Displacement：ＭＳＤ）は、粒子の拡散の大きさを示す指標である。平均二乗変位は、上記式（３）で定義される。このような平均二乗変位により、各粒子モデル７の各時間幅についての拡散の大きさを把握することができる。なお、上記式（３）において、「ｒ（０）」は、任意の時刻における、ある粒子モデルの座標値として定義される。また、「ｒ（ｔ'）」は、任意の時刻から時間（時間幅ｔ'）後における同じ粒子モデルの座標値として定義される。
【０１０２】
計算時間を短縮するために、平均二乗変位は、全原子モデル３を構成する全ての粒子モデル７のうち、予め定められた粒子モデル７について計算されるのが望ましい。図３に示されるように、平均二乗変位は、全原子モデル３の主鎖方向の中央側に配置される一部の粒子モデル７を対象に計算されるのが望ましい。これにより、全原子モデル３のうち、粒子モデル７の絡み合いが起こりにくい末端側の拡散の影響を除去できるため、転移域での平均二乗変位を精度よく計算することができる。
【０１０３】
平均二乗変位が計算される粒子モデル７の個数については、適宜設定することができる。平均二乗変位が計算される粒子モデル７の個数は、平均ニ乗変位が計算されるビーズ１６の個数と同様の観点より、全原子モデル３を構成する全ての粒子モデル７の個数の３０〜７０％に設定されるのが望ましい。図２０は、全原子モデル３の平均二乗変位と時間幅との関係を示すグラフである。なお、図２０のグラフは、全原子モデル３の設定される温度、圧力等の計算条件に応じて異なる結果が得られる。このため、粗視化モデル４のビーズ数と高分子鎖２のモノマー数との比の計算、及び、時間単位の換算を精度良く行うには、モノマー数換算工程Ｓ１及び時間換算工程Ｓ２において、全原子モデルの計算条件を統一する必要がある。
【０１０４】
次に、本実施形態の第２計算工程Ｓ２２では、コンピュータ１が、全原子モデル３の平均二乗変位が、転移域の拡散に対応しているか否かが判断される（工程Ｓ２２４）。全原子モデル３の平均二乗変位が、転移域の拡散に対応しているか否かの判断は、時間幅のおよそ１／２乗に比例する領域（即ち、転移域の拡散に対応する領域）を統計的に有意に含んでいるか否かに基づいて判断される。従って、第２計算工程Ｓ２２では、全原子モデル３の配置が大きく変化する前に運動計算が終了される。
【０１０５】
工程Ｓ２２４において、全原子モデル３の平均二乗変位が、転移域の拡散に対応していると判断された場合（工程Ｓ２２４で、「Ｙ」）、次の工程Ｓ２３（図１９に示す）が実施される。他方、全原子モデル３の平均二乗変位が、転移域の拡散に対応していないと判断された場合（工程Ｓ２２４で、「Ｎ」）、単位時間を（例えば、工程Ｓ２２２で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ２２５）、工程Ｓ２２２及び工程Ｓ２２３が実施される。これにより、第２計算工程Ｓ２２では、転移域の拡散に対応する全原子モデル３の平均二乗変位が計算されるため、後述の工程Ｓ２３において、粗視化モデル４の平均二乗変位（図１８に示す）と、全原子モデル３の平均二乗変位（図２０に示す）とを重複させることができる。
【０１０６】
次に、本実施形態の時間換算工程Ｓ２では、コンピュータ１が、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算する（工程Ｓ２３）。工程Ｓ２３では、粗視化モデル４の平均二乗変位と、全原子モデル３の平均二乗変位とに基づいて、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算している。
【０１０７】
工程Ｓ２３では、先ず、図１８に示した粗視化モデル４の平均二乗変位と時間との関係、及び、図２０に示した全原子モデル３の平均二乗変位と時間との関係が、統計誤差を考慮した非線形最小二乗法に基づいてフィッティングされる。図２１は、全原子モデル３と粗視化モデル４とをフィッティングした平均二乗変位と時間幅との関係を示すグラフである。
【０１０８】
上記フィッティングは、先ず、粗視化モデル４の平均二乗変位について、全原子モデル３の平均二乗変位と重なる（即ち、転移域の拡散に対応する）時間スケール（時間幅）の範囲を滑らかな曲線（図示省略）でフィッティングされる。次に、この曲線と、全原子モデル３の平均二乗変位とが重なるように、粗視化モデル４の時間単位（時間幅）と全原子モデル３の時間単位（時間幅）との比についてフィッティングされる。フィッティングの手法としては、例えば、論文３に記載されるレーベンバーグ・マーカート法（Levenberg-Marquardt Method）が用いられ、例えば、上記グラフ作成ソフトウェアにより、容易に行うことができる。
【０１０９】
そして、上記フィッティングの結果から、粗視化モデル４の時間単位（時間幅τ）と、全原子モデル３の時間単位（時間幅ns）とを対応づけることにより、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算することができる。換算された粗視化モデル４の時間単位は、コンピュータ１に入力される。
【０１１０】
このように、本実施形態の時間換算工程Ｓ２では、上記特許文献１の換算方法のように、全原子モデル３の配置を十分に緩和するまで（即ち、全原子モデル３同士の絡み合いがほどけるまで）緩和計算を行うことなく、粗視化モデル４の時間単位を、高分子鎖の時間単位に換算することができる。上記特許文献１の緩和計算は、本実施形態の時間換算工程Ｓ２に比べて、多くの時間を要する。従って、時間換算工程Ｓ２では、粗視化モデル４の時間単位を、高分子鎖の時間単位に短時間で換算することができる。
【０１１１】
しかも、本実施形態では、時間幅の広い範囲（即ち、転移域の拡散に対応する時間スケール）において、平均二乗変位をフィッティングすることができるため、末端間ベクトルの相関がなくなる時間として定義される時間のパラメータをピンポイントでフィッティングしていた前記特許文献１に比べて、精度よく換算することができる。
【０１１２】
上述したように、モノマー数換算工程Ｓ１では、全原子モデル３のRouseのビーズ数と、粗視化モデル４のRouseのビーズ数とをフィッティングさせることにより、粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を求めている。これにより、粗視化モデル４には、高分子鎖２のモノマーに対応する複数のビーズ１６が慣性半径の鎖長依存性を考慮して精度よく設定されるため、鎖長依存性を考慮しない上記特許文献１の換算方法に比べて、粗視化モデル４を用いた分子動力学計算の計算結果と、全原子モデル３を用いた分子動力学計算の計算結果とを精度よく対応させることができる。従って、時間換算工程Ｓ２では、粗視化モデル４の時間単位を、高分子鎖の時間単位に精度よく換算することができる。
【０１１３】
次に、図５に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、粗視化モデル４を用いた様々な条件のシミュレーションを実施する（工程Ｓ３）。工程Ｓ３のシミュレーションの一例としては、例えば、緩和状態の高分子材料モデル（図示省略）を用いた変形計算等を含んでいる。
【０１１４】
工程Ｓ３では、モノマー数換算工程Ｓ１で求められた粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比に基づいて、粗視化モデル４が設定される。従って、工程Ｓ３では、粗視化モデル４を用いたシミュレーションにより、高分子鎖２に近似した結果を得ることができる。
【０１１５】
次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル（図示省略）を用いたシミュレーション結果が良好である（即ち、所望の性能を有する）か否かを判断する（工程Ｓ４）。工程Ｓ４では、時間換算工程Ｓ２で求められた結果に基づいて、工程Ｓ３で計算された粗視化モデル４（図４に示す）の時間単位が、高分子鎖２（図２に示す）の時間単位に換算される。これにより、粗視化モデル４を用いた様々な条件の分子動力学計算の計算結果を、実際の高分子鎖２の運動として取り扱うことができる。
【０１１６】
工程Ｓ４において、シミュレーション結果が良好であると判断された場合（工程Ｓ４で、「Ｙ」）、シミュレーションされた高分子材料モデル（図示省略）に基づいて、高分子材料が製造される（工程Ｓ５）。他方、シミュレーション結果が良好でないと判断された場合（工程Ｓ４で、「Ｎ」）、高分子鎖２（図２に示す）の構造や条件等を変更して（工程Ｓ６）、工程Ｓ１〜工程Ｓ４が再度実施される。このように、本実施形態のシミュレーション方法では、所望の性能を有する未知の高分子材料を開発することができる。しかも、本実施形態のシミュレーション方法は、工程Ｓ４において、全原子モデル３を用いたシミュレーションを行わなくても、粗視化モデル４（図４に示す）を用いたシミュレーションにより、高い精度のシミュレーション結果を得ることができるため、高分子材料の開発コストを低減することができる。
【０１１７】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【０１１８】
図５に示した処理手順に従って、粗視化モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を求めるモノマー数換算工程と、粗視化モデルの時間単位を高分子鎖の時間単位に換算する時間換算工程が実施された（実施例）。
【０１１９】
実施例のモノマー数換算工程では、図６に示した処理手順に従って、全原子モデルのRouseのビーズ数を計算する第３計算工程と、粗視化モデルのRouseのビーズ数を計算する第４計算工程と、粗視化モデルのビーズを高分子鎖のモノマー数に換算する工程とが実施された。
【０１２０】
実施例の第３計算工程では、図７に示した処理手順に従って、モノマー数が異なる全原子モデル（本実施形態では、第１全原子モデル〜第３全原子モデル）が空間にそれぞれ配置され、全原子モデル毎にRouseのビーズ数が計算された。図１０は、全原子モデルのRouseのビーズ数とモノマー数との関係を示すグラフである。なお、分子動力学計算において、空間の初期密度が０．００１g／cm3に設定され、温度が３６０Ｋに設定され、圧力が１atmに設定された。また、各全原子モデルのモノマー数（鎖長）、及び、空間に配置された全原子モデルの本数等については、次のとおりである。
第１全原子モデル：
モノマー数：２０、空間内の本数：５０本、緩和時間：２０ns
第２全原子モデル：
モノマー数：３０、空間内の本数：７５本、緩和時間：５０ns
第３全原子モデル：
モノマー数：４０、空間内の本数：７５本、緩和時間：１００ns
【０１２１】
実施例の第４計算工程では、図１１に示した処理手順に従って、ビーズ数が異なる複数の粗視化モデル（本実施形態では、第１粗視化モデル〜第９粗視化モデル）が空間にそれぞれ配置され、粗視化モデルのRouseのビーズ数が計算された。図１４は、粗視化モデルのRouseのビーズ数とビーズ数との関係を示すグラフである。各粗視化モデルのビーズ数（鎖長）、及び、空間に配置された粗視化モデルの本数等については、次のとおりである。
第１粗視化モデル：
ビーズ数：１０、空間内の本数：１０００本、緩和時間：１５０００τ
第２粗視化モデル：
ビーズ数：１３、空間内の本数：１０００本、緩和時間：１５０００τ
第３粗視化モデル：
ビーズ数：２０、空間内の本数：１０００本、緩和時間：１５０００τ
第４粗視化モデル：
ビーズ数：２７、空間内の本数：１０００本、緩和時間：１５０００τ
第５粗視化モデル：
ビーズ数：５０、空間内の本数：２００本、緩和時間：１５０００００τ
第６粗視化モデル：
ビーズ数：８０、空間内の本数：１００本、緩和時間：３００００００τ
第７粗視化モデル：
ビーズ数：１００、空間内の本数：１００本、緩和時間：３００００００τ
第８粗視化モデル：
ビーズ数：２００、空間内の本数：１００本、緩和時間：４５０００００τ
第９粗視化モデル：
ビーズ数：５００、空間内の本数：１５０本、緩和時間：３００００００τ
【０１２２】
実施例の粗視化モデルのビーズを高分子鎖のモノマー数に換算する工程では、非線形最小二乗法に基づいて、全原子モデルのRouseのビーズ数とモノマー数との関係と、粗視化モデルのRouseのビーズ数とビーズ数との関係とがフィッティングされ、粗視化モデルのビーズ数と、高分子鎖２のモノマー数との比（即ち、１個のビーズ１６あたりのモノマー数）が計算された。１個のビーズあたりのモノマー数は、１．４８３±０．０１４であった。
【０１２３】
実施例の時間換算工程では、図１６に示した処理手順に従って、粗視化モデルの拡散の大きさを示す平均二乗変位を計算する第１計算工程と、全原子モデルの拡散の大きさを示す平均二乗変位を計算する第２計算工程と、高分子鎖の時間単位に換算する工程とが実施された。
【０１２４】
実施例の第１計算工程では、図１７に示した処理手順に従って、初期配置された粗視化モデルの緩和状態を計算する工程と、緩和状態の粗視化モデルを用いて運動計算を行い、平均二乗変位を計算する工程とが実施された。なお、粗視化モデルのビーズ数（鎖長）は、上記の１個のビーズあたりのモノマー数から８０と決定された。また、緩和時間は、１５０００τであった。図１８は、粗視化モデルの平均二乗変位と時間幅との関係を示すグラフである。
【０１２５】
実施例の第２計算工程では、図１９に示した処理手順に従って、第１計算工程で得られた緩和状態の粗視化モデルの配置に基づいて緩和状態の全原子モデルを作成する工程と、緩和状態の全原子モデルを用いて運動計算を行い、平均二乗変位を計算する工程が実施された。運動計算（分子動力学計算）での温度は、３６０Ｋに設定された。図２０は、全原子モデルの平均二乗変位と時間幅との関係を示すグラフである。
【０１２６】
実施例の高分子鎖の時間単位に換算する工程では、非線形最小二乗法に基づいて、粗視化モデルの平均二乗変位と時間との関係、及び、全原子モデルの平均二乗変位と時間との関係が、フィッティングされ、粗視化モデルの時間単位が、高分子鎖の時間単位に換算された。その結果、粗視化モデルの時間単位１τが、高分子鎖の時間単位４４．９psに換算された。
【０１２７】
比較のために、上記特許文献１の方法に基づいて、粗視化モデルの時間単位が、高分子鎖の時間単位に換算された（比較例）。比較例では、全原子モデルの末端間ベクトルの長さの時間相関関数を求める際に、全原子モデルを主鎖方向に強制的に引き伸ばす変形計算が行われた。比較例では、全原子モデルの末端間ベクトルの長さの時間相関関数を求める際に、全原子モデルの配置が十分に緩和するまで平衡計算が行われた。なお、換算方法の詳細については、上記特許文献１のとおりである。また、実施例及び比較例で使用されたコンピュータは、２４ＣＰＵ（２ノード）である。
【０１２８】
テストの結果、実施例では、粗視化モデルの時間単位を、高分子鎖の時間単位に換算するのに要した日数は、１５日であった。他方、比較例の方法について、換算するのに要する日数を概算したところ、数年であった。従って、本発明のシミュレーション方法は、粗視化モデルの時間単位を高分子鎖の時間単位に短時間で換算することができた。
【符号の説明】
【０１２９】
Ｓ２１      第１計算工程
Ｓ２２      第２計算工程
Ｓ２３      粗視化モデルの時間単位を、高分子鎖の時間単位に換算する工程
【要約】
【課題】  粗視化モデルの時間単位を高分子鎖の時間単位に、短時間で換算する。
【解決手段】  コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法である。本発明のシミュレーション方法は、高分子鎖の粗視化モデルを用いて運動計算を行い、粗視化モデルの拡散の大きさを示す平均二乗変位を計算する第１計算工程Ｓ２１、高分子鎖の全原子モデルを用いて運動計算を行い、全原子モデルの拡散の大きさを示す平均二乗変位を計算する第２計算工程Ｓ２２、及び、粗視化モデルの平均二乗変位と、全原子モデルの平均二乗変位とに基づいて、粗視化モデルの時間単位を、高分子鎖の時間単位に換算する工程Ｓ２３を含む。
【選択図】図１６【図１】
0006050903.tif 000006
【図２】
0006050903.tif 000007
【図３】
0006050903.tif 000008
【図４】
0006050903.tif 000009
【図５】
0006050903.tif 000010
【図６】
0006050903.tif 000011
【図７】
0006050903.tif 000012
【図８】
0006050903.tif 000013
【図９】
0006050903.tif 000014
【図１０】
0006050903.tif 000015
【図１１】
0006050903.tif 000016
【図１２】
0006050903.tif 000017
【図１３】
0006050903.tif 000018
【図１４】
0006050903.tif 000019
【図１５】
0006050903.tif 000020
【図１６】
0006050903.tif 000021
【図１７】
0006050903.tif 000022
【図１８】
0006050903.tif 000023
【図１９】
0006050903.tif 000024
【図２０】
0006050903.tif 000025
【図２１】
0006050903.tif 000026