JP2019159683A - 高分子材料のシミュレーション方法 - Google Patents

Abstract

【課題】 計算精度を向上させる。【解決手段】 コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法である。この方法は、高分子鎖を、高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、コンピュータに入力する工程Ｓ１と、分子動力学に基づく構造緩和の計算時において、粗視化分子モデルのKuhn長をPacking長で除した値である実効的な密度が、高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかのKuhn長をPacking長で除した値である実効的な密度に近づくように、粗視化分子モデルの隣り合う粒子間に相互作用を定義する工程Ｓ２と、コンピュータが、相互作用が定義された粗視化分子モデルを対象に、構造緩和を計算する工程Ｓ３とを含む。【選択図】図５

Description

本発明は、高分子材料を解析するためのシミュレーション方法に関する。

近年、ゴム材料等の高分子材料の反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。下記特許文献１及び２では、高分子鎖を複数のビーズでモデル化した粗視化モデルを作り、その後、このモデルの空間の単位及び時間の単位を、全原子モデルの空間の単位及び時間の単位に対応付けている。

特許第６０９７１３０号公報 特許第６０５０９０３号公報

上記特許文献２では、粗視化モデルの長さを調節することで、全原子モデルに対応付けている。発明者らは、鋭意研究を重ねた結果、上記特許文献２において、粗視化モデルのKuhn長をPacking長で除した値である実効的な密度が、実際の高分子鎖のそれから乖離する場合があり、熱平衡状態の粗視化モデルの座標を全原子モデルへ長さの単位換算を行うことで取得できる全原子モデルの座標の密度が、全原子モデルの平衡密度から乖離しうることを見出した。したがって、計算精度の向上には、さらなる改善の余地がある。

本発明は、以上のような実状に鑑み案出されたもので、計算精度を向上しうる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、前記高分子鎖を、前記高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、前記コンピュータに入力する工程と、分子動力学に基づく構造緩和の計算時において、前記粗視化分子モデルのKuhn長をPacking長で除した値である実効的な密度が、前記高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかのKuhn長をPacking長で除した値である実効的な密度に近づくように、前記粗視化分子モデルの隣り合う前記粒子間に相互作用を定義する工程と、前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルを対象に、前記構造緩和を計算する工程とを含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用は、前記粗視化分子モデルの屈曲性に影響を与えてもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用は、下記式（１）で定義されてもよい。

ここで、
Ｅ：相互作用ポテンシャル関数
Ｋ：相互作用パラメータ
θ：隣り合う３つの粒子がなす角度

本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用を定義する工程は、前記相互作用の相互作用パラメータが異なる複数の前記粗視化分子モデルを定義して、前記高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかの実効的な密度に近似する実効的な密度を有する前記粗視化分子モデルの前記相互作用を決定する工程を含んでもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルの構造緩和後の第１屈曲度合と、前記全原子モデル又は前記ユナイテッドアトムモデルの構造緩和後の第２屈曲度合とに基づいて、前記粗視化分子モデルの粒子数と、前記高分子鎖のモノマー数との比を計算する工程を含んでもよい。

本発明の高分子材料のシミュレーション方法は、分子動力学に基づく構造緩和の計算時において、粗視化分子モデルの実効的な密度が、高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかの実効的な密度に近づくように、粗視化分子モデルの隣り合う前記粒子間に相互作用を定義する工程を含んでいる。これにより、前記高分子鎖、前記全原子モデル、又は、前記ユナイテッドアトムモデルの空間（長さ）の単位を、前記粗視化分子モデルの空間（長さ）の単位に高精度に対応付けることができ、粗視化モデルを用いたシミュレーションの計算精度を向上させることができる。

高分子材料のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。 ポリブタジエンの構造式である。 粗視化分子モデルの一例を示す概念図である。 理想鎖に近似したときの粗視化分子モデルの部分拡大図である。 シミュレーション方法の一例を示すフローチャートである。 屈曲ポテンシャルの一例を説明する図である。 相互作用定義工程の処理手順の一例を示すフローチャートである。 全原子モデルの一例を示す概念図である。 第１計算工程の処理手順の一例を示すフローチャートである。 全原子モデルが配置された高分子材料モデルの一例を示す概念図である。 各全原子モデルの実効的な密度ｌ_Ｋ／ｐと、モノマー数との関係を示すグラフである。 第２計算工程の一例を示すフローチャートである。 粗視化分子モデルが配置された高分子材料モデルの一例を示す概念図である。 各粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐ、屈曲ポテンシャルの相互作用パラメータＫ、及び、粒子数Ｎの関係を示すグラフの一例である。 全原子モデルのKuhnセグメント数とモノマー数との関係と、粗視化分子モデルのKuhnセグメント数と粗視化粒子数との関係とをフィッティングした結果のグラフである。 全原子モデルのKuhnセグメント数とモノマー数との関係を示すグラフである。

以下、本発明の実施の一形態が図面に基づき説明される。
高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法である。高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。

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

高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、高分子材料として、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある。）が例示される。図２は、ポリブタジエンの構造式である。このポリブタジエンを構成する高分子鎖４は、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とからなるモノマー５｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基（−ＣＨ_２）に替えて、メチル基（−ＣＨ_３）が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。

本実施形態のシミュレーション方法では、高分子鎖４を表現した粗視化分子モデルが用いられる。図３は、粗視化分子モデルの一例を示す概念図である。

粗視化分子モデル６は、複数の粒子７と、粒子７、７間を結合する結合鎖８とを含んで構成されている。粗視化分子モデル６は、高分子鎖４を、高分子鎖４を構成する原子の数よりも少ない複数の粒子７を用いて表現されている。本実施形態の粗視化分子モデル６は、Kremer-Grestモデルである場合が例示されるが、特に限定されるわけではなく、例えば、ＤＰＤ（散逸粒子動力学法）に基づくモデル等であってもよい。粒子７は、例えば、図２に示した高分子鎖４のモノマー５に対応している。

粗視化分子モデル６は、例えば、高分子鎖４の実際の構造に基づいて表現された全原子モデル１１（図８に示す）や、ユナイテッドアトムモデル（以下、単に「ＵＡモデル」ということがある。）に比べて、大きな時空間を扱うことができる。一方、粗視化分子モデル６は、現実の高分子鎖４（図２に示す）とは異なる長さの単位が用いられている。このため、粗視化分子モデル６の分子動力学計算の計算結果を、実際の高分子鎖４の運動として取り扱うためには、粗視化分子モデル６の長さの単位を、高分子鎖４の長さの単位に精度よく換算することが重要である。

しかしながら、Kremer-Grestモデルや、ＤＰＤ（散逸粒子動力学法）に基づくモデルなどのように、個々の原子の座標を扱わない粗視化分子モデル６では、高分子鎖４の太さが考慮されていないため、上記特許文献２のように、高分子鎖４の屈曲性を一致させるように空間（長さ）の単位を換算しても、粗視化分子モデル６の実効的な密度が、高分子鎖
４、全原子モデル１１、又は、ＵＡモデルのそれと必ずしも一致しない。

ここで、「実効的な密度」とは、下記式（２）で定義されるKuhn長ｌ_Ｋを、下記式（３）で定義されるPacking長ｐで除した値ｌ_Ｋ／ｐとして表される量である。これは、ゴム域の１本の高分子鎖４が動くことのできるチューブ状の空間の範囲（体積）あたりの高分子鎖４の本数に比例する量として近似的に扱うことができる。

