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

Abstract

【課題】 ビーズでモデル化された分子鎖モデルの長さの単位を、高分子鎖の長さの単位に精度よく換算する。【解決手段】 コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法である。このシミュレーション方法は、高分子鎖の全原子モデル又はユナイテッドアトムモデルからなる第１分子鎖モデルの屈曲度合を示す第１パラメータを計算する工程Ｓ３、高分子鎖を複数のビーズで単純化した第２分子鎖モデルの屈曲度合を示す第２パラメータを計算する工程Ｓ６、第１パラメータと第２パラメータとに基づいて、第２分子鎖モデルのビーズ数と、高分子鎖のモノマー数との比を計算する第３計算工程Ｓ７、及び、第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比に基づいて、第２分子鎖モデルの長さの単位を、前記高分子鎖の長さの単位に換算する工程Ｓ８を含む。【選択図】図５

Description

本発明は、高分子材料のシミュレーション方法に関し、詳しくは、ビーズでモデル化された分子鎖モデルの長さの単位を、高分子鎖の長さの単位に精度よく換算することができる高分子材料のシミュレーション方法に関する。

近年、ゴム材料等の高分子材料の反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。この種のシミュレーション方法では、例えば、高分子鎖を複数のビーズでモデル化した分子鎖モデルが、コンピュータに入力される。そして、分子鎖モデルが予め定められた空間に配置され、コンピュータによる分子動力学計算が行われる。

分子鎖モデルを用いた分子動力学計算では、現実の高分子鎖とは異なる単位系が用いられている。このため、下記特許文献１は、分子鎖モデルの単位系を、高分子鎖の単位系に換算する方法を提案している。

下記特許文献１の換算方法は、先ず、分子鎖モデル及び全原子モデルの理想鎖のパラメータ（Rouseのパラメータ）をそれぞれ求めている。そして、下記特許文献１の換算方法は、分子鎖モデルの理想鎖のパラメータと、全原子モデルの理想鎖のパラメータとが同一となるとみなして、分子鎖モデルの単位系を、高分子鎖の単位系に換算している。

特開２０１４−２０６９１３号公報

上記特許文献１では、長さの単位の換算精度が十分ではなく、さらなる改善の余地があった。

本発明は、以上のような実状に鑑み案出されたもので、ビーズでモデル化された分子鎖モデルの長さの単位を、高分子鎖の長さの単位に精度よく換算できる高分子材料のシミュレーション方法を提供することを主たる目的としている。

本発明は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、前記高分子鎖の全原子モデル又はユナイテッドアトムモデルからなる第１分子鎖モデルを、前記コンピュータに入力する第１モデル設定工程、前記コンピュータが、分子動力学計算に基づいて、前記第１分子鎖モデルの緩和した構造を計算する工程、前記コンピュータが、構造緩和後の前記第１分子鎖モデルの屈曲度合を示す第１パラメータを計算する第１計算工程、前記高分子鎖を複数のビーズで単純化した第２分子鎖モデルを、前記コンピュータに入力する第２モデル設定工程、前記コンピュータが、分子動力学計算に基づいて、前記第２分子鎖モデルの緩和した構造を計算する工程、前記コンピュータが、構造緩和後の前記第２分子鎖モデルの屈曲度合を示す第２パラメータを計算する第２計算工程、前記コンピュータが、前記第１パラメータと前記第２パラメータとに基づいて、前記第２分子鎖モデルのビーズ数と、前記高分子鎖のモノマー数との比を計算する第３計算工程、並びに、前記第２分子鎖モデルの前記ビーズ数と前記高分子鎖の前記モノマー数との比に基づいて、前記第２分子鎖モデルの長さの単位を、前記高分子鎖の長さの単位に換算する工程を含むことを特徴とする。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１モデル設定工程は、前記モノマー数が異なる前記第１分子鎖モデルを設定する工程を含み、前記第１計算工程は、前記モノマー数が異なる前記第１分子鎖モデル毎に、前記第１パラメータを計算する工程を含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２モデル設定工程は、前記ビーズ数が異なる前記第２分子鎖モデルを設定する工程を含み、前記第２計算工程は、前記ビーズ数が異なる前記第２分子鎖モデル毎に、前記第２パラメータを計算する工程を含むのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第１パラメータ及び前記第２パラメータは、下記式（１）で計算されるのが望ましい。

ここで、
Ｎ_Kuhn：Kuhnのセグメント数
Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
Ｒ_g：第１分子鎖モデル又は第２分子鎖モデルの慣性半径

本発明に係る前記高分子材料のシミュレーション方法において、前記第１パラメータ及び前記第２パラメータは、下記式（２）で計算されるのが望ましい。

ここで、
Ｎ_Kuhn：Kuhnのセグメント数
Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
Ｒ：第１分子鎖モデル又は第２分子鎖モデルの末端間ベクトルの長さ

本発明に係る前記高分子材料のシミュレーション方法において、前記第１分子鎖モデルは、前記高分子鎖のモノマーをモデル化したモノマーモデルを、予め定められた前記モノマー数で連結して設定され、前記第１分子鎖モデルの全長Ｌは、前記各モノマーモデルの長さを平均した平均モノマー長と、前記モノマー数とを乗じることにより求められるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２分子鎖モデルの全長Ｌは、隣接する前記ビーズの平衡核間距離と、前記ビーズ数とを乗じることによって求められるのが望ましい。

本発明に係る前記高分子材料のシミュレーション方法において、前記第２計算工程は、前記第２パラメータと前記ビーズ数との第２関係を求める工程を含み、前記第３計算工程は、前記第２関係を、下記式（３）で定義される曲線にフィッティングする工程を含むのが望ましい。

ここで、
Ｎ_kuhn：Kuhnのセグメント数
Ｎ_CG：第２分子鎖モデルのビーズ数
ａ、ｂ、ｃ：フィッティングパラメータ

本発明に係る前記高分子材料のシミュレーション方法において、前記第１計算工程は、前記第１パラメータと前記モノマー数との第１関係を求める工程を含み、前記第３計算工程は、下記式（４）に基づいて、前記第１関係を、前記曲線にフィッティングする工程を含むのが望ましい。

ここで、
Ｎ_FA：第１分子鎖モデルのモノマー数
Ｎ_CG：第２分子鎖モデルのビーズ数
ｄ：フィッティングパラメータ

本発明の高分子材料のシミュレーション方法は、高分子鎖の全原子モデル又はユナイテッドアトムモデルからなる第１分子鎖モデルについて、屈曲度合を示す第１パラメータを求める。また、本発明の高分子材料のシミュレーション方法は、高分子鎖を複数のビーズで単純化した第２分子鎖モデルについて、屈曲度合を示す第２パラメータを求める。そして、本発明の高分子材料のシミュレーション方法は、第１パラメータと第２パラメータとに基づいて、第２分子鎖モデルのビーズ数と、高分子鎖のモノマー数との比を計算し、第２分子鎖モデルの長さの単位を、高分子鎖の長さの単位に換算している。

