(19)【発行国】日本国特許庁(JP)
(12)【公報種別】特許公報(B2)
(11)【特許番号】
(24)【登録日】2023-09-11
(45)【発行日】2023-09-20
(54)【発明の名称】高分子材料のシミュレーション方法
(51)【国際特許分類】
G06F 30/10 20200101AFI20230912BHJP
G06F 30/20 20200101ALI20230912BHJP
【FI】
G06F30/10
G06F30/20
(21)【出願番号】P 2019208059
(22)【出願日】2019-11-18
【審査請求日】2022-09-16
(73)【特許権者】
【識別番号】000183233
【氏名又は名称】住友ゴム工業株式会社
(74)【代理人】
【識別番号】100104134
【氏名又は名称】住友 慎太郎
(74)【代理人】
【識別番号】100156225
【氏名又は名称】浦 重剛
(74)【代理人】
【識別番号】100168549
【氏名又は名称】苗村 潤
(74)【代理人】
【識別番号】100200403
【氏名又は名称】石原 幸信
(74)【代理人】
【識別番号】100206586
【氏名又は名称】市田 哲
(72)【発明者】
【氏名】馬場 昭典
【審査官】合田 幸裕
(56)【参考文献】
【文献】特開2019-159683(JP,A)
【文献】特開2017-049807(JP,A)
【文献】特開2017-220167(JP,A)
【文献】特開2006-338449(JP,A)
【文献】特開2005-032058(JP,A)
(58)【調査した分野】(Int.Cl.,DB名)
G06F 30/10
G06F 30/20
IEEE Xplore
JSTPlus(JDreamIII)
(57)【特許請求の範囲】
【請求項1】
コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、
前記高分子鎖を、前記高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、前記コンピュータに入力する工程と、
前記粗視化分子モデルの隣り合う前記粒子間に、Kuhn長とPacking長との比に影響する相互作用パラメータを含む相互作用を定義する工程と、
前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルを用いて、分子動力学に基づく構造緩和を計算する工程とを含み、
前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記高分子鎖の全原子モデル又は前記高分子鎖のユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記Kuhn長、前記Packing長、及び、前記Kuhn長と前記Packing長との比のいずれかに関する少なくとも2つの物性値がそれぞれ近づくように、前記高分子鎖と前記粗視化分子モデルとの間の長さの単位換算定数、及び、前記相互作用パラメータを定義する工程を含む、
高分子材料のシミュレーション方法。
【請求項2】
前記相互作用は、下記式(1)で定義される、請求項1記載の高分子材料のシミュレーション方法。
【数1】
ここで、
E:相互作用ポテンシャル関数
K:相互作用パラメータ
θ:隣り合う3つの粒子がなす角度
【請求項3】
前記単位換算定数、及び、前記相互作用パラメータを定義する工程は、前記鎖長及び前記相互作用パラメータがそれぞれ異なる複数の前記粗視化分子モデルを定義する工程と、
複数の前記粗視化分子モデルの前記物性値を、前記鎖長と前記相互作用パラメータとの関数に近似させる工程と、
前記鎖長と前記相互作用パラメータとの関数を用いて、前記単位換算定数、及び、前記相互作用パラメータを決定する工程とを含む、請求項1又は2記載の高分子材料のシミュレーション方法。
【請求項4】
前記鎖長と前記相互作用パラメータとの関数は、前記Kuhn長と前記Packing長との比に近似するものであり、
前記鎖長と前記相互作用パラメータとの関数は、下記式(2)で定義される、請求項3記載の高分子材料のシミュレーション方法。
【数2】
ここで、
f
lk/p:Kuhn長とPacking長との比の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
c
0~c
5:f
lk/pのフィッティングパラメータ
【請求項5】
前記鎖長と前記相互作用パラメータとの関数は、前記粗視化分子モデルの全長と前記Kuhn長との比であるKuhnセグメント数に近似するものであり、
前記鎖長と前記相互作用パラメータとの関数は、下記式(3)で定義される、請求項3又は4記載の高分子材料のシミュレーション方法。
【数3】
ここで、
f
Nk:Kuhnセグメント数の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ:粗視化分子モデルの数密度
f
lk/p:Kuhn長とPacking長との比の近似関数
【請求項6】
前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の一方がそれぞれ近づくように、前記相互作用パラメータを固定して、前記単位換算定数を求める第1工程と、
鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の他方がそれぞれ近づくように、前記単位換算定数を固定して、前記相互作用パラメータを求める第2工程と、
前記単位換算定数と前記相互作用パラメータとが収束するまで、前記第1工程と前記第2工程とを交互に実施させる工程とを含む、請求項1ないし5のいずれかに記載の高分子材料のシミュレーション方法。
【請求項7】
前記相互作用を定義する工程は、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかの密度の鎖長依存性を考慮して、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかのモノマー数と、前記粗視化分子モデルの粒子数とを対応付ける工程を含む、請求項1ないし6のいずれかに記載の高分子材料のシミュレーション方法。
【請求項8】
前記密度の鎖長依存性は、下記式(4)で定義され、
前記モノマー数に対応する前記粒子数は、下記式(5)で定義される、請求項7記載の高分子材料のシミュレーション方法。
【数4】
ここで、
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ρ
0、ρ
1:ρ
FAのフィッティングパラメータ
【数5】
ここで、
N
CG:モノマー数がN
FAの分子鎖に対応する粗視化分子モデルの分子鎖1本の粒子数
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
l
FA:モノマー1つあたりの長さ
C
l:長さの単位換算定数
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
ρ
0:フィッティングパラメータ
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、コンピュータを用いて、高分子材料を解析するための方法に関する。
【背景技術】
【0002】
下記特許文献1は、コンピュータを用いて、高分子材料を解析するための方法を提案している。この種の方法では、先ず、高分子鎖を複数のビーズで粗視化した粗視化モデルが、コンピュータに入力される。そして、分子動力学計算に基づいて、粗視化モデルの緩和した構造が計算される。
【0003】
一般に、上記の分子動力学計算では、現実の高分子鎖とは異なる単位系が用いられている。例えば、上記の方法では、粗視化分子モデルの全長と、Kuhn長との比であるKuhnセグメント数に基づいて、粗視化モデルの長さの単位と、全原子モデルの長さの単位とが関連付けられている。さらに、上記の方法では、平均二乗変位に基づいて、粗視化モデルの時間の単位と、高分子鎖の時間の単位とが関連付けられている。
【先行技術文献】
【特許文献】
【0004】
【発明の概要】
【発明が解決しようとする課題】
【0005】
特許文献1では、粗視化モデルの単位と、高分子鎖の単位とを関連付けることができるものの、その精度については、さらなる改善の余地があった。
【0006】
粗視化モデルの単位と、高分子鎖の単位とを精度よく関連付けるには、例えば、粗視化モデルと高分子鎖との間において、Kuhn長lkとPacking長pとの比(以下、単に「実効的な密度」ということがある。)lk/pなどの物性値を、互いに近づけることが重要である。このような物性値を精度良く求めるために、例えば、粗視化モデルを用いた分子動力学計算において、十分に長い鎖長のモデルを用いて、最長緩和時間よりも十分に長い緩和計算を行うと、膨大な計算時間が必要となるため、現実的ではない。一方、上記の物理量を現実的な時間で計算できるような十分に短い鎖長の粗視化モデルを用いると、粗視化モデルの鎖長に応じて上記の物性値が変化する傾向がある。したがって、粗視化モデルと高分子鎖との間において、上記の物性値を互いに近づけることは困難であり、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な方法が強く求められていた。
【0007】
本発明は、以上のような実状に鑑み案出されたもので、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な高分子材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【0008】
本発明は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法であって、前記高分子鎖を、前記高分子鎖を構成する原子の数よりも少ない複数の粒子を用いて表現した粗視化分子モデルを、前記コンピュータに入力する工程と、前記粗視化分子モデルの隣り合う前記粒子間に、Kuhn長とPacking長との比に影響する相互作用パラメータを含む相互作用を定義する工程と、前記コンピュータが、前記相互作用が定義された前記粗視化分子モデルを用いて、分子動力学に基づく構造緩和を計算する工程とを含み、前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記高分子鎖の全原子モデル又は前記高分子鎖のユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記Kuhn長、前記Packing長、及び、前記Kuhn長と前記Packing長との比のいずれかに関する少なくとも2つの物性値がそれぞれ近づくように、前記高分子鎖と前記粗視化分子モデルとの間の長さの単位換算定数、及び、前記相互作用パラメータを定義する工程を含むことを特徴とする。
【0009】
本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用は、下記式(1)で定義されてもよい。
【数1】
ここで、
E:相互作用ポテンシャル関数
K:相互作用パラメータ
θ:隣り合う3つの粒子がなす角度
【0010】
本発明に係る前記高分子材料のシミュレーション方法において、前記単位換算定数、及び、前記相互作用パラメータを定義する工程は、前記鎖長及び前記相互作用パラメータがそれぞれ異なる複数の前記粗視化分子モデルを定義する工程と、複数の前記粗視化分子モデルの前記物性値を、前記鎖長と前記相互作用パラメータとの関数に近似させる工程と、 前記鎖長と前記相互作用パラメータとの関数を用いて、前記単位換算定数、及び、前記相互作用パラメータを決定する工程とを含んでもよい。
【0011】
本発明に係る前記高分子材料のシミュレーション方法において、前記鎖長と前記相互作用パラメータとの関数は、前記Kuhn長と前記Packing長との比に近似するものであり、前記鎖長と前記相互作用パラメータとの関数は、下記式(2)で定義されてもよい。
【数2】
ここで、
f
lk/p:Kuhn長とPacking長との比の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
c
0~c
5:f
lk/pのフィッティングパラメータ
【0012】
本発明に係る前記高分子材料のシミュレーション方法において、前記鎖長と前記相互作用パラメータとの関数は、前記粗視化分子モデルの全長と前記Kuhn長との比であるKuhnセグメント数に近似するものであり、前記鎖長と前記相互作用パラメータとの関数は、下記式(3)で定義されてもよい。
【数3】
ここで、
f
Nk:Kuhnセグメント数の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ:粗視化分子モデルの数密度
f
lk/p:Kuhn長とPacking長との比の近似関数
【0013】
本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用を定義する工程は、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の一方がそれぞれ近づくように、前記単位換算定数を固定して、前記相互作用パラメータを求める第1工程と、鎖長が異なる複数の前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかと、鎖長が異なる複数の前記粗視化分子モデルとの間において、前記物性値の他方がそれぞれ近づくように、前記相互作用パラメータを固定して、前記単位換算定数を求める第2工程と、前記単位換算定数と前記相互作用パラメータとが収束するまで、前記第1工程と前記第2工程とを交互に実施させる工程とを含
んでもよい。
【0014】
本発明に係る前記高分子材料のシミュレーション方法において、前記相互作用を定義する工程は、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかの密度の鎖長依存性を考慮して、前記高分子鎖、前記全原子モデル又は前記ユナイテッドアトムモデルのいずれかのモノマー数と、前記粗視化分子モデルの粒子数とを対応付ける工程を含んでもよい。
【0015】
本発明に係る前記高分子材料のシミュレーション方法において、前記密度の鎖長依存性は、下記式(4)で定義され、前記モノマー数に対応する前記粒子数は、下記式(5)で定義されてもよい。
【数4】
ここで、
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ρ
0、ρ
1:ρ
FAのフィッティングパラメータ
【数5】
ここで、
N
CG:モノマー数がN
FAの分子鎖に対応する粗視化分子モデルの分子鎖1本の粒子数
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
l
FA:モノマー1つあたりの長さ
C
l:長さの単位換算定数
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
ρ
0:フィッティングパラメータ
【発明の効果】
【0016】
本発明は、上記の方法を採用することにより、現実的な計算時間で、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることが可能な高分子材料のシミュレーション方法を提供することができる。
【図面の簡単な説明】
【0017】
【
図1】本発明の高分子材料のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。
【
図3】高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。
【
図4】粗視化分子モデルの一例を示す概念図である。
【
図5】理想鎖に近似したときの粗視化分子モデルの部分拡大図である。
【
図6】相互作用(角度ポテンシャル)の一例を説明する図である。
【
図7】相互作用定義工程の処理手順の一例を示すフローチャートである。
【
図8】第1準備工程の処理手順の一例を示すフローチャートである。
【
図10】全原子モデルが配置された高分子材料モデルの一例を示す概念図である。
【
図11】各全原子モデルの実効的な密度l
k/pと、モノマー数N
FAとの関係を示すグラフである。
【
図12】第2準備工程の処理手順の一例を示すフローチャートである。
【
図13】粗視化分子モデルが配置された高分子材料モデルの一例を示す概念図である。
【
図14】パラメータ決定工程の処理手順の一例を示すフローチャートである。
【
図15】各粗視化分子モデルの実効的な密度l
k/p、相互作用パラメータK、及び、粒子数N
CGの関係の一例を示すグラフである。
【
図16】粗視化分子モデルのKuhnセグメント数N
kと粒子数N
CGとの関係、及び、全原子モデルのKuhnセグメント数N
kとモノマー数N
FAとの関係を示すグラフである。
【
図17】粗視化分子モデルの実効的な密度l
k/pと粒子数N
CGとの関係、及び、全原子モデルの実効的な密度l
k/pとモノマー数N
FAとの関係を示すグラフである。
【
図18】全原子モデルの密度と、モノマー数との関係を示すグラフである。
【発明を実施するための形態】
【0018】
以下、本発明の実施の一形態が図面に基づき説明される。
高分子材料のシミュレーション方法(以下、単に「シミュレーション方法」ということがある)は、コンピュータを用いて、高分子鎖を有する高分子材料を解析するための方法である。高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。
【0019】
図1は、本発明の高分子材料のシミュレーション方法を実行するコンピュータの一例を示す斜視図である。コンピュータ1は、本体1a、キーボード1b、マウス1c及びディスプレイ装置1dを含んで構成されている。この本体1aには、例えば、演算処理装置(CPU)、ROM、作業用メモリ、磁気ディスクなどの記憶装置、及び、ディスクドライブ装置1a1、1a2が設けられている。また、記憶装置には、本実施形態のシミュレーション方法を実行するためのソフトウェア等が予め記憶されている。
【0020】
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、高分子材料として、cis-1,4ポリブタジエン(以下、単に「ポリブタジエン」ということがある。)が例示される。
図2は、ポリブタジエンの構造式である。
【0021】
ポリブタジエンを構成する高分子鎖4は、メチレン基(-CH
2-)とメチン基(-CH-)とからなるモノマー5{-[CH
2-CH=CH-CH
2]-}が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基(-CH
2)に替えて、メチル基(-CH
3)が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。
図3は、高分子材料のシミュレーション方法の処理手順の一例を示すフローチャートである。
【0022】
本実施形態のシミュレーション方法では、先ず、高分子鎖4(
図2に示す)を表現した粗視化分子モデル6が、コンピュータ1に入力される(工程S1)。
図4は、粗視化分子モデルの一例を示す概念図である。
【0023】
粗視化分子モデル6は、高分子鎖4(
図2に示す)を構成する原子の数よりも少ない複数の粒子7を用いて表現されている。これらの粒子7、7間には、結合鎖8が結合されている。本実施形態の粗視化分子モデル6は、Kremer-Grestモデルである場合が例示されるが、特に限定されるわけではなく、例えば、DPD(散逸粒子動力学法)に基づくモデル等であってもよい。
【0024】
本実施形態の各粒子7は、例えば、高分子鎖4のモノマー5(
図2に示す)に対応している。高分子鎖4(
図2に示す)がポリブタジエンである場合には、例えば、1.10個分のモノマー5を構造単位として、該構造単位が1個の粒子7に置換される。これは、下記の論文1に基づくKremer-Grestモデルに、後述の式(1)で定義される相互作用を追加した粗視化分子モデル6について、この粗視化分子モデル6に対応付ける後述の全原子モデルの力場として、下記の論文2に示されるL-OPLS力場を用いた場合に、温度360Kの分子動力学計算に対応させるためである。また、相互作用パラメータKは、例えば、1.12に設定される。これにより、粗視化分子モデル6には、複数(例えば、10~5000個)の粒子7が設定される。工程S1において、高分子鎖4のモノマー数(鎖長)N
FAと、粗視化分子モデル6の粒子数(鎖長)N
CGとの比C
m(1個の粒子7あたりのモノマー数N
FA)は、例えば、1.10である。これらの比C
m、及び、相互作用パラメータKは、後述の相互作用定義工程S2において、より適切な値に決定(修正)される。
論文1: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
論文2: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
【0025】
粒子7は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、粒子7には、質量、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ1に記憶される。
【0026】
結合鎖8は、粒子7、7間に平衡長を定義した結合ポテンシャルとして構成される。ここで、「平衡長」とは、結合ポテンシャルの値が最も小さくなるような粒子7、7間の結合距離として定義される。また、結合鎖8の結合ポテンシャルは、例えば、上記の論文1に基づいて設定することができる。このような粗視化分子モデル6は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ1に入力される。
【0027】
ところで、粗視化分子モデル6は、例えば、高分子鎖4(
図2に示す)、全原子モデル(
図9に示す)又はユナイテッドアトムモデル(以下、これらをまとめて「高分子鎖4等」ということがある。)に比べて、大きな時空間を扱うことができるという利点がある。一方、粗視化分子モデル6には、現実の高分子鎖4とは異なる単位系が用いられている。このため、粗視化分子モデル6の計算結果を、高分子鎖4の運動として取り扱うためには、粗視化分子モデル6の単位を、高分子鎖4の単位に精度よく関連付けることが重要である。
【0028】
Kremer-Grestモデルや、DPD(散逸粒子動力学法)に基づくモデル等のように、個々の原子の座標を扱わない粗視化分子モデル6では、高分子鎖4(
図2に示す)の個性(例えば、太さや曲がりやすさ等)が十分に考慮されていない。このため、例えば、上記の特許文献1のように、高分子鎖4(
図2に示す)のモノマー数(鎖長)N
FAと、粗視化分子モデル6の粒子数(鎖長)N
CGとの比C
mに基づいて、粗視化分子モデル6の長さの単位と、高分子鎖4の長さの単位とを関連付けても、粗視化分子モデル6の密度が、高分子鎖4等の密度と必ずしも一致しないことがある。
【0029】
粗視化分子モデル6と高分子鎖4との間で長さの単位を関連付ける際に、密度も同時に関連付けるには、粗視化分子モデル6と高分子鎖4との間で、実効的な密度l
k/pが精度よく一致するように、粗視化分子モデル6が作成されるのが望ましい。ここで、実効的な密度l
k/pとは、下記の式(6)で定義されるKuhn長l
kと、下記の式(7)で定義されるPacking長pとの比で表される量である。これは、ゴム域の1本の高分子鎖4(
図2に示す)が動くことのできるチューブ状の空間の範囲(体積)あたりの高分子鎖4の本数に比例する量として近似的に扱うことができる。
【0030】
【数6】
ここで、
L:粗視化分子モデル、又は、高分子鎖(全原子モデル又はユナイテッドアトムモデル)の全長
R
g:粗視化分子モデル、又は、高分子鎖(全原子モデル又はユナイテッドアトムモデル)の慣性半径
【0031】
【数7】
ここで、
M:粗視化分子モデル、又は、高分子鎖(全原子モデル又はユナイテッドアトムモデル)の質量
R
g:粗視化分子モデル、又は、高分子鎖(全原子モデル又はユナイテッドアトムモデル)の慣性半径
ρ:粗視化分子モデル、又は、高分子鎖(全原子モデル又はユナイテッドアトムモデル)の密度
N
A:アボガドロ数
【0032】
図5は、理想鎖に近似したときの粗視化分子モデルの部分拡大図である。上記の式(6)のKuhn長l
kは、長さの単位を持つ量であり、高分子鎖4の剛直性(即ち、曲がりにくさ)を表すパラメータである。このKuhn長l
kは、高分子物理で知られている理想鎖6’に、高分子鎖4が近似したときにおいて、接続された粒子7’(Kuhnセグメント)の間の距離に対応する。なお、「高分子物理で知られている理想鎖6’に高分子鎖4が近似する」とは、慣性半径が再現されるように、粒子7’が一定の間隔でランダムな向きに(理想鎖の結合鎖8’で)接続された直鎖に、粗視化分子モデル6が近似することを意味している。
【0033】
一方、上記の式(7)のPacking長pは、長さの単位を持つ量であり、高分子鎖4の1本の占める体積と、末端間距離の二乗平均との比として定義される。このPacking長pは、ゴム域の1本の高分子鎖4が動くことのできるチューブ状の空間15の太さに比例する量として近似的に扱うことができる。上記の式(6)及び上記の式(7)では、高分子鎖4の末端間距離の二乗平均が、慣性半径の二乗平均の6倍として見積もられている。
【0034】
粗視化分子モデル6と高分子鎖4等との間で、実効的な密度lk/pの乖離(ズレ)が大きくなると、粗視化分子モデル6の長さの単位と、高分子鎖4の長さの単位との間で関連付け(換算)を行った際に、密度の乖離が大きくなる。例えば、上記の特許文献1の手法を用いて、粗視化分子モデル6と高分子鎖4との時間単位の換算定数を求める際に、粗視化分子モデル6の構造緩和した座標に重なるように、全原子モデル11の座標を設定して、全原子モデル11の構造緩和した座標を作成すると、得られた座標の密度が、平衡密度から乖離することがある。このような粗視化分子モデル6と、高分子鎖4等との組み合わせが用いられた場合、例えば、高分子鎖4の拡散運動が、本来の正比例の関係から外れてしまい、時間単位の換算定数を精度よく求めることが困難となるおそれがある。このような事象は、高分子鎖4の太さや曲がりやすさの個性うち、曲がりやすさ(Kuhn長に関する物性値)のみを考慮して単位換算すると生じやすい。
【0035】
上記のような問題を解決するには、例えば、高分子鎖の太さ(Packing長に関する物性値)と、高分子鎖4の曲がりやすさとの双方を考慮することが有効である。粗視化分子モデル6のKuhn長及びPacking長は、粗視化分子モデル6に定義される相互作用等に依存するため、相互作用を適切に定義することが重要である。とりわけ、Kuhn長とPacking長との比は無次元量であり、長さの単位換算に影響を受けない。このため、粗視化分子モデル6のKuhn長とPacking長との比を、高分子鎖4等のKuhn長とPacking長との比と一致させることが望ましい。
【0036】
相互作用については、粗視化分子モデル6のKuhn長とPacking長との比に影響を与えることができるものであれば、特に限定されない。本実施形態のシミュレーション方法では、粗視化分子モデル6に、隣り合う(連続する)3つの粒子がなす角度(結合角)θに応じた相互作用(角度ポテンシャル)が定義される(相互作用定義工程S2)。
図6は、相互作用(角度ポテンシャル)E(θ)の一例を説明する図である。
【0037】
相互作用(角度ポテンシャル)E(θ)は、角度(結合角)θと、相互作用ポテンシャル関数E(θ)の強度を示す相互作用パラメータKとの関数であり、例えば、下記の式(1)で定義される。
【0038】
【数1】
ここで、
E:相互作用ポテンシャル関数
K:相互作用パラメータ
θ:隣り合う3つの粒子がなす角度
【0039】
上記の式(1)において、相互作用パラメータKが正の場合には、結合角θが小さくなるほど、相互作用ポテンシャル関数E(θ)が大きくなり、粗視化分子モデル6が曲がりにくくなる。一方、相互作用パラメータKが負の場合、相互作用ポテンシャル関数E(θ)は、相互作用パラメータKが正の場合の逆の傾向となり、粗視化分子モデル6が曲がりやすくなる。このように、相互作用ポテンシャル関数E(θ)の大きさは、相互作用パラメータKの値に応じて適宜設定される。
【0040】
長さの単位に換算する際に、粗視化分子モデル6の密度を、高分子鎖4(
図2に示す)等のいずれかの密度に近づけるためには、相互作用パラメータKを適切に決定することが重要である。
【0041】
本実施形態の相互作用定義工程S2では、構造緩和の計算時において、鎖長が異なる複数の高分子鎖4等と、高分子鎖4等の鎖長に対応する複数の粗視化分子モデル6との間において、Kuhn長、Packing長、及び、Kuhn長と前記Packing長との比のいずれかに関する2つの物性値(以下、これらの2つの物性値を、「第1物性値」及び「第2物性値」として区別する場合がある。)が近づくように、相互作用(相互作用パラメータK)、及び、下記式(8)で定義される長さの単位換算定数Clが設定される。
【0042】
【数8】
ここで、
N
FA:前記高分子鎖、前記高分子鎖の全原子モデル又は前記高分子鎖のユナイテッドアトムモデルのいずれかの分子鎖1本のモノマー数
l
FA:高分子鎖内のモノマー1つあたりの長さ
N
CG:モノマー数がN
FAの分子鎖に対応する粗視化分子モデルの分子鎖1本の粒子数
l
CG:粗視化分子モデルの結合の平衡長
【0043】
2つの物性値(第1物性値及び第2物性値)については、Kuhn長lk、Packing長p、及び、Kuhn長とPacking長との比のいずれかに関するものであれば、適宜設定することができる。本実施形態の第1物性値は、Kuhn長とPacking長との比として得られる実効的な密度lk/pである場合が例示される。一方、本実施形態の第2物性値は、Kuhn長lkに関連する物理量であるKuhnセグメント数Nkである場合が例示される。このKuhnセグメント数Nkは、高分子鎖4や粗視化分子モデル6の全長LとKuhn長lkの比(即ち、Nk=L/lk)として定義される。
【0044】
粗視化分子モデル6と高分子鎖4等との間において、各物性値や鎖長の対応関係は、相互作用パラメータK、及び、長さの単位換算定数C
lの双方のパラメータに依存する。このため、各物性値や鎖長の対応関係を、自己無撞着に求めることが重要である。なお、高分子鎖4(
図2に示す)のモノマー数(鎖長)N
FAと、粗視化分子モデル6の粒子数(鎖長)N
CGとの比C
m(1個の粒子7あたりのモノマー数N
FA)については、上記の式(8)で定義される長さの単位換算定数C
lを用いて、C
m=N
FA/N
CG=C
l・l
CG/l
FAで定義することができる。
図7は、相互作用定義工程S2の処理手順の一例を示すフローチャートである。
【0045】
本実施形態の相互作用定義工程S2では、先ず、鎖長の異なる複数の高分子鎖4等のいずれか(本例では、全原子モデル11)について、第1物性値(本例では、実効的な密度lk/p)、及び、第2物性値(本例では、Kuhnセグメント数Nk)が取得される(第1準備工程S21)。なお、第1物性値及び第2物性値が既知である場合には、第1準備工程S21が省略されてもよい。
【0046】
本実施形態の第1準備工程S21では、高分子鎖4(
図2に示す)等のいずれかのうち、全原子モデル11の第1物性値(本例では、実効的な密度l
k/p)、及び、第2物性値(本例では、Kuhnセグメント数N
k)が取得される場合が説明される。
図8は、第1準備工程S21の処理手順の一例を示すフローチャートである。
【0047】
本実施形態の第1準備工程S21では、先ず、鎖長(即ち、重合度(モノマー数))N
FAが異なる複数の全原子モデル11(
図9に示す)が設定される(工程S31)。
図9は、全原子モデル11の一例を示す概念図である。
【0048】
全原子モデル11は、高分子鎖4の実際の構造に基づいて、原子をモデル化した原子モデル12で表現したものである。全原子モデル11を用いた分子動力学計算では、現実の高分子鎖4に基づいた長さの単位が用いられている。
【0049】
全原子モデル11は、複数の原子モデル12と、原子モデル12、12間を結合するボンド13とを含んで構成されている。全原子モデル11では、
図2に示した高分子鎖4のモノマー5を表す単位構造に基づいて、原子モデル12がボンド13で連結されることにより、モノマーモデル14が設定される。このモノマーモデル14が、予め定められた複数のモノマー数(鎖長)N
FA毎に連結される。これにより、モノマー数N
FAが異なる複数の全原子モデル11が設定される。
【0050】
モノマー数(鎖長)NFAについては、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ1の性能等に基づいて、構造緩和計算が現実的な計算時間で完了しうる範囲内で設定されるのが望ましい。なお、モノマー数NFAは、計算精度を維持するために、極端に小さい値を除外するのが望ましい。本実施形態では、例えば、10~60からモノマー数NFAが選択されるが、このような態様に限定されない。
【0051】
原子モデル12は、後述の分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、原子モデル12には、質量、直径、電荷、又は、初期座標などのパラメータが定義される。本実施形態の原子モデル12は、高分子鎖4の炭素原子をモデル化した炭素原子モデル12C、及び、高分子鎖4の水素原子をモデル化した水素原子モデル12Hを含んでいる。
【0052】
ボンド13は、原子モデル12、12間を拘束するためのものである。本実施形態のボンド13は、炭素原子モデル12C、12Cを連結する主鎖13A、及び、炭素原子モデル12Cと水素原子モデル12Hとの間を連結する側鎖13Bとを含んでいる。これらの主鎖13A及び側鎖13Bは、例えば、平衡長とバネ定数とが定義されたバネとして取り扱われる。
【0053】
全原子モデル11は、各原子モデル12、12間の結合長さである結合長、ボンド13を介して連続する3つの原子モデル12がなす角度である結合角、及び、ボンド13を介して連続する4つの原子モデル12において、隣り合う3つの原子モデル12が作る二面角などが定義される。これにより、全原子モデル11は、三次元構造を有する。全原子モデル11は、慣例に従い、外力又は内力を受けることによって、結合長、結合角及び二面角が変化する。これにより、全原子モデル11は、その三次元構造を変化させることができる。
【0054】
結合長、結合角及び二面角についてのポテンシャルは、例えば、上記の論文2や論文3(J. Comput. Chem. 25, 1157-1174 (2004))に示されるGAFFや、L-OPLS等に基づいて定義されうる。ポテンシャルは、高分子鎖4の構造に応じて設定されるのが望ましい。このような全原子モデル11は、材料物性シミュレーションソフトウェア(例えば、(株)JSOL社製のJ-OCTA)を用いて作成することができる。各全原子モデル11は、コンピュータ1で取り扱い可能な数値データであり、コンピュータ1に入力される。
【0055】
次に、本実施形態の第1準備工程S21では、モノマー数(鎖長)N
FAが異なる各全原子モデル11の初期配置が決定される(工程S32)。工程S32では、予め定められた空間に、全原子モデル11が配置された高分子材料モデルが設定される。
図10は、全原子モデル11が配置された高分子材料モデル18の一例を示す概念図である。
【0056】
本実施形態の工程S32では、モノマー数(鎖長)NFAが異なる各全原子モデル11が、独立して設けられた空間16にそれぞれ配置される。これにより、工程S32では、各全原子モデル11の初期配置が決定された高分子材料モデル18がそれぞれ定義される。各高分子材料モデル18には、同一のモノマー数NFAの全原子モデル11が、複数本配置されている。全原子モデル11は、例えば、モンテカルロ法に基づいて、空間16内に配置されるのが望ましい。
【0057】
空間16は、解析対象の高分子材料の微小構造部分に相当する。本実施形態の空間16は、互いに向き合う三対の平面17、17を有する立方体として定義されている。各平面17には、周期境界条件が定義されている。したがって、一方の平面17と、反対側の平面17とが連続している(繋がっている)ものとして取り扱うことができる。
【0058】
各空間16に配置される全原子モデル11の本数については、モノマー数(鎖長)NFAや、空間16の大きさに基づいて適宜設定することができる。全原子モデル11の本数としては、周期境界条件を介した同一分子鎖のコピー同士の絡まりを防ぐ観点より、平衡時の周期境界長が、平衡時の慣性半径の3倍以上に長くなりうる本数に設定されるのが望ましい。本実施形態の本数は、好ましくは20本以上であり、また、好ましくは200本以下である。また、空間16の一辺の長さLaは、例えば、系内の原子モデル12の初期密度が、例えば0.001g/cm3となるように設定されるのが望ましい。
【0059】
次に、本実施形態の第1準備工程S21では、全原子モデル11の初期配置が決定された高分子材料モデル18において、隣接する全原子モデル11、11の原子モデル12、12間に、相互作用ポテンシャルP1が定義される(工程S33)。
【0060】
相互作用ポテンシャルP1は、例えば、結合長、結合角及び二面角のポテンシャルの定義と対応するように、GAFFやL-OPLS等に基づいて定義されうる。本実施形態の相互作用ポテンシャルP1は、LJ(Lennard-Jones)ポテンシャルULJ(rij)であり、下記の式(9)で定義される。このような相互作用ポテンシャルP1は、原子モデル12、12間の距離rijに応じて、斥力及び引力を定義することができる。
【0061】
【数9】
ここで、各定数及び変数は、 LJポテンシャルのパラメータであり、次のとおりである。
r
ij:原子モデル間の距離
r
c:カットオフ距離
ε:原子モデル間に定義されるLJポテンシャルの強度
σ:原子モデルの直径に相当
なお、距離r
ij及びカットオフ距離r
cは、各原子モデル12、12の中心間の距離について定義される。
【0062】
相互作用ポテンシャルP1は、第1ポテンシャル、第2ポテンシャル及び第3ポテンシャルを含んでいる。第1ポテンシャルは、炭素原子モデル12C、12C(
図9に示す)間に設定される。第2ポテンシャルは、水素原子モデル12H、12H(
図9に示す)間に設定される。第3ポテンシャルは、炭素原子モデル12Cと水素原子モデル12Hとの間に設定される。なお、上記の式(9)中の各定数は、例えば、上記の論文2又は論文3に基づいて、適宜設定することができる。これらの全原子モデル11の初期配置が決定された高分子材料モデル18は、コンピュータ1に入力される。
【0063】
次に、本実施形態の第1準備工程S21では、初期配置された全原子モデル11の構造緩和が計算される(工程S34)。工程S34では、全原子モデル11の初期配置が決定された各高分子材料モデル18(
図10に示す)について、分子動力学計算が実施される。分子動力学計算は、例えば、
図10に示した空間16について所定の時間、配置した全ての全原子モデル11が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全ての原子モデル12の動きが追跡され、コンピュータ1に記憶される。また、分子動力学計算の条件は、例えば、系内の原子モデル12の個数、体積及び温度は一定で行われる。このような分子動力学計算は、例えば、分子動力学計算プログラムLAMMPSを用いて行うことができる。
【0064】
工程S34では、全原子モデル11の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、全原子モデル11の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。本実施形態では、例えば、1atmの条件下で体積がほぼ一定となり、かつ、その後のシミュレーション時間が最長緩和時間の3倍以上となったときに、全原子モデル11の初期配置が十分に構造緩和されたと判断される。最長緩和時間は、例えば、末端間ベクトルV1(
図9に示す)の自己相関関数を、指数減衰関数に近似させたときの指数緩和時間として見積ることができる。
【0065】
次に、本実施形態の第1準備工程S21では、全ての全原子モデル(即ち、モノマー数(鎖長)N
FAが異なる全ての全原子モデル)11について、構造緩和が計算されたか否かが判断される(工程S35)。工程S35において、全ての全原子モデル11の構造緩和が計算されたと判断された場合(工程S35で、「Y」)、次の工程S36が実施される。他方、全ての全原子モデル11の構造緩和計算が終了していないと判断された場合(工程S35で、「N」)、構造緩和が計算されていない全原子モデル11(即ち、
図10に示した高分子材料モデル18)が選択される(工程S37)。そして、第1準備工程S21では、工程S34及び工程S35が実施される。これにより、全ての全原子モデル11の熱平衡状態(構造緩和状態)がそれぞれ計算される。
【0066】
次に、本実施形態の第1準備工程S21では、
図9に示した全原子モデル11の第1物性値及び第2物性値が取得される(工程S36)。
【0067】
第1物性値及び第2物性値については、Kuhn長lk、Packing長p、及び、Kuhn長とPacking長との比(実効的な密度lk/p)のいずれかに関するものであれば、それぞれ適宜設定することができる。なお、Kuhn長lkは、上記の式(6)で定義される。Packing長pは、上記式(7)で定義される。上述したように、本実施形態では、第1物性値が、Kuhn長とPacking長との比として得られる実効的な密度lk/pであり、第2物性値が、Kuhn長lkに関連する物理量であるKuhnセグメント数Nkである。上述のとおり、Kuhnセグメント数Nkは、全長LとKuhn長lkとの比(即ち、Nk=L/lk)として定義される。
【0068】
上記の式(6)の全長Lには、全原子モデル11の全長L(図示省略)が代入される。この全長Lは、例えば、上記の特許文献1の方法に基づいて、全原子モデル11を主鎖方向に強制的に引き伸ばす変形計算によって求めることができる。本実施形態では、特許文献(特開2018-010525号公報)に基づいて、平均モノマー長lFAに、モノマー数(鎖長)NFAを乗じることにより、全長Lが計算される(即ち、L=lFA・NFA)。各全原子モデル11の平均モノマー長lFAは、構造緩和後の全原子モデル11において、一対のモノマー間を共有結合の中点で区切り、隣接する一対の中点の間の距離であるモノマー長を平均することで計算される。なお、末端モノマーの長さについては、平均モノマー長の計算から除外されている。構造緩和した全ての全原子モデル(即ち、モノマー数NFAが異なる全ての全原子モデル)11について、それらの全長Lがそれぞれ計算される。
【0069】
上記の式(6)及び(7)の慣性半径Rgには、全原子モデル11の慣性半径Rgが代入される。この慣性半径Rgは、例えば、上記の論文1に記載されている式2.5に基づいて、工程S34の分子動力学計算で求められた全原子モデル11の原子モデル12の座標値から計算される。本実施形態では、構造緩和した全ての全原子モデル(即ち、モノマー数(鎖長)NFAが異なる全ての全原子モデル)11について、それらの慣性半径Rgがそれぞれ計算される。
【0070】
上記の式(7)の質量Mには、1本の全原子モデル11を構成する原子モデル12の合計質量が代入される。また、上記の式(7)の密度ρには、全原子モデル11の密度ρが代入される。これらの計算により、工程S36では、第1物性値(本例では、実効的な密度l
k/p)、及び、第2物性値(本例では、Kuhnセグメント数N
k)を取得することができる。第1物性値及び第2物性値は、コンピュータ1に記憶される。
図11は、各全原子モデル11の実効的な密度l
k/pと、モノマー数(鎖長)N
FAとの関係を示すグラフである。
図11では、実効的な密度l
k/pとモノマー数N
FAとの関係が、点23でプロットされている。
【0071】
次に、本実施形態の相互作用定義工程S2では、鎖長、及び、屈曲性に影響を与える相互作用がそれぞれ異なる複数の粗視化分子モデルについて、第1物性値及び第2物性値が取得される(第2準備工程S22)。
図12は、第2準備工程S22の処理手順の一例を示すフローチャートである。
【0072】
本実施形態の第2準備工程S22では、先ず、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6が設定される(工程S41)。本実施形態の工程S41では、複数の異なる粒子数NCGに基づいて、粒子7が結合鎖8で連結されることにより、粒子数NCGが異なる複数の粗視化分子モデル6が設定される。
【0073】
複数の異なる粒子数(鎖長)NCGについては、高分子材料の種類や、後述の分子動力学計算を実施するコンピュータ1の性能等に基づいて、適宜設定することができる。本実施形態の粒子数NCGは、例えば、5~1000から選択されているが、このような態様に限定されない。これらの粗視化分子モデル6は、高分子材料を分子動力学計算で取り扱うための数値データであり、コンピュータ1に入力される。
【0074】
次に、本実施形態の第2準備工程S22では、相互作用パラメータ(例えば、上記の式(1)の相互作用パラメータK)が互いに異なる複数の粗視化分子モデル6が定義される(工程S42)。
図6に示されるように、本実施形態の工程S42では、例えば、互いに異なる複数の相互作用パラメータK毎に、粗視化分子モデル6の隣り合う粒子7、7間に相互作用(角度ポテンシャル)E(θ)がそれぞれ定義される。
【0075】
互いに異なる相互作用パラメータKについては、例えば、高分子鎖4(
図2に示す)がとり得る範囲の実効的な密度l
k/pに基づいて、適宜設定することができる。本実施形態の相互作用パラメータKは、例えば、高分子鎖4等(本例では、全原子モデル11)の実効的な密度l
k/pが2.5~7.0程度の場合、0~1.5の範囲から選択される。この相互作用パラメータKは、その値が小さいほど、全原子モデル11が曲がりやすくなる。一方、相互作用パラメータKの値が大きいほど、全原子モデル11が曲がりにくくなる(剛直性が高くなる)。
【0076】
本実施形態の工程S42では、工程S41で設定された粒子数(鎖長)N
CGが異なる各粗視化分子モデル6について、相互作用パラメータKが異なる複数の粗視化分子モデル6がそれぞれ設定される。即ち、本実施形態の工程S42では、粒子数N
CGの種類(例えば、5種類(後述の
図15に示したN
CG=10、20、50、100、150))に、相互作用パラメータKの種類(例えば、4種類(
図15に示す))を乗じた種類(例えば、20種類)の粗視化分子モデル6が設定される。粒子数N
CG及び相互作用パラメータKが異なる各粗視化分子モデル6は、コンピュータ1に記憶される。
【0077】
次に、本実施形態の第2準備工程S22では、粒子数(鎖長)N
CG及び相互作用パラメータKが異なる各粗視化分子モデル6の初期配置が決定される(工程S43)。本実施形態の工程S43では、予め定められた空間に、各粗視化分子モデル6がそれぞれ配置された高分子材料モデルが設定される。
図13は、粗視化分子モデル6が配置された高分子材料モデル20の一例を示す概念図である。
【0078】
本実施形態の工程S43では、粒子数(鎖長)NCG及び相互作用パラメータKが異なる各粗視化分子モデル6が、独立して設けられた空間16にそれぞれ配置される。これにより、工程S43では、各粗視化分子モデル6の初期配置が決定された高分子材料モデル20がそれぞれ定義される。各高分子材料モデル20において、空間16には、同一の粒子数NCG、かつ、同一の相互作用パラメータKの粗視化分子モデル6が複数本配置されている。粗視化分子モデル6は、モンテカルロ法に基づいて、空間16内に配置されるのが望ましい。
【0079】
空間16は、解析対象の高分子材料の微小構造部分に相当し、
図10に示した空間16と同様に設定される。各空間16に配置される粗視化分子モデル6の本数については、粒子数(鎖長)N
CGや、空間16の大きさに基づいて適宜設定することができる。空間16に配置される粗視化分子モデル6の本数は、好ましくは50本以上であり、また、好ましくは2000本以下である。
【0080】
本実施形態の工程S43では、各高分子材料モデル20において、隣接する粗視化分子モデル6、6の粒子7、7間に、相互作用ポテンシャルP2がそれぞれ定義される。本実施形態の相互作用ポテンシャルP2は、LJポテンシャルであり、上記の式(9)で定義される。なお、上記の式(9)において、「原子モデル間」は、「粒子間」に置き換えて適用される。
【0081】
相互作用ポテンシャルP2の強度ε、相互作用ポテンシャルP2が作用する距離σ、及び、カットオフ距離rcは、例えば、上記の論文2又は論文3に基づいて設定されるのが望ましい。これらの粗視化分子モデル6の初期配置が決定された高分子材料モデル20は、コンピュータ1に入力される。
【0082】
次に、本実施形態の第2準備工程S22では、初期配置された粗視化分子モデル6の構造緩和が計算される(工程S44)。本実施形態の工程S44では、粗視化分子モデル6の初期配置が決定された各高分子材料モデル20(
図13に示す)について、分子動力学計算が実施される。工程S44では、各時刻での全ての粒子7の動きが追跡され、コンピュータ1に記憶される。
【0083】
本実施形態の工程S44では、粗視化分子モデル6の初期配置が十分に構造緩和されるまで、分子動力学計算が実施される。なお、初期配置の構造緩和の判断基準については、粗視化分子モデル6の人為的な初期配置が十分に排除されたとみなせる基準であれば、適宜設定されうる。本実施形態では、シミュレーション時間が、最長緩和時間の3倍以上となったときに、粗視化分子モデル6の初期配置が十分に構造緩和されたと判断している。
【0084】
次に、本実施形態の第2準備工程S22では、全ての粗視化分子モデル(即ち、粒子数(鎖長)NCG及び相互作用パラメータKが異なる全ての粗視化分子モデル)6について、構造緩和が計算されたか否かが判断される(工程S45)。工程S45において、全ての粗視化分子モデル6の構造緩和が計算されたと判断された場合(工程S45で、「Y」)、次の工程S46が実施される。
【0085】
一方、工程S45において、全ての粗視化分子モデル6の構造緩和計算が終了していないと判断された場合(工程S45で、「N」)、構造緩和が計算されていない粗視化分子モデル6(即ち、
図13に示した高分子材料モデル20)が選択され(工程S47)、工程S44及び工程S45が実施される。これにより、第2準備工程S22では、全ての粗視化分子モデル6の熱平衡状態(構造緩和状態)がそれぞれ計算される。
【0086】
次に、本実施形態の第2準備工程S22では、全ての粗視化分子モデル(即ち、粒子数(鎖長)NCG及び相互作用パラメータKが異なる全ての粗視化分子モデル)6について、第1物性値及び第2物性値が取得される(工程S46)。
【0087】
第1物性値及び第2物性値については、Kuhn長lk、Packing長p、及び、Kuhn長とPacking長との比(実効的な密度lk/p)のいずれかに関するものであれば、それぞれ適宜設定することができる。なお、Kuhn長lkは、上記式(6)で定義される。Packing長pは、上記式(7)で定義される。上述したように、本実施形態では、第1物性値が、Kuhn長とPacking長との比として得られる実効的な密度lk/pであり、第2物性値が、Kuhn長lkに関連する物理量であるKuhnセグメント数Nkである。上述のとおり、Kuhnセグメント数Nkは、全長LとKuhn長lkの比Nk=L/lkとして定義される。
【0088】
上記の式(6)の全長Lには、粗視化分子モデル6の全長L(図示省略)が代入される。この全長Lは、結合の平衡長lCGに、粒子数(鎖長)NCGを乗じることによって計算することができる(即ち、L=lCG・NCG)。上記の結合の平衡長lCGの値としては、FENEポテンシャルの平衡長(0.961σ)の代わりに、温度の影響等を考慮して補正した値(例えば、0.965σ)が用いられるのが望ましい。本実施形態では、構造緩和した全ての粗視化分子モデル(即ち、粒子数NCGが異なる全ての粗視化分子モデル)6について、それらの全長Lがそれぞれ計算される。
【0089】
上記の式(6)及び(7)の慣性半径Rgには、粗視化分子モデル6の慣性半径Rgが代入される。この慣性半径Rgは、例えば、上記の論文1に記載されている式2.5に基づいて、工程S44の分子動力学計算で求められた粗視化分子モデル6の粒子7の座標値から計算される。本実施形態では、構造緩和した全ての粗視化分子モデル(即ち、粒子数(鎖長)NCGが異なる全ての粗視化分子モデル)6について、それらの慣性半径Rgがそれぞれ計算される。
【0090】
上記の式(7)の質量Mには、1本の粗視化分子モデル6を構成する粒子7の合計質量が代入される。また、上記の式(7)の密度ρには、粗視化分子モデル6の密度ρ(例えば、0.85σ-3)が代入される。
【0091】
次に、
図7に示されるように、本実施形態の相互作用定義工程S2では、鎖長が異なる複数の全原子モデル11と、鎖長が異なる複数の粗視化分子モデル6との間において、第1物性値(本例では、実効的な密度l
k/p)及び第2物性値(本例では、Kuhnセグメント数N
k)が近づくように、長さの単位換算定数C
lと、相互作用パラメータKとが決定される(パラメータ決定工程S23)。本実施形態のパラメータ決定工程S23では、例えば、第1準備工程S21の結果と、第2準備工程S22の結果とを比較することによって、長さの単位換算定数C
lと相互作用パラメータKとを決定している。
【0092】
パラメータ決定工程S23の処理手順の一例としては、先ず、相互作用パラメータKについて、複数の候補値が決定される。これらの候補値には、第2準備工程S22で用いられた互いに異なる相互作用パラメータK(本例では、上述の0~1.5)が用いられる。
【0093】
次に、長さの単位換算定数Clについて、複数の候補値が決定される。これらの候補値は、互いに異なる長さの単位換算定数Clが適宜選択される。候補値が選択される数値範囲の一例としては、0.3~1.0である。
【0094】
次に、長さの単位換算定数Clの候補値について、第1準備工程S21で用いた各鎖長(モノマー数)NFAに最も近似する第2準備工程S22の鎖長(粒子数)NCGが、上記の式(8)に基づいて、それぞれ対応付けられる(即ち、NCG=NFA・lFA/(lCG・Cl))。なお、粒子数NCGが整数である必要がある場合には、小数点以下が四捨五入される。
【0095】
そして、長さの単位換算定数Clの候補値と、相互作用パラメータKの候補値とのペアそれぞれについて、例えば、対応する高分子鎖等(本例では、全原子モデル11)と、粗視化分子モデル6との第1物性値(本例では、実効的な密度lk/p)のズレの最小二乗和と、第2物性値(本例では、Kuhnセグメント数Nk)のズレの最小二乗和との積をそれぞれ算出し、算出した積が最も小さくなるような候補値のペアが、長さの単位換算定数Clと相互作用パラメータKとして決定される。
【0096】
このような方法は、単純な処理手順で、長さの単位換算定数Clと相互作用パラメータKとを決定することができるが、精度の高い結果を得るためには、第2準備工程S22において、鎖長(粒子数NCG)及び相互作用パラメータKが異なる多くの粗視化分子モデル6を用いた計算が必要となる。このため、計算コスト(即ち、計算時間及び計算費用)が増大しやすい欠点がある。
【0097】
上記のような欠点を克服可能な好ましい方法として、本実施形態のパラメータ決定工程S23の処理手順の一例について説明する。
図14は、パラメータ決定工程S23の処理手順の一例を示すフローチャートである。
【0098】
本実施形態のパラメータ決定工程S23では、先ず、長さの単位換算定数Clと相互作用パラメータKとを決定する前に、先ず、第2準備工程S22で得られた鎖長(粒子数)NCG及び相互作用パラメータKが異なる複数の粗視化分子モデル6について、第1物性値(本例では、実効的な密度lk/p)及び第2物性値(本例では、Kuhnセグメント数Nk)を、鎖長(粒子数)NCGと相互作用パラメータKとの関数に近似させる(工程S51)。これにより、本実施形態では、鎖長(粒子数)NCG及び相互作用パラメータKが異なる多くの粗視化分子モデル6を用いなくても、任意の鎖長(粒子数)NCG及び相互作用パラメータKから、第1物性値及び第2物性値を精度よく推定することができる。したがって、本実施形態の方法では、計算コストを削減することが可能となる。
【0099】
鎖長(粒子数)NCGと相互作用パラメータKの関数については、適宜定義することができる。本実施形態では、粗視化分子モデル6のKuhn長とPacking長との比(即ち、実効的な密度lk/p)を近似させる関数、及び、Kuhnセグメント数Nkを近似させる関数が用いられるのが望ましい。鎖長と相互作用パラメータKの関数は、下記の式(2)及び式(3)で定義されるのが好ましい。なお、下記の式(2)で定義される関数は、実効的な密度lk/pに近似させるものである。下記の式(3)で定義される関数は、Kuhnセグメント数Nkを近似させるものである。
【0100】
【数2】
ここで、
f
lk/p:Kuhn長とPacking長との比(実効的な密度l
k/p)の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
c
0~c
5:f
lk/pのフィッティングパラメータ
【0101】
【数3】
ここで、
f
Nk:Kuhnセグメント数の近似関数
N
CG:粗視化分子モデルの鎖長
K:相互作用パラメータ
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ:粗視化分子モデルの数密度
f
lk/p:Kuhn長とPacking長との比の近似関数
【0102】
上記の式(2)及び式(3)のフィッティングパラメータc
0~c
5は、第2準備工程S22で得られた各粗視化分子モデル6(即ち、粒子数(鎖長)N
CG及び相互作用パラメータKが異なる各粗視化分子モデル6)の物性値(実効的な密度l
k/p及びKuhnセグメント数N
k)に、上記の式(2)及び式(3)がフィッティングするように調節することで得られる。フィッティングには、例えば、Levenberg-Marquardt法を用いることができる。
図15は、各粗視化分子モデルの実効的な密度l
k/p、相互作用パラメータK、及び、粒子数N
CGの関係の一例を示すグラフである。
図15では、上記の式(2)のフィッティングの例を示している。フィッティングパラメータc
0~c
5は、コンピュータ1に記憶される。
【0103】
次に、本実施形態のパラメータ決定工程S23では、次の工程S52~S56において、鎖長(モノマー数)NFAが異なる高分子鎖4等(本例では、全原子モデル11)と、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間において、第1物性値と第2物性値の鎖長依存性が近づくように、工程S51で得られた上記の式(2)及び式(3)の近似関数を用いて、長さの単位換算定数Clと相互作用パラメータKが決定される。
【0104】
上記の2つのパラメータ(本例では、長さの単位換算定数Cl及び相互作用パラメータK)を決定するには、例えば、2つのパラメータの初期値(推定値)に基づいて、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等の第1物性値及び第2物性値と、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6の上記の近似関数で推定される第1物性値及び第2物性値とを比較し、各鎖長で第1物性値及び第2物性値が近くなるように、それらの第1物性値のズレの最小二乗和と、第2物性値のズレの最小二乗和との積とが小さくなるような2つのパラメータを探索する方法が考えられる。なお、粗視化分子モデル6の鎖長(粒子数)NCGは、長さの単位換算定数Clを用いて、高分子鎖4の鎖長(モノマー数NFA)に対応付けられる。また、上述の探索には、例えば、差分進化法などの遺伝的アルゴリズムを用いることができる。この方法では、第1物性値及び第2物性値や、これらの近似関数に依存することなく汎用的に用いることができる。しかしながら、この方法では、パラメータの探索範囲を事前に指定する必要があるなどの制限がある。
【0105】
本実施形態のパラメータ決定工程S23は、上述の差分進化法などの遺伝的アルゴリズムを用いずに、次の工程S52~工程S56において、2つのパラメータ(本例では、長さの単位換算定数Cl及び相互作用パラメータK)を、1つずつ交互に、Levenberg-Marquardt法などの非線形最小二乗法で更新しながら、2つのパラメータの双方が収束するまで繰り返し実施される。
【0106】
本実施形態のパラメータ決定工程S23では、高分子鎖4等と、粗視化分子モデル6との間で、第1物性値(本例では、実効的な密度lk/p)の長鎖長極限(極限値)が一致するように、相互作用パラメータKが選択される(工程S52)。
【0107】
本実施形態の工程S52では、先ず、第1準備工程S21で得られた高分子鎖4等の各鎖長(モノマー数)NFAでの第1物性値(以下、「第1物性値AFA」ということがある。)に、高分子鎖4等の第1物性値AFAと、高分子鎖4等の鎖長(モノマー数)NFAとの関係を示す関数AFA(NFA)をフィッティングさせる。関数AFA(NFA)については、適宜定義することができる。本実施形態の関数AFA(NFA)は、AFA(NFA)=A0+A1/NFAで定義される。A0及びA1は、フィッティングパラメータである。このような関数AFA(NFA)は、鎖長依存性を有する第1物性値AFAを、鎖長NFAから特定することができる。
【0108】
上記のフィッティングは、フィッティングパラメータA0及びA1を調節することで行うことができる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。フィッティングパラメータA0及びA1は、コンピュータ1に記憶される。
【0109】
次に、本実施形態の工程S52では、上記の関数A
FA(N
FA)について、鎖長(本例では、モノマー数)N
FAを無限大(即ち、長鎖長極限)とする第1物性値A
FAの極限値が求められる(本例では、
図11の実効的な密度の長鎖長極限l
k/p
∞が第1物性値A
FAの極限値A
0となる)。なお、第1物性値A
FAの極限値が発散する場合には、鎖長N
FAを無限大にしたときに、第1物性値A
FAの極限値が一定値に近づくように規格化した関数A’
FA(N
FA)を定義して、第1物性値A
FAの極限値に代えて、関数A’
FA(N
FA)の極限値が求められるのが望ましい。
【0110】
関数A’FA(NFA)は、上記の規格化のために適宜定義した関数A’’FA(NFA)で、関数AFA(NFA)を除することによって定義することができる(即ち、A’FA(NFA)=AFA(NFA)/A’’FA(NFA))。例えば、第1物性値AFAが本実施形態のような実効的な密度lk/pではなく、Kuhnセグメント数Nkの場合には、関数AFA(NFA)が鎖長(モノマー数NFA)に比例して増加する傾向があるため、関数A’’FA(NFA)がモノマー数NFAと等しい(即ち、A’’FA(NFA)=NFA)とみなして、関数A’FA(NFA)は、AFA(NFA)/NFAで定義することができる。
【0111】
次に、本実施形態の工程S52では、粗視化分子モデル6の第1物性値(以下、「第1物性値ACG」ということがある。)について、粗視化分子モデル6の鎖長(本例では、粒子数)NCG及び相互作用パラメータKとの関係を示す関数ACG(NCG,K)を定義する。関数ACG(NCG,K)については、適宜定義することができる。本実施形態の関数ACG(NCG,K)は、上記の式(2)に基づいて定義することができる(即ち、ACG(NCG,K)=flk/p(NCG,K))。このような関数ACG(NCG,K)は、鎖長依存性及び相互作用パラメータ依存性を有する第1物性値ACGを、鎖長NCG及び相互作用パラメータKから特定することができる。
【0112】
そして、本実施形態の工程S52では、関数ACG(NCG,K)に、第2準備工程S22で得られた粗視化分子モデル6の各鎖長NCG及び相互作用パラメータKでの第1物性値ACGをフィッティングさせる。本実施形態のフィッティングは、上記の式(2)のフィッティングパラメータc0~c5を調節することで行うことができる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。フィッティングパラメータ(本例では、c0~c5)は、コンピュータ1に記憶される。
【0113】
次に、本実施形態の工程S52では、関数A
CG(N
CG,K)について、鎖長(粒子数)N
CGを無限大(即ち、長鎖長極限)とする第1物性値A
CGの極限値の関数lim
NCG→∞A
CG(N
CG,K)が定義される。この第1物性値A
CGの極限値の関数は、
図15の近似曲線28で例示される。なお、第1物性値A
CGの極限値が発散する場合には、鎖長N
CGを無限大にしたときに、第1物性値A
CGの極限値の関数が有限となるように規格化した関数A’
CG(N
CG,K)を定義して、第1物性値A
CGの極限値の関数の代わりに、A’
CG(N
CG,K)の極限値の関数lim
NCG→∞A’
CG(N
CG,K)が定義されるのが望ましい。
【0114】
関数A’CG(NCG,K)は、上述の高分子鎖4等の関数A’ ’FA(NFA)を用いて、鎖長(モノマー数)NFAを、鎖長(粒子数)NCGで対応付けたもの(NCG=NFA・lFA/(Cl・lCG))によって置き換えることで(即ち、A’CG(NCG,K)=ACG(NCG,K)/A’’FA(NCG・(Cl・lCG)/lFA))、定義することができる。なお、上記の対応付けには、長さの単位換算定数Cl、並びに、上記の式(8)に示される粗視化分子モデルの結合の平衡長lCG、及び、高分子鎖内のモノマー1つあたりの長さlFAが用いられる。
【0115】
次に、本実施形態の工程S52では、高分子鎖4等の第1物性値A
FAの極限値(本例では、
図11の実効的な密度の長鎖長極限l
k/p
∞)と、粗視化分子モデル6の第1物性値A
CGの極限値の関数lim
NCG→∞A
CG(N
CG,K)(
図15の近似曲線28に例示される)が一致するような相互作用パラメータKが選択される。第1物性値A
FA及びA
CGが長鎖長極限で発散する場合には、第1物性値A
FAの極限値と、第1物性値A
CGの極限値の関数lim
NCG→∞A
CG(N
CG,K)の代わりに、関数A’
FA(N
FA)の極限値と、A’
CG(N
CG,K)の極限値の関数lim
NCG→∞A’
CG(N
CG,K)とが一致するような相互作用パラメータKが選択される。相互作用パラメータKの選択には、例えば、2分法などの求根アルゴリズムを用いることができる。なお、上記の極限値が一致するような相互作用パラメータKが存在しない場合には、第2準備工程S22において、鎖長(粒子数)N
CG及び相互作用パラメータKの範囲を大きくして、これまでの手順がやり直されてもよい。
【0116】
上述したように、本実施形態の第1物性値AFA及びACGは、実効的な密度lk/pである。この実効的な密度lk/pは、無次元量であり、長さの単位換算定数Clに依存しない。このため、高分子鎖4等と粗視化分子モデル6との間で、長さの単位換算定数Clを用いた単位換算を行う必要がなく、高分子鎖4等の第1物性値AFAの極限値と、粗視化分子モデル6の第1物性値ACGの極限値とを直接比較することができる。
【0117】
なお、第1物性値AFA、ACGが無次元量ではない物性値(例えば、Kuhn長lk)である場合には、高分子鎖4等と粗視化分子モデル6との間で、長さの単位換算定数Clを用いた単位換算を行う必要がある。この長さの単位換算には、長さの単位換算定数Clの初期値として、例えば、0.5nm/σが用いられる。これは、1σ=0.5nmを意味する。第1物性値AFA及びACGとしてKuhn長lkが用いられた場合、粗視化分子モデル6のKuhn長(単位は「σ」)は、高分子鎖4等のKuhn長(単位は「nm」)から、長さの単位換算定数Clを除することで換算することができる。
【0118】
本実施形態では、上述のとおり、相互作用が上記の式(1)で定義され、かつ、第1物性値が実効的な密度lk/pであり、かつ、粗視化分子モデル6の実効的な密度lk/pを上記の式(2)に近似させている。この場合、粗視化分子モデル6の第1物性値(実効的な密度lk/p)の極限値は、上記の式(2)の近似関数から、lim NCG→∞flk/p(NCG,K)=c0・exp[c3K(1+c5K)]として特定することができる。
【0119】
一方、高分子鎖4等の第1物性値(実効的な密度lk/p)を、上記の関数AFA(NFA)=A0+A1/NFAに近似させた場合、高分子鎖4等の第1物性値(実効的な密度lk/p)の極限値は、フィッティングパラメータA0と同一の値となる。
【0120】
高分子鎖4等の第1物性値の極限値(フィッティングパラメータA0)と、粗視化分子モデル6の第1物性値の極限値(c0・exp[c3K(1+c5K)])とが一致するような相互作用パラメータKは、例えば求根アルゴリズムを用いなくても、解析的に求めることができる。本実施形態のように、例えば、相互作用パラメータK≧0と仮定した場合、相互作用パラメータKは、K=-1+sqrt[1+4(c5/c3)ln(A0/c0)]/(2c5)となる。このように、実施形態によっては、求根アルゴリズムを用いることなく、相互作用パラメータKを解析的に求めることができる場合がある。
【0121】
本実施形態の工程S52では、高分子鎖4等と、粗視化分子モデル6との間で、第1物性値(本例では、実効的な密度lk/p)の長鎖長極限(極限値)が一致するように、相互作用パラメータKが選択されたが、このような態様に限定されるわけではなく、例えば、任意の相互作用パラメータKが指定されてもよい。選択された相互作用パラメータKは、コンピュータ1に記憶される。
【0122】
次に、本実施形態のパラメータ決定工程S23では、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等と、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間において、物性値の一方(本例では、第2物性値)がそれぞれ近づくように、相互作用パラメータKを固定して、長さの単位換算定数Clが求められる(第1工程S53)。本実施形態の第1工程S53では、工程S51で得られた第2物性値の近似関数(本例では、上記の式(3))を用いて、工程S52で選択された相互作用パラメータKで固定して、第2物性値(Kuhnセグメント数Nk)がそれぞれ近づくように(本例では、一致するように)、長さの単位換算定数Clが求められる。
【0123】
本実施形態の第1工程S53では、第1準備工程S21で得られた高分子鎖4等の各鎖長(モノマー数)NFAでの第2物性値(以下、「第2物性値BFA」ということがある)に、粗視化分子モデル6の第2物性値BCGと、粗視化分子モデル6の鎖長(粒子数)NCG及び相互作用パラメータKとの関係を示す関数BCG(NCG,K)をフィッティングさせる。関数BCG(NCG,K)については、適宜定義することができる。本実施形態の関数BCG(NCG,K)は、上記の式(3)に基づいて定義することができる(即ち、BCG(NCG,K)=fNk(NCG,K))。上記のフィッティングは、相互作用パラメータKを固定した状態で、上記の式(8)で定義される長さの単位換算定数C1のみをフィッティングパラメータとして調節することで行われる。フィッティングには、例えば、Levenberg-Marquardt法を用いることができる。
【0124】
上述したように、本実施形態の第2物性値は、Kuhnセグメント数Nkである。このKuhnセグメント数Nkは無次元量であり、長さの単位換算定数Clに依存しない。このため、高分子鎖4等と粗視化分子モデル6との間で、長さの単位換算を行う必要がなく、高分子鎖4等のKuhnセグメント数Nkと、粗視化分子モデル6のKuhnセグメント数Nkとを直接比較することができる。
【0125】
図16は、粗視化分子モデルのKuhnセグメント数N
kと粒子数(鎖長)N
CGとの関係、及び、全原子モデルのKuhnセグメント数N
kとモノマー数(鎖長)N
FAとの関係を示すグラフである。第2物性値がKuhnセグメント数N
kである場合、第1工程S53での長さの単位換算定数C
lは、相互作用パラメータKを考慮する点を除き、上記特許文献1に記載の長さの単位換算定数の算出手順に基づいて求めることができる。なお、第2物性値が無次元量でない物性値である場合には、高分子鎖4等と粗視化分子モデル6との間で、長さの単位換算を行う必要がある。単位換算の方法については、上述の第1物性値A
FA、A
CGと同様の方法が用いられる。
【0126】
本実施形態では、相互作用が上記の式(1)で定義され、粗視化分子モデル6の実効的な密度lk/pを上記式(2)に近似させており、粗視化分子モデルのKuhnセグメント数Nkを上記式(3)に近似させている。
【0127】
上記式(2)の近似関数では、粗視化分子モデル6の鎖長(粒子数)NCG、及び、相互作用パラメータKの広い範囲で、粗視化分子モデル6の物性値を高精度に近似させることができ、かつ、パラメータ数が少なく簡素になるように設計されている。このため、上記の式(2)の近似関数では、粗視化分子モデル6の鎖長NCG、及び、相互作用パラメータKの広い範囲において、実効的な密度lk/pを高精度かつ短い計算時間で近似させることができる。
【0128】
一方、上記の式(3)は、上記の式(2)、式(6)及び式(7)を全て満たすように算出されているため、上記の式(2)との整合性が良好である。また、上記の式(2)及び式(3)で近似させている物性値(本例では、実効的な密度lk/p及びKuhnセグメント数Nk)は、共に無次元量であるため、高分子鎖4等と粗視化分子モデル6との物性値のフィッティングを行う際の単位換算が不要である。このため、本実施形態では、無次元量ではない物性値を用いる場合に比べ、計算時間を短縮することができる。したがって、パラメータ決定工程S23において、長さの単位換算定数Cl、及び、相互作用パラメータKを、高い精度かつ効率的に推定できる。
【0129】
このように、本実施形態の方法では、上記の式(2)及び(3)を共に用いることで、長さの単位換算定数Clおよび相互作用パラメータKを高い精度で、かつ効率的に推定できる。求められた長さの単位換算定数Clは、他の工程で求められる長さの単位換算定数Clとは独立して、コンピュータ1に記憶される。
【0130】
次に、本実施形態のパラメータ決定工程S23では、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等と、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間において、物性値の他方(本例では、第1物性値)がそれぞれ近づくように、長さの単位換算定数Clを固定して、相互作用パラメータKが求められる(第2工程S54)。本実施形態の第2工程S54では、工程S51で得られた第1物性値の近似関数(本例では、上記の式(2))を用いて、第1工程S53で求められた長さの単位換算定数Clで固定して、第1物性値(実効的な密度lk/p)が近づくように(本例では、一致するように)、相互作用パラメータKが求められる。
【0131】
本実施形態の第2工程S54では、第1準備工程S21で得られた高分子鎖4等の各鎖長(モノマー数)N
FAでの第1物性値A
FAに、粗視化分子モデルの第1物性値A
CGと、粗視化分子モデル6の鎖長(粒子数)N
CG及び相互作用パラメータKとの近似関数(本例では、上述の関数A
CG(N
CG,K))をフィッティングさせる。フィッティングは、長さの単位換算定数C
lを固定した状態で、相互作用パラメータKをフィッティングパラメータとして調節することで行われる。このようなフィッティングには、例えば、Levenberg-Marquardt法を用いることができる。
図17は、粗視化分子モデルの実効的な密度l
k/pと粒子数N
CGとの関係、及び、全原子モデルの実効的な密度l
k/pとモノマー数N
FAとの関係を示すグラフである。
【0132】
上述したように、本実施形態の第1物性値は、実効的な密度l
k/pである。上記のような相互作用パラメータKを求めることは、
図17に示されるように、粗視化分子モデル6実効的な密度l
k/pと粒子数(鎖長)N
CGとの関係を示し、かつ、相互作用パラメータKが異なる曲線33のうち、
図11で示される高分子鎖4等の実効的な密度l
k/p(
図17の点23)に最も近似する曲線33を選択することに対応している。求められた相互作用パラメータKは、他の工程で求められる相互作用パラメータKとは独立して、コンピュータ1に記憶される。
【0133】
次に、本実施形態の相互作用定義工程S2では、長さの単位換算定数Clと相互作用パラメータKとがそれぞれ収束したか否かが判断される(工程S55)。収束したか否かの判断については、適宜行うことができる。本実施形態において、複数回実施された第1工程S53及び第2工程S54の前後において、長さの単位換算定数Cl及び相互作用パラメータKがそれぞれ数値誤差の範囲内(例えば、0.01%以内)に収まっている場合に、収束したと判断することができる。
【0134】
工程S55において、長さの単位換算定数Clと相互作用パラメータKとがそれぞれ収束したと判断された場合(工程S55において、「Y」)、次の工程S56が実施される。この工程S56では、収束した長さの単位換算定数Clが、長さの単位換算定数Clとして決定される。さらに、工程S56では、収束した相互作用パラメータKが、相互作用パラメータKとして決定される。このような長さの単位換算定数Cl及び相互作用パラメータKは、第1工程S53及び第2工程S54において、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等と鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間で、第1物性値(実効的な密度lk/p)及び第2物性値(Kuhnセグメント数)が近づくように、長さの単位換算定数Clと相互作用パラメータKとの双方を収束させて求められたものである。したがって、本実施形態の相互作用定義工程S2では、粗視化分子モデル6の単位と、高分子鎖4の単位とを精度よく対応づけることが可能な長さの単位換算定数Cl及び相互作用パラメータKを求めることができる。
【0135】
一方、工程S55において、長さの単位換算定数Clと相互作用パラメータKとがそれぞれ収束していないと判断された場合(工程S55で、「N」)、第1工程S53及び第2工程S54が再度実施される。なお、再度実施される第1工程S53では、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等と、鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間において、物性値の一方(本例では、第2物性値)がそれぞれ近づくように、先の第2工程S54で求められた相互作用パラメータKで固定して、長さの単位換算定数Clが再度求められる。これにより、パラメータ決定工程S23では、長さの単位換算定数Clと相互作用パラメータKとが収束するまで、第1工程S53と第2工程S54とが交互に実施されるため、長さの単位換算定数Clと相互作用パラメータKとを確実に収束させることができる。
【0136】
このように、本実施形態の相互作用定義工程S2では、鎖長(モノマー数)NFAが異なる複数の高分子鎖4等と鎖長(粒子数)NCGが異なる複数の粗視化分子モデル6との間で、第1物性値及び第2物性値が近づくように、互いに依存する関係にある長さの単位換算定数Clと相互作用パラメータKとをそれぞれ収束させて求めることができる。このため、相互作用定義工程S2では、粗視化分子モデル6の単位と、高分子鎖4の単位とを精度よく対応づけることが可能な長さの単位換算定数Cl及び相互作用パラメータKを確実に求めることができる。
【0137】
相互作用定義工程S2では、高分子鎖4等の密度の鎖長依存性を考慮した補正が行われてもよい。これは、論文1に基づく粗視化分子モデル6の密度が一定であるのに対して、高分子鎖4等の密度は鎖長(モノマー数)N
FAに依存するため、このような鎖長依存性が考慮されないと、長さの単位換算が行われた際に、密度に乖離が生じるおそれがあるためである。このような依存性を考慮した補正を行う場合には、例えば、高分子鎖4等の密度の鎖長依存性を、下記の式(4)で定義することができる(
図18(全原子モデルの密度と、モノマー数N
FAとの関係を示すグラフ)に示す)。
【0138】
【数4】
ここで、
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
ρ
0、ρ
1:ρ
FAのフィッティングパラメータ
【0139】
なお、上記の式(4)を用いる場合、高分子鎖4等のモノマー数(鎖長)NFAに対応する粗視化分子モデルの粒子数(鎖長)NCGを求めるには、上述の式(8)を用いて表されるNCG=NFA・lFA/(Cl・lCG)に代えて、下記の式(5)が用いられるのが望ましい。
【0140】
【数5】
ここで、
N
CG:モノマー数がN
FAの分子鎖に対応する粗視化分子モデルの分子鎖1本の粒子数
N
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの鎖長
l
FA:モノマー1つあたりの長さ
C
l:長さの単位換算定数
l
CG:粗視化分子モデルの粒子間の結合の平衡長
ρ
FA:高分子鎖、全原子モデル又はユナイテッドアトムモデルのいずれかの密度
ρ
0:フィッティングパラメータ
【0141】
上述の式(8)に代えて、上記の式(5)が用いられることにより、この実施形態では、高分子鎖4等の密度の鎖長依存性を補正しながら、高分子鎖4等と粗視化分子モデル6の間で鎖長を対応付けて、第1物性値(実効的な密度l
k/p)及び第2物性値(Kuhnセグメント数N
k)を比較することができる。なお、上記の式(5)は、パラメータ決定工程S23の第1工程S53及び第2工程S54において、高分子鎖4等の鎖長(モノマー数)N
FAと対応する粗視化分子モデル6の鎖長(粒子数)N
CGを算出する際に用いられてよい。なお、高分子鎖4(
図2に示す)のモノマー数(鎖長)N
FAと、粗視化分子モデル6の粒子数(鎖長)N
CGとの比C
m(1個の粒子7あたりのモノマー数N
FA)については、上記の式(5)で定義される長さの単位換算定数C
lを用いて、C
m=N
FA/N
CG(N
FA)=C
l・l
CG/(l
FA・((ρ
0+ρ
1/N
FA)/ρ
0)
1/3)で定義することができる。
【0142】
次に、
図3に示されるように、本実施形態のシミュレーション方法では、コンピュータ1が、相互作用パラメータKを含む相互作用(角度ポテンシャル)E(θ)が定義された粗視化分子モデル6(
図4及び
図13に示す)を用いて、分子動力学に基づく構造緩和を計算する(工程S3)。本実施形態の工程S3では、先ず、相互作用定義工程S2で求められた相互作用パラメータKに基づいて設定される相互作用(角度ポテンシャル)E(θ)が、粗視化分子モデル6に定義される。次に、本実施形態の工程S3では、上記の相互作用が定義された複数の粗視化分子モデル6を、空間16に配置した高分子材料モデル20(
図13に示す)が設定される。粗視化分子モデル6の本数については、適宜設定することができる。
【0143】
次に、本実施形態の工程S3では、
図13に示した高分子材料モデル20において、隣接する粗視化分子モデル6、6の粒子7、7間に、相互作用ポテンシャルP2がそれぞれ定義される。相互作用ポテンシャルP2は、上述の手順に基づいて定義される。そして、本実施形態の工程S3では、上記の相互作用(角度ポテンシャル)E(θ)が定義された粗視化分子モデル6を対象に、分子動力学に基づく構造緩和が計算される。
【0144】
本実施形態のシミュレーション方法では、構造緩和の計算時において、粗視化分子モデル6の粒子7、7間に定義された上記の相互作用(角度ポテンシャル)E(θ)により、粗視化分子モデル6の実効的な密度lk/pを、高分子鎖4の実効的な密度lk/pに近づけることができる。また、本実施形態のシミュレーション方法では、構造緩和計算後の高分子材料モデル20の座標と、上述の長さの単位換算定数Clとに基づいて、粗視化分子モデル6の粒子7を、高分子鎖4のモノマーに置き換えることができる。このような置き換えによって得られる全原子モデル11が配置された高分子材料モデル18の密度は、実際の高分子鎖4の密度に近似する。
【0145】
本実施形態のシミュレーション方法では、粗視化分子モデル6が配置された高分子材料モデル20を用いて、実際の高分子鎖4に近似する構造緩和が可能となる。このため、本実施形態のシミュレーション方法では、構造緩和計算後の高分子材料モデル18の密度を、実際の高分子鎖4の密度に近似させることができるため、計算精度を向上させることができる。
【0146】
本実施形態の粗視化分子モデル6は、上述の長さの単位換算定数Clに基づいて、高分子鎖4のモノマーCm個が1個の粒子7に置換されているため、粗視化分子モデル6の長さの単位を、高分子鎖4の長さの単位に精度よく換算することができる。したがって、本実施形態のシミュレーション方法は、計算精度を向上させることができる。
【0147】
上記の相互作用(角度ポテンシャル)E(θ)が定義された粗視化分子モデル6は、例えば、特許文献1の時間換算工程に用いられてもよい。これにより、上述の長さの単位換算定数Clに基づいて、粗視化分子モデル6の粒子7が高分子鎖4のモノマーに置き換えられて、高分子材料モデル18が作成されることにより、時間単位の換算に必要な長鎖長の全原子モデルの構造緩和した座標を高い精度で作成できる。したがって、粗視化分子モデル6の時間単位を、高分子鎖4の時間単位により精度よく換算することができる。
【0148】
次に、本実施形態のシミュレーション方法では、
図3に示されるように、コンピュータ1が、高分子材料モデル20(
図13に示す)を用いたシミュレーション結果が良好である(即ち、所望の性能を有する)か否かを判断する(工程S4)。工程S4では、例えば、工程S3の構造緩和後の高分子材料モデル20を用いた変形計算後に、粗視化分子モデル6(
図3に示す)の単位を、高分子鎖4(
図2に示す)の単位に換算した結果に基づいて、高分子鎖4の物性等が評価される。
【0149】
工程S4において、シミュレーション結果が良好であると判断された場合(工程S4で、「Y」)、シミュレーションされた高分子材料モデル20に基づいて、高分子材料が製造される(工程S5)。一方、工程S4において、シミュレーション結果が良好でないと判断された場合(工程S4で、「N」)、高分子鎖4の構造や条件等を変更して(工程S6)、工程S1~工程S4が再度実施される。これにより、本実施形態のシミュレーション方法では、実際に高分子材料を試作しなくても、所望の性能を有する高分子材料を選択することができる。
【0150】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【0151】
図3に示した処理手順に従って、
図2に示した高分子鎖(cis-1.4ポリブタジエン)をモデル化した粗視化分子モデルがコンピュータに入力された。そして、粗視化分子モデルに、上記式(1)で定義される相互作用が定義された(実施例1、実施例2及び実施例3)。
【0152】
実施例1~3の相互作用を定義する工程では、
図7、
図8、
図12及び
図14に示した処理手順に基づいて、鎖長が異なる複数の全原子モデルと、鎖長が異なる複数の粗視化分子モデルとの間において、第1物性値(実効的な密度l
k/p)及び第2物性値(Kuhnセグメント数N
k)が近づくように、高分子鎖と粗視化分子モデルとの間の長さの単位換算定数C
l、及び、相互作用パラメータKが求められた。
【0153】
なお、実施例1~3において、粗視化分子モデルの実効的な密度lk/p及び粗視化分子モデルのKuhnセグメント数Nkは、上記の式(2)及び上記の式(3)の関数にそれぞれ近似させた。実施例1~2では、長さの単位換算定数Cl及び相互作用パラメータKを1つずつ交互に非線形最小二乗法で更新しながら、2つのパラメータの双方が収束するまで繰り返し実施された。一方、実施例3では、差分進化法に基づいて、2つのパラメータが探索された。
【0154】
また、実施例2では、高分子鎖4等の密度の鎖長依存性を、上記の式(4)で定義し、また、モノマー数に対応するビーズ数を、下記式(5)で定義した。全原子モデル、及び、粗視化モデルについては、次のとおりである。
全原子モデル:
第1全原子モデル:
モノマー数:20、空間内の本数:50本
第2全原子モデル:
モノマー数:30、空間内の本数:75本
第3全原子モデル:
モノマー数:40、空間内の本数:75本
粗視化分子モデル:
第1粗視化分子モデル~第5粗視化分子モデル:
粒子数:10、空間内の本数:10000本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
第6粗視化分子モデル~第10粗視化分子モデル:
粒子数:20、空間内の本数:5000本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
第11粗視化分子モデル~第15粗視化分子モデル:
粒子数:50、空間内の本数:2000本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
第16粗視化分子モデル~第20粗視化分子モデル:
粒子数:100、空間内の本数:1000本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
第21粗視化分子モデル~第25粗視化分子モデル:
粒子数:150、空間内の本数:677本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
第26粗視化分子モデル~第30粗視化分子モデル:
粒子数:200、空間内の本数:500本
相互作用パラメータ:0、0.5、1.0、1.5及び2.0で互いに異なる。
【0155】
実施例1~3では、250モノマーに相当する鎖長を有し、かつ、相互作用が定義された粗視化分子モデルを、空間内に100本配置された高分子材料モデルが定義された。そして、実施例1~3では、空間内に配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。
【0156】
比較のために、
図2に示した高分子鎖(cis-1.4ポリブタジエン)をモデル化した粗視化分子モデルがコンピュータに入力された(比較例)。比較例の粗視化分子モデルには、上記式(1)で定義される相互作用が定義されたが、
図14に示したパラメータ決定工程S23において、工程S51から第1工程S53までの工程が実施された後、第2工程S54以下の処理が省略された。即ち、比較例では、鎖長が異なる複数の全原子モデルと、鎖長が異なる複数の粗視化分子モデルとの間において、長鎖長極限のみを考慮して求められた長さの単位換算定数C
lと相互作用パラメータKが用いられており、実施例1~3のように、第1物性値(実効的な密度l
k/p)及び第2物性値(Kuhnセグメント数N
k)が近づくように、高分子鎖と粗視化分子モデルとの間の長さの単位換算定数C
l、及び、相互作用パラメータKが求められていない。そして、比較例の粗視化分子モデルが空間内に100本配置された高分子材料モデルが定義された。そして、比較例では、空間内に配置された粗視化分子モデルを対象に構造緩和が計算され、構造緩和計算後の高分子材料モデルの密度が計算された。
【0157】
実施例1~3で計算された密度、及び、比較例で計算された密度が、
図2に示した高分子鎖(cis-1.4ポリブタジエン)をモデル化した全原子モデルの平衡密度と比較された。全原子モデルは、構造緩和された粗視化分子モデルから、高分子鎖4のモノマーC
m個が1個の粒子7に置換されることで作成された。この平衡密度は、下記の仕様に基づいて、空間内に配置された全原子モデルを対象に構造緩和が計算されることで計算された。
全原子モデルの鎖長:250モノマー
空間内の本数:100本
温度:360K
力場:L-OPLS
【0158】
テストの結果、比較例の粗視化分子モデルの密度は、全原子モデルの平衡密度から10%乖離した。一方、実施例1~3の粗視化分子モデルの密度は、全原子モデルの平衡密度からの乖離を1%未満に抑えることができた。したがって、実施例1~3は、比較例に比べて、粗視化分子モデルの単位と、高分子鎖の単位とを精度良く関連付けることができた。
【0159】
また、実施例2では、上記の式(4)及び(5)に基づいて、高分子鎖の密度の鎖長依存性を考慮した補正が行われたため、実施例1及び実施例3に比べて、全原子モデルの平衡密度からの乖離を小さくすることができた。
【符号の説明】
【0160】
S1 粗視化分子モデルを入力する工程
S2 相互作用を定義する工程
S3 構造緩和を計算する工程