JP7275653B2 - 粗視化分子モデルの作成方法 - Google Patents

Description

本発明は、高分子材料の分子鎖を複数の粗視化粒子で表現した粗視化分子モデルを作成するための方法に関する。

近年、ゴム配合の開発のために、高分子材料の反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。この種のシミュレーション方法では、例えば、高分子材料の高分子鎖を、複数の粗視化粒子で表現した粗視化分子モデルが用いられている。これらの粗視化粒子の初期座標は、例えば、モンテカルロ法に基づいて、空間内にランダムに配置される。

特許第５４２１６４９号公報

熱平衡状態の高精度なシミュレーションを行うには、粗視化分子モデルが、熱平衡状態の慣性半径を再現していることが重要である。しかしながら、熱平衡状態の慣性半径を再現せずに粗視化粒子の初期座標を決定した場合、熱平衡状態の慣性半径を再現するためには、分子動力学計算又はモンテカルロ法に基づく粗視化分子モデルの緩和計算を、長時間行う必要がある。

本発明は、以上のような実状に鑑み案出されたもので、熱平衡状態の慣性半径を再現しうる粗視化分子モデルを、緩和計算を長時間行うことなく短時間で作成することができる方法を提供することを主たる目的としている。

本発明は、高分子材料の分子鎖を複数の粗視化粒子で表現した粗視化分子モデルを、コンピュータを用いて作成するための方法であって、前記粗視化分子モデルの一端を構成する１番目の粗視化粒子から前記粗視化分子モデルの他端を構成するｍ番目（ｍ＝５以上の整数）の粗視化粒子までの各位置を、前記コンピュータの仮想空間内に順番に決定しながら配置する工程を含み、前記配置する工程は、ｉ番目（ｉは、正の整数）以降の粗視化粒子の位置を、前記ｉ番目の粗視化粒子と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用のみを考慮して仮決定する工程と、前記考慮した相互作用を用いて、前記高分子材料の熱平衡状態へと緩和される前記ｉ番目の粗視化粒子の位置を計算する工程と、前記計算された位置を前記ｉ番目の粗視化粒子の位置として決定する工程とを含むことを特徴とする。

本発明に係る前記粗視化分子モデルの作成方法において、前記相互作用は、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、及び、前記ｉ－２番目の粗視化粒子の間の結合角に関する第１ポテンシャルを含んでもよい。

本発明に係る前記粗視化分子モデルの作成方法において、前記相互作用は、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、前記ｉ－２番目の粗視化粒子、及び、前記ｉ－３番目の粗視化粒子の間の二面角に関する第２ポテンシャルを含んでもよい。

本発明に係る前記粗視化分子モデルの作成方法において、前記計算する工程に先立ち、下記式（１）に基づいて、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、及び、前記ｉ－２番目の粗視化粒子の間の結合角に関する第３ポテンシャルを定義する工程をさらに含んでもよい。

ここで、
Ｅ_i（θ）：第３ポテンシャル
Ｐ_i（θ）：熱平衡状態の結合角の頻度分布
Ｐ_i,local（θ）：ｉ番目の粗視化粒子と、ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布
θ：結合角
k_B：ボルツマン定数
Ｔ：温度

本発明に係る前記粗視化分子モデルの作成方法において、前記配置する工程の後に、前記粗視化粒子の重なりを許容するソフトポテンシャルを、前記粗視化粒子間に定義して、前記粗視化分子モデルの構造緩和を計算する工程をさらに含んでもよい。

本発明の粗視化分子モデルの作成方法は、ｉ番目（ｉは、正の整数）以降の粗視化粒子の位置を、前記ｉ番目の粗視化粒子と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用のみを考慮して仮決定する工程と、前記考慮した相互作用を用いて、前記高分子材料の熱平衡状態へと緩和される前記ｉ番目の粗視化粒子の位置を計算する工程と、前記計算された位置を前記ｉ番目の粗視化粒子の位置として決定する工程とを含んでいる。

本発明の作成方法では、前記ｉ番目の粗視化粒子と主鎖上で近傍にある前記粗視化粒子間の相互作用（即ち、前記ｉ番目の粗視化粒子と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用）のみを考慮して、前記ｉ番目の粗視化粒子の位置が決定されるため、高分子鎖の熱平衡状態の慣性半径を決定づける主鎖の曲がりやすさを再現することができる。これにより、全ての粗視化粒子が配置されることで、熱平衡状態の慣性半径を再現しうる粗視化分子モデルを作成することができる。したがって、本発明では、従来の方法のように、分子動力学計算に基づく粗視化分子モデルの緩和計算を長時間行う必要がないため、熱平衡状態の慣性半径を再現しうる粗視化分子モデルを、短時間で作成できる。

粗視化分子モデルの作成方法を実施するためのコンピュータの斜視図である。 ポリブタジエンの構造式である。 粗視化分子モデルが配置された仮想空間の一例を示す概念図である。 粗視化分子モデルの一例を示す概念図である。 粗視化分子モデルの部分拡大図である。 （ａ）、（ｂ）は、ポテンシャルが定義される粗視化粒子モデルを示す図である。 本実施形態の粗視化分子モデルの作成方法の処理手順の一例を示すフローチャートである。 初期工程の処理手順の一例を示すフローチャートである。 第１頻度分布計算工程の処理手順の一例を示すフローチャートである。 第２頻度分布計算工程の処理手順の一例を示すフローチャートである。 配置工程の処理手順の一例を示すフローチャートである。 ４番目の粗視化粒子の配置の一例を説明する図である。 計算工程の処理手順の一例を示すフローチャートである。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の粗視化分子モデルの作成方法（以下、単に「作成方法」ということがある。）は、コンピュータを用いて、高分子材料の分子鎖を複数の粗視化粒子で表現した粗視化分子モデルが作成される。

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

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

ポリブタジエンを構成する高分子鎖３は、メチレン基（－ＣＨ_２－）とメチン基（－ＣＨ－）とからなるモノマー５｛－［ＣＨ_２－ＣＨ＝ＣＨ－ＣＨ_２］－｝が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基（－ＣＨ_２）に替えて、メチル基（－ＣＨ_３）が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。

本実施形態の作成方法では、高分子鎖３を表現した粗視化分子モデルが用いられる。図３は、粗視化分子モデル６が配置された仮想空間８の一例を示す概念図である。図４は、粗視化分子モデルの一例を示す概念図である。

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

図５は、粗視化分子モデル６の部分拡大図である。粗視化分子モデル６の粗視化粒子９、９間には、平衡長が設定された第４ポテンシャルＰ４が定義される。本実施形態の第４ポテンシャルＰ４は、粗視化分子モデル６がKremer-Grestモデルである場合、下記式（２）で定義されるポテンシャル（以下、「ＬＪポテンシャルＵ_LJ（ｒ_ij）」ということがある。）と、下記式（３）で定義される結合ポテンシャルＵ_FENEとの和で定義される。