第１分子鎖モデル及び第２分子鎖モデルは、それらの屈曲度合に応じて、第１分子鎖モデル及び第２分子鎖モデルの慣性半径や、末端間ベクトルの長さ等の空間的な広がりが変化する。このため、第１パラメータ及び第２パラメータを用いることにより、屈曲度合の鎖長依存性を考慮して、同一鎖長で空間的な広がりが同じとなるような、第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比を精度よく求めることができる。

第２分子鎖モデルの長さの単位を、高分子鎖の長さの単位に換算するには、第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比（１ビーズあたりのモノマー数）に、モノマー１つあたりの長さを乗じ、さらに、ビーズ１つあたりの長さで除することで行うことができる。これにより、第２分子鎖モデルには、第１分子鎖モデルの特定の原子集団を特定のビーズに１対１で対応させた、所謂ユナイテッドアトムモデルを用いなくても、第１分子鎖モデルとの長さの対応が自明でないKremer Grestモデル等の抽象的なバネビーズモデルを用いることができる。

また、本発明の高分子材料のシミュレーション方法では、屈曲度合の鎖長依存性を考慮して、ビーズ数とモノマー数との比を計算しているため、鎖長の短い（即ち、モノマー数が小さい）第１分子鎖モデルを用いることができる。これにより、第１分子鎖モデルの緩和した構造を短時間で得ることができるため、換算に要する時間を短縮することができる。

本発明のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。 ポリブタジエンの構造式である。 第１分子鎖モデルの一例を示す概念図である。 第２分子鎖モデルの一例を示す概念図である。 本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。 第１緩和計算工程の処理手順の一例を示すフローチャートである。 第１分子鎖モデルが配置された高分子材料モデルの一例を示す概念図である。 第１分子鎖モデルのポテンシャルを説明する概念図である。 第１計算工程の処理手順の一例を示すフローチャートである。 第１分子鎖モデルの第１パラメータ（Kuhnのセグメント数）とモノマー数との第１関係を示すグラフである。 第２緩和計算工程の処理手順の一例を示すフローチャートである。 第２分子鎖モデルが配置された高分子材料モデルの一例を示す概念図である。 第２分子鎖モデルのポテンシャルを説明する概念図である。 第２計算工程の処理手順の一例を示すフローチャートである。 第２分子鎖モデルの第２パラメータ（Kuhnのセグメント数）とビーズ数との第２関係を示すグラフである。 第１分子鎖モデルの第１パラメータと、第２分子鎖モデルの第２パラメータとをフィッティングしたグラフである。

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

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

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

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

第１分子鎖モデル３は、全原子モデルとして構成される。全原子モデルは、図２に示した高分子鎖２の実際の構造に基づいて、原子をモデル化した粒子モデル７で表現したものである。なお、第１分子鎖モデル３は、ユナイテッドアトムモデル（図示省略）として構成されてもよい。ユナイテッドアトムモデルは、図２に示した高分子鎖２において、炭素原子と、炭素原子に結合した水素原子とを、一つの粒子モデル（図示省略）として扱うものである。このような第１分子鎖モデル３は、第２分子鎖モデル４に比べて、微細な構造や運動を取り扱うことができ、シミュレーション精度を向上しうる。

図４に示されるように、第２分子鎖モデル４は、高分子鎖２を複数のビーズ１６で単純化したものである。各ビーズ１６は、直径を持った球として表現される。第２分子鎖モデルの一例としては、Kremer-Grestモデルや、ＤＰＤ（散逸粒子動力学法）に基づくモデル等が採用されうる。本実施形態の第２分子鎖モデルとしては、Kremer-Grestモデルとして構成される。このような第２分子鎖モデル４は、第１分子鎖モデル３に比べて、大きな時空間を扱うことができる利点がある。

第１分子鎖モデル３（図３に示す）を用いた分子動力学計算では、現実の高分子鎖２（図２に示す）に基づいた長さの単位が用いられている。他方、第２分子鎖モデル４（図４に示す）を用いた分子動力学計算では、現実の高分子鎖２とは異なる長さの単位が用いられている。このため、第２分子鎖モデル４の分子動力学計算の計算結果を、実際の高分子鎖２の運動として取り扱うためには、第２分子鎖モデル４の長さの単位を、高分子鎖２の長さの単位に精度よく換算することが重要である。

本実施形態のシミュレーション方法では、第２分子鎖モデル４（図４に示す）を用いたシミュレーションに先立ち、第２分子鎖モデル４の長さの単位を、高分子鎖２（図２に示す）の長さの単位に換算している。図５は、本実施形態のシミュレーション方法の処理手順の一例を示すフローチャートである。

本実施形態のシミュレーション方法では、先ず、第１分子鎖モデル３を、コンピュータ１に入力する（第１モデル設定工程Ｓ１）。図３に示されるように、第１分子鎖モデル３は、複数の粒子モデル７と、粒子モデル７、７間を結合するボンド８とを含んで構成されている。

粒子モデル７及びボンド８は、図２に示した高分子鎖２のモノマー５を表す単位構造に基づいて連結されることにより、モノマーモデル１１が設定される。このモノマーモデル１１が、予め定められたモノマー数（即ち、分子量（重合度））に基づいて連結されることにより、第１分子鎖モデル３が設定される。

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

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

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

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

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

本実施形態の第１モデル設定工程Ｓ１は、モノマー数が異なる第１分子鎖モデル３が設定される。即ち、第１モデル設定工程Ｓ１では、複数のモノマー数に基づいて、モノマーモデル１１がそれぞれ連結される。これにより、モノマー数が異なる複数の第１分子鎖モデル３が設定される。また、モノマー数が異なる第１分子鎖モデル３の種類については、適宜設定されうる。モノマー数が異なる第１分子鎖モデル３は、例えば、２〜８種類程度設定されるのが望ましい。本実施形態では、モノマー数が２０、３０及び４０に設定された３種類の第１分子鎖モデル３が設定されている。これらの第１分子鎖モデル３は、コンピュータ１で取り扱い可能な数値データであり、コンピュータ１に入力される。

次に、本実施形態のシミュレーション方法では、コンピュータが、分子動力学計算に基づいて、第１分子鎖モデルの緩和した構造を計算する（第１緩和計算工程Ｓ２）。図６は、第１緩和計算工程Ｓ２の処理手順の一例を示すフローチャートである。

本実施形態の第１緩和計算工程Ｓ２では、先ず、コンピュータ１が、モノマー数が異なる各第１分子鎖モデル３の初期配置を決定する（工程Ｓ２１）。工程Ｓ２１では、予め定められた空間に、第１分子鎖モデル３が配置された高分子材料モデルが設定される。図７は、第１分子鎖モデル３が配置された高分子材料モデル１０の一例を示す概念図である。

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