ここで、
Ｌ：粗視化分子モデル又は高分子鎖（全原子モデル又はＵＡモデル）の全長
Ｒ_g：粗視化分子モデル又は高分子鎖（全原子モデル又はＵＡモデル）の慣性半径

ここで、
Ｍ：粗視化分子モデル又は高分子鎖（全原子モデル又はＵＡモデル）の質量
Ｒ_g：粗視化分子モデル又は高分子鎖（全原子モデル又はＵＡモデル）の慣性半径
ρ：粗視化分子モデル又は高分子鎖（全原子モデル又はＵＡモデル）の密度
Ｎ_A：アボガドロ数

図４は、理想鎖に近似したときの粗視化分子モデルの部分拡大図である。上記式（１）のKuhn長ｌ_Ｋは、長さの単位を持つ量であり、高分子鎖４（図２に示す）の剛直性（即ち、屈曲性とは逆の概念である曲がりにくさ）を表すパラメータである。図４に示されるように、Kuhn長ｌ_Ｋは、高分子物理で知られている理想鎖６’に高分子鎖４が近似したとき、即ち、慣性半径が再現されるように、粒子７’（Kuhnセグメントと呼ぶ）が一定の距離の間隔でランダムな向きに（理想鎖の結合鎖８’で）接続された直鎖として、粗視化分子モデル６が近似したときにおいて、接続された粒子７’（Kuhnセグメント）の間の距離に対応する。

一方、上記式（２）のPacking長ｐは、長さの単位を持つ量であり、高分子鎖１本の占める体積と、末端間距離の二乗平均との比として定義される。このPacking長ｐは、ゴム域の１本の高分子鎖４が動くことのできるチューブ状の空間１５の太さに比例する量として近似的に扱うことができる。上記式（２）及び上記式（３）では、高分子鎖４の末端間距離の二乗平均を、慣性半径の二乗平均の６倍として見積もっている。

粗視化分子モデル６と、高分子鎖４、全原子モデル１１又はＵＡモデルとの間で、実効的な密度ｌ_Ｋ／ｐの乖離（ズレ）が大きくなると、粗視化分子モデル６の長さの単位と、高分子鎖４（図２に示す）の長さの単位との間で換算を行った際に、密度のズレ（即ち、平衡密度からのズレ）が大きくなる。例えば、上記の特許文献２の手法を用いて、粗視化分子モデル６と高分子鎖４との時間単位の換算定数を求める際に、粗視化モデルの構造緩和した座標に重なるように全原子モデルの座標を設定して、全原子モデルの構造緩和した座標を作成すると、得られた座標の密度が、平衡密度から乖離する（即ち、構造緩和が不十分となる）。このような粗視化分子モデル６と、高分子鎖４、全原子モデル１１又はＵＡモデルとの組み合わせが用いられた場合、時間単位の換算定数を精度よく求めることができなくなる。

本実施形態のシミュレーション方法では、構造緩和の計算時において、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが、高分子鎖４（図２に示す）、高分子鎖４の全原子モデル１１（図８に示す）、又は、高分子鎖４のＵＡモデル（図示省略）の実効的な密度ｌ_Ｋ／ｐに近づくように、粗視化分子モデル６の隣り合う粒子７、７間に相互作用を定義して、構造緩和が計算されている。図５は、シミュレーション方法の一例を示すフローチャートである。なお、本実施形態のシミュレーション方法では、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが、全原子モデル１１（図８に示す）の実効的な密度ｌ_Ｋ／ｐに近づくように、粗視化分子モデル６の構造緩和が計算される態様が例示されるが、このような態様に限定されない。例えば、ＵＡモデルの実効的な密度ｌ_Ｋ／ｐに近づくように粗視化分子モデル６の構造緩和が計算されてもよいし、高分子鎖４の実効的な密度ｌ_Ｋ／ｐが既知の場合には、この既知の実効的な密度ｌ_Ｋ／ｐに近づくように、粗視化分子モデル６の構造緩和が計算されてもよい。

本実施形態のシミュレーション方法では、先ず、図３に示した粗視化分子モデル６が、コンピュータ１に入力される（工程Ｓ１）。

高分子鎖４（図２に示す）がポリブタジエンである場合には、論文１（ 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 ）に基づいて、例えば１．５５個分のモノマー５を構造単位として、該構造単位を１個の粒子７に置換される。これにより、粗視化分子モデル６には、複数（例えば、１０〜５０００個）の粒子７が設定される。

粒子７は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粒子７には、質量、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ１に記憶される。

結合鎖８は、粒子７、７間に平衡長を定義した結合ポテンシャルとして構成される。ここで、「平衡長」とは、粒子７、７間の結合距離である。平衡長は、隣り合う粒子７、７の中心７ａ、７ａ間の距離として定義される。また、結合鎖８の結合ポテンシャルには、例えば、上記論文１に基づいて設定されるのが望ましい。このような粗視化分子モデル６は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ１に入力される。

次に、本実施形態のシミュレーション方法では、コンピュータ１に、粗視化分子モデル６の隣り合う粒子間に、相互作用が定義される（相互作用定義工程Ｓ２）。上述したように、相互作用は、分子動力学に基づく構造緩和の計算時において、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが、高分子鎖４（図２に示す）の実効的な密度ｌ_Ｋ／ｐに近づけるためのものである。

相互作用としては、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが、高分子鎖４、全原子モデル１１又はＵＡモデル（本実施形態では、全原子モデル１１）の実効的な密度ｌ_Ｋ／ｐに近づくものであれば、特に限定されない。本実施形態の相互作用としては、粗視化分子モデル６の屈曲性に影響を与える相互作用（以下、単に「屈曲ポテンシャル」ということがある。）が採用される。

屈曲ポテンシャルとしては、粗視化分子モデル６の屈曲性に影響を与えるものであれば、特に限定されない。図６は、屈曲ポテンシャルの一例を説明する図である。本実施形態の屈曲ポテンシャルは、下記式（１）で定義される。

ここで、
Ｅ：相互作用ポテンシャル関数（屈曲ポテンシャル）
Ｋ：相互作用パラメータ
θ：隣り合う３つの粒子がなす角度

上記式（１）において、相互作用パラメータＫが正の場合、隣り合う３つの粒子７がなす角度θが小さくなるほど、屈曲ポテンシャルＥ（θ）が大きくなる。一方、相互作用パラメータＫが負の場合には、正の場合の逆の傾向となる。このように、相互作用パラメータＫの値を調節することで、屈曲ポテンシャルＥ（θ）を適宜設定することができる。このような屈曲ポテンシャルＥ（θ）により、曲がりやすい粗視化分子モデル６や、曲がりにくい粗視化分子モデル６を設定できる。

粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐを、高分子鎖４、全原子モデル１１又はＵＡモデル（本実施形態では、全原子モデル１１）の実効的な密度ｌ_Ｋ／ｐに効果的に近づけるためには、相互作用パラメータＫを適切に決定することが重要である。本実施形態の相互作用定義工程Ｓ２は、相互作用パラメータＫが異なる複数の粗視化分子モデルを定義して、これらの粗視化分子モデルのうち、高分子鎖４、全原子モデル１１又はＵＡモデル（本実施形態では、全原子モデル１１）の実効的な密度によく近似する実効的な密度を有する粗視化分子モデル６の相互作用１０を決定している。

さらに、本実施形態の相互作用定義工程Ｓ２では、上述の相互作用１０を決定するとともに、粗視化分子モデル６（図３に示す）の粒子数と、高分子鎖４（図２に示す）のモノマー数との比も計算している。図７は、相互作用定義工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態の相互作用定義工程Ｓ２では、先ず、全原子モデル１１の実効的な密度ｌ_Ｋ／ｐが取得される（第１計算工程Ｓ２１）。図８は、全原子モデルの一例を示す概念図である。

