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

Description

本発明は、コンピュータを用いて、高分子材料を解析するための方法に関する。

下記特許文献１は、コンピュータを用いて、高分子材料を解析するための方法を提案している。この種の方法では、先ず、高分子鎖を複数のビーズで粗視化した粗視化モデルが、コンピュータに入力される。そして、分子動力学計算に基づいて、粗視化モデルの緩和した構造が計算される。

一般に、上記の分子動力学計算では、現実の高分子鎖とは異なる単位系が用いられている。例えば、上記の方法では、粗視化分子モデルの全長と、Kuhn長との比であるKuhnセグメント数に基づいて、粗視化モデルの長さの単位と、全原子モデルの長さの単位とが関連付けられている。さらに、上記の方法では、平均二乗変位に基づいて、粗視化モデルの時間の単位と、高分子鎖の時間の単位とが関連付けられている。

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

特許文献１では、粗視化モデルの単位と、高分子鎖の単位とを関連付けることができるものの、その精度については、さらなる改善の余地があった。

粗視化モデルの単位と、高分子鎖の単位とを精度よく関連付けるには、例えば、粗視化モデルと高分子鎖との間において、Kuhn長ｌ_kとPacking長ｐとの比（以下、単に「実効的な密度」ということがある。）ｌ_k／ｐなどの物性値を、互いに近づけることが重要である。このような物性値を精度良く求めるために、例えば、粗視化モデルを用いた分子動力学計算において、十分に長い鎖長のモデルを用いて、最長緩和時間よりも十分に長い緩和計算を行うと、膨大な計算時間が必要となるため、現実的ではない。一方、上記の物理量を現実的な時間で計算できるような十分に短い鎖長の粗視化モデルを用いると、粗視化モデルの鎖長に応じて上記の物性値が変化する傾向がある。したがって、粗視化モデルと高分子鎖との間において、上記の物性値を互いに近づけることは困難であり、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な方法が強く求められていた。

本発明は、以上のような実状に鑑み案出されたもので、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、前記高分子鎖を、前記高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、前記コンピュータに入力する工程と、前記粗視化分子モデルの隣り合う前記粒子間に、Kuhn長とPacking長との比に影響する相互作用パラメータを含む相互作用を定義する工程と、前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルを用いて、分子動力学に基づく構造緩和を計算する工程とを含み、前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記高分子鎖の全原子モデル又は前記高分子鎖のユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記Kuhn長、前記Packing長、及び、前記Kuhn長と前記Packing長との比のいずれかに関する少なくとも２つの物性値がそれぞれ近づくように、前記高分子鎖と前記粗視化分子モデルとの間の長さの単位換算定数、及び、前記相互作用パラメータを定義する工程を含むことを特徴とする。

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

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

本発明に係る前記高分子材料のシミュレーション方法において、前記単位換算定数、及び、前記相互作用パラメータを定義する工程は、前記鎖長及び前記相互作用パラメータがそれぞれ異なる複数の前記粗視化分子モデルを定義する工程と、複数の前記粗視化分子モデルの前記物性値を、前記鎖長と前記相互作用パラメータとの関数に近似させる工程と、 前記鎖長と前記相互作用パラメータとの関数を用いて、前記単位換算定数、及び、前記相互作用パラメータを決定する工程とを含んでもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記鎖長と前記相互作用パラメータとの関数は、前記Kuhn長と前記Packing長との比に近似するものであり、前記鎖長と前記相互作用パラメータとの関数は、下記式（２）で定義されてもよい。

ここで、
ｆ_lk/p：Kuhn長とPacking長との比の近似関数
Ｎ_CG：粗視化分子モデルの鎖長
Ｋ：相互作用パラメータ
ｃ₀～ｃ₅：ｆ_lk/pのフィッティングパラメータ

本発明に係る前記高分子材料のシミュレーション方法において、前記鎖長と前記相互作用パラメータとの関数は、前記粗視化分子モデルの全長と前記Kuhn長との比であるKuhnセグメント数に近似するものであり、前記鎖長と前記相互作用パラメータとの関数は、下記式（３）で定義されてもよい。

ここで、
ｆ_Nk：Kuhnセグメント数の近似関数
Ｎ_CG：粗視化分子モデルの鎖長
Ｋ：相互作用パラメータ
ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
ρ：粗視化分子モデルの数密度
ｆ_lk/p：Kuhn長とPacking長との比の近似関数

本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の一方がそれぞれ近づくように、前記単位換算定数を固定して、前記相互作用パラメータを求める第１工程と、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の他方がそれぞれ近づくように、前記相互作用パラメータを固定して、前記単位換算定数を求める第２工程と、前記単位換算定数と前記相互作用パラメータとが収束するまで、前記第１工程と前記第２工程とを交互に実施させる工程とを含
んでもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用を定義する工程は、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかの密度の鎖長依存性を考慮して、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかのモノマー数と、前記粗視化分子モデルの粒子数とを対応付ける工程を含んでもよい。

本発明に係る前記高分子材料のシミュレーション方法において、前記密度の鎖長依存性は、下記式（４）で定義され、前記モノマー数に対応する前記粒子数は、下記式（５）で定義されてもよい。

ここで、
ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ρ_0、ρ₁：ρ_FAのフィッティングパラメータ

ここで、
Ｎ_CG：モノマー数がＮ_FAの分子鎖に対応する粗視化分子モデルの分子鎖１本の粒子数
Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ｌ_FA：モノマー１つあたりの長さ
Ｃ_l：長さの単位換算定数
ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
ρ₀：フィッティングパラメータ

本発明は、上記の方法を採用することにより、現実的な計算時間で、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な高分子材料のシミュレーション方法を提供することができる。

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

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

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

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

ポリブタジエンを構成する高分子鎖４は、メチレン基（－ＣＨ_２－）とメチン基（－ＣＨ－）とからなるモノマー５｛－［ＣＨ_２－ＣＨ＝ＣＨ－ＣＨ_２］－｝が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基（－ＣＨ_２）に替えて、メチル基（－ＣＨ_３）が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。図３は、高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。

本実施形態のシミュレーション方法では、先ず、高分子鎖４（図２に示す）を表現した粗視化分子モデル６が、コンピュータ１に入力される（工程Ｓ１）。図４は、粗視化分子モデルの一例を示す概念図である。

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

本実施形態の各粒子７は、例えば、高分子鎖４のモノマー５（図２に示す）に対応している。高分子鎖４（図２に示す）がポリブタジエンである場合には、例えば、１．１０個分のモノマー５を構造単位として、該構造単位が１個の粒子７に置換される。これは、下記の論文１に基づくKremer-Grestモデルに、後述の式（１）で定義される相互作用を追加した粗視化分子モデル６について、この粗視化分子モデル６に対応付ける後述の全原子モデルの力場として、下記の論文２に示されるＬ－ＯＰＬＳ力場を用いた場合に、温度３６０Kの分子動力学計算に対応させるためである。また、相互作用パラメータＫは、例えば、１．１２に設定される。これにより、粗視化分子モデル６には、複数（例えば、１０～５０００個）の粒子７が設定される。工程Ｓ１において、高分子鎖４のモノマー数（鎖長）Ｎ_FAと、粗視化分子モデル６の粒子数（鎖長）Ｎ_CGとの比Ｃ_m（１個の粒子７あたりのモノマー数Ｎ_FA）は、例えば、１．１０である。これらの比Ｃ_m、及び、相互作用パラメータＫは、後述の相互作用定義工程Ｓ２において、より適切な値に決定（修正）される。
論文１：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
論文２：S. W. I. Siuら著「Optimization of the OPLS-AA Force Field for Long Hydrocarbons」、J. Chem. Theory. Comput. vol.8, 1459-1470, 4 August 2012

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

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

ところで、粗視化分子モデル６は、例えば、高分子鎖４（図２に示す）、全原子モデル（図９に示す）又はユナイテッドアトムモデル（以下、これらをまとめて「高分子鎖４等」ということがある。）に比べて、大きな時空間を扱うことができるという利点がある。一方、粗視化分子モデル６には、現実の高分子鎖４とは異なる単位系が用いられている。このため、粗視化分子モデル６の計算結果を、高分子鎖４の運動として取り扱うためには、粗視化分子モデル６の単位を、高分子鎖４の単位に精度よく関連付けることが重要である。

Kremer-Grestモデルや、ＤＰＤ（散逸粒子動力学法）に基づくモデル等のように、個々の原子の座標を扱わない粗視化分子モデル６では、高分子鎖４（図２に示す）の個性（例えば、太さや曲がりやすさ等）が十分に考慮されていない。このため、例えば、上記の特許文献１のように、高分子鎖４（図２に示す）のモノマー数（鎖長）Ｎ_FAと、粗視化分子モデル６の粒子数（鎖長）Ｎ_CGとの比Ｃ_mに基づいて、粗視化分子モデル６の長さの単位と、高分子鎖４の長さの単位とを関連付けても、粗視化分子モデル６の密度が、高分子鎖４等の密度と必ずしも一致しないことがある。