空間９は、解析対象の高分子材料の微小構造部分に相当する。本実施形態の空間９は、互いに向き合う三対の平面９Ｓ、９Ｓを有する立方体として定義されている。各平面９Ｓには、周期境界条件が定義されている。これにより、空間９では、例えば、一方の平面９Ｓａから出て行った第１分子鎖モデル３の一部が、反対側の平面９Ｓｂから入ってくるように計算することができる。従って、一方の平面９Ｓａと、反対側の平面９Ｓｂとが連続している（繋がっている）ものとして取り扱うことができる。

各空間９に配置される第１分子鎖モデル３の本数については、モノマー数や、空間９の大きさに基づいて適宜設定される。なお、第１分子鎖モデル３の本数が少ないと、後述の分子動力学計算において、一方の平面９Ｓａから出て行きかつ反対側の平面９Ｓｂから入ってきた第１分子鎖モデル３の一端と、平面９Ｓｂ側に配置されている第１分子鎖モデル３の他端とが絡まりやすくなり、計算落ちするおそれがある。さらに、第１分子鎖モデル３の慣性半径を精度良く計算するためには、平衡時において第１分子鎖モデル３同士の絡まりを可能な限り防ぐ必要がある。このため、第１分子鎖モデル３の本数としては、平衡時の周期境界長が平衡時の慣性半径の３倍以上に長くなりうる本数、本実施形態では、２０本以上に設定されるのが望ましい。

逆に、第１分子鎖モデル３の本数が多くても、運動方程式の質点として取り扱われる粒子モデル７の粒子数が増大し、多くの計算時間を要するおそれがある。このような観点より、空間９に配置される第１分子鎖モデル３の本数は、２００本以下に設定されるのが望ましい。

空間９の一辺の長さＬａは、適宜設定することができる。本実施形態の長さＬａは、系内の粒子モデル７の初期密度が、例えば０．００１g／cm³となるように設定されるのが望ましい。また、空間９の大きさは、例えば１atmで安定な体積に設定されるのが望ましい。このような空間９では、後述の分子動力学計算において、第１分子鎖モデル３の回転運動をスムーズに計算することができる。

図８は、第１分子鎖モデル３のポテンシャルを説明する概念図である。工程Ｓ２１では、第１分子鎖モデル３の初期配置が決定された各高分子材料モデル１０（即ち、モノマー数が異なる第１分子鎖モデル３がそれぞれ配置された高分子材料モデル１０）において、隣接する第１分子鎖モデル３、３の粒子モデル７、７間に、相互作用ポテンシャルＰ１が定義される。本実施形態の相互作用ポテンシャルＰ１は、ＬＪポテンシャルＵ_LJ（ｒ_ij）であり、下記式（５）で定義される。

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

相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒ_ijがσよりも小さくなるほど、粒子モデル７、７間に作用する斥力が大きくなる。また、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒ_ijがσになるときに最小となり、粒子モデル７、７間に斥力や引力は働かない。さらに、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒ_ijがσよりも大になるほど、粒子モデル７、７間に作用する引力が働く。このように、相互作用ポテンシャルＰ１は、粒子モデル間の距離ｒ_ijに応じて、斥力及び引力を定義することができる。

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

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

次に、本実施形態の第１緩和計算工程Ｓ２では、コンピュータ１が、第１分子鎖モデル３の初期配置を十分に緩和できたか否かを判断する（工程Ｓ２３）。なお、初期配置の緩和の判断基準については、第１分子鎖モデル３の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。

工程Ｓ２３において、第１分子鎖モデル３の初期配置を十分に緩和できたと判断された場合（工程Ｓ２３で、「Ｙ」）、次の工程Ｓ２４が実施される。他方、第１分子鎖モデル３の初期配置を十分に緩和できていないと判断された場合（工程Ｓ２３で、「Ｎ」）は、単位時間（例えば、工程Ｓ２２で実施済みのシミュレーション時間程度まとめて）を進めて（工程Ｓ２５）、工程Ｓ２２（分子動力学計算）及び工程Ｓ２３が再度実施される。これにより、第１分子鎖モデル３の初期配置を、確実に緩和することができる。

次に、本実施形態の第１緩和計算工程Ｓ２では、コンピュータ１が、全ての第１分子鎖モデル（本実施形態では、モノマー数が異なる３種類の第１分子鎖モデル）３の緩和状態が計算されたか否かを判断する（工程Ｓ２４）。工程Ｓ２４において、全ての第１分子鎖モデル３の緩和計算が計算されたと判断された場合（工程Ｓ２４で、「Ｙ」）、第１緩和計算工程Ｓ２の一連の処理が終了し、次の第１計算工程Ｓ３（図５に示す）が実施される。他方、全ての第１分子鎖モデル３の緩和計算が終了していないと判断された場合（工程Ｓ２４で、「Ｎ」）、緩和状態が計算されていない他の第１分子鎖モデル３（即ち、図７に示した高分子材料モデル１０）が選択され（工程Ｓ２６）、工程Ｓ２２〜工程Ｓ２４が実施される。これにより、第１緩和計算工程Ｓ２では、モノマー数が異なる全ての第１分子鎖モデル３の平衡状態（緩和状態）が計算される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、構造緩和後の第１分子鎖モデル３の屈曲度合を示す第１パラメータを計算する（第１計算工程Ｓ３）。本実施形態の第１パラメータは、Kuhnのセグメント数Ｎ_Kuhnであり、下記式（１）又は下記式（２）で計算される。

ここで、
Ｎ_Kuhn：Kuhnのセグメント数
Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
Ｒ_g：第１分子鎖モデル又は第２分子鎖モデルの慣性半径

ここで、
Ｎ_Kuhn：Kuhnのセグメント数
Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
Ｒ：第１分子鎖モデル又は第２分子鎖モデルの末端間ベクトルの長さ

Kuhnのセグメント数Ｎ_Kuhnは、高分子物理で知られている理想鎖のパラメータの一つである。上記式（１）又は上記式（２）において、Kuhnのセグメント数の値が大きいほど、第１分子鎖モデル３（図３に示す）及び第２分子鎖モデル４（図４に示す）は、コンパクトに折り畳んだ形状となる。他方、Kuhnのセグメント数の値が小さいほど、第１分子鎖モデル３及び第２分子鎖モデル４の折れ曲がりが小さくなり、直線状に連続してのびる形状となる。従って、Kuhnのセグメント数は、第１分子鎖モデル３の屈曲度合を示す第１パラメータ、又は、後述の第２分子鎖モデル４の屈曲度合を示す第２パラメータとして用いられる。第１分子鎖モデル３及び第２分子鎖モデル４は、それらの屈曲度合に応じて、第１分子鎖モデル又は第２分子鎖モデルの慣性半径Ｒ_gや全長Ｌが変化する。