ここで、各定数及び変数は、Lennard-Jones及びFENEの各ポテンシャルのパラメータであり、次のとおりである。
ｒ_ij：粗視化粒子間の距離
ｒ_c：カットオフ距離
ｋ：粗視化粒子間のばね定数
ε：粗視化粒子間に定義されるＬＪポテンシャルの強度
σ：粗視化分子モデルの長さの単位
Ｒ₀：伸びきり長
なお、距離ｒ_ij、カットオフ距離ｒ_c、及び、伸びきり長Ｒ₀は、各粗視化粒子９の中心９ｏ、９ｏ間の距離について定義される。

このような第４ポテンシャルＰ４は、粗視化粒子９、９間を、平衡長を設定して連結するボンド１０として定義される。これにより、複数の粗視化粒子９がボンド１０によって伸縮自在に拘束された直鎖状の粗視化分子モデル６を設定することができる。

ＬＪポテンシャルＵ_LJ（ｒ_ij）及び結合ポテンシャルＵ_FENEの各定数及び各変数の値としては、適宜設定することができる。本実施形態では、論文１（ 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）に基づいて適宜設定される。

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

仮想空間８の一辺の長さＬａは、適宜設定することができる。本実施形態の長さＬａは、粗視化分子モデル６の拡がりを示す量である慣性半径（図示省略）の３倍以上が望ましい。これにより、仮想空間８は、後述の粗視化分子動力学計算において、周期境界条件のために周期的に並ぶ自分自身（イメージ）との衝突の発生を防いで、粗視化分子モデル６の空間的拡がりを適切に計算することができる。また、粗視化粒子９の密度は、解析対象の高分子材料に対応するように適宜設定される。例えばKremer-Grestモデルにおいて、粗視化粒子９の密度は、０．８５個／σ³である。粗視化分子モデル６の本数、及び、一辺の長さＬａについては、一辺の長さＬａが慣性半径より十分大きく（３倍以上）なるように、かつ、粗視化粒子９の密度が適切となるように調整されるのが望ましい。これにより、仮想空間８は、解析対象の高分子材料の少なくとも一部の空間領域を定義することができる。このような仮想空間８に、粗視化分子モデル６が配置されることにより、高分子材料をモデル化した高分子材料モデル１２が設定される。仮想空間８は、コンピュータ１に記憶される。

図５に示されるように、本実施形態の作成方法では、粗視化分子モデル６の一端を構成する１番目の粗視化粒子９ａから粗視化分子モデル６の他端を構成するｍ番目（ｍ＝５以上の整数）の粗視化粒子９ｍまでの各位置を、コンピュータ１の仮想空間８（図３に示す）内に順番に決定しながら配置している。粗視化粒子９の個数ｍについては、例えば、高分子材料の構造に応じて適宜設定される。本実施形態では、ｍを大きな整数に設定することで、大きな鎖長を有する粗視化分子モデル６を作成することができる。

ところで、従来の作成方法では、粗視化分子モデル６を作成する際に、熱平衡状態の慣性半径を定量的に考慮せずに、粗視化粒子９を、例えばランダムに配置していた。従って、熱平衡状態の慣性半径を再現した粗視化分子モデル６を作成するためには、粗視化粒子９を配置した後に、分子動力学計算に基づく緩和計算を、長時間行う必要があった。

発明者らは、鋭意研究を重ねた結果、ｉ番目（ｉは、正の整数）の粗視化粒子９_iについて、主鎖上の近傍（その粗視化粒子９_iの直前）のｋ個の粗視化粒子９（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用のみを考慮して、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置を計算して配置することで、熱平衡状態の慣性半径を再現しうる粗視化分子モデル６を短時間で作成しうることを知見した。

本実施形態の作成方法では、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置を計算するために、上述した第４ポテンシャルＰ４が定義される。さらに、必要に応じて、ｉ番目の粗視化粒子９_iに対して、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、及び、第３ポテンシャルＥ_i（θ）が定義される。なお、第１ポテンシャルＵ_angle（θ）、及び、第２ポテンシャルＵ_torsion（φ）は、例えば、通常のKremer-Grestモデルを使用する場合に省略することができる。なお、第３ポテンシャルＥ_i（θ）の定義が省略される場合については、後述する。図６（ａ）、（ｂ）は、ポテンシャルが定義される粗視化粒子モデルを示す図である。

図６（ａ）に示されるように、第１ポテンシャルＵ_angle（θ）は、ｉ番目の粗視化粒子９_i、ｉ－１番目の粗視化粒子９_i-1、及び、ｉ－２番目の粗視化粒子９_i-2の間の結合角θに関するものである。図６（ｂ）に示されるように、第２ポテンシャルＵ_torsion（φ）は、ｉ番目の粗視化粒子９_i、ｉ－１番目の粗視化粒子９_i-1、ｉ－２番目の粗視化粒子９_i-2、及び、ｉ－３番目の粗視化粒子９_i-3の間の二面角φに関するものである。第１ポテンシャルＵ_angle（θ）及び第２ポテンシャルＵ_torsion（φ）は、シミュレーションに用いられる粗視化粒子９、及び、粗視化分子モデル６に応じて適宜設定される。

上記の第１ポテンシャルＵ_angle（θ）及び第２ポテンシャルＵ_torsion（φ）は、図２に示した高分子鎖３の立体配座に基づくものである。このような第１ポテンシャルＵ_angle（θ）及び第２ポテンシャルＵ_torsion（φ）がｉ番目の粗視化粒子９_iに定義されることにより、高分子鎖３の立体配座を考慮した配置が可能となる。従って、第１ポテンシャルＵ_angle（θ）及び第２ポテンシャルＵ_torsion（φ）は、熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置の計算に役立つ。

第３ポテンシャルＥ_i（θ）は、下記式（１）に基づいて定義される。この第３ポテンシャルＥ_i（θ）は、平均場近似ポテンシャルである。平均場近似ポテンシャルは、ｉ番目の粗視化粒子と、ｉ番目の粗視化粒子９_iと主鎖上で近傍にあるｋ個の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）のみとの相互作用を、ｉ番目の粗視化粒子９_iに作用する全ての相互作用から除いた相互作用に基づくものである。結合角θは、図６（ａ）に示されるように、ｉ番目の粗視化粒子９_i、ｉ－１番目の粗視化粒子９_i-1、及び、ｉ－２番目の粗視化粒子９_i-2の間の角度として定義される。この第３ポテンシャルＥ_i（θ）は、ｉ番目の粗視化粒子９_iと主鎖上で近傍の相互作用を除いた周囲の相互作用（溶媒和など）の全体が、結合角θへ与える影響を表しており、高分子鎖の置かれている環境（例えば、密度、温度など）に依存して変化し得る。第３ポテンシャルＥ_i（θ）の寄与が常に十分小さいと考えられる場合（例えば、粗視化粒子間の非結合相互作用が斥力（または、ほぼ斥力）で、密度が希薄であるような系）については、第３ポテンシャルＥ_i（θ）の定義を省略（Ｅ_i（θ）＝０）してもよい。