粗視化分子モデル６と高分子鎖４との間で長さの単位を関連付ける際に、密度も同時に関連付けるには、粗視化分子モデル６と高分子鎖４との間で、実効的な密度ｌ_k／ｐが精度よく一致するように、粗視化分子モデル６が作成されるのが望ましい。ここで、実効的な密度ｌ_k／ｐとは、下記の式（６）で定義されるKuhn長ｌ_kと、下記の式（７）で定義されるPacking長ｐとの比で表される量である。これは、ゴム域の１本の高分子鎖４（図２に示す）が動くことのできるチューブ状の空間の範囲（体積）あたりの高分子鎖４の本数に比例する量として近似的に扱うことができる。

ここで、
Ｌ：粗視化分子モデル、又は、高分子鎖（全原子モデル又はユナイテッドアトムモデル）の全長
Ｒ_g：粗視化分子モデル、又は、高分子鎖（全原子モデル又はユナイテッドアトムモデル）の慣性半径

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

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

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

粗視化分子モデル６と高分子鎖４等との間で、実効的な密度ｌ_k／ｐの乖離（ズレ）が大きくなると、粗視化分子モデル６の長さの単位と、高分子鎖４の長さの単位との間で関連付け（換算）を行った際に、密度の乖離が大きくなる。例えば、上記の特許文献１の手法を用いて、粗視化分子モデル６と高分子鎖４との時間単位の換算定数を求める際に、粗視化分子モデル６の構造緩和した座標に重なるように、全原子モデル１１の座標を設定して、全原子モデル１１の構造緩和した座標を作成すると、得られた座標の密度が、平衡密度から乖離することがある。このような粗視化分子モデル６と、高分子鎖４等との組み合わせが用いられた場合、例えば、高分子鎖４の拡散運動が、本来の正比例の関係から外れてしまい、時間単位の換算定数を精度よく求めることが困難となるおそれがある。このような事象は、高分子鎖４の太さや曲がりやすさの個性うち、曲がりやすさ（Kuhn長に関する物性値）のみを考慮して単位換算すると生じやすい。

上記のような問題を解決するには、例えば、高分子鎖の太さ（Packing長に関する物性値）と、高分子鎖４の曲がりやすさとの双方を考慮することが有効である。粗視化分子モデル６のKuhn長及びPacking長は、粗視化分子モデル６に定義される相互作用等に依存するため、相互作用を適切に定義することが重要である。とりわけ、Kuhn長とPacking長との比は無次元量であり、長さの単位換算に影響を受けない。このため、粗視化分子モデル６のKuhn長とPacking長との比を、高分子鎖４等のKuhn長とPacking長との比と一致させることが望ましい。

相互作用については、粗視化分子モデル６のKuhn長とPacking長との比に影響を与えることができるものであれば、特に限定されない。本実施形態のシミュレーション方法では、粗視化分子モデル６に、隣り合う（連続する）３つの粒子がなす角度（結合角）θに応じた相互作用（角度ポテンシャル）が定義される（相互作用定義工程Ｓ２）。図６は、相互作用（角度ポテンシャル）Ｅ（θ）の一例を説明する図である。

相互作用（角度ポテンシャル）Ｅ（θ）は、角度（結合角）θと、相互作用ポテンシャル関数Ｅ（θ）の強度を示す相互作用パラメータＫとの関数であり、例えば、下記の式（１）で定義される。

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

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

長さの単位に換算する際に、粗視化分子モデル６の密度を、高分子鎖４（図２に示す）等のいずれかの密度に近づけるためには、相互作用パラメータＫを適切に決定することが重要である。

本実施形態の相互作用定義工程Ｓ２では、構造緩和の計算時において、鎖長が異なる複数の高分子鎖４等と、高分子鎖４等の鎖長に対応する複数の粗視化分子モデル６との間において、Kuhn長、Packing長、及び、Kuhn長と前記Packing長との比のいずれかに関する２つの物性値（以下、これらの２つの物性値を、「第１物性値」及び「第２物性値」として区別する場合がある。）が近づくように、相互作用（相互作用パラメータＫ）、及び、下記式（８）で定義される長さの単位換算定数Ｃ_lが設定される。

ここで、
Ｎ_FA：前記高分子鎖、前記高分子鎖の全原子モデル又は前記高分子鎖のユナイテッドアトムモデルのいずれかの分子鎖１本のモノマー数
ｌ_FA：高分子鎖内のモノマー１つあたりの長さ
Ｎ_CG：モノマー数がＮ_FAの分子鎖に対応する粗視化分子モデルの分子鎖１本の粒子数
ｌ_CG：粗視化分子モデルの結合の平衡長

２つの物性値（第１物性値及び第２物性値）については、Kuhn長ｌ_k、Packing長ｐ、及び、Kuhn長とPacking長との比のいずれかに関するものであれば、適宜設定することができる。本実施形態の第１物性値は、Kuhn長とPacking長との比として得られる実効的な密度ｌ_k／ｐである場合が例示される。一方、本実施形態の第２物性値は、Kuhn長ｌ_kに関連する物理量であるKuhnセグメント数Ｎ_kである場合が例示される。このKuhnセグメント数Ｎ_kは、高分子鎖４や粗視化分子モデル６の全長ＬとKuhn長ｌ_kの比（即ち、Ｎ_k＝Ｌ／ｌ_k）として定義される。

粗視化分子モデル６と高分子鎖４等との間において、各物性値や鎖長の対応関係は、相互作用パラメータＫ、及び、長さの単位換算定数Ｃ_lの双方のパラメータに依存する。このため、各物性値や鎖長の対応関係を、自己無撞着に求めることが重要である。なお、高分子鎖４（図２に示す）のモノマー数（鎖長）Ｎ_FAと、粗視化分子モデル６の粒子数（鎖長）Ｎ_CGとの比Ｃ_m（１個の粒子７あたりのモノマー数Ｎ_FA）については、上記の式（８）で定義される長さの単位換算定数Ｃ_lを用いて、Ｃ_m＝Ｎ_FA／Ｎ_CG＝Ｃ_l・ｌ_CG／ｌ_FAで定義することができる。図７は、相互作用定義工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態の相互作用定義工程Ｓ２では、先ず、鎖長の異なる複数の高分子鎖４等のいずれか（本例では、全原子モデル１１）について、第１物性値（本例では、実効的な密度ｌ_k／ｐ）、及び、第２物性値（本例では、Kuhnセグメント数Ｎ_k）が取得される（第１準備工程Ｓ２１）。なお、第１物性値及び第２物性値が既知である場合には、第１準備工程Ｓ２１が省略されてもよい。

本実施形態の第１準備工程Ｓ２１では、高分子鎖４（図２に示す）等のいずれかのうち、全原子モデル１１の第１物性値（本例では、実効的な密度ｌ_k／ｐ）、及び、第２物性値（本例では、Kuhnセグメント数Ｎ_k）が取得される場合が説明される。図８は、第１準備工程Ｓ２１の処理手順の一例を示すフローチャートである。

本実施形態の第１準備工程Ｓ２１では、先ず、鎖長（即ち、重合度（モノマー数））Ｎ_FAが異なる複数の全原子モデル１１（図９に示す）が設定される（工程Ｓ３１）。図９は、全原子モデル１１の一例を示す概念図である。

全原子モデル１１は、高分子鎖４の実際の構造に基づいて、原子をモデル化した原子モデル１２で表現したものである。全原子モデル１１を用いた分子動力学計算では、現実の高分子鎖４に基づいた長さの単位が用いられている。

全原子モデル１１は、複数の原子モデル１２と、原子モデル１２、１２間を結合するボンド１３とを含んで構成されている。全原子モデル１１では、図２に示した高分子鎖４のモノマー５を表す単位構造に基づいて、原子モデル１２がボンド１３で連結されることにより、モノマーモデル１４が設定される。このモノマーモデル１４が、予め定められた複数のモノマー数（鎖長）Ｎ_FA毎に連結される。これにより、モノマー数Ｎ_FAが異なる複数の全原子モデル１１が設定される。

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

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

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

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

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

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

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

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

各空間１６に配置される全原子モデル１１の本数については、モノマー数（鎖長）Ｎ_FAや、空間１６の大きさに基づいて適宜設定することができる。全原子モデル１１の本数としては、周期境界条件を介した同一分子鎖のコピー同士の絡まりを防ぐ観点より、平衡時の周期境界長が、平衡時の慣性半径の３倍以上に長くなりうる本数に設定されるのが望ましい。本実施形態の本数は、好ましくは２０本以上であり、また、好ましくは２００本以下である。また、空間１６の一辺の長さＬａは、例えば、系内の原子モデル１２の初期密度が、例えば０．００１g／cm³となるように設定されるのが望ましい。

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