上記式（１）では、例えば、第１分子鎖モデル３の末端間ベクトルＶ１（図３に示す）の長さのアンサンブル平均を用いる上記式（２）や上記特許文献１とは異なり、慣性半径Ｒ_gを用いて、Kuhnのセグメント数が求められる。これにより、上記式（１）は、第１分子鎖モデル３の末端側の粒子モデル７の高い運動性の影響を小さくできるため、Kuhnのセグメント数を精度よく求めることができる。本実施形態では、第１パラメータ及び第２パラメータが、上記式（１）を用いて求められる態様について説明する。

本実施形態の第１計算工程Ｓ３は、モノマー数が異なる第１分子鎖モデル３毎に、第１パラメータを計算している。図９は、第１計算工程Ｓ３の処理手順の一例を示すフローチャートである。

本実施形態の第１計算工程Ｓ３では、先ず、コンピュータ１が、第１分子鎖モデル３の全長Ｌ（図示省略）を計算する（工程Ｓ３１）。本実施形態の工程Ｓ３１では、先ず、各モノマーモデル１１の長さＬａ（図３に示す）を平均した平均モノマー長が求められる。長さＬａは、隣接するモノマーモデル１１の共有結合の中点１２、１２間の距離である。また、平均モノマー長は、第１分子鎖モデル３の末端に配置されたモノマーモデル１１、１１を除いたモノマーモデル１１の長さＬａの平均値である。そして、工程３１では、平均モノマー長と、モノマー数とを乗じることにより、第１分子鎖モデル３の全長Ｌが求められる。このような計算方法により、例えば、第１分子鎖モデル３を主鎖方向に強制的に引き伸ばす変形計算が実施されていた上記特許文献１の計算方法に比べて、作業効率を高めつつ、計算コストを大幅に低減しうる。また、複数の第１分子鎖モデル３の全長Ｌを容易に求めることができるため、複数の第１分子鎖モデル３の全長Ｌをさらに平均することで、統計誤差を減らすこともできる。第１分子鎖モデル３の全長Ｌは、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ３では、コンピュータ１が、第１分子鎖モデル３の慣性半径Ｒ_g（図示省略）を計算する（工程Ｓ３２）。第１分子鎖モデル３の慣性半径Ｒ_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 ）に記載されている式２．５に基づいて、工程Ｓ２２の分子動力学計算で求められた第１分子鎖モデル３の粒子モデル７の座標値から計算される。第１分子鎖モデル３の慣性半径Ｒ_gは、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ３では、コンピュータ１が、第１パラメータを計算する（工程Ｓ３３）。上述したように、第１パラメータ（即ち、第１分子鎖モデル３のKuhnのセグメント数）は、上記式（１）を用いて、第１分子鎖モデル３の全長Ｌ（図示省略）、及び、第１分子鎖モデル３の慣性半径Ｒ_g（図示省略）から計算される。第１パラメータは、コンピュータ１に入力される。

次に、本実施形態の第１計算工程Ｓ３では、コンピュータ１が、モノマー数が異なる全ての第１分子鎖モデル（本実施形態では、モノマー数が異なる３種類の第１分子鎖モデル）３について、第１パラメータ（即ち、Kuhnのセグメント数）が計算されたか否かを判断する（工程Ｓ３４）。工程Ｓ３４において、全ての第１分子鎖モデル３の第１パラメータが計算されたと判断された場合（工程Ｓ３４において、「Ｙ」）、第１計算工程Ｓ３の一連の処理が終了し、次の第２モデル設定工程Ｓ４が実施される。他方、全ての第１分子鎖モデル３の第１パラメータが計算されていないと判断された場合（工程Ｓ３４において、「Ｎ」）、第１パラメータが計算されていない他の第１分子鎖モデル３が選択され（工程Ｓ３５）、工程Ｓ３１〜Ｓ３４が実施される。これにより、第１計算工程Ｓ３では、モノマー数が異なる第１分子鎖モデル３毎に、第１パラメータが計算される。そして、第１計算工程Ｓ３では、第１パラメータとモノマー数との第１関係が求められる。

図１０は、第１分子鎖モデル３の第１パラメータ（Kuhnのセグメント数）とモノマー数との第１関係を示すグラフである。本実施形態の第１関係は、第１パラメータ（Kuhnのセグメント数）とモノマー数との比Ｎ_kuhn／Ｎ_FAと、モノマー数Ｎ_FAとの関係が示されている。図１０の例では、モノマー数が２０、３０及び４０の第１分子鎖モデル３の第１パラメータが示されている。

次に、本実施形態のシミュレーション方法は、第２分子鎖モデル４を、コンピュータ１に入力する（第２モデル設定工程Ｓ４）。図４に示されるように、第２分子鎖モデル４は、複数のビーズ１６と、ビーズ１６、１６間を結合する結合鎖１７とを含んで構成されている。本実施形態では、予め定められたビーズ数に基づいて、ビーズ１６が連結されることにより、第２分子鎖モデル４が設定される。

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

結合鎖１７は、ビーズ１６、１６間に平衡長を定義した結合ポテンシャルとして構成される。ここで、「平衡長」とは、ビーズ１６、１６間の結合距離である。この結合距離が変化した場合は、結合鎖１７によって、元の平衡長が保持される。これにより、第２分子鎖モデル４は、直鎖状の三次元構造を維持することができる。なお、平衡長は、隣り合うビーズ１６、１６の中心１６ａ、１６ａ間の距離として定義される。また、結合鎖１７の結合ポテンシャルには、例えば、上記論文２に基づいて設定されるのが望ましい。

ビーズ数については、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ１の性能等に基づいて、適宜設定することができる。ビーズ数は、例えば、５〜１０００から選択されるのが望ましい。なお、このようなビーズ数に限定されるわけではない。

本実施形態の第２モデル設定工程Ｓ４は、ビーズ数が異なる第２分子鎖モデル４が設定される。即ち、第２モデル設定工程Ｓ４は、複数のビーズ数に基づいて、ビーズ１６が連結されることにより、ビーズ数が異なる複数の第２分子鎖モデル４が設定される。また、ビーズ数が異なる第２分子鎖モデル４の種類については、適宜設定されうる。ビーズ数が異なる複数の第２分子鎖モデル４は、例えば、５〜１５種類程度が設定されるのが望ましい。本実施形態では、ビーズ数が１０、１３、２０、２７、５０、８０、１００、２００及び５００に設定された９種類の第２分子鎖モデル４が設定されている。