ここで、
Ｅ_i（θ）：第３ポテンシャル
Ｐ_i（θ）：熱平衡状態の結合角の頻度分布
Ｐ_i,local（θ）：ｉ番目の粗視化粒子と、ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布
θ：結合角
k_B：ボルツマン定数
Ｔ：温度

上記式（１）の第３ポテンシャルＥ_i（θ）は、ボルツマン因子exp（－Ｅ_i（θ）／k_BＴ）が、２つの結合角の頻度分布の比Ｐ_i（θ）／Ｐ_i,local（θ）と等しくなるように定義される。上記式（１）において、温度Ｔは、仮想空間８に設定される温度である。

熱平衡状態の結合角の頻度分布（以下、単に「平衡頻度分布」ということがある。）Ｐ_i（θ）は、適宜設定された温度及び密度に基づいて、高分子材料の熱平衡状態へと緩和された粗視化分子モデル６について、ボンド１０を介して連続する３つの粗視化粒子９がなす結合角θ（図６（ａ）に示す）と、その頻度との関係を示すものである。

一方、ｉ番目の粗視化粒子９_iと、ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９_i-k～９_i-1のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布（以下、単に「部分的平衡頻度分布」ということがある。）Ｐ_i,local（θ）は、ｉ番目の粗視化粒子９_iと、ｉ番目の粗視化粒子９_iの主鎖上の近傍の粗視化粒子９_i-k～９_i-1との間の相互作用（第４ポテンシャルＰ４、並びに、必要に応じて定義される第１ポテンシャル、及び、第２ポテンシャルを含む）に基づいて熱平衡状態に緩和するように配置された粗視化分子モデルについて、ボンド１０を介して連続する３つの粗視化粒子９がなす結合角θ（図６（ａ）に示す）と、その頻度との関係を示すものである。

平衡頻度分布Ｐ_i（θ）を部分的平衡頻度分布Ｐ_i,local（θ）で除した比Ｐ_i（θ）／Ｐ_i,local（θ）は、部分的に熱平衡状態にある結合角θの頻度分布を、熱平衡状態の結合角θの頻度分布にするための重みを示している。比Ｐ_i（θ）／Ｐ_i,local（θ）は、結合角θごとに計算され、結合角θの区間は、適宜設定される。この比Ｐ_i（θ）／Ｐ_i,local（θ）の自然対数が、ボルツマン定数k_B及び温度Ｔに乗じられることで、第３ポテンシャルＥ_i（θ）が定義される。このような第３ポテンシャルＥ_i（θ）は、図６（ａ）に示したｉ番目の粗視化粒子９_iが形成する結合角θを、部分的な熱平衡状態から熱平衡状態に調整することができる。

第３ポテンシャルＥ_i（θ）は、粗視化分子モデル６の相互作用の定義（第４ポテンシャルＰ４、第１ポテンシャル、及び、第２ポテンシャルなど）、及び、作成すべき粗視化分子モデル６が置かれる仮想空間８の環境（温度、密度など）に依存する定数関数である。従って、本実施形態の作成方法では、粗視化粒子９を仮想空間８に配置する工程に先立ち、第３ポテンシャルＥ_i（θ）が求められる。図７は、本実施形態の粗視化分子モデルの作成方法の処理手順の一例を示すフローチャートである。

本実施形態の作成方法では、第３ポテンシャルＥ_i（θ）が求められる（初期工程Ｓ１）。上述したように、粗視化粒子間の非結合相互作用が斥力（または、ほぼ斥力）で、密度が希薄であるような系である場合には、初期工程Ｓ１が省略（第３ポテンシャルＥ_i（θ）＝０と定義）されてもよい。本実施形態の初期工程Ｓ１では、第３ポテンシャルＥ_i（θ）を求めるために、未知の定数である平衡頻度分布Ｐ_i（θ）、及び、部分的平衡頻度分布Ｐ_i,local（θ）がそれぞれ取得される。図８は、初期工程Ｓ１の処理手順の一例を示すフローチャートである。

本実施形態の初期工程Ｓ１では、先ず、熱平衡状態の結合角の頻度分布Ｐ_i（θ）が求められる（第１頻度分布計算工程Ｓ１１）。図９は、第１頻度分布計算工程Ｓ１１の処理手順の一例を示すフローチャートである。

本実施形態の第１頻度分布計算工程Ｓ１１では、先ず、図３に示されるように、粗視化分子モデル６が、仮想空間８に配置される（工程Ｓ６１）。工程Ｓ６１では、粗視化分子モデル６の一端を構成する１番目の粗視化粒子９ａから、粗視化分子モデル６の他端を構成するｍ番目の粗視化粒子９ｍまでの各位置を、仮想空間８内に順番に決定しながら配置している。

なお、第１頻度分布計算工程Ｓ１１で定義される粗視化分子モデル６及び仮想空間８は、熱平衡状態の結合角の頻度分布Ｐ_i（θ）を求めるためだけに用いられるものである。従って、粗視化分子モデル６及び仮想空間８は、本実施形態の作成方法で求められる粗視化分子モデル６及び仮想空間８とは異なるものであり、独立して定義される。

第１頻度分布計算工程Ｓ１１で定義される粗視化粒子９の各位置は、他の粗視化粒子９との衝突が避けられることなく、モンテカルロ法に基づいてランダムに決定されてもよい。なお、粗視化粒子９の各位置は、他の粗視化粒子９との衝突を避けて配置されてもよい。また、粗視化粒子９の各位置をランダムに配置するのに代えて、例えば、図１１に示す後述の配置工程Ｓ２の処理手順のうち、第３ポテンシャルを設定する工程Ｓ２６を省略（即ち、Ｅ_i（θ）=０とみなす）した手順に基づいて、粗視化粒子９が配置されてもよい。

図５に示されるように、隣接する粗視化粒子９、９の間は、第４ポテンシャルＰ４に基づくボンド１０で結合されている。各粗視化粒子９には、第３ポテンシャルＥ_i（θ）は定義されない。本実施形態の工程Ｓ６１では、熱平衡状態の結合角の頻度分布Ｐ_i（θ）を求めるために作成すべき粗視化分子モデル６の環境（温度、密度など）に応じて、系（仮想空間８）の体積、及び、粗視化分子モデル６の本数が適宜設定される。仮想空間８に配置された１本の粗視化分子モデル６は、コンピュータ１に記憶される。

次に、本実施形態の第１頻度分布計算工程Ｓ１１では、図３に示した仮想空間８に配置された粗視化分子モデル６を対象に、分子動力学（ Molecular Dynamics : ＭＤ ）計算に基づいた構造緩和が計算される（工程Ｓ６２）。本実施形態の分子動力学計算では、例えば、仮想空間８について所定の時間、粗視化分子モデル６が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粗視化粒子９の動きが、単位時間毎に追跡される。