全原子モデル１１は、高分子鎖４（図２に示す）の実際の構造に基づいて、原子をモデル化した原子モデル１２で表現したものである。全原子モデル１１を用いた分子動力学計算では、現実の高分子鎖４に基づいた長さの単位が用いられている。従って、全原子モデル１１の実効的な密度ｌ_Ｋ／ｐを求めることで、高分子鎖４の実効的な密度ｌ_Ｋ／ｐを求めることができる。

本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の実効的な密度ｌ_Ｋ／ｐを求めるとともに、全原子モデル１１のKuhnセグメント数が求められている。なお、Kuhnセグメント数は、粗視化分子モデル６（図３に示す）の粒子数と、高分子鎖４（図２に示す）のモノマー数との比を計算するのに用いられる。図９は、第１計算工程Ｓ２１の処理手順の一例を示すフローチャートである。

本実施形態の第１計算工程Ｓ２１では、先ず、モノマー数が異なる複数の全原子モデル１１（図８に示す）が設定される（工程Ｓ４１）。図８に示されるように、全原子モデル１１は、複数の原子モデル１２と、原子モデル１２、１２間を結合するボンド１３とを含んで構成されている。図２に示した高分子鎖４のモノマー５を表す単位構造に基づいて、原子モデル１２がボンド１３で連結されることにより、モノマーモデル１４が設定される。このモノマーモデル１４が、予め定められた複数のモノマー数（即ち、分子量（重合度））に基づいて連結される。これにより、モノマー数が異なる複数の全原子モデル１１が設定される。

モノマー数については、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ１の性能等に基づいて、構造緩和計算が現実的な計算時間で完了しうる範囲内で設定されるのが望ましい。なお、モノマー数は、計算精度を維持するために、極端に小さい値を除外するのが望ましい。モノマー数の一例としては、１０〜６０から選択されうる。なお、このモノマー数に限定されるわけではない。

原子モデル１２は、後述の分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、原子モデル１２には、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の原子モデル１２は、高分子鎖４（図２に示す）の炭素原子をモデル化した炭素原子モデル１２Ｃ、及び、高分子鎖４の水素原子をモデル化した水素原子モデル１２Ｈを含んでいる。

ボンド１３は、原子モデル１２、１２間を拘束するものである。本実施形態のボンド１３は、炭素原子モデル１２Ｃ、１２Ｃを連結する主鎖１３Ａ、及び、炭素原子モデル１２Ｃと水素原子モデル１２Ｈとの間を連結する側鎖１３Ｂとを含んでいる。これらの主鎖１３Ａ及び側鎖１３Ｂは、例えば、平衡長とバネ定数とが定義されたバネとして取り扱われる。

全原子モデル１１は、各原子モデル１２、１２間の結合長さである結合長、ボンド１３を介して連続する３つの原子モデル１２がなす角度である結合角、及び、ボンド１３を介して連続する４つの原子モデル１２において、隣り合う３つの原子モデル１２が作る二面角などが定義される。これにより、全原子モデル１１は、三次元構造を有する。全原子モデル１１は、慣例に従い、外力又は内力を受けることによって、結合長、結合角及び二面角が変化する。これにより、全原子モデル１１は、その三次元構造を変化させることができる。

結合長、結合角及び二面角は、例えば、論文２（J. Comput. Chem. 25, 1157-1174 (2004)）に基づいて設定されるポテンシャル（ＧＡＦＦ）によって定義されうる。ポテンシャルは、高分子鎖４の構造に応じて設定されるのが望ましい。このような全原子モデル１１は、材料物性シミュレーションソフトウェア（例えば、（株）ＪＳＯＬ社製のＪ−ＯＣＴＡ）を用いて作成することができる。各全原子モデル１１は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ２１では、モノマー数が異なる各全原子モデル１１の初期配置が決定される（工程Ｓ４２）。工程Ｓ４２では、予め定められた空間に、全原子モデル１１が配置された高分子材料モデルが設定される。図１０は、全原子モデル１１が配置された高分子材料モデル１８の一例を示す概念図である。

本実施形態の工程Ｓ４２では、モノマー数が異なる各全原子モデル１１が、独立して設けられた空間１６にそれぞれ配置される。これにより、工程Ｓ４２では、各全原子モデル１１の初期配置が決定された高分子材料モデル１８がそれぞれ定義される。各高分子材料モデル１８には、同一のモノマー数の全原子モデル１１が複数本配置されている。全原子モデル１１は、モンテカルロ法に基づいて、空間１６内に配置されるのが望ましい。

空間１６は、解析対象の高分子材料の微小構造部分に相当する。本実施形態の空間１６は、互いに向き合う三対の平面１７、１７を有する立方体として定義されている。各平面１７には、周期境界条件が定義されている。従って、一方の平面１７と、反対側の平面１７とが連続している（繋がっている）ものとして取り扱うことができる。

各空間１６に配置される全原子モデル１１の本数については、モノマー数や、空間１６の大きさに基づいて適宜設定される。全原子モデル１１の本数としては、全原子モデル同士の絡まりを防ぐ観点より、平衡時の周期境界長が平衡時の慣性半径の３倍以上に長くなりうる本数が選択されるのが望ましい。本実施形態では、好ましくは２０本以上であり、また、好ましくは２００本以下である。

空間１６の一辺の長さＬａは、系内の原子モデル１２の初期密度が、例えば０．００１g／cm³となるように設定されるのが望ましい。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の初期配置が決定された高分子材料モデル１８において、隣接する全原子モデル１１、１１の原子モデル１２、１２間に、相互作用ポテンシャルＰ１が定義される（工程Ｓ４３）。

本実施形態の相互作用ポテンシャルＰ１は、ＬＪポテンシャルＵ_LJ（ｒ_ij）であり、下記式（４）で定義される。このような相互作用ポテンシャルは、原子モデル１２、１２間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

ここで、各定数及び変数は、Lennard-Jones ポテンシャルのパラメータであり、次のとおりである。
ｒ_ij：原子モデル間の距離
ｒ_c：カットオフ距離
ε：原子モデル間に定義されるＬＪポテンシャルの強度
σ：原子モデルの直径に相当
なお、距離ｒ_ij及びカットオフ距離ｒ_cは、各原子モデル１２、１２の中心間の距離として定義される。

相互作用ポテンシャルは、炭素原子モデル１２Ｃ、１２Ｃ（図８に示す）間に設定される第１ポテンシャル、水素原子モデル１２Ｈ、１２Ｈ（図８に示す）間に設定される第２ポテンシャル、及び、炭素原子モデル１２Ｃと水素原子モデル１２Ｈとの間に設定される第３ポテンシャルを含んでいる。なお、上記式（４）中の各定数は、上記論文１に基づいて、適宜設定することができる。これらの全原子モデル１１の初期配置が決定された高分子材料モデル１８は、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ２１では、初期配置された全原子モデル１１の構造緩和が計算される（工程Ｓ４４）。工程Ｓ４４では、全原子モデル１１の初期配置が決定された各高分子材料モデル１８（図１０に示す）について、分子動力学計算が実施される。分子動力学計算は、例えば、図１０に示した空間１６について所定の時間、配置した全ての全原子モデル１１が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全ての原子モデル１２の動きが追跡され、コンピュータ１に記憶される。また、分子動力学計算の条件は、例えば、系内の原子モデル１２の個数、体積及び温度は一定で行われる。このような分子動力学計算は、例えば、分子動力学計算プログラムLAMMPSを用いて行うことができる。