このように、本実施形態のシミュレーション方法では、第１分子鎖モデル３のモノマー数（本実施形態では、２０、３０及び４０）と、第２分子鎖モデルのビーズ数（本実施形態では、１０、１３、２０、２７、５０、８０、１００、２００及び５００）とが、１対１で対応していない（即ち、モノマー数とビーズ数とが同一でない）。さらに、本実施形態の第２モデル設定工程Ｓ４は、第１分子鎖モデル３のモノマー数よりも大きいビーズ数に設定された鎖長の大きい第２分子鎖モデル４を含んでいる。これらの第２分子鎖モデル４は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ１に入力される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、分子動力学計算に基づいて、第２分子鎖モデル４の緩和した構造を計算する（第２緩和計算工程Ｓ５）。図１１は、第２緩和計算工程Ｓ５の処理手順の一例を示すフローチャートである。

第２緩和計算工程Ｓ５は、先ず、コンピュータ１が、ビーズ数が異なる各第２分子鎖モデル４の初期配置を決定する（工程Ｓ５１）。工程Ｓ５１では、予め定められた空間（セル）に、第２分子鎖モデル４が配置された高分子材料モデルが設定される。図１２は、第２分子鎖モデル４が配置された高分子材料モデル１８の一例を示す概念図である。

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

空間９は、解析対象の高分子材料の微小構造部分に相当し、図８に示した空間９と同様に設定される。各空間９に配置される第２分子鎖モデル４の本数については、ビーズ数や、空間９の大きさに基づいて適宜設定される。なお、第２分子鎖モデル４の本数が少ないと、第１分子鎖モデル３（図７に示す）と同様に、第２分子鎖モデル４の一端と他端とが絡まりやすくなり、計算落ちするおそれがある。さらに、第２分子鎖モデル４の慣性半径を精度良く計算するためには、平衡時において第２分子鎖モデル４同士の絡まりを可能な限り防ぐ必要がある。このため、第２分子鎖モデル４の本数としては、平衡時の周期境界長が平衡時の慣性半径の３倍以上に長くなりうる本数が選択されるのが望ましい。逆に、第２分子鎖モデル４の本数が多くても、運動方程式の質点として取り扱われるビーズ１６の粒子数が増大し、多くの計算時間を要するおそれがある。このような観点より、空間９に配置される第２分子鎖モデル４の本数は、好ましくは５０本以上であり、また、好ましくは２０００本以下である。

図１３は、第２分子鎖モデル４のポテンシャルを説明する概念図である。工程Ｓ５１では、各高分子材料モデル１８において、隣接する第２分子鎖モデル４、４のビーズ１６、１６間に、相互作用ポテンシャルＰ２がそれぞれ定義される。本実施形態の相互作用ポテンシャルＰ２は、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（２）で定義される。なお、上記式（３）において、「粒子モデル間」は、「ビーズ間」に置き換えて適用される。

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

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

次に、本実施形態の第２緩和計算工程Ｓ５では、コンピュータ１が、第２分子鎖モデル４の初期配置を十分に緩和できたか否かを判断する（工程Ｓ５３）。なお、初期配置の緩和の判断基準については、第２分子鎖モデル４の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。
＄
工程Ｓ５３において、第２分子鎖モデル４の初期配置を十分に緩和できたと判断された場合（工程Ｓ５３で、「Ｙ」）、次の工程Ｓ５４が実施される。他方、第２分子鎖モデル４の初期配置を十分に緩和できていないと判断された場合（工程Ｓ５３で、「Ｎ」）、単位時間を（例えば、工程Ｓ５２で実施済みのシミュレーション時間程度まとめて）進めて（工程Ｓ５５）、工程Ｓ５２（分子動力学計算）及び工程Ｓ５３が再度実施される。これにより、第２分子鎖モデル４の初期配置を、確実に緩和することができる。

次に、本実施形態の第２緩和計算工程Ｓ５では、コンピュータ１が、全ての第２分子鎖モデル（本実施形態では、ビーズ数が異なる９種類の第２分子鎖モデル）４の緩和状態が計算されたか否かを判断する（工程Ｓ５４）。工程Ｓ５４において、全ての第２分子鎖モデル４の緩和状態が計算されたと判断された場合（工程Ｓ５４で、「Ｙ」）、第２緩和計算工程Ｓ５の一連の処理が終了し、次の第２計算工程Ｓ６が実施される。他方、全ての第２分子鎖モデル４の緩和状態が計算されていないと判断された場合（工程Ｓ５４で、「Ｎ」）、緩和状態が計算されていない他の第２分子鎖モデル４（即ち、図１２に示した高分子材料モデル１８）が選択され（工程Ｓ５６）、工程Ｓ５２〜工程Ｓ５４が実施される。これにより、ビーズ数が異なる全ての第２分子鎖モデル４の平衡状態（緩和状態）が計算される。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、構造緩和後の第２分子鎖モデル４の屈曲度合を示す第２パラメータを計算する（第２計算工程Ｓ６）。本実施形態の第２パラメータは、上記式（１）又は上記式（２）で計算されるKuhnのセグメント数Ｎ_Kuhnである。本実施形態の第２計算工程Ｓ６は、上記式（１）を用いて、ビーズ数が異なる第２分子鎖モデル４毎に、第２パラメータを計算している。図１４は、第２計算工程Ｓ６の処理手順の一例を示すフローチャートである。

本実施形態の第２計算工程Ｓ６では、先ず、コンピュータ１が、第２分子鎖モデル４の全長Ｌ（図示省略）を計算する（工程Ｓ６１）。第２分子鎖モデル４の全長Ｌは、第２分子鎖モデル４内の隣接するビーズ１６、１６の平衡核間距離と、ビーズ数とを乗じることによって求めることができる。なお、平衡核間距離は、振動の中心となるビーズ１６、１６間の距離であり、上記の分子動力学計算の結果に基づいて求められる。第２分子鎖モデル４の全長Ｌは、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ６では、コンピュータ１が、第２分子鎖モデル４の慣性半径Ｒ_g（図示省略）を計算する（工程Ｓ６２）。第２分子鎖モデル４の慣性半径Ｒ_gは、例えば、上記論文２に記載されている式２．５に基づいて、工程Ｓ５２の分子動力学計算で求められた第２分子鎖モデル４のビーズ１６の座標値から計算される。第２分子鎖モデル４の慣性半径Ｒ_gは、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ６では、コンピュータ１が、第２分子鎖モデル４の第２パラメータを計算する（工程Ｓ６３）。上述したように、第２パラメータ（即ち、第２分子鎖モデル４のKuhnのセグメント数）は、上記式（１）を用いて、第２分子鎖モデル４の全長Ｌ（図示省略）、及び、第２分子鎖モデル４の慣性半径Ｒ_g（図示省略）から計算される。第２パラメータは、コンピュータ１に入力される。