本実施形態の構造緩和の計算は、仮想空間８において、圧力及び温度が一定、又は、体積及び温度が一定に保たれる。工程Ｓ６２では、粗視化分子モデル６の初期配置を十分に緩和されるまで、シミュレーションの単位時間毎に計算される。これにより、工程Ｓ６２では、実際の高分子材料の分子運動に近似させて、粗視化分子モデル６の初期配置を精度よく緩和することができ、熱平衡状態の粗視化分子モデル６を計算することができる。

本実施形態の工程Ｓ６２では、末端の影響が無視できる程度（典型的には絡み合い点間距離の２倍程度）の鎖長の粗視化分子モデル６の構造緩和が行われる。これにより、第１頻度分布計算工程Ｓ１１では、実際の高分子材料に対応するような長い鎖長の粗視化分子モデル６を配置する必要がない。さらに、そのような長い鎖長の粗視化分子モデル６について、慣性半径の３倍程度の周期境界長とするために、非常に多数の粗視化分子モデル６を配置する必要もない。このため、本実施形態の第１頻度分布計算工程Ｓ１１では、計算時間を大幅に短縮しうる。上記の構造緩和の計算は、例えば（株）ＪＳＯＬ社製のソフトマテリアル総合シミュレーター（Ｊ－ＯＣＴＡ）に含まれるＣＯＧＮＡＣを用いて処理することができる。なお、工程Ｓ６２では、分子動力学計算に代えて、メトロポリス法などのモンテカルロ法を用いて、熱平衡状態の粗視化分子モデル６が求められてもよい。熱平衡状態の粗視化分子モデル６は、コンピュータ１に記憶される。

次に、本実施形態の第１頻度分布計算工程Ｓ１１では、熱平衡状態の結合角θの頻度分布が求められる（工程Ｓ６３）。工程Ｓ６３では、先ず、工程Ｓ６２で求められた熱平衡状態の粗視化分子モデル６について、ボンド１０を介して連続する３つの粗視化粒子９の全ての組み合わせ（ペア）の結合角θ（図５に示す）が計算される。そして、各結合角θについて、適宜分割された区間ごとの出現頻度が計算されることにより、熱平衡状態の結合角と、各区間の頻度との関係を示す頻度分布Ｐ_i（θ）が求められる。熱平衡状態の結合角の頻度分布Ｐ_i（θ）は、コンピュータ１に記憶される。

次に、本実施形態の初期工程Ｓ１では、部分的に熱平衡状態にある結合角の頻度分布が求められる（第２頻度分布計算工程Ｓ１２）。図１０は、第２頻度分布計算工程Ｓ１２の処理手順の一例を示すフローチャートである。

本実施形態の第２頻度分布計算工程Ｓ１２では、先ず、図３に示されるように、粗視化分子モデル６が、仮想空間８に配置される（工程Ｓ６６）。工程Ｓ６６では、第１頻度分布計算工程Ｓ１１の工程Ｓ６１と同様に、粗視化分子モデル６の一端を構成する１番目の粗視化粒子９ａから粗視化分子モデル６の他端を構成するｍ番目の粗視化粒子９ｍまでの各位置を、仮想空間８内に順番に決定しながら配置している。

なお、第２頻度分布計算工程Ｓ１２で定義される粗視化分子モデル６及び仮想空間８は、部分的に熱平衡状態にある結合角の頻度分布（部分的平衡頻度分布）Ｐ_i,local（θ）を求めるためだけに用いられるものである。従って、粗視化分子モデル６及び仮想空間８は、第１頻度分布計算工程Ｓ１１で定義される粗視化分子モデル６及び仮想空間８や、本実施形態の作成方法で求められる粗視化分子モデル６及び仮想空間８とは異なるものであり、独立して定義される。

本実施形態の粗視化粒子９の各位置は、例えば、図１１に示す後述の配置工程Ｓ２の処理手順のうち、第３ポテンシャルを設定する工程Ｓ２６を省略（即ち、Ｅ_i（θ）=０とみなして）した手順に基づいて決定することができる。これにより、部分的に熱平衡状態の粗視化分子モデル６を作成することができる。図４に示されるように、粗視化粒子９、９間は、第４ポテンシャルＰ４に基づくボンド１０で結合されている。

第２頻度分布計算工程Ｓ１２で定義される粗視化粒子９には、計算に用いられる粗視化分子モデル６に基づき、必要に応じて、第１ポテンシャルＵ_angle（θ）、及び、第２ポテンシャルＵ_torsion（φ）が定義される。第１ポテンシャルＵ_angle（θ）、及び、第２ポテンシャルＵ_torsion（φ）の定義を省略できる場合については、上述のとおりである。一方、第３ポテンシャルＥ_i（θ）は定義されてない（即ち、Ｅ_i（θ）=０として定義される）。本実施形態の工程Ｓ６６では、１本の粗視化分子モデル６が、仮想空間８に配置されている。粗視化分子モデル６には、分子動力学計算に基づく構造緩和が実施されない。本実施形態では、ｉ番目（ｉは、正の整数）以降の粗視化粒子９iについて、ｉ番目の粗視化粒子９_iの主鎖上の近傍のｋ個（ｋ＝３以上の整数）の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９_i-k～９_i-1）のみとの相互作用（第４ポテンシャルＰ４、及び、必要に応じて定義される第１ポテンシャル、第２ポテンシャルを含む）に基づいて、熱平衡状態に緩和するように配置される。これにより、粗視化分子モデル６を、部分的な相互作用に基づく熱平衡状態にあるものとして扱うことができる。

次に、本実施形態の第２頻度分布計算工程Ｓ１２では、部分的に熱平衡状態にある結合角の頻度分布（部分的平衡頻度分布）Ｐ_i,local（θ）が求められる（工程Ｓ６７）。本実施形態の工程Ｓ６７では、先ず、部分的に熱平衡状態にある粗視化分子モデル６について、ボンド１０を介して連続する３つの粗視化粒子９の全ての組み合わせ（ペア）の結合角θが計算される。そして、各結合角θについて、適宜分割された区間ごとの出現頻度が計算されることにより、部分的に熱平衡状態にある結合角の頻度分布Ｐ_i,local（θ）が求められる。部分的平衡頻度分布Ｐ_i,local（θ）は、コンピュータ１に記憶される。

次に、本実施形態の初期工程Ｓ１では、熱平衡状態の結合角の頻度分布Ｐ_i（θ）、及び、部分的に熱平衡状態にある結合角の頻度分布Ｐ_i,local（θ）に基づいて、第３ポテンシャルＥ_i（θ）が求められる（工程Ｓ１３）。工程Ｓ１３では、第３ポテンシャルＥ_i（θ）を定義する上記式（１）に、熱平衡状態の結合角の頻度分布Ｐ_i（θ）、及び、部分的に熱平衡状態にある結合角の頻度分布Ｐ_i,local（θ）が代入される。これにより、初期工程Ｓ１では、第３ポテンシャルＥ_i（θ）を求めることができる。第３ポテンシャルＥ_i（θ）は、コンピュータ１に記憶される。