工程Ｓ４４では、全原子モデル１１の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、全原子モデル１１の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。工程Ｓ４４では、構造緩和計算によって、空間１６の最終的な大きさが、例えば1atmの平衡体積に設定されるのが望ましい。

次に、本実施形態の第１計算工程Ｓ２１では、全ての全原子モデル（即ち、モノマー数が異なる全ての全原子モデル）１１の構造緩和が計算されたか否かが判断される（工程Ｓ４５）。工程Ｓ４５において、全ての全原子モデル１１の構造緩和が計算されたと判断された場合（工程Ｓ４５で、「Ｙ」）、次の工程Ｓ４６が実施される。他方、全ての全原子モデル１１の構造緩和計算が終了していないと判断された場合（工程Ｓ４５で、「Ｎ」）、構造緩和が計算されていない全原子モデル１１（即ち、図１０に示した高分子材料モデル１８）が選択され（工程Ｓ４７）、工程Ｓ４４〜工程Ｓ４５が実施される。これにより、全ての全原子モデル１１の熱平衡状態（構造緩和状態）がそれぞれ計算される。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の全長Ｌ（図示省略）が計算される（工程Ｓ４６）。工程Ｓ４６では、構造緩和した全ての全原子モデル（即ち、モノマー数が異なる全ての全原子モデル）１１について、その全長Ｌがそれぞれ計算される。全長Ｌは、例えば、上記特許文献１の方法に基づいて、全原子モデル１１を主鎖方向に強制的に引き伸ばす変形計算によって求めることができる。各全原子モデル１１の全長Ｌは、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の慣性半径Ｒ_gが計算される（工程Ｓ４８）。工程Ｓ４８では、構造緩和した全ての全原子モデル（即ち、モノマー数が異なる全ての全原子モデル）１１について、その慣性半径Ｒ_gがそれぞれ計算される。全原子モデル１１の慣性半径Ｒ_gは、例えば、上記論文１に記載されている式２．５に基づいて、工程Ｓ４４の分子動力学計算で求められた全原子モデル１１の原子モデル１２の座標値から計算される。各全原子モデル１１の慣性半径は、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の密度ρが計算される（工程Ｓ４９）。密度ρについては、適宜計算することができる。

本実施形態の工程Ｓ４９では、重合度が異なる（即ち、モノマー数が異なる）全原子モデル１１について、それぞれ構造緩和された高分子材料モデル１８（図１０に示す）の密度ρが計算される。計算された密度ρは、コンピュータ１に記憶される。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の第２屈曲度合が計算される（工程Ｓ５０）。本実施形態の第２屈曲度合は、構造緩和後の全原子モデル１１のKuhnセグメント数である。工程Ｓ５０では、構造緩和した全ての全原子モデル（即ち、モノマー数が異なる全ての全原子モデル）１１について、そのKuhnセグメント数がそれぞれ計算される。Kuhnセグメント数は、高分子物理で知られている理想鎖のパラメータの一つであり、下記式（５）で計算される。

ここで、
Ｎ_kuhn：Kuhnセグメント数
Ｌ：全原子モデルの全長
Ｒ_g：全原子モデルの慣性半径

上記式（５）では、例えば、全原子モデル１１の末端間ベクトルＶ１（図８に示す）の長さのアンサンブル平均を用いる上記特許文献１とは異なり、工程Ｓ４８で求められた全原子モデルの慣性半径Ｒ_gと、工程Ｓ４６で求められた全原子モデルの全長Ｌとを用いて、Kuhnセグメント数Ｎ_kuhnが求められる。これにより、本実施形態では、全原子モデル１１の末端側の原子モデル１２の高い運動性の影響を小さくできるため、Kuhnセグメント数Ｎ_kuhnを精度よく求めることができる。

上記式（５）において、Kuhnセグメント数Ｎ_kuhnの値が大きいほど、全原子モデル１１は、コンパクトに折り畳んだ形状となる。他方、Kuhnセグメント数Ｎ_kuhnの値が小さいほど、全原子モデル１１の折れ曲がりが小さくなり、直線状に連続してのびる形状となる。従って、Kuhnセグメント数Ｎ_kuhnは、全原子モデルの屈曲度合を示すパラメータとして用いられる。Kuhnセグメント数Ｎ_kuhnは、コンピュータ１に記憶される。

次に、本実施形態の第１計算工程Ｓ２１では、全原子モデル１１の実効的な密度ｌ_Ｋ／ｐが計算される（工程Ｓ５１）。本実施形態では、重合度依存性の少ない（長鎖長極限の）全原子モデル１１の実効的な密度ｌ_Ｋ／ｐが計算される。

本実施形態の工程Ｓ５１では、先ず、重合度が異なる（即ち、モノマー数が異なる）構造緩和した全原子モデル１１の実効的な密度ｌ_Ｋ／ｐがそれぞれ計算される。上述したように、実効的な密度ｌ_Ｋ／ｐは、上記式（２）で定義されるKuhn長ｌ_Ｋを、上記式（３）で定義されるPacking長ｐで除した値ｌ_Ｋ／ｐで取得することができる。

上記式（２）の全原子モデル１１の全長Ｌには、工程Ｓ４６で求められた全長Ｌが代入される。上記式（２）、（３）の全原子モデル１１の慣性半径Ｒ_gには、工程Ｓ４８で求められた慣性半径Ｒ_gが代入される。上記式（３）の全原子モデル１１の質量Ｍには、１本の全原子モデル１１を構成する原子モデル１２の合計質量が代入される。上記式（３）の全原子モデル１１の密度ρには、工程Ｓ４９で求められた全原子モデル１１の密度ρが代入される。これにより、工程Ｓ５１では、各全原子モデル１１の実効的な密度ｌ_Ｋ／ｐを計算することができる。図１１は、各全原子モデルの実効的な密度ｌ_Ｋ／ｐと、モノマー数との関係を示すグラフである。

次に、工程Ｓ５１では、各全原子モデルの実効的な密度ｌ_Ｋ／ｐの近似曲線２１が求められる。近似曲線２１としては、例えば、ｌ_Ｋ／ｐ＝ｃ１×（１−Ｍ／ｃ２）を用いることができる。なお、変数ｃ１及びｃ２は、フィッティングパラメータである。次に、工程Ｓ５１では、近似曲線２１の分子量∞の値（上記近似曲線２１の変数ｃ１の値）として、実効的な密度の長鎖長極限ｌ_Ｋ／ｐ_∞が求められる。実効的な密度の長鎖長極限ｌ_Ｋ／ｐ_∞は、コンピュータ１に記憶される。

次に、本実施形態の相互作用定義工程Ｓ２では、相互作用パラメータＫが異なる複数の粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが取得される（第２計算工程Ｓ２２）。本実施形態の第２計算工程Ｓ２２では、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐを求めるとともに、粗視化分子モデル６のKuhnセグメント数が求められている。図１２は、第２計算工程Ｓ２２の一例を示すフローチャートである。

本実施形態の第２計算工程Ｓ２２は、先ず、粒子数（鎖長）が異なる複数の粗視化分子モデル６が設定される（工程Ｓ６１）。図３に示されるように、粗視化分子モデル６は、複数の粒子７と、粒子７、７間を結合する結合鎖８とを含んで構成されている。