次に、本実施形態の第２計算工程Ｓ６では、コンピュータ１が、ビーズ数が異なる全ての第２分子鎖モデル（本実施形態では、ビーズ数が異なる９種類の第２分子鎖モデル）４について、第２パラメータが計算されたか否かを判断する（工程Ｓ６４）。工程Ｓ６４において、全ての第２分子鎖モデル４の第２パラメータが計算されたと判断された場合（工程Ｓ６４において、「Ｙ」）、第２計算工程Ｓ６の一連の処理が終了し、次の第３計算工程Ｓ７が実施される。他方、全ての第２分子鎖モデル４の第２パラメータが計算されていないと判断された場合（工程Ｓ６４において、「Ｎ」）、第２パラメータが計算されていない他の第２分子鎖モデル４が選択され（工程Ｓ６５）、工程Ｓ６１〜Ｓ６４が実施される。これにより、第２計算工程Ｓ６では、ビーズ数が異なる第２分子鎖モデル４毎に、第２パラメータが計算される。そして、第２計算工程Ｓ６では、第２パラメータとビーズ数との第２関係が求められる。

図１５は、第２分子鎖モデル４の第２パラメータ（Kuhnのセグメント数）とビーズ数との第２関係を示すグラフである。本実施形態の第２関係は、第２パラメータ（Kuhnのセグメント数）とビーズ数との比Ｎ_kuhn／Ｎ_CGと、ビーズ数Ｎ_CGとの関係が示されている。図１５の例では、ビーズ数（鎖長）が１０、１３、２０、２７、５０、８０、１００、２００及び５００の第２分子鎖モデル４のKuhnのセグメント数が示されている。

次に、本実施形態のシミュレーション方法は、コンピュータ１が、第２分子鎖モデル４（図２に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比を計算する（第３計算工程Ｓ７）。第３計算工程Ｓ７では、第１分子鎖モデル３（図３に示す）の第１パラメータと、第２分子鎖モデル４（図４に示す）の第２パラメータとに基づいて、第２分子鎖モデル４の１個のビーズ１６あたりの高分子鎖２のモノマー数を求めている。

第３計算工程Ｓ７では、先ず、統計誤差を考慮した非線形最小二乗法に基づいて、図１０に示した第１分子鎖モデル３の第１パラメータとモノマー数との第１関係と、図１５に示した第２分子鎖モデル４の第２パラメータとビーズ数との第２関係とがフィッティングされる。図１６は、第１分子鎖モデル３の第１パラメータと、第２分子鎖モデル４の第２パラメータとをフィッティングしたグラフである。

第３計算工程Ｓ７では、第２分子鎖モデル４の第２パラメータとビーズ数との第２関係が、図１６に示した滑らかな曲線２０（破線で示している）でフィッティングされる。曲線は、下記式（３）で定義される。このような曲線２０は、３つのフィッティングパラメータａ、ｂ、ｃが用いられるため、末端の運動性が高い影響を考慮して、第２関係に精度よくフィッティングすることができる。

ここで、
Ｎ_kuhn：Kuhnのセグメント数
Ｎ_CG：第２分子鎖モデルのビーズ数
ａ、ｂ、ｃ：フィッティングパラメータ

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

ここで、
Ｎ_FA：第１分子鎖モデルのモノマー数
Ｎ_CG：第２分子鎖モデルのビーズ数
ｄ：フィッティングパラメータ

フィッティングの計算方法としては、例えば、論文３（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）により、容易に行うことができる。

フィッティングに用いる統計誤差としては、第１分子鎖モデル３の慣性半径の統計誤差、及び、第２分子鎖モデル４の慣性半径の統計誤差が計算される。第１分子鎖モデル３の慣性半径の統計誤差は、モノマー数が異なる第１分子鎖モデル３毎に、分子鎖毎の慣性半径の時間平均についての分散を、空間９に配置された第１分子鎖モデル３の本数で除した値を求め、その値の平方根から求めることができる。また、第２分子鎖モデル４の慣性半径の統計誤差は、ビーズ数が異なる第２分子鎖モデル４毎に、分子鎖毎の慣性半径の時間平均についての分散を、空間９に配置された第２分子鎖モデル４の本数で除した値を求め、その値の平方根から求めることができる。

第３計算工程Ｓ７では、図１６に示したグラフに基づいて、第２分子鎖モデル４（図４に示す）のビーズ数と、高分子鎖２（図２に示す）のモノマー数との比が計算される。なお、図１６に示したグラフの例において、第２分子鎖モデル４のビーズ数と、高分子鎖２のモノマー数との比（即ち、１個のビーズ１６あたりのモノマー数）は、１．５である。このような第２分子鎖モデル４のビーズ数と、高分子鎖２のモノマー数との比は、コンピュータ１に入力される。

このように、本実施形態のシミュレーション方法では、第１分子鎖モデル３（図３に示す）の第１パラメータ（Kuhnのセグメント数）と、第２分子鎖モデル４の第２パラメータ（Kuhnのセグメント数）とをフィッティングさせることにより、屈曲度合の鎖長依存性を考慮して、同一鎖長で空間的な広がりが同じとなるような、第２分子鎖モデル４のビーズ数と、高分子鎖２のモノマー数との比を精度よく計算することができる。

また、本実施形態では、第１分子鎖モデル３の第１パラメータ、及び、第２分子鎖モデル４の第２パラメータがそれぞれ複数求められているため、精度良くフィッティングすることができる。なお、第１分子鎖モデル３の第１パラメータ、及び、第２分子鎖モデル４の第２パラメータを確実にフィッティングさせるために、第１分子鎖モデル３のモノマー数が設定される範囲と、第２分子鎖モデル４のビーズ数が設定される範囲とが重複しているのが望ましい。

本実施形態では、上記フィッティングにより、第１分子鎖モデル３のモノマー数（本実施形態では、２０、３０及び４０）と、第２分子鎖モデル４のビーズ数（本実施形態では、１０、１３、２０、２７、５０、８０、１００、２００及び５００）とが１対１で対応していなくても、鎖長依存性を考慮したフィッティング（内挿及び外挿）により、屈曲度合の鎖長依存性を考慮して、第１分子鎖モデル３と第２分子鎖モデル４とを対応付けることができる。これにより、第２分子鎖モデル４のビーズ数に比べて、モノマー数が小さい（即ち、鎖長の短い）第１分子鎖モデルを用いて内挿及び外挿することができるため、換算に要する時間を短縮することができる。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、第２分子鎖モデル４のビーズ数と高分子鎖２のモノマー数との比に基づいて、第２分子鎖モデル４の長さの単位を、高分子鎖２の長さの単位に換算する（工程Ｓ８）。第２分子鎖モデル４の長さの単位は、σ（シグマ）である。高分子鎖２の長さの単位は、ｍ（メートル）である。