相互作用ポテンシャルＰ１は、例えば、結合長、結合角及び二面角のポテンシャルの定義と対応するように、ＧＡＦＦやＬ－ＯＰＬＳ等に基づいて定義されうる。本実施形態の相互作用ポテンシャルＰ１は、ＬＪ（Lennard-Jones）ポテンシャルＵ_LJ（ｒ_ij）であり、下記の式（９）で定義される。このような相互作用ポテンシャルＰ１は、原子モデル１２、１２間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

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

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

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

工程Ｓ３４では、全原子モデル１１の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、全原子モデル１１の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。本実施形態では、例えば、１atmの条件下で体積がほぼ一定となり、かつ、その後のシミュレーション時間が最長緩和時間の3倍以上となったときに、全原子モデル１１の初期配置が十分に構造緩和されたと判断される。最長緩和時間は、例えば、末端間ベクトルＶ１（図９に示す）の自己相関関数を、指数減衰関数に近似させたときの指数緩和時間として見積ることができる。

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

次に、本実施形態の第１準備工程Ｓ２１では、図９に示した全原子モデル１１の第１物性値及び第２物性値が取得される（工程Ｓ３６）。

第１物性値及び第２物性値については、Kuhn長ｌ_k、Packing長ｐ、及び、Kuhn長とPacking長との比（実効的な密度ｌ_k／ｐ）のいずれかに関するものであれば、それぞれ適宜設定することができる。なお、Kuhn長ｌ_kは、上記の式（６）で定義される。Packing長ｐは、上記式（７）で定義される。上述したように、本実施形態では、第１物性値が、Kuhn長とPacking長との比として得られる実効的な密度ｌ_k／ｐであり、第２物性値が、Kuhn長ｌ_kに関連する物理量であるKuhnセグメント数Ｎ_kである。上述のとおり、Kuhnセグメント数Ｎ_kは、全長ＬとKuhn長ｌ_kとの比（即ち、Ｎ_k＝Ｌ／ｌ_k）として定義される。

上記の式（６）の全長Ｌには、全原子モデル１１の全長Ｌ（図示省略）が代入される。この全長Ｌは、例えば、上記の特許文献１の方法に基づいて、全原子モデル１１を主鎖方向に強制的に引き伸ばす変形計算によって求めることができる。本実施形態では、特許文献（特開２０１８－０１０５２５号公報）に基づいて、平均モノマー長ｌ_FAに、モノマー数（鎖長）Ｎ_FAを乗じることにより、全長Ｌが計算される（即ち、Ｌ＝ｌ_FA・Ｎ_FA）。各全原子モデル１１の平均モノマー長ｌ_FAは、構造緩和後の全原子モデル１１において、一対のモノマー間を共有結合の中点で区切り、隣接する一対の中点の間の距離であるモノマー長を平均することで計算される。なお、末端モノマーの長さについては、平均モノマー長の計算から除外されている。構造緩和した全ての全原子モデル（即ち、モノマー数Ｎ_FAが異なる全ての全原子モデル）１１について、それらの全長Ｌがそれぞれ計算される。

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

上記の式（７）の質量Ｍには、１本の全原子モデル１１を構成する原子モデル１２の合計質量が代入される。また、上記の式（７）の密度ρには、全原子モデル１１の密度ρが代入される。これらの計算により、工程Ｓ３６では、第１物性値（本例では、実効的な密度ｌ_k／ｐ）、及び、第２物性値（本例では、Kuhnセグメント数Ｎ_k）を取得することができる。第１物性値及び第２物性値は、コンピュータ１に記憶される。図１１は、各全原子モデル１１の実効的な密度ｌ_k／ｐと、モノマー数（鎖長）Ｎ_FAとの関係を示すグラフである。図１１では、実効的な密度ｌ_k／ｐとモノマー数Ｎ_FAとの関係が、点２３でプロットされている。

次に、本実施形態の相互作用定義工程Ｓ２では、鎖長、及び、屈曲性に影響を与える相互作用がそれぞれ異なる複数の粗視化分子モデルについて、第１物性値及び第２物性値が取得される（第２準備工程Ｓ２２）。図１２は、第２準備工程Ｓ２２の処理手順の一例を示すフローチャートである。

本実施形態の第２準備工程Ｓ２２では、先ず、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６が設定される（工程Ｓ４１）。本実施形態の工程Ｓ４１では、複数の異なる粒子数Ｎ_CGに基づいて、粒子７が結合鎖８で連結されることにより、粒子数Ｎ_CGが異なる複数の粗視化分子モデル６が設定される。

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

次に、本実施形態の第２準備工程Ｓ２２では、相互作用パラメータ（例えば、上記の式（１）の相互作用パラメータＫ）が互いに異なる複数の粗視化分子モデル６が定義される（工程Ｓ４２）。図６に示されるように、本実施形態の工程Ｓ４２では、例えば、互いに異なる複数の相互作用パラメータＫ毎に、粗視化分子モデル６の隣り合う粒子７、７間に相互作用（角度ポテンシャル）Ｅ（θ）がそれぞれ定義される。

互いに異なる相互作用パラメータＫについては、例えば、高分子鎖４（図２に示す）がとり得る範囲の実効的な密度ｌ_k／ｐに基づいて、適宜設定することができる。本実施形態の相互作用パラメータＫは、例えば、高分子鎖４等（本例では、全原子モデル１１）の実効的な密度ｌ_k／ｐが２．５～７．０程度の場合、０～１．５の範囲から選択される。この相互作用パラメータＫは、その値が小さいほど、全原子モデル１１が曲がりやすくなる。一方、相互作用パラメータＫの値が大きいほど、全原子モデル１１が曲がりにくくなる（剛直性が高くなる）。

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

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

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

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

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

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

次に、本実施形態の第２準備工程Ｓ２２では、初期配置された粗視化分子モデル６の構造緩和が計算される（工程Ｓ４４）。本実施形態の工程Ｓ４４では、粗視化分子モデル６の初期配置が決定された各高分子材料モデル２０（図１３に示す）について、分子動力学計算が実施される。工程Ｓ４４では、各時刻での全ての粒子７の動きが追跡され、コンピュータ１に記憶される。

本実施形態の工程Ｓ４４では、粗視化分子モデル６の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、粗視化分子モデル６の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。本実施形態では、シミュレーション時間が、最長緩和時間の3倍以上となったときに、粗視化分子モデル６の初期配置が十分に構造緩和されたと判断している。

次に、本実施形態の第２準備工程Ｓ２２では、全ての粗視化分子モデル（即ち、粒子数（鎖長）Ｎ_CG及び相互作用パラメータＫが異なる全ての粗視化分子モデル）６について、構造緩和が計算されたか否かが判断される（工程Ｓ４５）。工程Ｓ４５において、全ての粗視化分子モデル６の構造緩和が計算されたと判断された場合（工程Ｓ４５で、「Ｙ」）、次の工程Ｓ４６が実施される。

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

次に、本実施形態の第２準備工程Ｓ２２では、全ての粗視化分子モデル（即ち、粒子数（鎖長）Ｎ_CG及び相互作用パラメータＫが異なる全ての粗視化分子モデル）６について、第１物性値及び第２物性値が取得される（工程Ｓ４６）。

第１物性値及び第２物性値については、Kuhn長ｌ_k、Packing長ｐ、及び、Kuhn長とPacking長との比（実効的な密度ｌ_k／ｐ）のいずれかに関するものであれば、それぞれ適宜設定することができる。なお、Kuhn長ｌ_kは、上記式（６）で定義される。Packing長ｐは、上記式（７）で定義される。上述したように、本実施形態では、第１物性値が、Kuhn長とPacking長との比として得られる実効的な密度ｌ_k／ｐであり、第２物性値が、Kuhn長ｌ_kに関連する物理量であるKuhnセグメント数Ｎ_kである。上述のとおり、Kuhnセグメント数Ｎ_kは、全長ＬとKuhn長ｌ_kの比Ｎ_k＝Ｌ／ｌ_kとして定義される。

上記の式（６）の全長Ｌには、粗視化分子モデル６の全長Ｌ（図示省略）が代入される。この全長Ｌは、結合の平衡長ｌ_CGに、粒子数（鎖長）Ｎ_CGを乗じることによって計算することができる（即ち、Ｌ＝ｌ_CG・Ｎ_CG）。上記の結合の平衡長ｌ_CGの値としては、FENEポテンシャルの平衡長（０．９６１σ）の代わりに、温度の影響等を考慮して補正した値（例えば、０．９６５σ）が用いられるのが望ましい。本実施形態では、構造緩和した全ての粗視化分子モデル（即ち、粒子数Ｎ_CGが異なる全ての粗視化分子モデル）６について、それらの全長Ｌがそれぞれ計算される。

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

上記の式（７）の質量Ｍには、１本の粗視化分子モデル６を構成する粒子７の合計質量が代入される。また、上記の式（７）の密度ρには、粗視化分子モデル６の密度ρ（例えば、０．８５σ^-3）が代入される。