本実施形態の第２計算工程Ｓ２２は、複数の異なる粒子数（鎖長）に基づいて、粒子７が連結されることにより、粒子数が異なる複数の粗視化分子モデル６が設定される。複数の異なる粒子数については、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ１の性能等に基づいて、適宜設定することができる。粒子数は、例えば、５〜１０００から選択されるのが望ましい。なお、このような粒子数に限定されるわけではない。これらの粗視化分子モデル６は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ２２では、相互作用パラメータＫが異なる複数の粗視化分子モデル６が定義される（工程Ｓ６２）。図６に示されるように、本実施形態の相互作用１０は、上記式（１）で定義される屈曲ポテンシャルＥ（θ）である。工程Ｓ６２では、複数の異なる相互作用パラメータＫに基づいて、粗視化分子モデル６の隣り合う粒子７、７間に屈曲ポテンシャルＥ（θ）が定義されることにより、相互作用パラメータＫが異なる複数の粗視化分子モデル６が定義される。

複数の異なる相互作用パラメータＫについては、例えば、高分子鎖４がとり得る範囲の実効的な密度ｌ_Ｋ／ｐに基づいて適宜設定することができる。本実施形態の相互作用パラメータＫは、例えば、高分子鎖４の実効的な密度ｌ_Ｋ／ｐが２．５〜７．０程度の場合、０〜１．５の範囲から選択されるのが望ましい。

工程Ｓ６２では、工程Ｓ６１で設定された粒子数が異なる各粗視化分子モデル６について、相互作用パラメータＫが異なる複数の粗視化分子モデル６がそれぞれ設定される。即ち、工程Ｓ６２では、粒子数Ｎの種類の数（例えば、２種類（図１４に示したＮ＝５０、１００））に、相互作用の種類の数（例えば、３種類（図１４に示したＫ＝０．８、１．２、１．６））を乗じた種類の数（例えば、６種類）の粗視化分子モデル６が設定される。粒子数及び相互作用が異なる各粗視化分子モデル６は、コンピュータ１に記憶される。

次に、本実施形態の第２計算工程Ｓ２２では、粒子数及び相互作用が異なる各粗視化分子モデル６の初期配置が決定される（工程Ｓ６３）。工程Ｓ６３では、予め定められた空間に、各粗視化分子モデル６がそれぞれ配置された高分子材料モデルが設定される。図１３は、粗視化分子モデル６が配置された高分子材料モデル２０の一例を示す概念図である。

本実施形態の工程Ｓ６３では、粒子数及び相互作用が異なる各粗視化分子モデル６が、独立して設けられた空間１６にそれぞれ配置される。これにより、工程Ｓ６３では、各粗視化分子モデル６の初期配置が決定された高分子材料モデル２０がそれぞれ定義される。各高分子材料モデル２０において、空間１６には、同一の粒子数かつ同一の相互作用の粗視化分子モデル６が複数本配置されている。粗視化分子モデル６は、モンテカルロ法に基づいて、空間１６内に配置されるのが望ましい。

空間１６は、解析対象の高分子材料の微小構造部分に相当し、図１０に示した空間１６と同様に設定される。各空間１６に配置される粗視化分子モデル６の本数については、粒子数や、空間１６の大きさに基づいて適宜設定される。空間１６に配置される粗視化分子モデル６の本数は、好ましくは５０本以上であり、また、好ましくは２０００本以下である。

工程Ｓ６３では、各高分子材料モデル２０において、隣接する粗視化分子モデル６、６の粒子７、７間に、相互作用ポテンシャルＰ２がそれぞれ定義される。本実施形態の相互作用ポテンシャルＰ２は、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（４）で定義される。なお、上記式（４）において、「原子モデル間」は、「粒子間」に置き換えて適用される。

相互作用ポテンシャルＰ２の強度ε、相互作用ポテンシャルＰ２が作用する距離σ、及び、カットオフ距離ｒ_cは、例えば、上記論文２に基づいて設定されるのが望ましい。これらの粗視化分子モデル６の初期配置が決定された高分子材料モデル２０は、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ２２では、初期配置された粗視化分子モデル６の構造緩和が計算される（工程Ｓ６４）。工程Ｓ６４では、粗視化分子モデル６の初期配置が決定された各高分子材料モデル２０（図１３に示す）について、分子動力学計算が実施される。工程Ｓ６４では、各時刻での全ての粒子７の動きが追跡され、コンピュータ１に記憶される。工程Ｓ６４では、粗視化分子モデル６の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、粗視化分子モデル６の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。

次に、本実施形態の第２計算工程Ｓ２２では、全ての粗視化分子モデル（即ち、粒子数及び相互作用が異なる全ての粗視化分子モデル）６の構造緩和が計算されたか否かが判断される（工程Ｓ６５）。工程Ｓ６５において、全ての粗視化分子モデル６の構造緩和が計算されたと判断された場合（工程Ｓ６５で、「Ｙ」）、次の工程Ｓ６６が実施される。他方、全ての粗視化分子モデル６の構造緩和計算が終了していないと判断された場合（工程Ｓ６５で、「Ｎ」）、構造緩和が計算されていない粗視化分子モデル６（即ち、図１３に示した高分子材料モデル２０）が選択され（工程Ｓ６７）、工程Ｓ６４〜工程Ｓ６５が実施される。これにより、全ての粗視化分子モデル６の熱平衡状態（構造緩和状態）がそれぞれ計算される。

次に、本実施形態の第２計算工程Ｓ２２では、粗視化分子モデル６の全長Ｌ（図示省略）が計算される（工程Ｓ６６）。工程Ｓ６６では、構造緩和した全ての粗視化分子モデル（即ち、粒子数及び相互作用が異なる全ての粗視化分子モデル）６について、その全長Ｌがそれぞれ計算される。本実施形態では、粗視化分子モデル６内の隣接粒子の平衡核間距離と、粒子数とを乗じることによって、全長Ｌが求められている。なお、平衡核間距離は、上記の分子動力学計算の結果に基づいて求められる。粗視化分子モデル６の全長Ｌは、コンピュータ１にそれぞれ入力される。

次に、本実施形態の第２計算工程Ｓ２２では、粗視化分子モデル６の慣性半径Ｒ_gが計算される（工程Ｓ６８）。工程Ｓ６８では、構造緩和した全ての粗視化分子モデル（即ち、粒子数及び相互作用が異なる全ての粗視化分子モデル）６について、その慣性半径Ｒ_gがそれぞれ計算される。全原子モデル１１の慣性半径は、例えば、上記論文１に記載されている式２．５に基づいて、工程Ｓ６４の分子動力学計算で求められた粗視化分子モデル６の粒子７の座標値から計算される。各粗視化分子モデル６の慣性半径は、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ２２では、粗視化分子モデル６の第１屈曲度合が求められる（工程Ｓ６９）。本実施形態の第１屈曲度合は、構造緩和後の粗視化分子モデル６のKuhnセグメント数である。工程Ｓ６９では、構造緩和した全ての粗視化分子モデル（即ち、粒子数及び相互作用が異なる全ての粗視化分子モデル）６のKuhnセグメント数が、それぞれ計算される。上述したように、Kuhnセグメント数は、粗視化分子モデルの慣性半径Ｒ_gと、粗視化分子モデルの全長Ｌとを用いて、上記式（５）で計算される。粗視化分子モデル６のKuhnセグメント数Ｎ_kuhnは、コンピュータ１に記憶される。

次に、本実施形態の第２計算工程Ｓ２２では、粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐが計算される（工程Ｓ７０）。本実施形態では、各相互作用（相互作用パラメータＫ）について、重合度依存性の少ない（即ち、長鎖長極限の）粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが計算される。