上述したように、本実施形態において、第２分子鎖モデル４のビーズ数と高分子鎖２のモノマー数との比は、１個のビーズ１６あたりのモノマー数である。従って、工程Ｓ８では、第２分子鎖モデル４のビーズ数と高分子鎖２のモノマー数との比（１つのビーズ１６あたりのモノマー数）に、高分子鎖２のモノマー１つあたりの長さ（即ち、第１分子鎖モデル３の全長Ｌ／モノマー数）を乗じ、さらに、第２分子鎖モデル４のビーズ１つあたりの長さ（第２分子鎖モデル４の全長Ｌ／ビーズ数）で除することで、第２分子鎖モデル４の長さの単位を、高分子鎖２の長さの単位に換算できる。

本実施形態において、第２分子鎖モデル４のビーズ数と、高分子鎖２のモノマー数との比（即ち、１個のビーズ１６あたりのモノマー数）が、例えば１．５であるとする。第２分子鎖モデル４のビーズ１つあたりの長さが０．９６σであり、高分子鎖２（第１分子鎖モデル３）のモノマー１つあたりの長さが４．１×１０⁻¹⁰ｍであるとする。この場合、第２分子鎖モデル４のビーズ数と高分子鎖２のモノマー数との比（１．５）に、高分子鎖２のモノマー１つあたりの長さ（４．１×１０⁻¹⁰ｍ）を乗じ、さらに、第２分子鎖モデル４のビーズ１つあたりの長さ（０．９６σ）で除する。これにより、第２分子鎖モデル４の長さ（１σ）を、高分子鎖２の長さ６．４×１０⁻¹⁰ｍに換算することができる。

本実施形態では、第２分子鎖モデル４のビーズ数と、高分子鎖２のモノマー数との比を用いることにより、ユナイテッドアトムモデルよりも粗視化されたKremer Grestモデル等のバネビーズモデルの長さの単位を、高分子鎖２（全原子モデル）の長さの単位に精度よく換算できる。これにより、第２分子鎖モデル４を用いたシミュレーション結果を、高分子鎖２の実空間に対応付けることができるため、第１分子鎖モデル３を用いたシミュレーションに比べて、よりマクロな物性の定量的予測が可能となる。

次に、図５に示されるように、本実施形態のシミュレーション方法では、コンピュータ１が、第２分子鎖モデル４を用いた様々な条件のシミュレーションを実施する（工程Ｓ９）。工程Ｓ９のシミュレーションの一例としては、例えば、緩和状態の高分子材料モデル（図示省略）を用いた変形計算等を含んでいる。

次に、本実施形態のシミュレーション方法では、コンピュータ１が、高分子材料モデル（図示省略）を用いたシミュレーション結果が良好である（即ち、所望の性能を有する）か否かを判断する（工程Ｓ１０）。工程Ｓ１０では、工程Ｓ８での換算結果に基づいて、工程Ｓ９のシミュレーション結果から高分子鎖の物性等が評価される。

工程Ｓ１０において、シミュレーション結果が良好であると判断された場合（工程Ｓ１０で、「Ｙ」）、シミュレーションされた高分子材料モデル（図示省略）に基づいて、高分子材料が製造される（工程Ｓ１１）。他方、シミュレーション結果が良好でないと判断された場合（工程Ｓ１０で、「Ｎ」）、高分子鎖２（図２に示す）の構造や条件等を変更して（工程Ｓ１２）、工程Ｓ１〜工程Ｓ１０が再度実施される。

このように、本実施形態のシミュレーション方法では、所望の性能を有する未知の高分子材料を開発することができる。しかも、本実施形態のシミュレーション方法は、工程Ｓ９において、第１分子鎖モデル３（図３に示す）を用いたシミュレーションを行わなくても、第２分子鎖モデル４（図４に示す）を用いたシミュレーションにより、高い精度のシミュレーション結果を得ることができるため、高分子材料の開発コストを低減することができる。

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

図５に示した処理手順に従って、高分子鎖（cis-１．４ポリブタジエン）をモデル化した第１分子鎖モデル（図３に示す）及び第２分子鎖モデル（図４に示す）が設定された。そして、第１分子鎖モデルの第１パラメータと、第２分子鎖モデルの第２パラメータとに基づいて、第２分子鎖モデルのビーズ数と、高分子鎖のモノマー数との比が計算された。そして、第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比に基づいて、第２分子鎖モデルの長さの単位が、高分子鎖の長さの単位に換算された（実施例）。

実施例において、第１パラメータは、モノマー数が異なる第１分子鎖モデル毎に計算された。図１０は、第１分子鎖モデルの第１パラメータ（Kuhnのセグメント数）とモノマー数との第１関係を示すグラフである。なお、分子動力学計算において、空間の初期密度が０．００１g／cm³に設定され、温度が３６０Ｋに設定され、圧力が１atmに設定された。また、各全原子モデルのモノマー数（鎖長）、及び、空間に配置された全原子モデルの本数等については、次のとおりである。
第１分子鎖モデル（１）：
モノマー数：２０、空間内の本数：１５０本、緩和時間：１００ns
第１分子鎖モデル（２）：
モノマー数：３０、空間内の本数：７５本、緩和時間：５０ns
第１分子鎖モデル（３）：
モノマー数：４０、空間内の本数：７５本、緩和時間：１００ns

実施例において、第２パラメータは、ビーズ数が異なる第２分子鎖モデル毎に計算された。図１５は、第２分子鎖モデルの第２パラメータ（Kuhnのセグメント数）とビーズ数との第２関係を示すグラフである。各粗視化モデルのビーズ数（鎖長）、及び、空間に配置された粗視化モデルの本数等については、次のとおりである。
第２分子鎖モデル（１）：
ビーズ数：１０、空間内の本数：１０００本、緩和時間：１５０００τ
第２分子鎖モデル（２）：
ビーズ数：１３、空間内の本数：１０００本、緩和時間：１５０００τ
第２分子鎖モデル（３）：
ビーズ数：２０、空間内の本数：１０００本、緩和時間：１５０００τ
第２分子鎖モデル（４）：
ビーズ数：２７、空間内の本数：１０００本、緩和時間：１５０００τ
第２分子鎖モデル（５）：
ビーズ数：５０、空間内の本数：２００本、緩和時間：１５０００００τ
第２分子鎖モデル（６）：
ビーズ数：８０、空間内の本数：１００本、緩和時間：３００００００τ
第２分子鎖モデル（７）：
ビーズ数：１００、空間内の本数：１００本、緩和時間：３００００００τ
第２分子鎖モデル（８）：
ビーズ数：２００、空間内の本数：１００本、緩和時間：４５０００００τ
第２分子鎖モデル（９）：
ビーズ数：５００、空間内の本数：１５０本、緩和時間：３００００００τ

実施例において、第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比（即ち、１個のビーズあたりのモノマー数）は、図１０の第１関係を示すグラフと、図１５の第２関係を示すグラフとを、上記式（３）及び上記式（４）を用いてフィッティングすることによって計算された。