次に、図７に示されるように、本実施形態の相互作用定義工程Ｓ２では、鎖長が異なる複数の全原子モデル１１と、鎖長が異なる複数の粗視化分子モデル６との間において、第１物性値（本例では、実効的な密度ｌ_k／ｐ）及び第２物性値（本例では、Kuhnセグメント数Ｎ_k）が近づくように、長さの単位換算定数Ｃ_lと、相互作用パラメータＫとが決定される（パラメータ決定工程Ｓ２３）。本実施形態のパラメータ決定工程Ｓ２３では、例えば、第１準備工程Ｓ２１の結果と、第２準備工程Ｓ２２の結果とを比較することによって、長さの単位換算定数Ｃ_lと相互作用パラメータＫとを決定している。

パラメータ決定工程Ｓ２３の処理手順の一例としては、先ず、相互作用パラメータＫについて、複数の候補値が決定される。これらの候補値には、第２準備工程Ｓ２２で用いられた互いに異なる相互作用パラメータＫ（本例では、上述の０～１．５）が用いられる。

次に、長さの単位換算定数Ｃ_lについて、複数の候補値が決定される。これらの候補値は、互いに異なる長さの単位換算定数Ｃ_lが適宜選択される。候補値が選択される数値範囲の一例としては、０．３～１．０である。

次に、長さの単位換算定数Ｃ_lの候補値について、第１準備工程Ｓ２１で用いた各鎖長（モノマー数）Ｎ_FAに最も近似する第２準備工程Ｓ２２の鎖長（粒子数）Ｎ_CGが、上記の式（８）に基づいて、それぞれ対応付けられる（即ち、Ｎ_CG＝Ｎ_FA・ｌ_FA／（ｌ_CG・Ｃ_l））。なお、粒子数Ｎ_CGが整数である必要がある場合には、小数点以下が四捨五入される。

そして、長さの単位換算定数Ｃ_lの候補値と、相互作用パラメータＫの候補値とのペアそれぞれについて、例えば、対応する高分子鎖等（本例では、全原子モデル１１）と、粗視化分子モデル６との第１物性値（本例では、実効的な密度ｌ_k／ｐ）のズレの最小二乗和と、第２物性値（本例では、Kuhnセグメント数Ｎ_k）のズレの最小二乗和との積をそれぞれ算出し、算出した積が最も小さくなるような候補値のペアが、長さの単位換算定数Ｃ_lと相互作用パラメータＫとして決定される。

このような方法は、単純な処理手順で、長さの単位換算定数Ｃ_lと相互作用パラメータＫとを決定することができるが、精度の高い結果を得るためには、第２準備工程Ｓ２２において、鎖長（粒子数Ｎ_CG）及び相互作用パラメータＫが異なる多くの粗視化分子モデル６を用いた計算が必要となる。このため、計算コスト（即ち、計算時間及び計算費用）が増大しやすい欠点がある。

上記のような欠点を克服可能な好ましい方法として、本実施形態のパラメータ決定工程Ｓ２３の処理手順の一例について説明する。図１４は、パラメータ決定工程Ｓ２３の処理手順の一例を示すフローチャートである。

本実施形態のパラメータ決定工程Ｓ２３では、先ず、長さの単位換算定数Ｃ_lと相互作用パラメータＫとを決定する前に、先ず、第２準備工程Ｓ２２で得られた鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫが異なる複数の粗視化分子モデル６について、第１物性値（本例では、実効的な密度ｌ_k／ｐ）及び第２物性値（本例では、Kuhnセグメント数Ｎ_k）を、鎖長（粒子数）Ｎ_CGと相互作用パラメータＫとの関数に近似させる（工程Ｓ５１）。これにより、本実施形態では、鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫが異なる多くの粗視化分子モデル６を用いなくても、任意の鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫから、第１物性値及び第２物性値を精度よく推定することができる。したがって、本実施形態の方法では、計算コストを削減することが可能となる。

鎖長（粒子数）Ｎ_CGと相互作用パラメータＫの関数については、適宜定義することができる。本実施形態では、粗視化分子モデル６のKuhn長とPacking長との比（即ち、実効的な密度ｌ_k／ｐ）を近似させる関数、及び、Kuhnセグメント数Ｎ_kを近似させる関数が用いられるのが望ましい。鎖長と相互作用パラメータＫの関数は、下記の式（２）及び式（３）で定義されるのが好ましい。なお、下記の式（２）で定義される関数は、実効的な密度ｌ_k／ｐに近似させるものである。下記の式（３）で定義される関数は、Kuhnセグメント数Ｎ_kを近似させるものである。

ここで、
ｆ_lk/p：Kuhn長とPacking長との比（実効的な密度ｌ_k／ｐ）の近似関数
Ｎ_CG：粗視化分子モデルの鎖長
Ｋ：相互作用パラメータ
ｃ₀～ｃ₅：ｆ_lk/pのフィッティングパラメータ

ここで、
ｆ_Nk：Kuhnセグメント数の近似関数
Ｎ_CG：粗視化分子モデルの鎖長
Ｋ：相互作用パラメータ
ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
ρ：粗視化分子モデルの数密度
ｆ_lk/p：Kuhn長とPacking長との比の近似関数

上記の式（２）及び式（３）のフィッティングパラメータｃ₀～ｃ₅は、第２準備工程Ｓ２２で得られた各粗視化分子モデル６（即ち、粒子数（鎖長）Ｎ_CG及び相互作用パラメータＫが異なる各粗視化分子モデル６）の物性値（実効的な密度ｌ_k／ｐ及びKuhnセグメント数Ｎ_k）に、上記の式（２）及び式（３）がフィッティングするように調節することで得られる。フィッティングには、例えば、Levenberg-Marquardt法を用いることができる。図１５は、各粗視化分子モデルの実効的な密度ｌ_k／ｐ、相互作用パラメータＫ、及び、粒子数Ｎ_CGの関係の一例を示すグラフである。図１５では、上記の式（２）のフィッティングの例を示している。フィッティングパラメータｃ₀～ｃ₅は、コンピュータ１に記憶される。

次に、本実施形態のパラメータ決定工程Ｓ２３では、次の工程Ｓ５２～Ｓ５６において、鎖長（モノマー数）Ｎ_FAが異なる高分子鎖４等（本例では、全原子モデル１１）と、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間において、第１物性値と第２物性値の鎖長依存性が近づくように、工程Ｓ５１で得られた上記の式（２）及び式（３）の近似関数を用いて、長さの単位換算定数Ｃ_lと相互作用パラメータＫが決定される。

上記の２つのパラメータ（本例では、長さの単位換算定数Ｃ_l及び相互作用パラメータＫ）を決定するには、例えば、２つのパラメータの初期値（推定値）に基づいて、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等の第１物性値及び第２物性値と、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６の上記の近似関数で推定される第１物性値及び第２物性値とを比較し、各鎖長で第１物性値及び第２物性値が近くなるように、それらの第１物性値のズレの最小二乗和と、第２物性値のズレの最小二乗和との積とが小さくなるような２つのパラメータを探索する方法が考えられる。なお、粗視化分子モデル６の鎖長（粒子数）Ｎ_CGは、長さの単位換算定数Ｃ_lを用いて、高分子鎖４の鎖長（モノマー数Ｎ_FA）に対応付けられる。また、上述の探索には、例えば、差分進化法などの遺伝的アルゴリズムを用いることができる。この方法では、第１物性値及び第２物性値や、これらの近似関数に依存することなく汎用的に用いることができる。しかしながら、この方法では、パラメータの探索範囲を事前に指定する必要があるなどの制限がある。

本実施形態のパラメータ決定工程Ｓ２３は、上述の差分進化法などの遺伝的アルゴリズムを用いずに、次の工程Ｓ５２～工程Ｓ５６において、２つのパラメータ（本例では、長さの単位換算定数Ｃ_l及び相互作用パラメータＫ）を、１つずつ交互に、Levenberg-Marquardt法などの非線形最小二乗法で更新しながら、２つのパラメータの双方が収束するまで繰り返し実施される。

本実施形態のパラメータ決定工程Ｓ２３では、高分子鎖４等と、粗視化分子モデル６との間で、第１物性値（本例では、実効的な密度ｌ_k／ｐ）の長鎖長極限（極限値）が一致するように、相互作用パラメータＫが選択される（工程Ｓ５２）。