本実施形態の工程Ｓ７０では、先ず、全ての粗視化分子モデル（即ち、粒子数及び相互作用が異なる全ての粗視化分子モデル）６の実効的な密度ｌ_Ｋ／ｐがそれぞれ計算される。上述したように、実効的な密度ｌ_Ｋ／ｐは、上記式（２）で定義されるKuhn長ｌ_Ｋを、上記式（３）で定義されるPacking長ｐで除した値ｌ_Ｋ／ｐで取得することができる。

上記式（２）の粗視化分子モデル６の全長Ｌには、工程Ｓ６６で求められた全長Ｌが代入される。上記式（２）、（３）の粗視化分子モデル６の慣性半径Ｒ_gには、工程Ｓ６８で求められた慣性半径Ｒ_gが代入される。上記式（３）の粗視化分子モデル６の質量Ｍには、１本の粗視化分子モデル６を構成する粒子７の合計質量が代入される。上記式（３）の粗視化分子モデル６の密度ρには、例えば、０．８５mass／σ³が代入される。これにより、工程Ｓ７０では、各粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐを計算することができる。

本実施形態の工程Ｓ７０では、相互作用（相互作用パラメータＫ）毎に、粒子数Ｎが異なる粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐが取得される。図１４は、各粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐ、屈曲ポテンシャルの相互作用パラメータＫ、及び、粒子数Ｎの関係を示すグラフの一例である。図１４では、相互作用パラメータＫが異なる３つの粗視化分子モデルについて、２種類の粒子数（鎖長）Ｎの実効的な密度ｌ_Ｋ／ｐ（粒子数Ｎ＝５０を一点鎖線、粒子数Ｎ＝１００を二点鎖線で示す）と、実効的な密度の長鎖長極限の推定値ｌ_Ｋ／ｐ（実線）を示している。

次に、工程Ｓ７０では、相互作用（相互作用パラメータＫ）毎に、粒子数Ｎが異なる粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐが取得されたグラフについて、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐの近似曲線２８がそれぞれ求められる。近似曲線２８としては、例えば、ｌ_Ｋ／ｐ＝ｃ１×（１−Ｍ／ｃ２）×ｅｘｐ（ｃ３×Ｋ）を用いることができる。なお、変数ｃ１、ｃ２及びｃ３は、フィッティングパラメータである。次に、工程Ｓ７０では、近似曲線２８の分子量∞の値（前記の式ではｃ１×ｅｘｐ（ｃ３×Ｋ））として、粗視化分子モデル６の長鎖長極限での実効的な密度ｌ_Ｋ／ｐ_∞が、相互作用（相互作用パラメータＫ）の関数として取得される。そして、取得された実効的な密度ｌ_Ｋ／ｐ_∞の近似曲線の全てのパラメータ値が、コンピュータ１に記憶される。

次に、本実施形態の相互作用定義工程Ｓ２では、高分子鎖（全原子モデル１１）の実効的な密度ｌ_Ｋ／ｐ（図１１に示したｌ_Ｋ／ｐ_∞）に近似する実効的な密度ｌ_Ｋ／ｐ（図１４に示したｌ_Ｋ／ｐ_∞）を有する粗視化分子モデル６の相互作用が決定される（工程Ｓ２３）。工程Ｓ２３では、先ず、図１１に示した全原子モデル１１の実効的な密度ｌ_Ｋ／ｐ_∞と、図１４に示した相互作用（相互作用パラメータＫ）の関数としての長鎖長極限の粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐ_∞とが比較される。そして、全原子モデル１１の実効的な密度ｌ_Ｋ／ｐに近似する密度ｌ_Ｋ／ｐを有する粗視化分子モデル６の相互作用（相互作用パラメータＫ）が決定される。

本実施形態の相互作用定義工程Ｓ２では、相互作用パラメータＫが異なる複数の粗視化分子モデル６を定義して、各粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐが計算されている。このため、高分子鎖４の実効的な密度ｌ_Ｋ／ｐに近似する実効的な密度ｌ_Ｋ／ｐを有する粗視化分子モデルの相互作用（相互作用パラメータＫ）を決定することができる。なお、相互作用パラメータＫの決定は、フィッティングに用いられた上記式ｌ_Ｋ／ｐ＝ｃ１×（１−Ｍ／ｃ２）×ｅｘｐ（ｃ３×Ｋ）において、全原子モデル１１のｌ_Ｋ／ｐ_∞を用いて、Ｋ＝ｌｎ（ｌ_Ｋ／ｐ_∞／ｃ１）／ｃ３ によって計算できる。

さらに、本実施形態の相互作用定義工程Ｓ２では、重合度依存性の少ない（長鎖長極限の）全原子モデル１１の実効的な密度ｌ_Ｋ／ｐ（図１１に示したｌ_Ｋ／ｐ_∞）と、重合度依存性の少ない（長鎖長極限の）粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐ（図１４に示したｌ_Ｋ／ｐ_∞）とが比較されている。これにより、相互作用定義工程Ｓ２では、全原子モデル１１の重合度や、粗視化分子モデル６の重合度に左右されることなく、粗視化分子モデル６の実効的な密度を、高分子鎖４の実効的な密度に近づけることができる相互作用パラメータＫを決定することができる。決定された粗視化分子モデル６の相互作用（相互作用パラメータＫ）は、コンピュータ１に記憶される。

次に、本実施形態の相互作用定義工程Ｓ２では、粗視化分子モデル６（図３に示す）の粒子数と、高分子鎖４（図２に示す）のモノマー数との比を計算する（工程Ｓ２４）。工程Ｓ２４では、相互作用が定義された粗視化分子モデル６の構造緩和後の第１屈曲度合と、全原子モデル１１又はユナイテッドアトムモデル（図示省略）の構造緩和後の第２屈曲度合とに基づいて、粗視化分子モデル３の粒子数と、高分子鎖４のモノマー数との比（粗視化分子モデル６の１個の粒子７あたりの高分子鎖４のモノマー数）を求めている。

第１屈曲度合は、工程Ｓ６９で取得された粗視化分子モデル６のKuhnセグメント数である。工程Ｓ２４では、先ず、全ての粗視化分子モデル６のうち、工程Ｓ２３で決定された相互作用パラメータＫを有する複数の粗視化分子モデル６のみを対象に、Kuhnセグメント数と粒子数との関係が取得される。

一方、第２屈曲度合は、工程Ｓ５０で取得された全原子モデル１１のKuhnセグメント数である。工程Ｓ２４では、全原子モデル１１のKuhnセグメント数とモノマー数との関係が取得される。

次に、工程Ｓ２４では、全原子モデル１１のKuhnセグメント数とモノマー数との関係と、粗視化分子モデル６のKuhnセグメント数と粒子数との関係とがフィッティングされる。図１５は、全原子モデル１１のKuhnセグメント数と、粗視化分子モデル６のKuhnセグメント数とをフィッティングした後の、Kuhnセグメント数とモノマー数との比（Ｎ_Kuhn／Ｎ_FA）およびKuhnセグメント数と粒子数との比（Ｎ_Kuhn／Ｎ_CG）の対応関係を示すグラフである。

本実施形態の工程Ｓ２４では、粗視化分子モデル６のKuhnのセグメント数と、粒子数Ｎ_CGとの関係、及び、全原子モデル１１のKuhnセグメント数と、モノマー数Ｎ_FAとの関係がフィッティングされる。これにより、鎖長（粒子数Ｎ_CG、モノマー数Ｎ_FA）の依存性を考慮して、同一鎖長で空間的な広がりが同じとなるような、粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比を計算することができる。