次に、本実施形態の作成方法では、コンピュータ１が、粗視化分子モデル６を構成する粗視化粒子９の各位置を、仮想空間８内に配置する（配置工程Ｓ２）。本実施形態の配置工程Ｓ２は、粗視化分子モデル６の一端を構成する１番目の粗視化粒子９ａから粗視化分子モデル６の他端を構成するｍ番目の粗視化粒子９ｍまでの各位置を、仮想空間８内に順番に決定しながら配置している。図１１は、配置工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態の配置工程Ｓ２は、先ず、図５に示した１番目の粗視化粒子９ａが、仮想空間８内に配置される（工程Ｓ２１）。本実施形態の工程Ｓ２１では、１番目の粗視化粒子９ａ位置が、仮想空間８内に既に配置されている他の粗視化粒子９との衝突が避けられることなくランダムに配置される。

なお、粗視化粒子９以外の粒子（例えば、高分子材料に含まれるフィラーをモデル化したフィラー粒子等）が仮想空間８に配置されている場合、その粒子との衝突を避けた位置に、１番目の粗視化粒子９ａの位置が決定されるのが望ましい。これにより、粗視化粒子９が他の粒子（図示省略）に拘束されるのを防ぐことができる。１番目の粗視化粒子９ａの位置は、コンピュータ１に記憶される。

次に、本実施形態の配置工程Ｓ２では、工程Ｓ２２以降において、ｉ番目以降（２番目以降）の粗視化粒子９が、仮想空間８に配置される。本実施形態の配置工程Ｓ２の工程Ｓ２２以降では、２番目の粗視化粒子９ｂからｍ番目の粗視化粒子９ｍまでの各位置を、仮想空間８内に順番に決定しながら配置している。図１２（ａ）～（ｃ）は、４番目の粗視化粒子９ｄの配置の一例を説明する図である。

本実施形態の配置工程Ｓ２では、先ず、図１２（ａ）に示されるように、ｉ番目の粗視化粒子９_i（例えば、図１２では、４番目の粗視化粒子９ｄ）の位置が仮決定される（工程Ｓ２２）。工程Ｓ２２では、先ず、粗視化粒子９に第４ポテンシャルＰ４が定義される。さらに、工程Ｓ２２では、必要に応じて、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、及び、第３ポテンシャルＥ_i（θ）が、ｉ番目の粗視化粒子９_iに対して、直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との間にのみ定義される。ただし、位置の仮決定の際には、計算の効率化のために、これらの相互作用をそれぞれ単純化した形で用いて良い。

第４ポテンシャルＰ４のＬＪポテンシャルＵ_LJ（ｒ_ij）については、例えば、距離r_min未満で∞、それ以上で０となる関数に単純化することができる。これにより、相互作用を考慮する２つの粗視化粒子間の距離が常に距離r_min以上となるような粗視化粒子９の位置を選ぶことができる。第４ポテンシャルＰ４の結合ポテンシャルＵ_FENEについては、例えば、平衡核間距離Ｌ１で０、それ以外で∞となる関数に単純化することができる。これにより、結合長が常に平衡核間距離Ｌ１となるような粗視化粒子９の位置を選ぶことができる。

第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、第３ポテンシャルＥ_i（θ）については、例えば、常に０となる関数に単純化することができる。これにより、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、第３ポテンシャルＥ_i（θ）を考慮せずに、粗視化粒子９の位置を仮決定することができる。

上記のような相互作用の単純化により、２番目以降のｉ番目の粗視化粒子９_iの位置は、ｉ－１番目の粗視化粒子の位置に対する方向をランダムに選び、距離Ｌ１進んだ位置として仮決定することができる。ただし、仮決定した位置が、主鎖上近傍のｋ個の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）９のいずれかと衝突する（距離がr_minよりも小さい）ときは、ｉ－１番目の粗視化粒子の位置に対する方向を選び直すのが望ましい。粗視化粒子９_iのランダムな位置については、平衡核間距離Ｌ１を半径とする球面上の点について、一様乱数となるように決定されるのが望ましい。なお、粗視化分子モデル６がKremer-Grestモデルの場合、一様乱数の範囲は、平衡核間距離Ｌ１～0.965σ、及び、距離r_min～0.9σに設定されるのが望ましい。

工程Ｓ２２の手順を具体的に説明すると、２番目の粗視化粒子９ｂの位置は、１番目の粗視化粒子９ａとの結合距離を平衡核間距離Ｌ１に固定して、１番目の粗視化粒子９との衝突を避けたランダムな方向に仮決定される。３番目の粗視化粒子９ｃの位置は、１番目の粗視化粒子９ａ及び２番目の粗視化粒子９ｂとの衝突を避けたランダムな方向に仮決定される。

４番目以降の粗視化粒子９ｄについては、ｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）の位置が、ｉ－２番目の粗視化粒子９_i-2（例えば、２番目の粗視化粒子９ｂ）と、ｉ－３番目の粗視化粒子９_i-3（例えば、１番目の粗視化粒子９ａ）との相互作用のみを考慮して、それらとの衝突を避けたランダムな方向に仮決定される。

ｉ番目の粗視化粒子９_iの位置が仮決定された後、ｉ番目の粗視化粒子９_iとｉ－１番目の粗視化粒子９_i-1との間は、ボンド１０で結合される。仮決定されたｉ番目の粗視化粒子９_iの位置は、コンピュータ１に記憶される。

このように、本実施形態の作成方法では、それぞれの粗視化粒子９について、その直前のｋ個の粗視化粒子９との相互作用のみを考慮して、それらとの衝突を避けた位置に、粗視化粒子９の位置を仮決定しているため、例えば、全ての粗視化粒子９との衝突を避けた位置に仮決定される場合に比べて、粗視化粒子９の位置を短時間で仮決定することができる。

また、粗視化粒子９以外の粒子（例えば、フィラー粒子（図示省略））が仮想空間８に配置されている場合、その粒子との衝突を避けた位置に、ｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）の位置が決定されるのが望ましい。なお、粗視化粒子９以外の粒子の配置によっては、例えば袋小路のような微小な密閉空間が存在すると、次の粗視化粒子９の位置が決定できなくなることがあり得る。そのような場合、配置工程Ｓ２が最初から実施されるのが望ましい。

次に、本実施形形態の配置工程Ｓ２では、仮決定されたｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９_i-k～９_i-1）との相互作用を考慮して、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置が計算される（計算工程Ｓ２３）。計算工程Ｓ２３では、メトロポリス法に基づいて、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置が決定される。図１３は、計算工程Ｓ２３の処理手順の一例を示すフローチャートである。

