【文献】
藤尾容正,タイヤ用ゴムの新材料創出にむけた階層的化学計算の応用 [online],2009年,[検索日 2018.01.04],インターネット,URL,http://www.nsc.riken.jp/sympo2009/poster-session/abstract/P-28Bito.pdf
【文献】
HARMANDARIS, V.A.; ADHIKARI, N.P.; VAN DER VEGT, N.F.A.; KREMER K. ,Hierarchical Modeling of Polystyrene: From Atomistic to Coase-Grained Simulations,Macromolecules [online],ACS,2006年 8月26日,Vol. 39, No.19,pp. 6708 - 6719,[検索日 2018.01.04],インターネット,URL,http://pubs.acs.org/doi/abs/10.1021/ma0606399
【文献】
物理学大事典,株式会社朝倉書店,2005年10月25日,初版,p. 573
(58)【調査した分野】(Int.Cl.,DB名)
【発明を実施するための形態】
【0015】
以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の全原子モデルの作成方法(以下、単に「作成方法」ということがある)は、コンピュータを用いて、高分子材料を構成する任意のホモポリマー高分子鎖に対して、モノマー単位での幾何異性体(シス・トランス異性体)の特徴を再現し、かつ、平衡状態に近似する全原子モデルを作成するための方法である。
【0016】
図1に示されるように、コンピュータ1は、本体1a、キーボード1b、マウス1c及びディスプレイ装置1dを含む。この本体1aには、演算処理装置(CPU)、ROM、作業用メモリー、磁気ディスクなどの記憶装置及びディスクドライブ装置1a1、1a2などが設けられる。なお、記憶装置には、本実施形態の作成方法を実行するための処理手順(プログラム)が予め記憶される。
【0017】
高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態の高分子材料は、
図2に示されるように、cis-1,4ポリブタジエン(以下、単に「ポリブタジエン」ということがある)である場合が例示される。このポリブタジエンは、ホモポリマー高分子鎖(単一重合体)から構成される。このホモポリマー高分子鎖は、メチレン基(−CH
2−)とメチン基(−CH−)とからなるモノマー{−[CH
2−CH=CH−CH
2]−}の繰り返し(重合度nで連結すること)によって、その大部分(本実施形態では、全体)が構成されている。このホモポリマー高分子鎖は、モノマー単位での幾何異性体に特徴を有している。ここで、「モノマー単位での幾何異性体に特徴を有する」とは、モノマー単位の幾何異性体(シス・トランス異性体)が存在する場合、幾何異性体の割合やホモポリマー高分子鎖中での配置等に何らかの規則性があることを意味している。本実施形態のホモポリマー高分子鎖は、モノマーが100%シス型の構造からなるという特徴を有している。なお、ホモポリマー高分子鎖は、モノマーが100%シス型ではない任意の幾何異性体を有するホモポリマー高分子鎖に対しても、本発明の作成方法は有効である。
【0018】
また、ホモポリマー高分子鎖は、後述するレプテーション理論、又は、レプテーション理論から派生した理論(以下、これらの理論を総括して「レプテーション理論」という。)が用いられることにより、物性を予測しうることが知られている。なお、これらの理論に基づくシミュレーションは、例えば、ソフトマテリアル統合シミュレーターOCTAのPASTA( Polymer rheology Analyzer with Slip-link model of enTAnglement )を用いて実施することができる。
【0019】
なお、高分子材料としては、ポリブタジエン以外の高分子材料が採用されてもよい。例えば、ホモポリマー高分子鎖が、100%シス等のモノマー単位での幾何異性体の特徴を有している天然ゴムや、モノマー単位での幾何異性体の特徴を人工的に付与した高分子材料が採用されてもよい。
【0020】
図3には、本実施形態の作成方法の処理手順の一例が示されている。本実施形態では、先ず、コンピュータ1に、ホモポリマー高分子鎖をモデル化した粗視化モデルが設定される(ステップS1)。なお、粗視化モデルとしては、全原子モデルの分子動力学計算よりも計算コストが小さいものであれば、任意の粗視化ポテンシャル(例えば、angle、torsionのポテンシャル形式)を含むものでもよい。
【0021】
ステップS1では、
図2及び
図4に示されるように、先ず、ホモポリマー高分子鎖のモノマー又はモノマーの一部分をなす構造単位2を、直径を持った球で表現されるビーズ(粒子)3で置換する。本実施形態では、ホモポリマー高分子鎖のチューブ半径R1に基づいて、1個のビーズ3に対応する構造単位2のモノマーの個数が決定される。このチューブ半径R1は、現実のホモポリマー高分子鎖の物性や運動を予測するのに用いられるものである。また、本実施形態のチューブ半径R1は、ホモポリマー高分子鎖に対して有効なレプテーション理論で定義される。
【0022】
レプテーション理論は、現実のホモポリマー高分子鎖の物性や運動を予測するのに用いられるものである。レプテーション理論では、一本の高分子鎖をチューブ半径の細長い“ひも”で粗視化して、物性挙動を解析することに成功している。この“ひも”は、内部で高分子鎖の局所的相関を維持し、他の高分子鎖を排除するチューブ半径程度のブロッブと呼ばれる“糸まり”の連結で近似することができる。本発明では、ブロッブの内部に全原子模型を適用することで、粗視化した高分子鎖のモデルから、分子的な視点で解析を行うことができる全原子模型の配置を効率良く求めている。
【0023】
レプテーション理論において、チューブ半径R1は、ホモポリマー高分子鎖がチューブ状にモデル化されるときの半径である。チューブ半径R1は、下記式(1)で定義される。
【0024】
【数1】
ここで、各変数は、次のとおりである。
<R
2>
0:ホモポリマー高分子鎖の末端間距離の自乗平均
M:ホモポリマー高分子鎖の分子量
ρ:ホモポリマー高分子鎖の密度(kg/m
3)
T:絶対温度(K)
G
e:平衡弾性率(Pa)
R:気体定数(m
2kgs
−2K
−1 mol
−1)
【0025】
上記式(1)において、ホモポリマー高分子鎖の密度ρ、絶対温度T、平衡弾性率G
e及び<R
2>
0/Mは、実験値である。特に、<R
2>
0/Mは、高分子材料(例えば、ポリブタジエン)に中性子ビームを照射する中性子小角散乱(SANS)を実施することによって求められる実験値である。なお、各変数は、例えば、論文( L,J.Fetters , D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings」、Physical Properties of Polymers Handbook Second Edition P445-P452 )に基づいて設定することもできる。
【0026】
例えば、ホモポリマー高分子鎖がポリブタジエンである場合、チューブ半径R1は、4.71nmである。このチューブ半径R1に基づいて、1.55個分のモノマーを構造単位2として、対応する一個のビーズ3で当該構造単位2を置換している。なお、1.55個分のモノマーを構造単位2としたのは、論文(L,J.Fetters ,D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings 」Physical Properties of Polymers Handbook Second Edition P448」)、及び、論文( 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を設定することができる。
【0027】
ビーズ3は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ビーズ3には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ1に記憶される。
【0028】
また、ビーズ3、3間には、平衡長が定義されたポテンシャルが設定される。これにより、ビーズ3、3間には、該ビーズ3、3を拘束する結合鎖4が定義される。このように、ステップS1では、ビーズ3及び結合鎖4が定義されることにより、粗視化モデル6を設定することができる。
【0029】
ここで、「平衡長」とは、ビーズ3、3間の結合距離である。この結合距離が変化した場合は、結合鎖4によって、元の平衡長が保持される。これにより、粗視化モデル6は、直鎖状の三次元構造を維持することができる。なお、平衡長は、隣り合うビーズ3、3の中心間の距離として定義される。また、結合鎖4のポテンシャルには、例えば、論文( 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 )に基づいて設定されるのが望ましい。
【0030】
次に、分子動力学計算に基づいて、粗視化モデルの平衡状態が計算される(平衡状態計算ステップS2)。
図5には、本実施形態の平衡状態計算ステップS2の処理手順の一例が示されている。
【0031】
本実施形態の平衡状態計算ステップS2では、先ず、
図6に示されるように、同一の粗視化モデル6に含まれるビーズ3、3間、及び異なる粗視化モデル6、6のビーズ3、3間に、ビーズ間の距離の関数である相互ポテンシャルP1がそれぞれ定義される(ステップS21)。本実施形態の相互ポテンシャルP1は、LJ(Lennard-Jones)ポテンシャルであり、下記式(2)で定義される。
【0032】
【数2】
ここで、各定数及び変数は、次のとおりである。
r
ij:各ビーズ間距離
r
c:カットオフ距離
ε:各ビーズ間に定義される相互ポテンシャルの強度
σ:各ビーズの直径に相当
なお、ビーズ間の距離r
ij及びカットオフ距離r
cは、各ビーズ3、3の中心3a、3a間の距離として定義される。
【0033】
本実施形態において、相互ポテンシャルP1の強度ε、相互ポテンシャルP1が作用する距離σ、及び、カットオフ距離r
cは、論文(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 に基づいて、下記のように設定される。
強度ε:1.0
距離σ:1.0
カットオフ距離r
c:2
1/6σ
【0034】
これにより、相互ポテンシャルP1は、ビーズ間距離r
ijが2
1/6σ未満になる場合にのみ、ビーズ3、3間に斥力を生じさせる。なお、ビーズ間距離r
ijが2
1/6σ以上になった場合には、相互ポテンシャルP1がゼロになり、ビーズ3、3間に斥力が生じない。このような相互ポテンシャルP1は、分子動力学計算において、実際の高分子材料の分子運動に近似させることができる。このような相互ポテンシャルP1は、コンピュータ1に記憶される。
【0035】
次に、
図7に示されるように、予め定められた体積を持つ仮想空間8内に、粗視化モデル6が配置される(ステップS22)。この仮想空間8は、解析対象の高分子材料の微小構造部分に相当する。
【0036】
本実施形態の仮想空間8は、1辺の長さLsが、例えば20〜100[σ]の立方体として定義される。また、仮想空間8には、例えば、粗視化モデル6が100〜800本程度配置される。このような仮想空間8は、コンピュータ1に記憶される。
【0037】
次に、コンピュータ1が、粗視化モデル6の分子動力学計算を行う(ステップS23)。本実施形態の分子動力学計算では、例えば、仮想空間8について所定の時間、配置した全ての粗視化モデル6が古典力学に従うものとして、ランジュバンの運動方程式が適用される。そして、各時刻での全てのビーズ3(
図6に示す)の動きが追跡される。また、分子動力学計算を行うに際しては、系内の粒子、体積及び温度は一定で行われる。
【0038】
次に、粗視化モデル6の初期配置を十分に緩和できたか否かが判断される(ステップS24)。このステップS24では、粗視化モデル6の初期配置を十分に緩和できたと判断された場合、次のリバースマッピングステップS3が実施される。一方、粗視化モデル6の初期配置を十分に緩和できていないと判断された場合は、単位ステップを進めて(ステップS25)、ステップS23(分子動力学計算)が再度実施される。これにより、平衡状態計算ステップS2では、粗視化モデル6の平衡状態(構造が緩和した状態)を確実に計算することができる。
【0039】
次に、コンピュータ1が、計算された平衡状態の粗視化モデル6の各ビーズ3に、モノマー又はモノマーの一部分をなすモノマーモデル10を割り当てて、平衡状態の全原子モデルを作成する(リバースマッピングステップS3)。このようなリバースマッピングステップS3は、例えば、全原子モデルを仮想空間8に配置して、始めから平衡状態を計算する方法に比べて、モノマー単位での幾何異性体の特徴が再現された平衡状態の全原子モデルの配置を短時間に作成するのに役立つ。
図8には、本実施形態のリバースマッピングステップS3の処理手順の一例が示されている。
【0040】
本実施形態のリバースマッピングステップS3は、
図9に示されるように、ホモポリマー高分子鎖に有効なレプテーション理論で定義されるホモポリマー高分子鎖のチューブ半径R1に基づいて、各ビーズ3にモノマーモデル10を割り当てる。なお、上述したように、ホモポリマー高分子鎖がポリブタジエンである場合、チューブ半径R1は、4.71nmである。また、チューブ半径R1の基準となる中心線13は、粗視化モデル6の結合鎖4の中心線4cと一致している。
【0041】
また、モノマーモデル10は、粒子で表現された少なくとも一つ、本実施形態では1つの構造単位2に含まれる炭素原子の個数分の炭素原子11を含んで構成される。即ち、一つのモノマーモデル10には、1.55個分のモノマーに含まれる炭素原子の個数(例えば、6個又は7個)分の炭素原子11が含まれる。なお、本実施形態のモノマーモデル10には、ポリブタジエンを構成する水素原子が省略されている。
【0042】
炭素原子11は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、炭素原子11には、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ1に記憶される。
【0043】
また、炭素原子11、11間には、平衡長が定義されたポテンシャルが設定される。これにより、炭素原子11、11を連結(拘束)する結合鎖12が定義される。なお、結合鎖12に定義されるポテンシャルは、例えば、論文( S. L. Mayo, B. D. Olafson & W. A. Goddard III 著、「DREIDING: A Generic Force Field for Molecular Simulations」、J. Phys. Chem. 1990, 94, 8897 )に基づいて、汎用パラメータである Dreiding が採用されるのが望ましい。なお、結合鎖12に定義されるポテンシャルとしては、汎用パラメータである Dreiding のみならず、例えば、量子化学計算や第一原理計算などに基づいて得られた精度の高いパラメータが採用されてもよい。
【0044】
また、結合鎖12には、ホモポリマー高分子鎖の構造(本実施形態では、シス構造)に基づいて、結合鎖12を介して連続する3つの炭素原子11、11がなす角度である結合角、及び、結合鎖12を介して連続する4つの炭素原子11において、隣り合う4つの粒子が作る二面角などがそれぞれ定義される。
【0045】
本実施形態のリバースマッピングステップS3は、先ず、モノマー単位での幾何異性体の特徴が再現されるように、最初の1個目のビーズである第1ビーズ3a、及び、最初から2個目のビーズである第2ビーズ3bのモノマー配置が決定される(初期マッピングステップS31)。ここで、「モノマー単位での幾何異性体の特徴を再現する」とは、例えば、100%シス型のモノマーを対象とする場合、シス−トランスの選択肢のあるすべての結合部において、シス型の構造に再現することを意味している。本実施形態では、結合鎖12に定義されている結合角及び二面角等を維持することにより、モノマー単位での幾何異性体の特徴を再現している。
図10には、本実施形態の初期マッピングステップS31の処理手順の一例が示されている。
【0046】
初期マッピングステップS31では、先ず、
図9に示されるように、第1ビーズ3aを構成する第1モノマーモデル10a、及び、第2ビーズ3bを構成する第2モノマーモデル10bを、モノマー単位での幾何異性体の特徴を再現して、第1ビーズ3a及び第2ビーズ3bに配置する(ステップS311)。
【0047】
第1モノマーモデル10a及び第2モノマーモデル10bは、モノマー単位での幾何異性体の構造(シス構造)を維持した状態で、結合鎖12を介して一体に連結されている。本実施形態のステップS311では、第1モノマーモデル10a及び第2モノマーモデル10bを、第1ビーズ3a及び第2ビーズ3bの方向ベクトルV1に沿って配置している。これにより、ステップS311では、第1モノマーモデル10a及び第2モノマーモデル10bを、モノマー単位での幾何異性体の特徴を再現しつつ、第1ビーズ3aと第2ビーズ3bとの間をのびる結合鎖4と略平行に配置することができる。なお、第1モノマーモデル10aの最初の一個目の炭素原子11sは、第1ビーズ3aの表面に固定されるのが望ましい。
【0048】
次に、第1ビーズ3aと第1モノマーモデル10aの炭素原子11との間、及び、第2ビーズ3bと第2モノマーモデル10bの炭素原子11との間に、引力が生じるポテンシャル(以下、単に「引力ポテンシャル」ということがある)P2を設定して分子動力学計算を行う(初期計算ステップS312)。本実施形態の引力ポテンシャルP2は、例えば、下記式(3)で定義される。
P2=−500/r … (3)
ここで、各定数及び変数は次のとおりである。
r:ビーズと炭素原子との距離
なお、ビーズと炭素原子との距離rは、ビーズ3の中心と、炭素原子11の中心との間の距離として定義される。
【0049】
引力ポテンシャルP2によって生じる引力の大きさは、ビーズ3と炭素原子11との距離rが小さいほど大きくなる。これにより、初期計算ステップS312では、
図11に示されるように、第1ビーズ3a及び第1モノマーモデル10aの炭素原子11、並びに、第2ビーズ3b及び第2モノマーモデル10bの炭素原子11を接近させることができる。
【0050】
従って、初期計算ステップS312では、第1モノマーモデル10a及び第2モノマーモデル10bを、平衡状態の粗視化モデル6の第1ビーズ3a及び第2ビーズ3bに沿って配置することができる。また、初期計算ステップS312では、各モノマーモデル10a、10bの全ての炭素原子11に、引力ポテンシャルP2を均等に作用させているため、例えば、人為的に配置される場合に比べて、各モノマーモデル10a、10bの構造が不自然になるのを防ぐことができる。従って、初期計算ステップS312では、各モノマーモデル10a、10bの炭素原子11の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に近づけることができる。
【0051】
上記のような作用を効果的に発揮させるために、分子動力学計算が実施される前に、モノマーモデル10の結合鎖12に定義されるポテンシャルのうち、炭素原子11、11間の二重結合に関する二面角のポテンシャルを大きくするのが望ましい。これにより、モノマーモデル10は、引力ポテンシャルP2の影響を受けて、二面角(炭素原子11、11間の二重結合)が回転するのを防ぐことができる。従って、モノマーモデル10は、ホモポリマー高分子鎖の構造(本実施形態では、シス構造)を維持することができ、炭素原子11の配列を、ホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。なお、初期計算ステップS312において、結合鎖12に定義される二面角のポテンシャルの大きさは、例えば、汎用ポテンシャル Dreiding の50倍〜150倍(本実施形態では、100倍)が望ましい。
【0052】
次に、第1ビーズ3aと第1モノマーモデル10aの結合鎖12との間の距離L1、及び、第2ビーズ3bと第2モノマーモデル10bの結合鎖12との間の距離L2が、チューブ半径R1以下であるか否かが判断される(ステップS313)。なお、各距離L1、L2は、結合鎖12の中心線12cと、ビーズ3の中心cとの最短距離とする。
【0053】
このステップS313では、全ての第1モノマーモデル10a及び全ての第2モノマーモデル10bを構成する結合鎖12において、各距離L1、L2がチューブ半径R1以下であると判断された場合、次のステップS314が実施される。
【0054】
一方、各距離L1、L2がチューブ半径R1以下でないと判断された場合は、単位ステップを進めて(ステップS315)、ステップS312〜ステップS313が再度実施される。これにより、初期マッピングステップS31では、第1モノマーモデル10a及び第2モノマーモデル10bを、各距離L1、L2がチューブ半径R1以下になるように確実に配置することができる。従って、初期マッピングステップS31では、レプテーション理論に基づいて、各モノマーモデル10a、10bの炭素原子11を配置することができるため、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。
【0055】
次に、初期計算ステップS312が終了した後に、第1モノマーモデル10aの位置を固定する(ステップS314)。これにより、第1モノマーモデル10aは、その位置の固定が解除されるまでの間、第1ビーズ3aと第1モノマーモデル10aの結合鎖12との間の距離L1を、チューブ半径R1以下に保持することができる。なお、第1モノマーモデル10aの位置を固定は、該第1モノマーモデル10aの位置情報を、変更不能に固定することにより設定することができる。
【0056】
また、第1モノマーモデル10aは、その位置が固定されても、以後の分子動力学計算の計算対象となり、計算時間が増大するおそれがある。このため、固定された第1モノマーモデル10aは、その位置情報が記憶された後に、仮想空間8(
図7に示す)から一時的に削除されてもよい。これにより、第1モノマーモデル10aは、計算対象外となるため、計算時間の増大を効果的に防ぐことができる。
【0057】
次に、第2ビーズ3bから最後のビーズである第Nビーズ3nまでのモノマー配置を順次決定する(主マッピングステップS32)。
図12には、本実施形態の主マッピングステップS32の処理手順の一例が示されている。
【0058】
主マッピングステップS32では、先ず、初期条件として、定数Mに2が設定される(ステップS321)。この定数Mは、コンピュータ1に記憶される。
【0059】
次に、
図13に示されるように、最初からM個目のビーズである第Mビーズ(例えば、第2ビーズ3b)を構成する第Mモノマーモデル(例えば、第2モノマーモデル10b)に1モノマー分延長して配される第M+1モノマーモデル(例えば、第3モノマーモデル10c)を設定する(ステップS322)。
【0060】
ステップS322では、第Mモノマーモデル(例えば、第2モノマーモデル10b)に、モノマー単位での幾何異性体の特徴を維持しつつ、その構造が安定する方向に、第M+1モノマーモデル(例えば、第3モノマーモデル10c)が延長して配される。ここで、構造が安定する方向とは、例えば、第2モノマーモデル10b及び第3モノマーモデル10cが、予め定められた平衡長、結合角及び二面角を維持できる方向を意味している。また、第Mモノマーモデル及び第M+1モノマーモデルは、結合鎖12を介して一体に連結されている。また、「モノマー単位での幾何異性体の特徴を維持する」とは、例えば、100%シス型のモノマーを対象とする場合、シス−トランスの選択肢のあるすべての結合部において、シス型の構造を維持することを意味している。
【0061】
次に、第Mビーズ(例えば、第2ビーズ3b)と第Mモノマーモデル(例えば、第2モノマーモデル10b)の炭素原子11との間、及び、第M+1ビーズ(例えば、第3ビーズ3c)と第M+1モノマーモデル(例えば、第3モノマーモデル10c)の炭素原子11との間に、引力ポテンシャルP2を設定して分子動力学計算を行う(主計算ステップS323)。引力ポテンシャルP2は、初期計算ステップS312と同様に、上記式(3)で定義されるのが望ましい。
【0062】
主計算ステップS323では、
図14に示されるように、引力ポテンシャルP2により、第2ビーズ3b及び第2モノマーモデル10bの炭素原子11、並びに、第3ビーズ3c及び第3モノマーモデル10cの炭素原子11を接近させることができる。従って、主計算ステップS323では、第2モノマーモデル10b及び第3モノマーモデル10cを、第2ビーズ3b及び第3ビーズ3cに沿って配置することができる。
【0063】
さらに、主計算ステップS323では、初期マッピングステップS31で分子動力学計算が行われた第2モノマーモデル10bの位置が固定されていないため、引力ポテンシャルP2による第3モノマーモデル10cの動きが阻害されない。このため、第3モノマーモデル10cは、シス構造を維持しつつ、第3ビーズ3cに円滑に接近することができる。従って、主計算ステップS323では、各モノマーモデル10b、10cの炭素原子11の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に近づけることができる。
【0064】
また、主計算ステップS323では、初期計算ステップS312と同様に、モノマーモデル10の結合鎖12に定義されるポテンシャルのうち、炭素原子11、11間の二重結合に関する二面角のポテンシャルを大きくするのが望ましい。これにより、モノマーモデル10は、ホモポリマー高分子鎖の構造(シス構造)を確実に維持することができる。
【0065】
次に、第Mビーズ(例えば、第2ビーズ3b)と第Mモノマーモデル(例えば、第2モノマーモデル10b)の結合鎖12との間の距離L2、及び、第M+1ビーズ(例えば、第3ビーズ3c)と第M+1モノマーモデル(例えば、第3モノマーモデル10c)の結合鎖12との間の距離L3が、チューブ半径R1以下であるか否かが判断される(ステップ324)。
【0066】
このステップS324では、例えば、第2モノマーモデル10b及び第3モノマーモデル10cを構成する全ての炭素原子11において、各距離L2、L3がチューブ半径R1以下であると判断された場合、次のステップS325が実施される。
【0067】
一方、各距離L2、L3がチューブ半径R1以下でないと判断された場合は、単位ステップを進めて(ステップS326)、ステップS323〜ステップS324が実施される。これにより、主マッピングステップS32では、各距離L2、L3をチューブ半径R1以下に確実にすることができる。従って、主マッピングステップS32では、レプテーション理論に基づいて、各モノマーモデル10b、10cの炭素原子11を配列することができるため、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。
【0068】
次に、主計算ステップS323が終了した後に、第Mモノマーモデル(例えば、第2モノマーモデル)を固定する(ステップS325)。これにより、第2モノマーモデル10bは、その位置の固定が解除されるまでの間、第2ビーズ3bと第2モノマーモデル10bの結合鎖12との間の距離L2を、チューブ半径R1以下に保持することができる。また、固定された第2モノマーモデル10bは、第1モノマーモデル10aと同様に、その位置情報が記憶された後に、一時的に削除してもよい。
【0069】
次に、第2ビーズ3bから第Nビーズ(図示省略)までのモノマー配置が決定されたか否か(定数MがN−1であるか否か)が判断される(ステップS327)。このステップS327では、第2ビーズ3bから第Nビーズまでのモノマー配置が決定されたと判断された場合、次のステップS33が実施される。
【0070】
一方、第2ビーズ3bから第Nビーズまでのモノマー配置が決定されていない(定数M<N−1)と判断された場合は、定数Mを1増加させて(ステップS328)、ステップS322〜S327が再度実施される。これにより、主マッピングステップS32では、第2モノマーモデル10b〜第Nモノマーモデル(図示省略)を、ビーズ3とモノマーモデル10の結合鎖12との間の距離Lがチューブ半径R1以下になるように配置することができる。従って、主マッピングステップS32では、ホモポリマー高分子鎖に対して有効なレプテーション理論に基づいて、第2モノマーモデル10b〜第Nモノマーモデルの炭素原子11の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。
【0071】
次に、主マッピングステップS32が終了した後に、第1モノマーモデル10a〜第Nモノマーモデル10nの位置の固定を解除して、粗視化モデル6を削除する(ステップS33)。これにより、
図15に示されるように、ステップS33では、第1モノマーモデル10a〜第Nモノマーモデル10nの炭素原子11のみからなる原子モデル16を設定することができる。なお、各モノマーモデル10が一時的に削除されている場合は、ステップS33に先立って、モノマーモデル10を表示させるのが望ましい。
【0072】
次に、原子モデル16の炭素原子11に、相互ポテンシャルを設定して、原子モデル16の分子動力学計算を行う(ステップS34)。
図16に示されるように、相互ポテンシャルP3は、LJ(Lennard-Jones)ポテンシャルであり、上記式(2)で定義される。また、相互ポテンシャルP3は、同一の原子モデル16に含まれる炭素原子11、11間、及び異なる原子モデル16、16の炭素原子11、11間にそれぞれ定義される。
【0073】
このステップS34では、相互ポテンシャルP3を設定して分子動力学計算が実施されることにより、原子モデル16を効果的に緩和することができる。なお、原子モデル16は、平衡状態の粗視化モデル6に基づいて設定されているため、非常に少ないステップ数(例えば、10000〜100000ステップ)で緩和することができる。
【0074】
なお、主マッピングステップS32が終了した直後の原子モデル16は、その炭素原子11が、隣接する原子モデル16の炭素原子11と重なって配置される場合がある。上記式(2)で定義される相互ポテンシャルP3は、粒子間距離が小さくなるほど大きくなる。このため、互いに重なる炭素原子11、11間では、相互ポテンシャルP3が非常に大きくなり、計算が強制的に中断されるおそれがある。
【0075】
従って、ステップS34では、炭素原子11の直径に相当する相互ポテンシャルP3のパラメータσを小さくした状態で、分子動力学計算が開始されるのが望ましい。このパラメータσには、例えば、最も隣接する炭素原子11、11の粒子間距離が設定されるのが望ましい。これにより、ステップS34では、相互ポテンシャルP3が大きくなることに起因する計算の中断を防ぐことができる。また、ステップS34では、炭素原子11の直径に相当する相互ポテンシャルP3のパラメータσを徐々に大きくして、分子動力学計算が実施されるのが望ましい。これにより、ステップS34では、平衡状態の原子モデル16を効率的に計算することができる。
【0076】
次に、
図17に示されるように、ホモポリマー高分子鎖の構造に基づいて、原子モデル16の炭素原子11に、水素原子14が付加される(ステップS35)。これにより、炭素原子11及び水素原子14を含む全原子モデル20を設定することができる。
【0077】
水素原子14は、炭素原子11と同様に、分子動力学計算において、運動方程式の質点として取り扱われる。このため、水素原子14にも、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ1に記憶される。
【0078】
また、水素原子14と炭素原子11との間には、平衡長が定義されたポテンシャルが設定される。これにより、水素原子14と炭素原子11との間を拘束する結合鎖18が定義される。なお、結合鎖18に定義されるポテンシャルは、炭素原子11、11間の結合鎖12に定義されるポテンシャルと同様のものが定義されるのが望ましい。また、結合鎖18には、ホモポリマー高分子鎖の構造に基づいて、結合角及び二面角が設定される。また、結合鎖18には、ホモポリマー高分子鎖の構造に基づいて、結合鎖18を介して連続する3つの炭素原子11又は水素原子14がなす角度である結合角、及び、結合鎖18を介して連続する4つの炭素原子11又は水素原子14が作る二面角などがそれぞれ定義される。
【0079】
なお、結合鎖18に当初の平衡長が定義された状態で水素原子14が付加されると、隣接する全原子モデル20、20の水素原子14、14同士、又は、水素原子14と炭素原子11とが衝突するおそれがある。このため、本実施形態のステップS35では、先ず、結合鎖18の平衡長を小さくした状態で、水素原子14を炭素原子11に付加している。なお、結合鎖18の平衡長には、例えば、当初の平衡長の1/30〜1/10程度に設定されるのが望ましい。
【0080】
次に、隣接する全原子モデル20、20の水素原子14、14間、及び、水素原子14と炭素原子との間に、相互ポテンシャル(図示省略)を設定する。この相互ポテンシャルは、LJ(Lennard-Jones)ポテンシャルであり、上記式(2)で定義される。そして、分子動力学計算を行いながら、結合鎖18の平衡長を徐々に大きくして、当初の平衡長を定義する。これにより、水素原子14は、隣接する全原子モデル20の水素原子14及び炭素原子11を離間させながら、結合鎖18に当初の平衡長を定義することができるため、水素原子14、14同士の衝突、又は水素原子14と炭素原子11との衝突を防ぐことができる。
【0081】
このように、リバースマッピングステップS3では、初期マッピングステップS31及び主マッピングステップS32において、ホモポリマー高分子鎖に対して有効なレプテーション理論に基づいて、各モノマーモデル10の炭素原子11の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。さらに、ステップ35では、炭素原子11の配列に基づいて、水素原子14を精度よく配置することができる。従って、本発明の作成方法では、現実のホモポリマー高分子鎖に精度よく近づけた全原子モデル20を確実に作成することができる。
【実施例】
【0082】
図3に示される手順に従って、ホモポリマー高分子鎖をビーズでモデル化した粗視化モデルを設定し、分子動力学計算に基づいて、粗視化モデルの平衡状態が計算された。また、平衡状態の粗視化モデルの各ビーズに、モノマーモデルを割り当てて平衡状態の全原子モデルが作成された。なお、各モノマーモデルは、各ビーズと各モノマーモデルの結合鎖との間の距離が、レプテーション理論で定義されるホモポリマー高分子鎖のチューブ半径以下になるように割り当てられた(実施例)。
【0083】
また、比較のために、ホモポリマー高分子鎖のチューブ半径を考慮することなく、平衡状態の粗視化モデルの各ビーズに、モノマーモデルを割り当てて平衡状態の全原子モデルが作成された(比較例1)。さらに、ホモポリマー高分子鎖をモデル化した全原子モデルを設定し、分子動力学計算に基づいて、平衡状態の全原子モデルが計算された(比較例2)。
【0084】
そして、実施例、比較例1及び比較例2の作成方法が夫々実施され、全原子モデルの緩和(平衡化)に要する時間が測定された。さらに、実施例1及び比較例1では、下記式(4)に示されるリバースマッピング率が計算された。
リバースマッピング率(%)=Rn/Ra×100…(4)
ここで、各変数は次のとおりである。
Rn:リバースマッピングに成功した粗視化モデルの個数
Ra:粗視化モデルの全個数
【0085】
なお、各ポテンシャルパラメータは、明細書中の記載通りであり、共通仕様は次のとおりである。
コンピュータ:SGI社製のAltix−340
CPUのコア数:8コア(8プロセス並列計算)
ホモポリマー高分子鎖
構造:cis-1,4ポリブタジエン(重合度:2070)
個数:124本
仮想空間:
1辺の長さLs:48.2σ
【0086】
テストの結果、実施例では、全原子モデルの緩和に要した時間が3週間であった。一方、比較例2では、最長緩和時間(全原子モデルが緩和する時間)まで計算するのに必要な時間を見積もった結果、約700万年であり、全原子モデルの平衡状態を計算できなかった。従って、実施例では、平衡状態の全原子モデルを短時間かつ確実に計算しうることが確認できた。
【0087】
また、実施例のリバースマッピング率は100%であったのに対し、比較例1のリバースマッピング率は、実施例の緩和計算が終了した3週間を経過した時点で、12%であった。従って、実施例は、比較例1に比べて、短時間かつ確実に計算しうることが確認できた。