本実施形態の工程Ｓ５２では、先ず、第１準備工程Ｓ２１で得られた高分子鎖４等の各鎖長（モノマー数）Ｎ_FAでの第１物性値（以下、「第１物性値Ａ_FA」ということがある。）に、高分子鎖４等の第１物性値Ａ_FAと、高分子鎖４等の鎖長（モノマー数）Ｎ_FAとの関係を示す関数Ａ_FA（Ｎ_FA）をフィッティングさせる。関数Ａ_FA（Ｎ_FA）については、適宜定義することができる。本実施形態の関数Ａ_FA（Ｎ_FA）は、Ａ_FA（Ｎ_FA）＝Ａ₀＋Ａ₁／Ｎ_FAで定義される。Ａ₀及びＡ₁は、フィッティングパラメータである。このような関数Ａ_FA（Ｎ_FA）は、鎖長依存性を有する第１物性値Ａ_FAを、鎖長Ｎ_FAから特定することができる。

上記のフィッティングは、フィッティングパラメータＡ₀及びＡ₁を調節することで行うことができる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。フィッティングパラメータＡ₀及びＡ₁は、コンピュータ１に記憶される。

次に、本実施形態の工程Ｓ５２では、上記の関数Ａ_FA（Ｎ_FA）について、鎖長（本例では、モノマー数）Ｎ_FAを無限大（即ち、長鎖長極限）とする第１物性値Ａ_FAの極限値が求められる（本例では、図１１の実効的な密度の長鎖長極限ｌ_k／ｐ_∞が第１物性値Ａ_FAの極限値Ａ₀となる）。なお、第１物性値Ａ_FAの極限値が発散する場合には、鎖長Ｎ_FAを無限大にしたときに、第１物性値Ａ_FAの極限値が一定値に近づくように規格化した関数Ａ’_FA（Ｎ_FA）を定義して、第１物性値Ａ_FAの極限値に代えて、関数Ａ’_FA（Ｎ_FA）の極限値が求められるのが望ましい。

関数Ａ’_FA（Ｎ_FA）は、上記の規格化のために適宜定義した関数Ａ’’_FA（Ｎ_FA）で、関数Ａ_FA（Ｎ_FA）を除することによって定義することができる（即ち、Ａ’_FA（Ｎ_FA）＝Ａ_FA（Ｎ_FA）／Ａ’’_FA（Ｎ_FA））。例えば、第１物性値Ａ_FAが本実施形態のような実効的な密度ｌ_k／ｐではなく、Kuhnセグメント数Ｎ_kの場合には、関数Ａ_FA（Ｎ_FA）が鎖長（モノマー数Ｎ_FA）に比例して増加する傾向があるため、関数Ａ’’_FA（Ｎ_FA）がモノマー数Ｎ_FAと等しい（即ち、Ａ’’_FA（Ｎ_FA）＝Ｎ_FA）とみなして、関数Ａ’_FA（Ｎ_FA）は、Ａ_FA（Ｎ_FA）／Ｎ_FAで定義することができる。

次に、本実施形態の工程Ｓ５２では、粗視化分子モデル６の第１物性値（以下、「第１物性値Ａ_CG」ということがある。）について、粗視化分子モデル６の鎖長（本例では、粒子数）Ｎ_CG及び相互作用パラメータＫとの関係を示す関数Ａ_CG（Ｎ_CG，Ｋ）を定義する。関数Ａ_CG（Ｎ_CG，Ｋ）については、適宜定義することができる。本実施形態の関数Ａ_CG（Ｎ_CG，Ｋ）は、上記の式（２）に基づいて定義することができる（即ち、Ａ_CG（Ｎ_CG，Ｋ）＝ｆ_ｌk/p（Ｎ_CG，Ｋ））。このような関数Ａ_CG（Ｎ_CG，Ｋ）は、鎖長依存性及び相互作用パラメータ依存性を有する第１物性値Ａ_CGを、鎖長Ｎ_CG及び相互作用パラメータＫから特定することができる。

そして、本実施形態の工程Ｓ５２では、関数Ａ_CG（Ｎ_CG，Ｋ）に、第２準備工程Ｓ２２で得られた粗視化分子モデル６の各鎖長Ｎ_CG及び相互作用パラメータＫでの第１物性値Ａ_CGをフィッティングさせる。本実施形態のフィッティングは、上記の式（２）のフィッティングパラメータｃ０～ｃ５を調節することで行うことができる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。フィッティングパラメータ（本例では、ｃ０～ｃ５）は、コンピュータ１に記憶される。

次に、本実施形態の工程Ｓ５２では、関数Ａ_CG（Ｎ_CG，Ｋ）について、鎖長（粒子数）Ｎ_CGを無限大（即ち、長鎖長極限）とする第１物性値Ａ_CGの極限値の関数ｌｉｍ _ＮCG→∞Ａ_CG（Ｎ_CG，Ｋ）が定義される。この第１物性値Ａ_CGの極限値の関数は、図１５の近似曲線２８で例示される。なお、第１物性値Ａ_CGの極限値が発散する場合には、鎖長Ｎ_CGを無限大にしたときに、第１物性値Ａ_CGの極限値の関数が有限となるように規格化した関数Ａ’_CG（Ｎ_CG，Ｋ）を定義して、第１物性値Ａ_CGの極限値の関数の代わりに、Ａ’_CG（Ｎ_CG，Ｋ）の極限値の関数ｌｉｍ_ＮCG→∞Ａ’_CG（Ｎ_CG，Ｋ）が定義されるのが望ましい。

関数Ａ’_CG（Ｎ_CG，Ｋ）は、上述の高分子鎖４等の関数Ａ’ ’_FA（Ｎ_FA）を用いて、鎖長（モノマー数）Ｎ_FAを、鎖長（粒子数）Ｎ_CGで対応付けたもの（Ｎ_CG＝Ｎ_FA・ｌ_FA／（Ｃ_l・ｌ_CG））によって置き換えることで（即ち、Ａ’_CG（Ｎ_CG，Ｋ）＝Ａ_CG（Ｎ_CG，Ｋ）／Ａ’’_FA（Ｎ_CG・（Ｃ_l・ｌ_CG）／ｌ_FA））、定義することができる。なお、上記の対応付けには、長さの単位換算定数Ｃ_l、並びに、上記の式（８）に示される粗視化分子モデルの結合の平衡長ｌ_CG、及び、高分子鎖内のモノマー１つあたりの長さｌ_FAが用いられる。

次に、本実施形態の工程Ｓ５２では、高分子鎖４等の第１物性値Ａ_FAの極限値（本例では、図１１の実効的な密度の長鎖長極限ｌ_k／ｐ_∞）と、粗視化分子モデル６の第１物性値Ａ_CGの極限値の関数ｌｉｍ _ＮCG→∞Ａ_CG（Ｎ_CG，Ｋ）（図１５の近似曲線２８に例示される）が一致するような相互作用パラメータＫが選択される。第１物性値Ａ_FA及びＡ_CGが長鎖長極限で発散する場合には、第１物性値Ａ_FAの極限値と、第１物性値Ａ_CGの極限値の関数ｌｉｍ _ＮCG→∞Ａ_CG（Ｎ_CG，Ｋ）の代わりに、関数Ａ’_FA（Ｎ_FA）の極限値と、Ａ’_CG（Ｎ_CG，Ｋ）の極限値の関数ｌｉｍ _ＮCG→∞Ａ’_CG（Ｎ_CG，Ｋ）とが一致するような相互作用パラメータＫが選択される。相互作用パラメータＫの選択には、例えば、２分法などの求根アルゴリズムを用いることができる。なお、上記の極限値が一致するような相互作用パラメータＫが存在しない場合には、第２準備工程Ｓ２２において、鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫの範囲を大きくして、これまでの手順がやり直されてもよい。

上述したように、本実施形態の第１物性値Ａ_FA及びＡ_CGは、実効的な密度ｌ_k／ｐである。この実効的な密度ｌ_k／ｐは、無次元量であり、長さの単位換算定数Ｃ_lに依存しない。このため、高分子鎖４等と粗視化分子モデル６との間で、長さの単位換算定数Ｃ_lを用いた単位換算を行う必要がなく、高分子鎖４等の第１物性値Ａ_FAの極限値と、粗視化分子モデル６の第１物性値Ａ_CGの極限値とを直接比較することができる。

なお、第１物性値Ａ_FA、Ａ_CGが無次元量ではない物性値（例えば、Kuhn長ｌ_k）である場合には、高分子鎖４等と粗視化分子モデル６との間で、長さの単位換算定数Ｃ_lを用いた単位換算を行う必要がある。この長さの単位換算には、長さの単位換算定数Ｃ_lの初期値として、例えば、０．５ｎｍ／σが用いられる。これは、１σ＝０．５ｎｍを意味する。第１物性値Ａ_FA及びＡ_CGとしてKuhn長ｌ_kが用いられた場合、粗視化分子モデル６のKuhn長（単位は「σ」）は、高分子鎖４等のKuhn長（単位は「ｎｍ」）から、長さの単位換算定数Ｃ_lを除することで換算することができる。