本実施形態の計算工程Ｓ２３では、先ず、仮決定されたｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）のポテンシャルエネルギーｅが計算される（工程Ｓ９１）。ポテンシャルエネルギーｅは、ｉ番目の粗視化粒子９_iと主鎖と、ｉ番目の粗視化粒子９_iの主鎖上の近傍のｋ個（ｋ＝３以上の整数）の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９i-k～９i-1）との第４ポテンシャルＰ４が含まれる。なお、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、及び、第３ポテンシャルＥ_i（θ）が定義されている場合には、第４ポテンシャルＰ４、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、及び、第３ポテンシャルＥ_i（θ）を含めたポテンシャルの和である。

第１ポテンシャルＵ_angle（θ）及び第３ポテンシャルＥ_i（θ）は、３番目以降の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）、ｉ－１番目の粗視化粒子９_i-1（例えば、３番目の粗視化粒子９ｃ）、及び、ｉ－２番目の粗視化粒子９_i-2（例えば、２番目の粗視化粒子９ｂ）の間（組）に設定されたものである。従って、工程Ｓ９１において、３番目以降の粗視化粒子９_iが含まれない組は、第１ポテンシャルＵ_angle（θ）及び第３ポテンシャルＥ_i（θ）が考慮されない。同様に、i番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）が含まれない組の粗視化粒子９については、第４ポテンシャルＰ４も考慮されない。

第２ポテンシャルＵ_torsion（φ）は、４番目以降の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）、ｉ－１番目の粗視化粒子９_i-1（例えば、３番目の粗視化粒子９ｃ）、ｉ－２番目の粗視化粒子９_i-2（例えば、２番目の粗視化粒子９ｂ）、及び、ｉ－３番目の粗視化粒子９_i-3（例えば、１番目の粗視化粒子９ａ）の間（組）に設定されたものである。従って、工程Ｓ９１において、４番目以降の粗視化粒子９_iが含まれない組は、第２ポテンシャルＵ_torsion（φ）が考慮されない。

このように、工程Ｓ９１では、ポテンシャルエネルギーｅに、第３ポテンシャルＥ_i（θ）を加算することができる。これにより、工程Ｓ９１では、ｉ番目の粗視化粒子９_iと、ｉ番目の粗視化粒子９_iの主鎖上の近傍のｋ個（ｋ＝３以上の整数）の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９i-k～９i-1）との相互作用を、第３ポテンシャルＥ_i（θ）によって補正して、結合角θの分布が熱平衡状態に近似するように調整されたポテンシャルエネルギーｅを計算することができる。ポテンシャルエネルギーｅは、コンピュータ１に記憶される。

次に、本実施形態の計算工程Ｓ２３では、図１２（ｂ）に示されるように、仮決定されたｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）とは別に、ｉ番目の粗視化粒子９_iの移動候補位置２０が設定される（工程Ｓ９２）。工程Ｓ９２では、仮決定する工程Ｓ２２と同様の手順に基づいて、ｉ番目の粗視化粒子９_iの移動候補位置２０が、ｉ－２番目の粗視化粒子９_i-2（例えば、２番目の粗視化粒子９ｂ）、及び、９_i-3（１番目の粗視化粒子９ａ）のみを考慮し、かつ、それらとの衝突を避けたランダムな位置に設定される。

なお、粗視化粒子９以外の粒子（例えば、フィラー粒子（図示省略））が仮想空間８に配置されている場合、その粒子との衝突を避けたランダムな位置に、ｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）の移動候補位置２０が決定されるのが望ましい。

次に、本実施形態の計算工程Ｓ２３では、２番目以降の粗視化粒子９について、移動候補位置２０に配置された場合のｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）のポテンシャルエネルギーｅ’が計算される（工程Ｓ９３）。このポテンシャルエネルギーｅ’は、移動候補位置２０に配置されたｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）と、ｉ番目の粗視化粒子９_iの主鎖上の近傍のｋ個（ｋ＝３以上の整数）の粗視化粒子９（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子９i-k～９i-1）の第４ポテンシャルＰ４が含まれる。なお、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）及び第３ポテンシャルＥ_i（θ）が定義されていれば、ポテンシャルエネルギーｅ’は、第４ポテンシャルＰ４、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）及び第３ポテンシャルＥ_i（θ）を含めたポテンシャルの和である。従って、工程Ｓ９３では、結合角θの分布が熱平衡状態に近似するように調整されたポテンシャルエネルギーｅ’を計算することができる。ポテンシャルエネルギーｅ’は、コンピュータ１に記憶される。

次に、本実施形態の計算工程Ｓ２３では、一様乱数Ｒを生成させる（工程Ｓ９４）。工程Ｓ９４では、一般的な疑似乱数生成アルゴリズム（pseudo random number generator）を用いて、０以上１未満の実数の一様乱数Ｒを生成させている。一様乱数Ｒは、コンピュータ１に記憶される。

次に、本実施形態の計算工程Ｓ２３では、メトロポリス法に基づいて、仮決定されたｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）のポテンシャルエネルギーｅと、移動候補位置２０のポテンシャルエネルギーｅ’との差（変化）についての確率exp［（ｅ－ｅ’）／k_BＴ］が、工程Ｓ９４で生成された一様乱数Ｒの値よりも大きいか否かが判断される（工程Ｓ９５）。確率exp［（ｅ－ｅ’）／k_BＴ］において、k_Bはボルツマン定数であり、Ｔは仮想空間８の温度である。

工程Ｓ９５において、確率exp［（ｅ－ｅ’）／k_BＴ］が、一様乱数Ｒの値よりも大きい場合（工程Ｓ９５において、「Ｙ」）、次の工程Ｓ９６において、ｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）の移動候補位置２０が、熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの仮決定された位置（以下、「仮決定位置」ということがある。）として更新（採択）される。移動候補位置２０のポテンシャルエネルギーｅ’は、ｉ番目の粗視化粒子９_iのポテンシャルエネルギーｅとして更新（採択）される。更新（採択）された仮決定位置２２、及び、更新されたポテンシャルエネルギーeは、コンピュータ１に記憶される。そして、次の工程Ｓ９７が行われる。

一方、工程Ｓ９５において、確率exp［（ｅ－ｅ’）／k_BＴ］が、一様乱数Ｒの値以下である場合（工程Ｓ９５において、「Ｎ」）、その移動候補位置２０は棄却され、次の工程Ｓ９７が実施される。

工程Ｓ９７では、工程Ｓ９２～工程Ｓ９６が、予め定められた回数繰り返されたか否かが判断される。工程Ｓ９２～工程Ｓ９６が繰り返される回数については、適宜設定することができる。本実施形態では、計算コストと計算精度とのバランスを考慮して、５～２０回程度に設定される。

工程Ｓ９７において、工程Ｓ９２～工程Ｓ９６が、予め定められた回数繰り返されたと判断された場合（工程Ｓ９７において、「Ｙ」）、次の工程Ｓ９８において、ｉ番目の粗視化粒子９i（例えば、４番目の粗視化粒子９ｄ）の仮決定位置が、熱平衡状態へと緩和されるｉ番目の粗視化粒子９iの位置として決定される。決定された位置２２は、コンピュータ１に記憶される。