工程Ｓ２４では、粗視化分子モデル６のKuhnのセグメント数と、粒子数Ｎ_CGとの関係が、図１５に示した滑らかな曲線２５（破線で示している）でフィッティングされる。図１５では、鎖長の依存性を視覚化するため、Kuhnのセグメント数と鎖長との比（Ｎ_Kuhn／Ｎ_CG、及び、Ｎ_Kuhn／Ｎ_FA）を縦軸としている。曲線は、下記式（６）で定義される。このような曲線２５は、３つのフィッティングパラメータａ、ｂ、ｃが用いられるため、鎖長の依存性（とりわけ、末端の運動性が高い影響）を考慮して精度よくフィッティングすることができる。

ここで、
Ｎ_kuhn：Kuhnのセグメント数
Ｎ_CG：粗視化分子モデルの粒子数
ａ、ｂ、ｃ：フィッティングパラメータ

次に、工程Ｓ２４では、下記式（７）に基づいて、全原子モデル１１のKuhnセグメント数Ｎ_kuhnと、モノマー数Ｎ_FAとの関係が、曲線２５にフィッティングされる。下記式（７）では、単独のフィッティングパラメータｄが用いられるため、粗視化分子モデル６に比べて統計誤差の大きい全原子モデル１１の関係を、曲線２５に精度よくフィッティングすることができる。

ここで、
Ｎ_FA：全原子モデルのモノマー数
Ｎ_CG：粗視化分子モデルの粒子数
ｄ：フィッティングパラメータ

フィッティングの方法やフィッティングに用いる統計誤差については、特許文献（特許第６０５０９０３号公報）に記載と同一の方法が用いられる。図１５に示したグラフの例において、粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比（即ち、１個の粒子７あたりのモノマー数）は、１．２である。このような粗視化分子モデル６（図３に示す）の粒子数と、高分子鎖４（図２に示す）のモノマー数との比は、コンピュータ１に入力される。

本実施形態の工程Ｓ２４では、全原子モデル１１（図５に示す）のKuhnセグメント数と、粗視化分子モデル６（図３に示す）のKuhnセグメント数とをフィッティングさせることにより、互いに略同一形状となる全原子モデル１１と粗視化分子モデル６とを対応づけることができるため、粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比を精度よく計算することができる。さらに、本実施形態では、全原子モデル１１のKuhnセグメント数、及び、粗視化分子モデル６のKuhnセグメント数がそれぞれ複数求められているため、精度良くフィッティングすることができる。

さらに、本実施形態の工程Ｓ２４では、工程Ｓ２３で決定された相互作用パラメータＫを有する複数の粗視化分子モデル６のみを対象に、Kuhnセグメント数と粒子数との関係が取得されている。これにより、実効的な密度ｌ_Ｋ／ｐが互いに近似する全原子モデル１１及び粗視化分子モデル６のKuhnセグメント数でフィッティングすることができるため、粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比を精度よく計算することができる。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、相互作用が定義された粗視化分子モデル６を対象に、構造緩和が計算される（工程Ｓ３）。工程Ｓ３では、先ず、図１３に示されるように、相互作用が定義された複数の粗視化分子モデル６を空間１６に配置した高分子材料モデル２０が設定される。粗視化分子モデル６の本数については、例えば、上記範囲内に設定されるのが望ましい。

次に、本実施形態の工程Ｓ３では、高分子材料モデル１８において、隣接する粗視化分子モデル６、６の粒子７、７間に、相互作用ポテンシャルＰ２がそれぞれ定義される。相互作用ポテンシャルＰ２については、上述のとおりである。そして、工程Ｓ３では、相互作用が定義された粗視化分子モデル６を対象に、分子動力学に基づく構造緩和が計算される。

このように、本実施形態のシミュレーション方法は、粗視化分子モデル６の隣り合う粒子７、７間に定義された相互作用に基づいて、粗視化分子モデル６の実効的な密度ｌ_Ｋ／ｐを、高分子鎖４（図２に示す）の実効的な密度ｌ_Ｋ／ｐに近づけることができる。これにより、構造緩和計算後の高分子材料モデル２０の座標から、工程Ｓ２４で求められた粗視化分子モデル６の粒子数と高分子鎖４のモノマー数との比に基づいて、粗視化分子モデル６の粒子７を、高分子鎖４のモノマーに置き換えて得られる全原子モデル１１が配置された高分子材料モデル１８の密度を、追加の構造緩和計算を行わなくても、実際の高分子鎖４の密度に近似させることができる。また、追加の構造緩和計算が不要になるため、全原子モデル１１の慣性半径も、実際の高分子鎖４の慣性半径に近似させることができる。すなわち、粗視化分子モデル６が配置された高分子材料モデル２０を用いて、実際の高分子鎖４に近似する構造緩和が可能となるため、構造緩和計算後の高分子材料モデル１８の密度を、実際の高分子鎖４の密度に近似させることができる。本実施形態のシミュレーション方法では、計算精度を向上させることができる。

さらに、本実施形態の粗視化分子モデル６は、工程Ｓ２４で求められた粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比（即ち、１個の粒子７あたりのモノマー数）に基づいて、高分子鎖４のモノマーを１個の粒子７に置換されているため、粗視化分子モデル６の長さの単位を、高分子鎖４の長さの単位に精度よく換算することができる。従って、本実施形態のシミュレーション方法は、計算精度を向上させることができる。

相互作用が定義された粗視化分子モデル６は、例えば、特許文献（特許第６０５０９０３号公報）の時間換算工程に用いられることにより、構造緩和計算後の高分子材料モデル２０から、工程Ｓ２４で求められた粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比に基づいて、粗視化分子モデル６の粒子を高分子鎖４のモノマーに置き換えて高分子材料モデル１８を作成することで、時間単位の換算に必要な長鎖長の全原子モデルの構造緩和した座標を、高い精度で作成できるため、粗視化分子モデル６の時間単位を、高分子鎖の時間単位により精度よく換算するのに役立つ。

次に、本実施形態のシミュレーション方法では、図５に示されるように、コンピュータ１が、高分子材料モデル（図示省略）を用いたシミュレーション結果が良好である（即ち、所望の性能を有する）か否かを判断する（工程Ｓ４）。工程Ｓ４では、例えば、工程Ｓ３の構造緩和後の高分子材料モデル２０を用いた変形計算後に、粗視化分子モデル６（図３に示す）の長さの単位を、高分子鎖４（図２に示す）の長さの単位に精度よく換算した結果に基づいて、高分子鎖４の物性等が評価される。

工程Ｓ４において、シミュレーション結果が良好であると判断された場合（工程Ｓ４で、「Ｙ」）、シミュレーションされた高分子材料モデル（図示省略）に基づいて、高分子材料が製造される（工程Ｓ５）。他方、シミュレーション結果が良好でないと判断された場合（工程Ｓ４で、「Ｎ」）、高分子鎖４（図２に示す）の構造や条件等を変更して（工程Ｓ６）、工程Ｓ１〜工程Ｓ４が再度実施される。これにより、本実施形態のシミュレーション方法では、実際に高分子材料を試作しなくても、所望の性能を有する高分子材料を作成することができる。