比較のために、上記特許文献１の方法に基づいて、第１分子鎖モデル及び第２分子鎖モデルの理想鎖のパラメータを求め、分子鎖モデルの理想鎖のパラメータと、全原子モデルの理想鎖のパラメータとが同一となるとみなして、分子鎖モデルの長さの単位を、高分子鎖の長さの単位に換算した（比較例）。なお、換算方法の詳細については、上記特許文献１のとおりである。また、比較例の第１分子鎖モデル及び第２分子鎖モデルについては、次のとおりである。
第１分子鎖モデル：
モノマー数：４０、空間内の本数：７５本、緩和時間：１００ns
第２分子鎖モデル：
ビーズ数：２７、空間内の本数：１０００本、緩和時間：１５０００τ

そして、実施例の換算精度及び計算時間と、比較例の換算精度及び計算時間とがそれぞれ比較された。実施例及び比較例で使用されたコンピュータは、２４コア（２ノード）である。

テストの結果、実施例の第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比（即ち、１個のビーズあたりのモノマー数）は、１．３６１±０．０１３であった。他方、比較例の第２分子鎖モデルのビーズ数と高分子鎖のモノマー数との比は、１．１９０±０．０２４であり、実施例と大きく異なった。これは、比較例が屈曲度合の鎖長依存性が考慮されていないためである。このように、実施例では、屈曲度合の鎖長依存性を考慮することができるため、精度よく換算できた。

実施例の計算時間が６３日（このうち、第１分子鎖モデルを用いた計算が１７日）であり、比較例の計算時間が８日（このうち、第１分子鎖モデルを用いた計算が７日）であった。また、比較例の方法で、実施例と同等の精度で換算するには、ほぼ鎖長依存性を無視できる長鎖長極限となる鎖長２５０モノマー程度（第２分子鎖モデルのビーズ数：１７０程度）での計算を実施する必要がある。この場合、モノマー数が７０の第１分子鎖モデルを用いて、最長緩和時間の３倍程度の時間の緩和計算を行うと、およそ３μs程度のシミュレーション時間が必要となり、実時間では、現在の一般的な計算機（２４コア）で、３年程度かかると推定される。従って、実施例は、同等の換算精度の比較において、比較例よりも短時間で換算することができた。

Ｓ３ 第１パラメータを計算する工程
Ｓ６ 第２パラメータを計算する工程
Ｓ７ 第２分子鎖モデルのビーズ数と、高分子鎖のモノマー数との比を計算する工程
Ｓ８ 第２分子鎖モデルの長さの単位を、高分子鎖の長さの単位に換算する工程

Claims (9)

1. コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、
   前記高分子鎖の全原子モデル又はユナイテッドアトムモデルからなる第１分子鎖モデルを、前記コンピュータに入力する第１モデル設定工程、
   前記コンピュータが、分子動力学計算に基づいて、前記第１分子鎖モデルの緩和した構造を計算する工程、
   前記コンピュータが、構造緩和後の前記第１分子鎖モデルの屈曲度合を示す第１パラメータを計算する第１計算工程、
   前記高分子鎖を複数のビーズで単純化した第２分子鎖モデルを、前記コンピュータに入力する第２モデル設定工程、
   前記コンピュータが、分子動力学計算に基づいて、前記第２分子鎖モデルの緩和した構造を計算する工程、
   前記コンピュータが、構造緩和後の前記第２分子鎖モデルの屈曲度合を示す第２パラメータを計算する第２計算工程、
   前記コンピュータが、前記第１パラメータと前記第２パラメータとに基づいて、前記第２分子鎖モデルのビーズ数と、前記高分子鎖のモノマー数との比を計算する第３計算工程、並びに
   前記第２分子鎖モデルの前記ビーズ数と前記高分子鎖の前記モノマー数との比に基づいて、前記第２分子鎖モデルの長さの単位を、前記高分子鎖の長さの単位に換算する工程を含むことを特徴とする高分子材料のシミュレーション方法。
2. 前記第１モデル設定工程は、前記モノマー数が異なる前記第１分子鎖モデルを設定する工程を含み、
   前記第１計算工程は、前記モノマー数が異なる前記第１分子鎖モデル毎に、前記第１パラメータを計算する工程を含む請求項１記載の高分子材料のシミュレーション方法。
3. 前記第２モデル設定工程は、前記ビーズ数が異なる前記第２分子鎖モデルを設定する工程を含み、
   前記第２計算工程は、前記ビーズ数が異なる前記第２分子鎖モデル毎に、前記第２パラメータを計算する工程を含む請求項１又は２記載の高分子材料のシミュレーション方法。
4. 前記第１パラメータ及び前記第２パラメータは、下記式（１）で計算される請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
   

   

   
   ここで、
   Ｎ_Kuhn：Kuhnのセグメント数
   Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
   Ｒ_g：第１分子鎖モデル又は第２分子鎖モデルの慣性半径
5. 前記第１パラメータ及び前記第２パラメータは、下記式（２）で計算される請求項１乃至３のいずれかに記載の高分子材料のシミュレーション方法。
   

   

   
   ここで、
   Ｎ_Kuhn：Kuhnのセグメント数
   Ｌ：第１分子鎖モデル又は第２分子鎖モデルの全長
   Ｒ：第１分子鎖モデル又は第２分子鎖モデルの末端間ベクトルの長さ
6. 前記第１分子鎖モデルは、前記高分子鎖のモノマーをモデル化したモノマーモデルを、予め定められた前記モノマー数で連結して設定され、
   前記第１分子鎖モデルの全長Ｌは、前記各モノマーモデルの長さを平均した平均モノマー長と、前記モノマー数とを乗じることにより求められる請求項４又は５に記載の高分子材料のシミュレーション方法。
7. 前記第２分子鎖モデルの全長Ｌは、隣接する前記ビーズの平衡核間距離と、前記ビーズ数とを乗じることによって求められる請求項４乃至６のいずれかに記載の高分子材料のシミュレーション方法。
8. 前記第２計算工程は、前記第２パラメータと前記ビーズ数との第２関係を求める工程を含み、
   前記第３計算工程は、前記第２関係を、下記式（３）で定義される曲線にフィッティングする工程を含む請求項４乃至７のいずれかに記載の高分子材料のシミュレーション方法。
   

   

   
   ここで、
   Ｎ_kuhn：Kuhnのセグメント数
   Ｎ_CG：第２分子鎖モデルのビーズ数
   ａ、ｂ、ｃ：フィッティングパラメータ
9. 前記第１計算工程は、前記第１パラメータと前記モノマー数との第１関係を求める工程を含み、
   前記第３計算工程は、下記式（４）に基づいて、前記第１関係を、前記曲線にフィッティングする工程を含む請求項８に記載の高分子材料のシミュレーション方法。
   

   

   
   ここで、
   Ｎ_FA：第１分子鎖モデルのモノマー数
   Ｎ_CG：第２分子鎖モデルのビーズ数
   ｄ：フィッティングパラメータ

Images