一方、工程Ｓ９７において、工程Ｓ９２～工程Ｓ９６が、予め定められた回数繰り返されていないと判断された場合（工程Ｓ９７において、「Ｎ」）、工程Ｓ９２～工程Ｓ９６が再度実施される。

このように、本実施形態の計算工程Ｓ２３では、メトロポリス法の棄却採択法に基づいて、ｉ番目の粗視化粒子９_iが形成する結合角θの熱平衡状態のアンサンブルを得ることができる。しかも、計算工程Ｓ２３では、予め定められた回数に達するまで、ｉ番目の粗視化粒子９iの熱平衡状態へと緩和される位置を棄却・採択する試行（工程Ｓ９２～工程Ｓ９６）が行われる。これにより、本実施形態の計算工程Ｓ２３では、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_i（例えば、４番目の粗視化粒子９ｄ）の位置２２を、精度良く計算することができる。

本実施形態の計算工程Ｓ２３では、ｉ番目の粗視化粒子９iの熱平衡状態へと緩和される位置を棄却・採択する試行（工程Ｓ９２～工程Ｓ９６）が、予め定められた回数繰り返されたが、このような態様に限定されない。棄却・採択する試行（工程Ｓ９２～工程Ｓ９６）が、予め定められた上限値を超えた場合に、１番目の粗視化粒子９ａの配置からやり直されてもよい。これにより、粗視化分子モデルの熱平衡状態への緩和が不十分になるのを防ぐことができる。

次に、本実施形態の配置工程Ｓ２では、図１２（ｃ）に示されるように、計算工程Ｓ２３で計算された位置（決定された位置）２２が、ｉ番目の粗視化粒子９_iの位置に設定される（工程Ｓ２４）。工程Ｓ２４では、仮決定されたｉ番目の粗視化粒子９_iの座標値が、計算された位置２２の座標値に変換される。変換された座標値は、コンピュータ１に記憶される。

次に、本実施形態の配置工程Ｓ２では、粗視化分子モデル６の他端を構成するｍ番目の粗視化粒子９ｍ（図５に示す）まで配置されたか否かが判断される（工程Ｓ２５）。工程Ｓ２５において、ｍ番目の粗視化粒子９ｍまで配置されたと判断された場合（工程Ｓ２５において、「Ｙ」）、配置工程Ｓ２の一連の処理が終了し、次の工程Ｓ３（図７に示す）が実施される。他方、工程Ｓ２５において、ｍ番目の粗視化粒子９ｍまで配置されていないと判断された場合（工程Ｓ２５において、「Ｎ」）、次の粗視化粒子９（即ち、ｉ←ｉ＋１）が選択され（工程Ｓ２６）、工程Ｓ２２～工程Ｓ２５が再度実施される。

これにより、配置工程Ｓ２では、図５に示されるように、粗視化分子モデル６の一端を構成する１番目の粗視化粒子９ａから粗視化分子モデルの他端を構成するｍ番目の粗視化粒子９ｍまでの各位置を、仮想空間８内に順番に決定しながら配置することができる。

次に、本実施形態の作成方法では、全ての粗視化分子モデル６（図３に示す）が仮想空間８に配置されたか否かが判断される（工程Ｓ３）。仮想空間８に配置される粗視化分子モデル６の本数については、高分子材料モデル１２（図３に示す）が用いられるシミュレーション等に応じて、適宜設定される。

工程Ｓ３において、全ての粗視化分子モデル６が仮想空間８に配置されたと判断された場合（工程Ｓ３において、「Ｙ」）、次の工程Ｓ４が実施される。他方、工程Ｓ３において、全ての粗視化分子モデル６が仮想空間８に配置されていないと判断された場合（工程Ｓ３において、「Ｎ」）、配置工程Ｓ２が再度実施されて、新たな粗視化分子モデル６が仮想空間８に配置される。これにより、本実施形態の作成方法では、全ての粗視化分子モデル６を仮想空間８に配置した高分子材料モデル１２を作成することができる。

このように、本実施形態の作成方法では、図１２に示されるように、仮決定されたｉ番目の粗視化粒子９_iと、ｉ－ｋ番目からｉ－１番目までの粗視化粒子９のみとの相互作用を考慮して、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子９_iの位置が計算されるため、高分子鎖の熱平衡状態の慣性半径を決定づける主鎖の曲がりやすさを再現することができる。これにより、本実施形態の作成方法では、熱平衡状態の慣性半径を再現しうる粗視化分子モデル６を作成することができる。したがって、本実施形態の作成方法では、従来の方法のように、分子動力学計算に基づく粗視化分子モデル６の緩和計算を長時間行う必要がないため、熱平衡状態の慣性半径を再現しうる粗視化分子モデル６を、短時間で作成できる。

次に、本実施形態の作成方法では、配置工程Ｓ２の後に、コンピュータ１が、図３に示した粗視化粒子９、９間に粒子の重なりを許容するソフトポテンシャルＵ（ｒ）を定義して、粗視化分子モデル６の構造緩和を計算する（工程Ｓ４）。工程Ｓ４では、仮想空間８に配置された全ての粗視化分子モデル６の粗視化粒子９、９間に、粒子の重なりを許容するソフトポテンシャルＵ（ｒ）が定義される。ソフトポテンシャルＵ（ｒ）は、上記論文１に基づいて、下記式（４）で定義される。

ここで、
ｒ：粗視化粒子間の距離
Ａ：粗視化粒子間に定義されるポテンシャルの強度
σ：粗視化分子モデルの長さの単位
なお、距離ｒは、各粗視化粒子９の中心９ｏ、９ｏ間の距離として定義される。

このようなソフトポテンシャルＵ（ｒ）は、斥力場だけを有しており、粗視化粒子９の重なりを許容することができる。なお、第４ポテンシャルＰ４は無効に設定される。

そして、工程Ｓ４では、仮想空間８に配置された全ての粗視化分子モデル６を対象に、分子動力学計算に基づいた構造緩和が計算される。構造緩和の計算は、仮想空間８に配置された全ての粗視化粒子９の重なりが除去されるまで行われる。なお、各粗視化分子モデル６は、熱平衡状態の慣性半径が再現されているため、粗視化分子モデルの緩和計算を長時間行う必要がない。

緩和計算後には、ソフトポテンシャルＵ（ｒ）が無効に設定されるとともに、第４ポテンシャルＰ４が有効に設定される。これにより、本実施形態の作成方法では、全ての粗視化粒子９の重なりが除去された高分子材料モデル１２を作成することができる。高分子材料モデル１２は、コンピュータ１に記憶される。

