以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の粗視化分子モデルの作成方法(以下、単に「作成方法」ということがある。)は、コンピュータを用いて、高分子材料の分子鎖を複数の粗視化粒子で表現した粗視化分子モデルが作成される。
図1は、粗視化分子モデルの作成方法を実施するためのコンピュータの斜視図である。コンピュータ1は、本体1a、キーボード1b、マウス1c及びディスプレイ装置1dを含んでいる。この本体1aには、演算処理装置(CPU)、ROM、作業用メモリ、磁気ディスクなどの記憶装置及びディスクドライブ装置1a1、1a2などが設けられている。なお、記憶装置には、本実施形態の作成方法を実行するための処理手順(プログラム)が予め記憶される。
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態では、高分子材料として、cis-1,4ポリブタジエン(以下、単に「ポリブタジエン」ということがある。)が例示される。図2は、ポリブタジエンの構造式である。
ポリブタジエンを構成する高分子鎖3は、メチレン基(−CH2−)とメチン基(−CH−)とからなるモノマー5{−[CH2−CH=CH−CH2]−}が、重合度で連結されて構成されている。また、高分子材料の末端には、メチレン基(−CH2)に替えて、メチル基(−CH3)が連結される。なお、高分子材料には、ポリブタジエン以外の高分子材料が用いられてもよい。
本実施形態の作成方法では、高分子鎖3を表現した粗視化分子モデルが用いられる。図3は、粗視化分子モデル6が配置された仮想空間8の一例を示す概念図である。図4は、粗視化分子モデルの一例を示す概念図である。
図4に示されるように、粗視化分子モデル6は、図2に示した高分子鎖3を、高分子鎖3を構成する原子の数よりも少ない複数の粗視化粒子9を用いて表現されている。粗視化粒子9は、例えば、図2に示した高分子鎖3のモノマー5に対応して設定されている。本実施形態の粗視化分子モデル6は、Kremer-Grestモデルである場合が例示されるが、特に限定されるわけではなく、例えば、DPD(散逸粒子動力学法)に基づくモデル等であってもよい。
図5は、粗視化分子モデル6の部分拡大図である。粗視化分子モデル6の粗視化粒子9、9間には、平衡長が設定された第4ポテンシャルP4が定義される。本実施形態の第4ポテンシャルP4は、粗視化分子モデル6がKremer-Grestモデルである場合、下記式(2)で定義されるポテンシャル(以下、「LJポテンシャルULJ(rij)」ということがある。)と、下記式(3)で定義される結合ポテンシャルUFENEとの和で定義される。
ここで、各定数及び変数は、Lennard-Jones及びFENEの各ポテンシャルのパラメータであり、次のとおりである。
r
ij:粗視化粒子間の距離
r
c:カットオフ距離
k:粗視化粒子間のばね定数
ε:粗視化粒子間に定義されるLJポテンシャルの強度
σ:粗視化分子モデルの長さの単位
R
0:伸びきり長
なお、距離r
ij、カットオフ距離r
c、及び、伸びきり長R
0は、各粗視化粒子9の中心9o、9o間の距離について定義される。
このような第4ポテンシャルP4は、粗視化粒子9、9間を、平衡長を設定して連結するボンド10として定義される。これにより、複数の粗視化粒子9がボンド10によって伸縮自在に拘束された直鎖状の粗視化分子モデル6を設定することができる。
LJポテンシャルULJ(rij)及び結合ポテンシャルUFENEの各定数及び各変数の値としては、適宜設定することができる。本実施形態では、論文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)に基づいて適宜設定される。
図3に示した仮想空間8は、解析対象の高分子材料の微小構造部分に相当する。本実施形態の仮想空間8は、互いに向き合う三対の平面11、11を有する立方体として定義されている。各平面11には、周期境界条件が定義されている。従って、一方の平面11と、反対側の平面11とが連続している(繋がっている)ものとして取り扱うことができる。
仮想空間8の一辺の長さLaは、適宜設定することができる。本実施形態の長さLaは、粗視化分子モデル6の拡がりを示す量である慣性半径(図示省略)の3倍以上が望ましい。これにより、仮想空間8は、後述の粗視化分子動力学計算において、周期境界条件のために周期的に並ぶ自分自身(イメージ)との衝突の発生を防いで、粗視化分子モデル6の空間的拡がりを適切に計算することができる。また、粗視化粒子9の密度は、解析対象の高分子材料に対応するように適宜設定される。例えばKremer-Grestモデルにおいて、粗視化粒子9の密度は、0.85個/σ3である。粗視化分子モデル6の本数、及び、一辺の長さLaについては、一辺の長さLaが慣性半径より十分大きく(3倍以上)なるように、かつ、粗視化粒子9の密度が適切となるように調整されるのが望ましい。これにより、仮想空間8は、解析対象の高分子材料の少なくとも一部の空間領域を定義することができる。このような仮想空間8に、粗視化分子モデル6が配置されることにより、高分子材料をモデル化した高分子材料モデル12が設定される。仮想空間8は、コンピュータ1に記憶される。
図5に示されるように、本実施形態の作成方法では、粗視化分子モデル6の一端を構成する1番目の粗視化粒子9aから粗視化分子モデル6の他端を構成するm番目(m=5以上の整数)の粗視化粒子9mまでの各位置を、コンピュータ1の仮想空間8(図3に示す)内に順番に決定しながら配置している。粗視化粒子9の個数mについては、例えば、高分子材料の構造に応じて適宜設定される。本実施形態では、mを大きな整数に設定することで、大きな鎖長を有する粗視化分子モデル6を作成することができる。
ところで、従来の作成方法では、粗視化分子モデル6を作成する際に、熱平衡状態の慣性半径を定量的に考慮せずに、粗視化粒子9を、例えばランダムに配置していた。従って、熱平衡状態の慣性半径を再現した粗視化分子モデル6を作成するためには、粗視化粒子9を配置した後に、分子動力学計算に基づく緩和計算を、長時間行う必要があった。
発明者らは、鋭意研究を重ねた結果、i番目(iは、正の整数)の粗視化粒子9iについて、主鎖上の近傍(その粗視化粒子9iの直前)のk個の粗視化粒子9(k=3以上の整数)の粗視化粒子(i−k番目からi−1番目までの全ての粗視化粒子)との相互作用のみを考慮して、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置を計算して配置することで、熱平衡状態の慣性半径を再現しうる粗視化分子モデル6を短時間で作成しうることを知見した。
本実施形態の作成方法では、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置を計算するために、上述した第4ポテンシャルP4が定義される。さらに、必要に応じて、i番目の粗視化粒子9iに対して、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、及び、第3ポテンシャルEi(θ)が定義される。なお、第1ポテンシャルUangle(θ)、及び、第2ポテンシャルUtorsion(φ)は、例えば、通常のKremer-Grestモデルを使用する場合に省略することができる。なお、第3ポテンシャルEi(θ)の定義が省略される場合については、後述する。図6(a)、(b)は、ポテンシャルが定義される粗視化粒子モデルを示す図である。
図6(a)に示されるように、第1ポテンシャルUangle(θ)は、i番目の粗視化粒子9i、i−1番目の粗視化粒子9i-1、及び、i−2番目の粗視化粒子9i-2の間の結合角θに関するものである。図6(b)に示されるように、第2ポテンシャルUtorsion(φ)は、i番目の粗視化粒子9i、i−1番目の粗視化粒子9i-1、i−2番目の粗視化粒子9i-2、及び、i−3番目の粗視化粒子9i-3の間の二面角φに関するものである。第1ポテンシャルUangle(θ)及び第2ポテンシャルUtorsion(φ)は、シミュレーションに用いられる粗視化粒子9、及び、粗視化分子モデル6に応じて適宜設定される。
上記の第1ポテンシャルUangle(θ)及び第2ポテンシャルUtorsion(φ)は、図2に示した高分子鎖3の立体配座に基づくものである。このような第1ポテンシャルUangle(θ)及び第2ポテンシャルUtorsion(φ)がi番目の粗視化粒子9iに定義されることにより、高分子鎖3の立体配座を考慮した配置が可能となる。従って、第1ポテンシャルUangle(θ)及び第2ポテンシャルUtorsion(φ)は、熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置の計算に役立つ。
第3ポテンシャルEi(θ)は、下記式(1)に基づいて定義される。この第3ポテンシャルEi(θ)は、平均場近似ポテンシャルである。平均場近似ポテンシャルは、i番目の粗視化粒子と、i番目の粗視化粒子9iと主鎖上で近傍にあるk個の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子)のみとの相互作用を、i番目の粗視化粒子9iに作用する全ての相互作用から除いた相互作用に基づくものである。結合角θは、図6(a)に示されるように、i番目の粗視化粒子9i、i−1番目の粗視化粒子9i-1、及び、i−2番目の粗視化粒子9i-2の間の角度として定義される。この第3ポテンシャルEi(θ)は、i番目の粗視化粒子9iと主鎖上で近傍の相互作用を除いた周囲の相互作用(溶媒和など)の全体が、結合角θへ与える影響を表しており、高分子鎖の置かれている環境(例えば、密度、温度など)に依存して変化し得る。第3ポテンシャルEi(θ)の寄与が常に十分小さいと考えられる場合(例えば、粗視化粒子間の非結合相互作用が斥力(または、ほぼ斥力)で、密度が希薄であるような系)については、第3ポテンシャルEi(θ)の定義を省略(Ei(θ)=0)してもよい。
ここで、
E
i(θ):第3ポテンシャル
P
i(θ):熱平衡状態の結合角の頻度分布
P
i,local(θ):i番目の粗視化粒子と、i−k番目からi−1番目までの全ての粗視化粒子のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布
θ:結合角
k
B:ボルツマン定数
T:温度
上記式(1)の第3ポテンシャルEi(θ)は、ボルツマン因子exp(−Ei(θ)/kBT)が、2つの結合角の頻度分布の比Pi(θ)/Pi,local(θ)と等しくなるように定義される。上記式(1)において、温度Tは、仮想空間8に設定される温度である。
熱平衡状態の結合角の頻度分布(以下、単に「平衡頻度分布」ということがある。)Pi(θ)は、適宜設定された温度及び密度に基づいて、高分子材料の熱平衡状態へと緩和された粗視化分子モデル6について、ボンド10を介して連続する3つの粗視化粒子9がなす結合角θ(図6(a)に示す)と、その頻度との関係を示すものである。
一方、i番目の粗視化粒子9iと、i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1のみとの相互作用について、部分的に熱平衡状態の結合角の頻度分布(以下、単に「部分的平衡頻度分布」ということがある。)Pi,local(θ)は、i番目の粗視化粒子9iと、i番目の粗視化粒子9iの主鎖上の近傍の粗視化粒子9i-k〜9i-1との間の相互作用(第4ポテンシャルP4、並びに、必要に応じて定義される第1ポテンシャル、及び、第2ポテンシャルを含む)に基づいて熱平衡状態に緩和するように配置された粗視化分子モデルについて、ボンド10を介して連続する3つの粗視化粒子9がなす結合角θ(図6(a)に示す)と、その頻度との関係を示すものである。
平衡頻度分布Pi(θ)を部分的平衡頻度分布Pi,local(θ)で除した比Pi(θ)/Pi,local(θ)は、部分的に熱平衡状態にある結合角θの頻度分布を、熱平衡状態の結合角θの頻度分布にするための重みを示している。比Pi(θ)/Pi,local(θ)は、結合角θごとに計算され、結合角θの区間は、適宜設定される。この比Pi(θ)/Pi,local(θ)の自然対数が、ボルツマン定数kB及び温度Tに乗じられることで、第3ポテンシャルEi(θ)が定義される。このような第3ポテンシャルEi(θ)は、図6(a)に示したi番目の粗視化粒子9iが形成する結合角θを、部分的な熱平衡状態から熱平衡状態に調整することができる。
第3ポテンシャルEi(θ)は、粗視化分子モデル6の相互作用の定義(第4ポテンシャルP4、第1ポテンシャル、及び、第2ポテンシャルなど)、及び、作成すべき粗視化分子モデル6が置かれる仮想空間8の環境(温度、密度など)に依存する定数関数である。従って、本実施形態の作成方法では、粗視化粒子9を仮想空間8に配置する工程に先立ち、第3ポテンシャルEi(θ)が求められる。図7は、本実施形態の粗視化分子モデルの作成方法の処理手順の一例を示すフローチャートである。
本実施形態の作成方法では、第3ポテンシャルEi(θ)が求められる(初期工程S1)。上述したように、粗視化粒子間の非結合相互作用が斥力(または、ほぼ斥力)で、密度が希薄であるような系である場合には、初期工程S1が省略(第3ポテンシャルEi(θ)=0と定義)されてもよい。本実施形態の初期工程S1では、第3ポテンシャルEi(θ)を求めるために、未知の定数である平衡頻度分布Pi(θ)、及び、部分的平衡頻度分布Pi,local(θ)がそれぞれ取得される。図8は、初期工程S1の処理手順の一例を示すフローチャートである。
本実施形態の初期工程S1では、先ず、熱平衡状態の結合角の頻度分布Pi(θ)が求められる(第1頻度分布計算工程S11)。図9は、第1頻度分布計算工程S11の処理手順の一例を示すフローチャートである。
本実施形態の第1頻度分布計算工程S11では、先ず、図3に示されるように、粗視化分子モデル6が、仮想空間8に配置される(工程S61)。工程S61では、粗視化分子モデル6の一端を構成する1番目の粗視化粒子9aから、粗視化分子モデル6の他端を構成するm番目の粗視化粒子9mまでの各位置を、仮想空間8内に順番に決定しながら配置している。
なお、第1頻度分布計算工程S11で定義される粗視化分子モデル6及び仮想空間8は、熱平衡状態の結合角の頻度分布Pi(θ)を求めるためだけに用いられるものである。従って、粗視化分子モデル6及び仮想空間8は、本実施形態の作成方法で求められる粗視化分子モデル6及び仮想空間8とは異なるものであり、独立して定義される。
第1頻度分布計算工程S11で定義される粗視化粒子9の各位置は、他の粗視化粒子9との衝突が避けられることなく、モンテカルロ法に基づいてランダムに決定されてもよい。なお、粗視化粒子9の各位置は、他の粗視化粒子9との衝突を避けて配置されてもよい。また、粗視化粒子9の各位置をランダムに配置するのに代えて、例えば、図11に示す後述の配置工程S2の処理手順のうち、第3ポテンシャルを設定する工程S26を省略(即ち、Ei(θ)=0とみなす)した手順に基づいて、粗視化粒子9が配置されてもよい。
図5に示されるように、隣接する粗視化粒子9、9の間は、第4ポテンシャルP4に基づくボンド10で結合されている。各粗視化粒子9には、第3ポテンシャルEi(θ)は定義されない。本実施形態の工程S61では、熱平衡状態の結合角の頻度分布Pi(θ)を求めるために作成すべき粗視化分子モデル6の環境(温度、密度など)に応じて、系(仮想空間8)の体積、及び、粗視化分子モデル6の本数が適宜設定される。仮想空間8に配置された1本の粗視化分子モデル6は、コンピュータ1に記憶される。
次に、本実施形態の第1頻度分布計算工程S11では、図3に示した仮想空間8に配置された粗視化分子モデル6を対象に、分子動力学( Molecular Dynamics : MD )計算に基づいた構造緩和が計算される(工程S62)。本実施形態の分子動力学計算では、例えば、仮想空間8について所定の時間、粗視化分子モデル6が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での粗視化粒子9の動きが、単位時間毎に追跡される。
本実施形態の構造緩和の計算は、仮想空間8において、圧力及び温度が一定、又は、体積及び温度が一定に保たれる。工程S62では、粗視化分子モデル6の初期配置を十分に緩和されるまで、シミュレーションの単位時間毎に計算される。これにより、工程S62では、実際の高分子材料の分子運動に近似させて、粗視化分子モデル6の初期配置を精度よく緩和することができ、熱平衡状態の粗視化分子モデル6を計算することができる。
本実施形態の工程S62では、末端の影響が無視できる程度(典型的には絡み合い点間距離の2倍程度)の鎖長の粗視化分子モデル6の構造緩和が行われる。これにより、第1頻度分布計算工程S11では、実際の高分子材料に対応するような長い鎖長の粗視化分子モデル6を配置する必要がない。さらに、そのような長い鎖長の粗視化分子モデル6について、慣性半径の3倍程度の周期境界長とするために、非常に多数の粗視化分子モデル6を配置する必要もない。このため、本実施形態の第1頻度分布計算工程S11では、計算時間を大幅に短縮しうる。上記の構造緩和の計算は、例えば(株)JSOL社製のソフトマテリアル総合シミュレーター(J−OCTA)に含まれるCOGNACを用いて処理することができる。なお、工程S62では、分子動力学計算に代えて、メトロポリス法などのモンテカルロ法を用いて、熱平衡状態の粗視化分子モデル6が求められてもよい。熱平衡状態の粗視化分子モデル6は、コンピュータ1に記憶される。
次に、本実施形態の第1頻度分布計算工程S11では、熱平衡状態の結合角θの頻度分布が求められる(工程S63)。工程S63では、先ず、工程S62で求められた熱平衡状態の粗視化分子モデル6について、ボンド10を介して連続する3つの粗視化粒子9の全ての組み合わせ(ペア)の結合角θ(図5に示す)が計算される。そして、各結合角θについて、適宜分割された区間ごとの出現頻度が計算されることにより、熱平衡状態の結合角と、各区間の頻度との関係を示す頻度分布Pi(θ)が求められる。熱平衡状態の結合角の頻度分布Pi(θ)は、コンピュータ1に記憶される。
次に、本実施形態の初期工程S1では、部分的に熱平衡状態にある結合角の頻度分布が求められる(第2頻度分布計算工程S12)。図10は、第2頻度分布計算工程S12の処理手順の一例を示すフローチャートである。
本実施形態の第2頻度分布計算工程S12では、先ず、図3に示されるように、粗視化分子モデル6が、仮想空間8に配置される(工程S66)。工程S66では、第1頻度分布計算工程S11の工程S61と同様に、粗視化分子モデル6の一端を構成する1番目の粗視化粒子9aから粗視化分子モデル6の他端を構成するm番目の粗視化粒子9mまでの各位置を、仮想空間8内に順番に決定しながら配置している。
なお、第2頻度分布計算工程S12で定義される粗視化分子モデル6及び仮想空間8は、部分的に熱平衡状態にある結合角の頻度分布(部分的平衡頻度分布)Pi,local(θ)を求めるためだけに用いられるものである。従って、粗視化分子モデル6及び仮想空間8は、第1頻度分布計算工程S11で定義される粗視化分子モデル6及び仮想空間8や、本実施形態の作成方法で求められる粗視化分子モデル6及び仮想空間8とは異なるものであり、独立して定義される。
本実施形態の粗視化粒子9の各位置は、例えば、図11に示す後述の配置工程S2の処理手順のうち、第3ポテンシャルを設定する工程S26を省略(即ち、Ei(θ)=0とみなして)した手順に基づいて決定することができる。これにより、部分的に熱平衡状態の粗視化分子モデル6を作成することができる。図4に示されるように、粗視化粒子9、9間は、第4ポテンシャルP4に基づくボンド10で結合されている。
第2頻度分布計算工程S12で定義される粗視化粒子9には、計算に用いられる粗視化分子モデル6に基づき、必要に応じて、第1ポテンシャルUangle(θ)、及び、第2ポテンシャルUtorsion(φ)が定義される。第1ポテンシャルUangle(θ)、及び、第2ポテンシャルUtorsion(φ)の定義を省略できる場合については、上述のとおりである。一方、第3ポテンシャルEi(θ)は定義されてない(即ち、Ei(θ)=0として定義される)。本実施形態の工程S66では、1本の粗視化分子モデル6が、仮想空間8に配置されている。粗視化分子モデル6には、分子動力学計算に基づく構造緩和が実施されない。本実施形態では、i番目(iは、正の整数)以降の粗視化粒子9iについて、i番目の粗視化粒子9iの主鎖上の近傍のk個(k=3以上の整数)の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1)のみとの相互作用(第4ポテンシャルP4、及び、必要に応じて定義される第1ポテンシャル、第2ポテンシャルを含む)に基づいて、熱平衡状態に緩和するように配置される。これにより、粗視化分子モデル6を、部分的な相互作用に基づく熱平衡状態にあるものとして扱うことができる。
次に、本実施形態の第2頻度分布計算工程S12では、部分的に熱平衡状態にある結合角の頻度分布(部分的平衡頻度分布)Pi,local(θ)が求められる(工程S67)。本実施形態の工程S67では、先ず、部分的に熱平衡状態にある粗視化分子モデル6について、ボンド10を介して連続する3つの粗視化粒子9の全ての組み合わせ(ペア)の結合角θが計算される。そして、各結合角θについて、適宜分割された区間ごとの出現頻度が計算されることにより、部分的に熱平衡状態にある結合角の頻度分布Pi,local(θ)が求められる。部分的平衡頻度分布Pi,local(θ)は、コンピュータ1に記憶される。
次に、本実施形態の初期工程S1では、熱平衡状態の結合角の頻度分布Pi(θ)、及び、部分的に熱平衡状態にある結合角の頻度分布Pi,local(θ)に基づいて、第3ポテンシャルEi(θ)が求められる(工程S13)。工程S13では、第3ポテンシャルEi(θ)を定義する上記式(1)に、熱平衡状態の結合角の頻度分布Pi(θ)、及び、部分的に熱平衡状態にある結合角の頻度分布Pi,local(θ)が代入される。これにより、初期工程S1では、第3ポテンシャルEi(θ)を求めることができる。第3ポテンシャルEi(θ)は、コンピュータ1に記憶される。
次に、本実施形態の作成方法では、コンピュータ1が、粗視化分子モデル6を構成する粗視化粒子9の各位置を、仮想空間8内に配置する(配置工程S2)。本実施形態の配置工程S2は、粗視化分子モデル6の一端を構成する1番目の粗視化粒子9aから粗視化分子モデル6の他端を構成するm番目の粗視化粒子9mまでの各位置を、仮想空間8内に順番に決定しながら配置している。図11は、配置工程S2の処理手順の一例を示すフローチャートである。
本実施形態の配置工程S2は、先ず、図5に示した1番目の粗視化粒子9aが、仮想空間8内に配置される(工程S21)。本実施形態の工程S21では、1番目の粗視化粒子9a位置が、仮想空間8内に既に配置されている他の粗視化粒子9との衝突が避けられることなくランダムに配置される。
なお、粗視化粒子9以外の粒子(例えば、高分子材料に含まれるフィラーをモデル化したフィラー粒子等)が仮想空間8に配置されている場合、その粒子との衝突を避けた位置に、1番目の粗視化粒子9aの位置が決定されるのが望ましい。これにより、粗視化粒子9が他の粒子(図示省略)に拘束されるのを防ぐことができる。1番目の粗視化粒子9aの位置は、コンピュータ1に記憶される。
次に、本実施形態の配置工程S2では、工程S22以降において、i番目以降(2番目以降)の粗視化粒子9が、仮想空間8に配置される。本実施形態の配置工程S2の工程S22以降では、2番目の粗視化粒子9bからm番目の粗視化粒子9mまでの各位置を、仮想空間8内に順番に決定しながら配置している。図12(a)〜(c)は、4番目の粗視化粒子9dの配置の一例を説明する図である。
本実施形態の配置工程S2では、先ず、図12(a)に示されるように、i番目の粗視化粒子9i(例えば、図12では、4番目の粗視化粒子9d)の位置が仮決定される(工程S22)。工程S22では、先ず、粗視化粒子9に第4ポテンシャルP4が定義される。さらに、工程S22では、必要に応じて、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、及び、第3ポテンシャルEi(θ)が、i番目の粗視化粒子9iに対して、直前のk個(k=3以上の整数)の粗視化粒子(i−k番目からi−1番目までの全ての粗視化粒子)との間にのみ定義される。ただし、位置の仮決定の際には、計算の効率化のために、これらの相互作用をそれぞれ単純化した形で用いて良い。
第4ポテンシャルP4のLJポテンシャルULJ(rij)については、例えば、距離rmin未満で∞、それ以上で0となる関数に単純化することができる。これにより、相互作用を考慮する2つの粗視化粒子間の距離が常に距離rmin以上となるような粗視化粒子9の位置を選ぶことができる。第4ポテンシャルP4の結合ポテンシャルUFENEについては、例えば、平衡核間距離L1で0、それ以外で∞となる関数に単純化することができる。これにより、結合長が常に平衡核間距離L1となるような粗視化粒子9の位置を選ぶことができる。
第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、第3ポテンシャルEi(θ)については、例えば、常に0となる関数に単純化することができる。これにより、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、第3ポテンシャルEi(θ)を考慮せずに、粗視化粒子9の位置を仮決定することができる。
上記のような相互作用の単純化により、2番目以降のi番目の粗視化粒子9iの位置は、i−1番目の粗視化粒子の位置に対する方向をランダムに選び、距離L1進んだ位置として仮決定することができる。ただし、仮決定した位置が、主鎖上近傍のk個の粗視化粒子(i−k番目からi−1番目までの全ての粗視化粒子)9のいずれかと衝突する(距離がrminよりも小さい)ときは、i−1番目の粗視化粒子の位置に対する方向を選び直すのが望ましい。粗視化粒子9iのランダムな位置については、平衡核間距離L1を半径とする球面上の点について、一様乱数となるように決定されるのが望ましい。なお、粗視化分子モデル6がKremer-Grestモデルの場合、一様乱数の範囲は、平衡核間距離L1〜0.965σ、及び、距離rmin〜0.9σに設定されるのが望ましい。
工程S22の手順を具体的に説明すると、2番目の粗視化粒子9bの位置は、1番目の粗視化粒子9aとの結合距離を平衡核間距離L1に固定して、1番目の粗視化粒子9との衝突を避けたランダムな方向に仮決定される。3番目の粗視化粒子9cの位置は、1番目の粗視化粒子9a及び2番目の粗視化粒子9bとの衝突を避けたランダムな方向に仮決定される。
4番目以降の粗視化粒子9dについては、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の位置が、i−2番目の粗視化粒子9i-2(例えば、2番目の粗視化粒子9b)と、i−3番目の粗視化粒子9i-3(例えば、1番目の粗視化粒子9a)との相互作用のみを考慮して、それらとの衝突を避けたランダムな方向に仮決定される。
i番目の粗視化粒子9iの位置が仮決定された後、i番目の粗視化粒子9iとi−1番目の粗視化粒子9i-1との間は、ボンド10で結合される。仮決定されたi番目の粗視化粒子9iの位置は、コンピュータ1に記憶される。
このように、本実施形態の作成方法では、それぞれの粗視化粒子9について、その直前のk個の粗視化粒子9との相互作用のみを考慮して、それらとの衝突を避けた位置に、粗視化粒子9の位置を仮決定しているため、例えば、全ての粗視化粒子9との衝突を避けた位置に仮決定される場合に比べて、粗視化粒子9の位置を短時間で仮決定することができる。
また、粗視化粒子9以外の粒子(例えば、フィラー粒子(図示省略))が仮想空間8に配置されている場合、その粒子との衝突を避けた位置に、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の位置が決定されるのが望ましい。なお、粗視化粒子9以外の粒子の配置によっては、例えば袋小路のような微小な密閉空間が存在すると、次の粗視化粒子9の位置が決定できなくなることがあり得る。そのような場合、配置工程S2が最初から実施されるのが望ましい。
次に、本実施形形態の配置工程S2では、仮決定されたi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)と、その直前のk個(k=3以上の整数)の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1)との相互作用を考慮して、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置が計算される(計算工程S23)。計算工程S23では、メトロポリス法に基づいて、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置が決定される。図13は、計算工程S23の処理手順の一例を示すフローチャートである。
本実施形態の計算工程S23では、先ず、仮決定されたi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)のポテンシャルエネルギーeが計算される(工程S91)。ポテンシャルエネルギーeは、i番目の粗視化粒子9iと主鎖と、i番目の粗視化粒子9iの主鎖上の近傍のk個(k=3以上の整数)の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1)との第4ポテンシャルP4が含まれる。なお、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、及び、第3ポテンシャルEi(θ)が定義されている場合には、第4ポテンシャルP4、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、及び、第3ポテンシャルEi(θ)を含めたポテンシャルの和である。
第1ポテンシャルUangle(θ)及び第3ポテンシャルEi(θ)は、3番目以降の粗視化粒子9i(例えば、4番目の粗視化粒子9d)、i−1番目の粗視化粒子9i-1(例えば、3番目の粗視化粒子9c)、及び、i−2番目の粗視化粒子9i-2(例えば、2番目の粗視化粒子9b)の間(組)に設定されたものである。従って、工程S91において、3番目以降の粗視化粒子9iが含まれない組は、第1ポテンシャルUangle(θ)及び第3ポテンシャルEi(θ)が考慮されない。同様に、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)が含まれない組の粗視化粒子9については、第4ポテンシャルP4も考慮されない。
第2ポテンシャルUtorsion(φ)は、4番目以降の粗視化粒子9i(例えば、4番目の粗視化粒子9d)、i−1番目の粗視化粒子9i-1(例えば、3番目の粗視化粒子9c)、i−2番目の粗視化粒子9i-2(例えば、2番目の粗視化粒子9b)、及び、i−3番目の粗視化粒子9i-3(例えば、1番目の粗視化粒子9a)の間(組)に設定されたものである。従って、工程S91において、4番目以降の粗視化粒子9iが含まれない組は、第2ポテンシャルUtorsion(φ)が考慮されない。
このように、工程S91では、ポテンシャルエネルギーeに、第3ポテンシャルEi(θ)を加算することができる。これにより、工程S91では、i番目の粗視化粒子9iと、i番目の粗視化粒子9iの主鎖上の近傍のk個(k=3以上の整数)の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1)との相互作用を、第3ポテンシャルEi(θ)によって補正して、結合角θの分布が熱平衡状態に近似するように調整されたポテンシャルエネルギーeを計算することができる。ポテンシャルエネルギーeは、コンピュータ1に記憶される。
次に、本実施形態の計算工程S23では、図12(b)に示されるように、仮決定されたi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)とは別に、i番目の粗視化粒子9iの移動候補位置20が設定される(工程S92)。工程S92では、仮決定する工程S22と同様の手順に基づいて、i番目の粗視化粒子9iの移動候補位置20が、i−2番目の粗視化粒子9i-2(例えば、2番目の粗視化粒子9b)、及び、9i-3(1番目の粗視化粒子9a)のみを考慮し、かつ、それらとの衝突を避けたランダムな位置に設定される。
なお、粗視化粒子9以外の粒子(例えば、フィラー粒子(図示省略))が仮想空間8に配置されている場合、その粒子との衝突を避けたランダムな位置に、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の移動候補位置20が決定されるのが望ましい。
次に、本実施形態の計算工程S23では、2番目以降の粗視化粒子9について、移動候補位置20に配置された場合のi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)のポテンシャルエネルギーe’が計算される(工程S93)。このポテンシャルエネルギーe’は、移動候補位置20に配置されたi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)と、i番目の粗視化粒子9iの主鎖上の近傍のk個(k=3以上の整数)の粗視化粒子9(i−k番目からi−1番目までの全ての粗視化粒子9i-k〜9i-1)の第4ポテンシャルP4が含まれる。なお、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)及び第3ポテンシャルEi(θ)が定義されていれば、ポテンシャルエネルギーe’は、第4ポテンシャルP4、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)及び第3ポテンシャルEi(θ)を含めたポテンシャルの和である。従って、工程S93では、結合角θの分布が熱平衡状態に近似するように調整されたポテンシャルエネルギーe’を計算することができる。ポテンシャルエネルギーe’は、コンピュータ1に記憶される。
次に、本実施形態の計算工程S23では、一様乱数Rを生成させる(工程S94)。工程S94では、一般的な疑似乱数生成アルゴリズム(pseudo random number generator)を用いて、0以上1未満の実数の一様乱数Rを生成させている。一様乱数Rは、コンピュータ1に記憶される。
次に、本実施形態の計算工程S23では、メトロポリス法に基づいて、仮決定されたi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)のポテンシャルエネルギーeと、移動候補位置20のポテンシャルエネルギーe’との差(変化)についての確率exp[(e−e’)/kBT]が、工程S94で生成された一様乱数Rの値よりも大きいか否かが判断される(工程S95)。確率exp[(e−e’)/kBT]において、kBはボルツマン定数であり、Tは仮想空間8の温度である。
工程S95において、確率exp[(e−e’)/kBT]が、一様乱数Rの値よりも大きい場合(工程S95において、「Y」)、次の工程S96において、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の移動候補位置20が、熱平衡状態へと緩和されるi番目の粗視化粒子9iの仮決定された位置(以下、「仮決定位置」ということがある。)として更新(採択)される。移動候補位置20のポテンシャルエネルギーe’は、i番目の粗視化粒子9iのポテンシャルエネルギーeとして更新(採択)される。更新(採択)された仮決定位置22、及び、更新されたポテンシャルエネルギーeは、コンピュータ1に記憶される。そして、次の工程S97が行われる。
一方、工程S95において、確率exp[(e−e’)/kBT]が、一様乱数Rの値以下である場合(工程S95において、「N」)、その移動候補位置20は棄却され、次の工程S97が実施される。
工程S97では、工程S92〜工程S96が、予め定められた回数繰り返されたか否かが判断される。工程S92〜工程S96が繰り返される回数については、適宜設定することができる。本実施形態では、計算コストと計算精度とのバランスを考慮して、5〜20回程度に設定される。
工程S97において、工程S92〜工程S96が、予め定められた回数繰り返されたと判断された場合(工程S97において、「Y」)、次の工程S98において、i番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の仮決定位置が、熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置として決定される。決定された位置22は、コンピュータ1に記憶される。
一方、工程S97において、工程S92〜工程S96が、予め定められた回数繰り返されていないと判断された場合(工程S97において、「N」)、工程S92〜工程S96が再度実施される。
このように、本実施形態の計算工程S23では、メトロポリス法の棄却採択法に基づいて、i番目の粗視化粒子9iが形成する結合角θの熱平衡状態のアンサンブルを得ることができる。しかも、計算工程S23では、予め定められた回数に達するまで、i番目の粗視化粒子9iの熱平衡状態へと緩和される位置を棄却・採択する試行(工程S92〜工程S96)が行われる。これにより、本実施形態の計算工程S23では、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9i(例えば、4番目の粗視化粒子9d)の位置22を、精度良く計算することができる。
本実施形態の計算工程S23では、i番目の粗視化粒子9iの熱平衡状態へと緩和される位置を棄却・採択する試行(工程S92〜工程S96)が、予め定められた回数繰り返されたが、このような態様に限定されない。棄却・採択する試行(工程S92〜工程S96)が、予め定められた上限値を超えた場合に、1番目の粗視化粒子9aの配置からやり直されてもよい。これにより、粗視化分子モデルの熱平衡状態への緩和が不十分になるのを防ぐことができる。
次に、本実施形態の配置工程S2では、図12(c)に示されるように、計算工程S23で計算された位置(決定された位置)22が、i番目の粗視化粒子9iの位置に設定される(工程S24)。工程S24では、仮決定されたi番目の粗視化粒子9iの座標値が、計算された位置22の座標値に変換される。変換された座標値は、コンピュータ1に記憶される。
次に、本実施形態の配置工程S2では、粗視化分子モデル6の他端を構成するm番目の粗視化粒子9m(図5に示す)まで配置されたか否かが判断される(工程S25)。工程S25において、m番目の粗視化粒子9mまで配置されたと判断された場合(工程S25において、「Y」)、配置工程S2の一連の処理が終了し、次の工程S3(図7に示す)が実施される。他方、工程S25において、m番目の粗視化粒子9mまで配置されていないと判断された場合(工程S25において、「N」)、次の粗視化粒子9(即ち、i←i+1)が選択され(工程S26)、工程S22〜工程S25が再度実施される。
これにより、配置工程S2では、図5に示されるように、粗視化分子モデル6の一端を構成する1番目の粗視化粒子9aから粗視化分子モデルの他端を構成するm番目の粗視化粒子9mまでの各位置を、仮想空間8内に順番に決定しながら配置することができる。
次に、本実施形態の作成方法では、全ての粗視化分子モデル6(図3に示す)が仮想空間8に配置されたか否かが判断される(工程S3)。仮想空間8に配置される粗視化分子モデル6の本数については、高分子材料モデル12(図3に示す)が用いられるシミュレーション等に応じて、適宜設定される。
工程S3において、全ての粗視化分子モデル6が仮想空間8に配置されたと判断された場合(工程S3において、「Y」)、次の工程S4が実施される。他方、工程S3において、全ての粗視化分子モデル6が仮想空間8に配置されていないと判断された場合(工程S3において、「N」)、配置工程S2が再度実施されて、新たな粗視化分子モデル6が仮想空間8に配置される。これにより、本実施形態の作成方法では、全ての粗視化分子モデル6を仮想空間8に配置した高分子材料モデル12を作成することができる。
このように、本実施形態の作成方法では、図12に示されるように、仮決定されたi番目の粗視化粒子9iと、i−k番目からi−1番目までの粗視化粒子9のみとの相互作用を考慮して、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子9iの位置が計算されるため、高分子鎖の熱平衡状態の慣性半径を決定づける主鎖の曲がりやすさを再現することができる。これにより、本実施形態の作成方法では、熱平衡状態の慣性半径を再現しうる粗視化分子モデル6を作成することができる。したがって、本実施形態の作成方法では、従来の方法のように、分子動力学計算に基づく粗視化分子モデル6の緩和計算を長時間行う必要がないため、熱平衡状態の慣性半径を再現しうる粗視化分子モデル6を、短時間で作成できる。
次に、本実施形態の作成方法では、配置工程S2の後に、コンピュータ1が、図3に示した粗視化粒子9、9間に粒子の重なりを許容するソフトポテンシャルU(r)を定義して、粗視化分子モデル6の構造緩和を計算する(工程S4)。工程S4では、仮想空間8に配置された全ての粗視化分子モデル6の粗視化粒子9、9間に、粒子の重なりを許容するソフトポテンシャルU(r)が定義される。ソフトポテンシャルU(r)は、上記論文1に基づいて、下記式(4)で定義される。
ここで、
r:粗視化粒子間の距離
A:粗視化粒子間に定義されるポテンシャルの強度
σ:粗視化分子モデルの長さの単位
なお、距離rは、各粗視化粒子9の中心9o、9o間の距離として定義される。
このようなソフトポテンシャルU(r)は、斥力場だけを有しており、粗視化粒子9の重なりを許容することができる。なお、第4ポテンシャルP4は無効に設定される。
そして、工程S4では、仮想空間8に配置された全ての粗視化分子モデル6を対象に、分子動力学計算に基づいた構造緩和が計算される。構造緩和の計算は、仮想空間8に配置された全ての粗視化粒子9の重なりが除去されるまで行われる。なお、各粗視化分子モデル6は、熱平衡状態の慣性半径が再現されているため、粗視化分子モデルの緩和計算を長時間行う必要がない。
緩和計算後には、ソフトポテンシャルU(r)が無効に設定されるとともに、第4ポテンシャルP4が有効に設定される。これにより、本実施形態の作成方法では、全ての粗視化粒子9の重なりが除去された高分子材料モデル12を作成することができる。高分子材料モデル12は、コンピュータ1に記憶される。
本実施形態の作成方法で作成された粗視化分子モデル6及び高分子材料モデル12は、例えば、一般的に行われている単軸引張り試験に基づいて、一方向に(例えば、X軸方向に0%〜20%)伸長させる変形シミュレーション等に用いることができる。従って、本実施形態の作成方法は、高分子材料の開発及び製造に役立つ。
本実施形態の配置工程S2では、計算工程S23に先立ち、第1ポテンシャルUangle(θ)、第2ポテンシャルUtorsion(φ)、及び、第3ポテンシャルEi(θ)を定義することができる。これらのポテンシャルの定義の有無や定義の内容については、i番目の粗視化粒子9iの位置の仮設定、及び、高分子材料の熱平衡状態へと緩和を計算できれば特に限定されない。例えば、粗視化分子モデル6がKremer-Grestモデルの場合、第1ポテンシャルUangle(θ)、及び、第2ポテンシャルUtorsion(φ)は省略される。さらに、気体状態のKGモデルなど、溶媒和の影響をほぼ無視できる場合には、第3ポテンシャルEi(θ)が0に近くなる。このため、第3ポテンシャルEi(θ)を省略することができる。
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
図7に示した処理手順に従って、粗視化分子モデルがコンピュータで作成された(実施例、比較例)。
実施例では、図8〜11に示した処理手順に従って、粗視化分子モデルの一端を構成する1番目の粗視化粒子から粗視化分子モデルの他端を構成するm番目の粗視化粒子までの各位置を、仮想空間内に順番に決定しながら配置する配置工程が実施された。実施例の配置工程では、図13に示した処理手順に従って、i番目(iは、正の整数)の粗視化粒子を、i番目の粗視化粒子と、その直前のk個(k=3以上の整数)の粗視化粒子(i−k番目からi−1番目までの全ての粗視化粒子)との相互作用のみを考慮して仮決定する工程と、考慮した相互作用を用いて、高分子材料の熱平衡状態へと緩和されるi番目の粗視化粒子の位置を計算する計算工程と、計算された位置をi番目の粗視化粒子の位置として決定する工程とが実施された。メトロポリス法に基づく工程S92〜S96の繰り返し回数は、各粗視化粒子について10回とした。
実施例では、第1ポテンシャル、第2ポテンシャル、第3ポテンシャル及び第4ポテンシャルが定義された。そして、実施例では、配置する工程の後に、粗視化粒子の重なりを許容するソフトポテンシャルを、粗視化粒子間に定義して、粗視化分子モデルの構造緩和が計算された。
比較例では、従来の方法に基づいて、粗視化分子モデルの一端を構成する1番目の粗視化粒子から粗視化分子モデルの他端を構成するm番目の粗視化粒子までの各位置が、モンテカルロ法に基づいて、仮想空間内にランダムに配置された。そして、比較例では、熱平衡状態の慣性半径を再現した粗視化分子モデルを作成するために、分子動力学計算に基づく粗視化分子モデルの緩和計算が実施された。
そして、実施例及び比較例の粗視化分子モデルについて、慣性半径の二乗平均が求められた。共通仕様は、次のとおりである。
粗視化分子モデル:
本数:1万本
1本の粗視化分子モデルに含まれる粗視化粒子の個数(鎖長)m:100個
テストの結果、実施例の慣性半径の二乗平均は37.5±0.6σであり、比較例の慣性半径の二乗平均は36.3±0.4σであった。これらの値は、実際の慣性半径の二乗平均に近似している。一方、実施例の作成時間は、比較例の作成時間の2万分の1に短縮できた。従って、実施例は、熱平衡状態の慣性半径を再現しうる粗視化分子モデルを、短時間で作成することができた。