本実施形態の相互作用定義工程Ｓ２では、粗視化分子モデル６の相互作用１０を定義するとともに、粗視化分子モデル６（図３に示す）の粒子数と、高分子鎖４（図２に示す）のモノマー数との比が計算されたが、このような態様に限定されない。例えば、相互作用定義工程Ｓ２では、粗視化分子モデル６の相互作用１０のみが定義されるものでもよい。このようなシミュレーション方法では、粗視化分子モデル６の粒子数と、高分子鎖４のモノマー数との比に基づいて、高分子鎖４のモノマーを１個の粒子７に置換した粗視化分子モデル６が用いられなくても、相互作用が定義された粗視化分子モデル６の構造緩和を計算する工程Ｓ３において、実際の高分子鎖４に近似する構造緩和が可能となる。従って、このようなシミュレーション方法では、計算精度を向上させつつ、計算時間を短縮しうる。

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

図５に示した処理手順に従って、高分子鎖（cis-１．４ポリブタジエン）をモデル化した粗視化分子モデルがコンピュータに入力され、分子動力学に基づく構造緩和の計算時において、粗視化分子モデルの実効的な密度が、高分子鎖の実効的な密度に近づくように、粗視化分子モデルの隣り合う粒子間に相互作用が定義された（実施例）。

実施例の相互作用を定義する工程では、図７に示した処理手順に基づいて、相互作用パラメータが異なる複数の粗視化分子モデルを定義して、高分子鎖の実効的な密度に近似する実効的な密度を有する粗視化分子モデルの相互作用が決定された。

実施例の相互作用を定義する工程において、全原子モデルの実効的な密度ｌ_Ｋ／ｐを取得する第１計算工程では、図９に示した処理手順に基づいて、モノマー数（分子量）が異なる複数の全原子モデル（本実施形態では、第１全原子モデル〜第３全原子モデル）が設定され、全原子モデルの実効的な密度ｌ_Ｋ／ｐが取得された。図１１は、各全原子モデルの実効的な密度ｌ_Ｋ／ｐと、モノマー数との関係を示すグラフである。図１６は、全原子モデルのKuhnセグメント数（Ｎ_Kuhn／Ｎ_FA）とモノマー数との関係を示すグラフである。なお、分子動力学計算において、空間の初期密度が０．００１g／cm³に設定され、温度が３６０Ｋに設定され、圧力が１atmに設定された。また、各全原子モデルのモノマー数、及び、空間に配置された全原子モデルの本数等については、次のとおりである。
第１全原子モデル：
モノマー数：２０、空間内の本数：５０本
第２全原子モデル：
モノマー数：３０、空間内の本数：７５本
第３全原子モデル：
モノマー数：４０、空間内の本数：７５本

実施例の相互作用を定義する工程において、粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐを取得する第２計算工程では、図１２に示した処理手順に基づいて、粒子数（鎖長）及び相互作用（相互作用パラメータＫ）が異なる複数の粗視化分子モデル（本実施形態では、第１粗視化モデル〜第６粗視化モデル）が設定され、全原子モデルの実効的な密度ｌ_Ｋ／ｐが取得された。図１４は、各粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐ、屈曲ポテンシャルの相互作用パラメータＫ、及び、粒子数Ｎの関係を示すグラフである。各粗視化モデルの粒子数（分子量）、相互作用（相互作用パラメータＫ）、及び、空間に配置された粗視化モデルの本数等については、次のとおりである。
第１粗視化モデル：
粒子数：５０、相互作用パラメータＫ：０．８、空間内の本数：２００本
第２粗視化モデル：
粒子数：５０、相互作用パラメータＫ：１．２、空間内の本数：２００本
第３粗視化モデル：
粒子数：５０、相互作用パラメータＫ：１．６、空間内の本数：２００本
第４粗視化モデル：
粒子数：１００、相互作用パラメータＫ：０．８、空間内の本数：１００本
第５粗視化モデル：
粒子数：１００、相互作用パラメータＫ：１．２、空間内の本数：１００本
第６粗視化モデル：
粒子数：２００、相互作用パラメータＫ：１．６、空間内の本数：１００本

実施例の相互作用を定義する工程では、相互作用（相互作用パラメータＫ）毎に取得された粗視化分子モデルの実効的な密度ｌ_Ｋ／ｐ（図１４に示した密度ｌ_Ｋ／ｐ_∞）のうち、全原子モデルの実効的な密度ｌ_Ｋ／ｐ_∞に近似する密度ｌ_Ｋ／ｐ（密度ｌ_Ｋ／ｐ_∞）を有する粗視化分子モデルの相互作用（相互作用パラメータＫ）が決定された。そして、決定された相互作用パラメータＫに基づいて、図３に示した粗視化分子モデルの隣り合う粒子間に相互作用が定義され、空間内に１００本配置された高分子材料モデルが定義された。そして、高分子材料モデルに配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。

比較のために、高分子鎖（cis-１．４ポリブタジエン）をモデル化した粗視化分子モデルがコンピュータに入力された（比較例）。比較例の粗視化分子モデルには、実施例のような屈曲性に影響を与える相互作用が定義されていない。そして、比較例の粗視化分子モデルが空間内に１００本配置された高分子材料モデルが定義された。そして、高分子材料モデルに配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。

比較例の高分子材料モデルの密度は、その座標から、密度が一定になるまで追加の構造緩和計算を行った後の密度に比べて、２倍小さかった。平衡化した後の慣性半径は、平衡時に期待される値より２割小さくなった。一方、実施例の高分子材料モデルの密度は、追加の構造緩和計算後の密度に比べて、その差が１％未満であった。このため、実施例では、追加の構造緩和計算は省略できることが確認できた。また、実施例の計算コストは、比較例の計算コストと略同一であった。従って、実施例は、比較例に比べて、計算コストを維持しつつ、計算精度を向上させることができた。

Ｓ１ 粗視化分子モデルを入力する工程
Ｓ２ 相互作用を定義する工程
Ｓ３ 構造緩和を計算する工程

Claims (5)

1. コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、
   前記高分子鎖を、前記高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、前記コンピュータに入力する工程と、
   分子動力学に基づく構造緩和の計算時において、前記粗視化分子モデルのKuhn長をPacking長で除した値である実効的な密度が、前記高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかのKuhn長をPacking長で除した値である実効的な密度に近づくように、前記粗視化分子モデルの隣り合う前記粒子間に相互作用を定義する工程と、
   前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルを対象に、前記構造緩和を計算する工程とを含む、
   高分子材料のシミュレーション方法。
2. 前記相互作用は、前記粗視化分子モデルの屈曲性に影響を与える請求項１記載の高分子材料のシミュレーション方法。
3. 前記相互作用は、下記式（１）で定義される請求項２記載の高分子材料のシミュレーション方法。
   

   

   ここで、
   Ｅ：相互作用ポテンシャル関数
   Ｋ：相互作用パラメータ
   θ：隣り合う３つの粒子がなす角度
4. 前記相互作用を定義する工程は、前記相互作用の相互作用パラメータが異なる複数の前記粗視化分子モデルを定義して、前記高分子鎖、前記高分子鎖の全原子モデル、又は、前記高分子鎖のユナイテッドアトムモデルのいずれかの実効的な密度に近似する実効的な密度を有する前記粗視化分子モデルの前記相互作用を決定する工程を含む請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
5. 前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルの構造緩和後の第１屈曲度合と、前記全原子モデル又は前記ユナイテッドアトムモデルの構造緩和後の第２屈曲度合とに基づいて、前記粗視化分子モデルの粒子数と、前記高分子鎖のモノマー数との比を計算する工程を含む請求項１乃至４のいずれかに記載の高分子材料のシミュレーション方法。

Images