本実施形態の作成方法で作成された粗視化分子モデル６及び高分子材料モデル１２は、例えば、一般的に行われている単軸引張り試験に基づいて、一方向に（例えば、Ｘ軸方向に０％～２０％）伸長させる変形シミュレーション等に用いることができる。従って、本実施形態の作成方法は、高分子材料の開発及び製造に役立つ。

本実施形態の配置工程Ｓ２では、計算工程Ｓ２３に先立ち、第１ポテンシャルＵ_angle（θ）、第２ポテンシャルＵ_torsion（φ）、及び、第３ポテンシャルＥ_i（θ）を定義することができる。これらのポテンシャルの定義の有無や定義の内容については、ｉ番目の粗視化粒子９_iの位置の仮設定、及び、高分子材料の熱平衡状態へと緩和を計算できれば特に限定されない。例えば、粗視化分子モデル６がKremer-Grestモデルの場合、第１ポテンシャルＵ_angle（θ）、及び、第２ポテンシャルＵ_torsion（φ）は省略される。さらに、気体状態のＫＧモデルなど、溶媒和の影響をほぼ無視できる場合には、第３ポテンシャルＥ_i（θ）が0に近くなる。このため、第３ポテンシャルＥ_i（θ）を省略することができる。

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

図７に示した処理手順に従って、粗視化分子モデルがコンピュータで作成された（実施例、比較例）。

実施例では、図８～１１に示した処理手順に従って、粗視化分子モデルの一端を構成する１番目の粗視化粒子から粗視化分子モデルの他端を構成するｍ番目の粗視化粒子までの各位置を、仮想空間内に順番に決定しながら配置する配置工程が実施された。実施例の配置工程では、図１３に示した処理手順に従って、ｉ番目（ｉは、正の整数）の粗視化粒子を、ｉ番目の粗視化粒子と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用のみを考慮して仮決定する工程と、考慮した相互作用を用いて、高分子材料の熱平衡状態へと緩和されるｉ番目の粗視化粒子の位置を計算する計算工程と、計算された位置をｉ番目の粗視化粒子の位置として決定する工程とが実施された。メトロポリス法に基づく工程Ｓ９２～Ｓ９６の繰り返し回数は、各粗視化粒子について１０回とした。

実施例では、第１ポテンシャル、第２ポテンシャル、第３ポテンシャル及び第４ポテンシャルが定義された。そして、実施例では、配置する工程の後に、粗視化粒子の重なりを許容するソフトポテンシャルを、粗視化粒子間に定義して、粗視化分子モデルの構造緩和が計算された。

比較例では、従来の方法に基づいて、粗視化分子モデルの一端を構成する１番目の粗視化粒子から粗視化分子モデルの他端を構成するｍ番目の粗視化粒子までの各位置が、モンテカルロ法に基づいて、仮想空間内にランダムに配置された。そして、比較例では、熱平衡状態の慣性半径を再現した粗視化分子モデルを作成するために、分子動力学計算に基づく粗視化分子モデルの緩和計算が実施された。

そして、実施例及び比較例の粗視化分子モデルについて、慣性半径の二乗平均が求められた。共通仕様は、次のとおりである。
粗視化分子モデル：
本数：１万本
１本の粗視化分子モデルに含まれる粗視化粒子の個数（鎖長）ｍ：１００個

テストの結果、実施例の慣性半径の二乗平均は３７．５±０．６σであり、比較例の慣性半径の二乗平均は３６．３±０．４σであった。これらの値は、実際の慣性半径の二乗平均に近似している。一方、実施例の作成時間は、比較例の作成時間の２万分の１に短縮できた。従って、実施例は、熱平衡状態の慣性半径を再現しうる粗視化分子モデルを、短時間で作成することができた。

Ｓ２２ ｉ番目以降の粗視化粒子の位置を仮決定する工程
Ｓ２３ ｉ番目の粗視化粒子の位置を計算する工程
Ｓ２４ 計算された位置をｉ番目の粗視化粒子の位置として決定する工程

Claims (5)

1. 高分子材料の分子鎖を複数の粗視化粒子で表現した粗視化分子モデルを、コンピュータを用いて作成するための方法であって、
   前記粗視化分子モデルの一端を構成する１番目の粗視化粒子から前記粗視化分子モデルの他端を構成するｍ番目（ｍ＝５以上の整数）の粗視化粒子までの各位置を、前記コンピュータの仮想空間内に順番に決定しながら配置する工程を含み、
   前記配置する工程は、ｉ番目（ｉ＝４以上の整数）以降の粗視化粒子の位置を、前記ｉ番目の粗視化粒子と、その直前のｋ個（ｋ＝３以上の整数）の粗視化粒子（ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子）との相互作用のみを考慮して仮決定する工程と、
   前記考慮した相互作用を用いて、前記高分子材料の熱平衡状態へと緩和される前記ｉ番目の粗視化粒子の位置を計算する工程と、
   前記計算された位置を前記ｉ番目の粗視化粒子の位置として決定する工程とを含み、
   前記計算する工程は、
   仮決定された前記ｉ番目の粗視化粒子とは別に、前記相互作用のみを考慮して、前記ｉ番目の粗視化粒子の移動候補位置を設定する工程と、
   前記移動候補位置を、前記熱平衡状態へと緩和される前記ｉ番目の粗視化粒子の仮決定された位置として更新するか否かを判断する工程と、
   前記判断する工程の後、前記仮決定された位置を、前記熱平衡状態へと緩和されるｉ番目の粗視化粒子の位置として決定する工程とを含む、
   粗視化分子モデルの作成方法。
2. 前記相互作用は、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、及び、前記ｉ－２番目の粗視化粒子の間の結合角に関する第１ポテンシャルを含む、請求項１記載の粗視化分子モデルの作成方法。
3. 前記相互作用は、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、前記ｉ－２番目の粗視化粒子、及び、前記ｉ－３番目の粗視化粒子の間の二面角に関する第２ポテンシャルを含む、請求項１又は２記載の粗視化分子モデルの作成方法。
4. 前記計算する工程に先立ち、下記式（１）に基づいて、前記ｉ番目の粗視化粒子、ｉ－１番目の粗視化粒子、及び、前記ｉ－２番目の粗視化粒子の間の結合角に関する第３ポテンシャルを定義する工程をさらに含む、請求項１乃至３のいずれかに記載の粗視化分子モデルの作成方法。
   

   

   ここで、
   Ｅ_i（θ）：第３ポテンシャル
   Ｐ_i（θ）：熱平衡状態の結合角の頻度分布
   Ｐ_i,local（θ）：ｉ番目の粗視化粒子と、ｉ－ｋ番目からｉ－１番目までの全ての粗視化粒子のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布
   θ：結合角
   k_B：ボルツマン定数
   Ｔ：温度
5. 前記配置する工程の後に、前記粗視化粒子の重なりを許容するソフトポテンシャルを、前記粗視化粒子間に定義して、前記粗視化分子モデルの構造緩和を計算する工程をさらに含む、請求項１乃至４のいずれかに記載の粗視化分子モデルの作成方法。

Images