本実施形態では、上述のとおり、相互作用が上記の式（１）で定義され、かつ、第１物性値が実効的な密度ｌ_k／ｐであり、かつ、粗視化分子モデル６の実効的な密度ｌ_k／ｐを上記の式（２）に近似させている。この場合、粗視化分子モデル６の第１物性値（実効的な密度ｌ_k／ｐ）の極限値は、上記の式（２）の近似関数から、ｌｉｍ _ＮCG→∞ｆ_lk/p（Ｎ_CG，Ｋ）＝ｃ₀・ｅｘｐ［ｃ₃Ｋ（１＋ｃ₅Ｋ）］として特定することができる。

一方、高分子鎖４等の第１物性値（実効的な密度ｌ_k／ｐ）を、上記の関数Ａ_FA（Ｎ_FA）＝Ａ₀＋Ａ₁／Ｎ_FAに近似させた場合、高分子鎖４等の第１物性値（実効的な密度ｌ_k／ｐ）の極限値は、フィッティングパラメータＡ₀と同一の値となる。

高分子鎖４等の第１物性値の極限値（フィッティングパラメータＡ₀）と、粗視化分子モデル６の第１物性値の極限値（ｃ₀・ｅｘｐ［ｃ₃Ｋ（１＋ｃ₅Ｋ）］）とが一致するような相互作用パラメータＫは、例えば求根アルゴリズムを用いなくても、解析的に求めることができる。本実施形態のように、例えば、相互作用パラメータＫ≧０と仮定した場合、相互作用パラメータＫは、Ｋ＝－１＋ｓｑｒｔ［１＋４（ｃ₅／ｃ₃）ｌｎ（Ａ₀／ｃ₀）］／（２ｃ₅）となる。このように、実施形態によっては、求根アルゴリズムを用いることなく、相互作用パラメータＫを解析的に求めることができる場合がある。

本実施形態の工程Ｓ５２では、高分子鎖４等と、粗視化分子モデル６との間で、第１物性値（本例では、実効的な密度ｌ_k／ｐ）の長鎖長極限（極限値）が一致するように、相互作用パラメータＫが選択されたが、このような態様に限定されるわけではなく、例えば、任意の相互作用パラメータＫが指定されてもよい。選択された相互作用パラメータＫは、コンピュータ１に記憶される。

次に、本実施形態のパラメータ決定工程Ｓ２３では、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等と、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間において、物性値の一方（本例では、第２物性値）がそれぞれ近づくように、相互作用パラメータＫを固定して、長さの単位換算定数Ｃ_lが求められる（第１工程Ｓ５３）。本実施形態の第１工程Ｓ５３では、工程Ｓ５１で得られた第２物性値の近似関数（本例では、上記の式（３））を用いて、工程Ｓ５２で選択された相互作用パラメータＫで固定して、第２物性値（Kuhnセグメント数Ｎ_k）がそれぞれ近づくように（本例では、一致するように）、長さの単位換算定数Ｃ_lが求められる。

本実施形態の第１工程Ｓ５３では、第１準備工程Ｓ２１で得られた高分子鎖４等の各鎖長（モノマー数）Ｎ_FAでの第２物性値（以下、「第２物性値Ｂ_FA」ということがある）に、粗視化分子モデル６の第２物性値Ｂ_CGと、粗視化分子モデル６の鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫとの関係を示す関数Ｂ_CG（Ｎ_CG，Ｋ）をフィッティングさせる。関数Ｂ_CG（Ｎ_CG，Ｋ）については、適宜定義することができる。本実施形態の関数Ｂ_CG（Ｎ_CG，Ｋ）は、上記の式（３）に基づいて定義することができる（即ち、Ｂ_CG（Ｎ_CG，Ｋ）＝ｆ_Nk（Ｎ_CG，Ｋ））。上記のフィッティングは、相互作用パラメータＫを固定した状態で、上記の式（８）で定義される長さの単位換算定数Ｃ₁のみをフィッティングパラメータとして調節することで行われる。フィッティングには、例えば、Levenberg-Marquardt法を用いることができる。

上述したように、本実施形態の第２物性値は、Kuhnセグメント数Ｎ_kである。このKuhnセグメント数Ｎ_kは無次元量であり、長さの単位換算定数Ｃ_lに依存しない。このため、高分子鎖４等と粗視化分子モデル６との間で、長さの単位換算を行う必要がなく、高分子鎖４等のKuhnセグメント数Ｎ_kと、粗視化分子モデル６のKuhnセグメント数Ｎ_kとを直接比較することができる。

図１６は、粗視化分子モデルのKuhnセグメント数Ｎ_kと粒子数（鎖長）Ｎ_CGとの関係、及び、全原子モデルのKuhnセグメント数Ｎ_kとモノマー数（鎖長）Ｎ_FAとの関係を示すグラフである。第２物性値がKuhnセグメント数Ｎ_kである場合、第１工程Ｓ５３での長さの単位換算定数Ｃ_lは、相互作用パラメータＫを考慮する点を除き、上記特許文献１に記載の長さの単位換算定数の算出手順に基づいて求めることができる。なお、第２物性値が無次元量でない物性値である場合には、高分子鎖４等と粗視化分子モデル６との間で、長さの単位換算を行う必要がある。単位換算の方法については、上述の第１物性値Ａ_FA、Ａ_CGと同様の方法が用いられる。

本実施形態では、相互作用が上記の式（１）で定義され、粗視化分子モデル６の実効的な密度ｌ_k／ｐを上記式（２）に近似させており、粗視化分子モデルのKuhnセグメント数Ｎ_kを上記式（３）に近似させている。

上記式（２）の近似関数では、粗視化分子モデル６の鎖長（粒子数）Ｎ_CG、及び、相互作用パラメータＫの広い範囲で、粗視化分子モデル６の物性値を高精度に近似させることができ、かつ、パラメータ数が少なく簡素になるように設計されている。このため、上記の式（２）の近似関数では、粗視化分子モデル６の鎖長Ｎ_CG、及び、相互作用パラメータＫの広い範囲において、実効的な密度ｌ_k／ｐを高精度かつ短い計算時間で近似させることができる。

一方、上記の式（３）は、上記の式（２）、式（６）及び式（７）を全て満たすように算出されているため、上記の式（２）との整合性が良好である。また、上記の式（２）及び式（３）で近似させている物性値（本例では、実効的な密度ｌ_k／ｐ及びKuhnセグメント数Ｎ_k）は、共に無次元量であるため、高分子鎖４等と粗視化分子モデル６との物性値のフィッティングを行う際の単位換算が不要である。このため、本実施形態では、無次元量ではない物性値を用いる場合に比べ、計算時間を短縮することができる。したがって、パラメータ決定工程Ｓ２３において、長さの単位換算定数Ｃ_l、及び、相互作用パラメータＫを、高い精度かつ効率的に推定できる。

このように、本実施形態の方法では、上記の式（２）及び（３）を共に用いることで、長さの単位換算定数Ｃ_lおよび相互作用パラメータＫを高い精度で、かつ効率的に推定できる。求められた長さの単位換算定数Ｃ_lは、他の工程で求められる長さの単位換算定数Ｃ_lとは独立して、コンピュータ１に記憶される。

次に、本実施形態のパラメータ決定工程Ｓ２３では、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等と、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間において、物性値の他方（本例では、第１物性値）がそれぞれ近づくように、長さの単位換算定数Ｃ_lを固定して、相互作用パラメータＫが求められる（第２工程Ｓ５４）。本実施形態の第２工程Ｓ５４では、工程Ｓ５１で得られた第１物性値の近似関数（本例では、上記の式（２））を用いて、第１工程Ｓ５３で求められた長さの単位換算定数Ｃ_lで固定して、第１物性値（実効的な密度ｌ_k／ｐ）が近づくように（本例では、一致するように）、相互作用パラメータＫが求められる。

本実施形態の第２工程Ｓ５４では、第１準備工程Ｓ２１で得られた高分子鎖４等の各鎖長（モノマー数）Ｎ_FAでの第１物性値Ａ_FAに、粗視化分子モデルの第１物性値Ａ_CGと、粗視化分子モデル６の鎖長（粒子数）Ｎ_CG及び相互作用パラメータＫとの近似関数（本例では、上述の関数Ａ_CG（Ｎ_CG，Ｋ））をフィッティングさせる。フィッティングは、長さの単位換算定数Ｃ_lを固定した状態で、相互作用パラメータＫをフィッティングパラメータとして調節することで行われる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。図１７は、粗視化分子モデルの実効的な密度ｌ_k／ｐと粒子数Ｎ_CGとの関係、及び、全原子モデルの実効的な密度ｌ_k／ｐとモノマー数Ｎ_FAとの関係を示すグラフである。

上述したように、本実施形態の第１物性値は、実効的な密度ｌ_k／ｐである。上記のような相互作用パラメータＫを求めることは、図１７に示されるように、粗視化分子モデル６実効的な密度ｌ_k／ｐと粒子数（鎖長）Ｎ_CGとの関係を示し、かつ、相互作用パラメータＫが異なる曲線３３のうち、図１１で示される高分子鎖４等の実効的な密度ｌ_k／ｐ（図１７の点２３）に最も近似する曲線３３を選択することに対応している。求められた相互作用パラメータＫは、他の工程で求められる相互作用パラメータＫとは独立して、コンピュータ１に記憶される。

次に、本実施形態の相互作用定義工程Ｓ２では、長さの単位換算定数Ｃ_lと相互作用パラメータＫとがそれぞれ収束したか否かが判断される（工程Ｓ５５）。収束したか否かの判断については、適宜行うことができる。本実施形態において、複数回実施された第１工程Ｓ５３及び第２工程Ｓ５４の前後において、長さの単位換算定数Ｃ_l及び相互作用パラメータＫがそれぞれ数値誤差の範囲内（例えば、０．０１％以内）に収まっている場合に、収束したと判断することができる。

工程Ｓ５５において、長さの単位換算定数Ｃ_lと相互作用パラメータＫとがそれぞれ収束したと判断された場合（工程Ｓ５５において、「Ｙ」）、次の工程Ｓ５６が実施される。この工程Ｓ５６では、収束した長さの単位換算定数Ｃ_lが、長さの単位換算定数Ｃ_lとして決定される。さらに、工程Ｓ５６では、収束した相互作用パラメータＫが、相互作用パラメータＫとして決定される。このような長さの単位換算定数Ｃ_l及び相互作用パラメータＫは、第１工程Ｓ５３及び第２工程Ｓ５４において、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等と鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間で、第１物性値（実効的な密度ｌ_k／ｐ）及び第２物性値（Kuhnセグメント数）が近づくように、長さの単位換算定数Ｃ_lと相互作用パラメータＫとの双方を収束させて求められたものである。したがって、本実施形態の相互作用定義工程Ｓ２では、粗視化分子モデル６の単位と、高分子鎖４の単位とを精度よく対応づけることが可能な長さの単位換算定数Ｃ_l及び相互作用パラメータＫを求めることができる。

一方、工程Ｓ５５において、長さの単位換算定数Ｃ_lと相互作用パラメータＫとがそれぞれ収束していないと判断された場合（工程Ｓ５５で、「Ｎ」）、第１工程Ｓ５３及び第２工程Ｓ５４が再度実施される。なお、再度実施される第１工程Ｓ５３では、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等と、鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間において、物性値の一方（本例では、第２物性値）がそれぞれ近づくように、先の第２工程Ｓ５４で求められた相互作用パラメータＫで固定して、長さの単位換算定数Ｃ_lが再度求められる。これにより、パラメータ決定工程Ｓ２３では、長さの単位換算定数Ｃ_lと相互作用パラメータＫとが収束するまで、第１工程Ｓ５３と第２工程Ｓ５４とが交互に実施されるため、長さの単位換算定数Ｃ_lと相互作用パラメータＫとを確実に収束させることができる。

このように、本実施形態の相互作用定義工程Ｓ２では、鎖長（モノマー数）Ｎ_FAが異なる複数の高分子鎖４等と鎖長（粒子数）Ｎ_CGが異なる複数の粗視化分子モデル６との間で、第１物性値及び第２物性値が近づくように、互いに依存する関係にある長さの単位換算定数Ｃ_lと相互作用パラメータＫとをそれぞれ収束させて求めることができる。このため、相互作用定義工程Ｓ２では、粗視化分子モデル６の単位と、高分子鎖４の単位とを精度よく対応づけることが可能な長さの単位換算定数Ｃ_l及び相互作用パラメータＫを確実に求めることができる。

相互作用定義工程Ｓ２では、高分子鎖４等の密度の鎖長依存性を考慮した補正が行われてもよい。これは、論文１に基づく粗視化分子モデル６の密度が一定であるのに対して、高分子鎖４等の密度は鎖長（モノマー数）Ｎ_FAに依存するため、このような鎖長依存性が考慮されないと、長さの単位換算が行われた際に、密度に乖離が生じるおそれがあるためである。このような依存性を考慮した補正を行う場合には、例えば、高分子鎖４等の密度の鎖長依存性を、下記の式（４）で定義することができる（図１８（全原子モデルの密度と、モノマー数Ｎ_FAとの関係を示すグラフ）に示す）。

ここで、
ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ρ_0、ρ₁：ρ_FAのフィッティングパラメータ

なお、上記の式（４）を用いる場合、高分子鎖４等のモノマー数（鎖長）Ｎ_FAに対応する粗視化分子モデルの粒子数（鎖長）Ｎ_CGを求めるには、上述の式（８）を用いて表されるＮ_CG＝Ｎ_FA・ｌ_FA／（Ｃ_l・ｌ_CG）に代えて、下記の式（５）が用いられるのが望ましい。

ここで、
Ｎ_CG：モノマー数がＮ_FAの分子鎖に対応する粗視化分子モデルの分子鎖１本の粒子数
Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ｌ_FA：モノマー１つあたりの長さ
Ｃ_l：長さの単位換算定数
ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
ρ₀：フィッティングパラメータ

上述の式（８）に代えて、上記の式（５）が用いられることにより、この実施形態では、高分子鎖４等の密度の鎖長依存性を補正しながら、高分子鎖４等と粗視化分子モデル６の間で鎖長を対応付けて、第１物性値（実効的な密度ｌ_k／ｐ）及び第２物性値（Kuhnセグメント数Ｎ_k）を比較することができる。なお、上記の式（５）は、パラメータ決定工程Ｓ２３の第１工程Ｓ５３及び第２工程Ｓ５４において、高分子鎖４等の鎖長（モノマー数）Ｎ_FAと対応する粗視化分子モデル６の鎖長（粒子数）Ｎ_CGを算出する際に用いられてよい。なお、高分子鎖４（図２に示す）のモノマー数（鎖長）Ｎ_FAと、粗視化分子モデル６の粒子数（鎖長）Ｎ_CGとの比Ｃ_m（１個の粒子７あたりのモノマー数Ｎ_FA）については、上記の式（５）で定義される長さの単位換算定数Ｃ_lを用いて、Ｃ_m＝Ｎ_FA／Ｎ_CG（Ｎ_FA）＝Ｃ_l・ｌ_CG／（ｌ_FA・（（ρ₀＋ρ₁／Ｎ_FA）／ρ₀）^１／３）で定義することができる。

次に、図３に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、相互作用パラメータＫを含む相互作用（角度ポテンシャル）Ｅ（θ）が定義された粗視化分子モデル６（図４及び図１３に示す）を用いて、分子動力学に基づく構造緩和を計算する（工程Ｓ３）。本実施形態の工程Ｓ３では、先ず、相互作用定義工程Ｓ２で求められた相互作用パラメータＫに基づいて設定される相互作用（角度ポテンシャル）Ｅ（θ）が、粗視化分子モデル６に定義される。次に、本実施形態の工程Ｓ３では、上記の相互作用が定義された複数の粗視化分子モデル６を、空間１６に配置した高分子材料モデル２０（図１３に示す）が設定される。粗視化分子モデル６の本数については、適宜設定することができる。

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

本実施形態のシミュレーション方法では、構造緩和の計算時において、粗視化分子モデル６の粒子７、７間に定義された上記の相互作用（角度ポテンシャル）Ｅ（θ）により、粗視化分子モデル６の実効的な密度ｌ_k／ｐを、高分子鎖４の実効的な密度ｌ_k／ｐに近づけることができる。また、本実施形態のシミュレーション方法では、構造緩和計算後の高分子材料モデル２０の座標と、上述の長さの単位換算定数Ｃ_lとに基づいて、粗視化分子モデル６の粒子７を、高分子鎖４のモノマーに置き換えることができる。このような置き換えによって得られる全原子モデル１１が配置された高分子材料モデル１８の密度は、実際の高分子鎖４の密度に近似する。

本実施形態のシミュレーション方法では、粗視化分子モデル６が配置された高分子材料モデル２０を用いて、実際の高分子鎖４に近似する構造緩和が可能となる。このため、本実施形態のシミュレーション方法では、構造緩和計算後の高分子材料モデル１８の密度を、実際の高分子鎖４の密度に近似させることができるため、計算精度を向上させることができる。

本実施形態の粗視化分子モデル６は、上述の長さの単位換算定数Ｃ_lに基づいて、高分子鎖４のモノマーＣ_m個が１個の粒子７に置換されているため、粗視化分子モデル６の長さの単位を、高分子鎖４の長さの単位に精度よく換算することができる。したがって、本実施形態のシミュレーション方法は、計算精度を向上させることができる。

上記の相互作用（角度ポテンシャル）Ｅ（θ）が定義された粗視化分子モデル６は、例えば、特許文献１の時間換算工程に用いられてもよい。これにより、上述の長さの単位換算定数Ｃ_lに基づいて、粗視化分子モデル６の粒子７が高分子鎖４のモノマーに置き換えられて、高分子材料モデル１８が作成されることにより、時間単位の換算に必要な長鎖長の全原子モデルの構造緩和した座標を高い精度で作成できる。したがって、粗視化分子モデル６の時間単位を、高分子鎖４の時間単位により精度よく換算することができる。

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

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

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

図３に示した処理手順に従って、図２に示した高分子鎖（cis-１．４ポリブタジエン）をモデル化した粗視化分子モデルがコンピュータに入力された。そして、粗視化分子モデルに、上記式（１）で定義される相互作用が定義された（実施例１、実施例２及び実施例３）。

実施例１～３の相互作用を定義する工程では、図７、図８、図１２及び図１４に示した処理手順に基づいて、鎖長が異なる複数の全原子モデルと、鎖長が異なる複数の粗視化分子モデルとの間において、第１物性値（実効的な密度ｌ_k／ｐ）及び第２物性値（Kuhnセグメント数Ｎ_k）が近づくように、高分子鎖と粗視化分子モデルとの間の長さの単位換算定数Ｃ_l、及び、相互作用パラメータＫが求められた。

なお、実施例１～３において、粗視化分子モデルの実効的な密度ｌ_k／ｐ及び粗視化分子モデルのKuhnセグメント数Ｎ_kは、上記の式（２）及び上記の式（３）の関数にそれぞれ近似させた。実施例１～２では、長さの単位換算定数Ｃ_l及び相互作用パラメータＫを１つずつ交互に非線形最小二乗法で更新しながら、２つのパラメータの双方が収束するまで繰り返し実施された。一方、実施例３では、差分進化法に基づいて、２つのパラメータが探索された。

また、実施例２では、高分子鎖４等の密度の鎖長依存性を、上記の式（４）で定義し、また、モノマー数に対応するビーズ数を、下記式（５）で定義した。全原子モデル、及び、粗視化モデルについては、次のとおりである。
全原子モデル：
第１全原子モデル：
モノマー数：２０、空間内の本数：５０本
第２全原子モデル：
モノマー数：３０、空間内の本数：７５本
第３全原子モデル：
モノマー数：４０、空間内の本数：７５本
粗視化分子モデル：
第１粗視化分子モデル～第５粗視化分子モデル：
粒子数：１０、空間内の本数：１００００本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。
第６粗視化分子モデル～第１０粗視化分子モデル：
粒子数：２０、空間内の本数：５０００本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。
第１１粗視化分子モデル～第１５粗視化分子モデル：
粒子数：５０、空間内の本数：２０００本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。
第１６粗視化分子モデル～第２０粗視化分子モデル：
粒子数：１００、空間内の本数：１０００本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。
第２１粗視化分子モデル～第２５粗視化分子モデル：
粒子数：１５０、空間内の本数：６７７本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。
第２６粗視化分子モデル～第３０粗視化分子モデル：
粒子数：２００、空間内の本数：５００本
相互作用パラメータ：０、０．５、１．０、１．５及び２．０で互いに異なる。

実施例１～３では、２５０モノマーに相当する鎖長を有し、かつ、相互作用が定義された粗視化分子モデルを、空間内に１００本配置された高分子材料モデルが定義された。そして、実施例１～３では、空間内に配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。

比較のために、図２に示した高分子鎖（cis-１．４ポリブタジエン）をモデル化した粗視化分子モデルがコンピュータに入力された（比較例）。比較例の粗視化分子モデルには、上記式（１）で定義される相互作用が定義されたが、図１４に示したパラメータ決定工程Ｓ２３において、工程Ｓ５１から第１工程Ｓ５３までの工程が実施された後、第２工程Ｓ５４以下の処理が省略された。即ち、比較例では、鎖長が異なる複数の全原子モデルと、鎖長が異なる複数の粗視化分子モデルとの間において、長鎖長極限のみを考慮して求められた長さの単位換算定数Ｃ_lと相互作用パラメータＫが用いられており、実施例１～３のように、第１物性値（実効的な密度ｌ_k／ｐ）及び第２物性値（Kuhnセグメント数Ｎ_k）が近づくように、高分子鎖と粗視化分子モデルとの間の長さの単位換算定数Ｃ_l、及び、相互作用パラメータＫが求められていない。そして、比較例の粗視化分子モデルが空間内に１００本配置された高分子材料モデルが定義された。そして、比較例では、空間内に配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。

実施例１～３で計算された密度、及び、比較例で計算された密度が、図２に示した高分子鎖（cis-１．４ポリブタジエン）をモデル化した全原子モデルの平衡密度と比較された。全原子モデルは、構造緩和された粗視化分子モデルから、高分子鎖４のモノマーＣ_m個が１個の粒子７に置換されることで作成された。この平衡密度は、下記の仕様に基づいて、空間内に配置された全原子モデルを対象に構造緩和が計算されることで計算された。
全原子モデルの鎖長：２５０モノマー
空間内の本数：１００本
温度：３６０Ｋ
力場：Ｌ－ＯＰＬＳ

テストの結果、比較例の粗視化分子モデルの密度は、全原子モデルの平衡密度から１０％乖離した。一方、実施例１～３の粗視化分子モデルの密度は、全原子モデルの平衡密度からの乖離を１％未満に抑えることができた。したがって、実施例１～３は、比較例に比べて、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることができた。

また、実施例２では、上記の式（４）及び（５）に基づいて、高分子鎖の密度の鎖長依存性を考慮した補正が行われたため、実施例１及び実施例３に比べて、全原子モデルの平衡密度からの乖離を小さくすることができた。

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

Claims (8)

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

   

   ここで、
   Ｅ：相互作用ポテンシャル関数
   Ｋ：相互作用パラメータ
   θ：隣り合う３つの粒子がなす角度
3. 前記単位換算定数、及び、前記相互作用パラメータを定義する工程は、前記鎖長及び前記相互作用パラメータがそれぞれ異なる複数の前記粗視化分子モデルを定義する工程と、
   複数の前記粗視化分子モデルの前記物性値を、前記鎖長と前記相互作用パラメータとの関数に近似させる工程と、
   前記鎖長と前記相互作用パラメータとの関数を用いて、前記単位換算定数、及び、前記相互作用パラメータを決定する工程とを含む、請求項１又は２記載の高分子材料のシミュレーション方法。
4. 前記鎖長と前記相互作用パラメータとの関数は、前記Kuhn長と前記Packing長との比に近似するものであり、
   前記鎖長と前記相互作用パラメータとの関数は、下記式（２）で定義される、請求項３記載の高分子材料のシミュレーション方法。
   

   

   ここで、
   ｆ_lk/p：Kuhn長とPacking長との比の近似関数
   Ｎ_CG：粗視化分子モデルの鎖長
   Ｋ：相互作用パラメータ
   ｃ₀～ｃ₅：ｆ_lk/pのフィッティングパラメータ
5. 前記鎖長と前記相互作用パラメータとの関数は、前記粗視化分子モデルの全長と前記Kuhn長との比であるKuhnセグメント数に近似するものであり、
   前記鎖長と前記相互作用パラメータとの関数は、下記式（３）で定義される、請求項３又は４記載の高分子材料のシミュレーション方法。
   

   

   ここで、
   ｆ_Nk：Kuhnセグメント数の近似関数
   Ｎ_CG：粗視化分子モデルの鎖長
   Ｋ：相互作用パラメータ
   ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
   ρ：粗視化分子モデルの数密度
   ｆ_lk/p：Kuhn長とPacking長との比の近似関数
6. 前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の一方がそれぞれ近づくように、前記相互作用パラメータを固定して、前記単位換算定数を求める第１工程と、
   鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の他方がそれぞれ近づくように、前記単位換算定数を固定して、前記相互作用パラメータを求める第２工程と、
   前記単位換算定数と前記相互作用パラメータとが収束するまで、前記第１工程と前記第２工程とを交互に実施させる工程とを含む、請求項１ないし５のいずれかに記載の高分子材料のシミュレーション方法。
7. 前記相互作用を定義する工程は、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかの密度の鎖長依存性を考慮して、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかのモノマー数と、前記粗視化分子モデルの粒子数とを対応付ける工程を含む、請求項１ないし６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記密度の鎖長依存性は、下記式（４）で定義され、
   前記モノマー数に対応する前記粒子数は、下記式（５）で定義される、請求項７記載の高分子材料のシミュレーション方法。
   

   

   ここで、
   ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
   Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
   ρ_0、ρ₁：ρ_FAのフィッティングパラメータ
   

   

   ここで、
   Ｎ_CG：モノマー数がＮ_FAの分子鎖に対応する粗視化分子モデルの分子鎖１本の粒子数
   Ｎ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
   ｌ_FA：モノマー１つあたりの長さ
   Ｃ_l：長さの単位換算定数
   ｌ_CG：粗視化分子モデルの粒子間の結合の平衡長
   ρ_FA：高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
   ρ₀：フィッティングパラメータ

Images