JP6414929B2

JP6414929B2 - 全原子モデルの作成方法

Info

Publication number: JP6414929B2
Application number: JP2014035775A
Authority: JP
Inventors: 容正尾藤; 多田　俊生; 俊生多田; 宏高野; 克美萩田
Original assignee: Sumitomo Rubber Industries Ltd; Keio University
Current assignee: Sumitomo Rubber Industries Ltd; Keio University
Priority date: 2013-04-23
Filing date: 2014-02-26
Publication date: 2018-10-31
Anticipated expiration: 2034-02-26
Also published as: JP2014225226A

Description

本発明は、現実のホモポリマー高分子鎖に近い全原子モデルの配置を効率的に作成する方法に関する。

高分子材料の反応を、コンピュータを用いて評価するために、発明者らは、下記非特許文献１において、高分子材料を構成する高分子鎖から平衡状態の全原子モデルを作成する方法を提案している。この作成方法では、先ず、高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルを設定する。次に、コンピュータが、分子動力学計算に基づいて、粗視化モデルの平衡状態を計算する。そして、コンピュータが、計算された平衡状態の粗視化モデルの各ビーズに、モノマー又はモノマーの一部分をなすモノマーモデルを割り当てて、平衡状態の全原子モデルを作成する（リバースマッピング）。

尾藤容正著、"タイヤ用ゴムの新材料創出に向けた階層的化学計算の応用"、［online］、理化学研究所次世代スーパーコンピュータ開発実施本部、［平成２４年１１月２７日検索］、インターネット＜URL:http://www.nsc.riken.jp/sympo2009/poster-session/abstract/P-28Bito.pdf＞ Macromolecules 2006, 39, 6708

上記のような方法では、例えば、全原子モデルを用いて平衡状態を計算する方法に比べて、平衡状態の全原子モデルを、短時間で作成することができた。しかしながら、作成された全原子モデルでは、現実の高分子鎖に十分に近づけることが難しいという問題があった。

また、現実の高分子鎖に近づける方法として、モノマー単位での幾何異性体の特徴を再現せずに全原子モデルを作成した後に、モノマー単位での幾何異性体の特徴を再現するための人為的なポテンシャル関数を設定して分子動力学計算を実施する方法がある。しかしながら、このような分子動力学計算では、複数本の鎖を含む系の扱うため、１本の鎖の構造を最適化する計算に比べて、多くの計算時間を要するという問題があった。しかも、１００％シス型の全原子モデルを構築する場合には、鎖の長さに応じて、級数的に計算時間が増加する。とりわけ、タイヤ材料ゴムなどの絡み合うほど長い高分子鎖では、非常に多くの計算時間を要するという問題があった。

発明者らは、鋭意研究を重ねた結果、ホモポリマー高分子鎖である場合、ホモポリマー高分子鎖に対して有効なレプテーション理論で定義される高分子鎖のチューブ半径に基づいて、各ビーズにモノマー又はモノマーの一部分を表すモノマーモデルを割り当てることが有効であるとの知見を得て、本発明を完成させるに至った。

本発明は、以上のような実状に鑑み案出されたもので、現実のホモポリマー高分子鎖に精度よく近づけうる全原子モデルの配置を作成する方法を提供することを主たる目的としている。

本発明は、任意のホモポリマー高分子鎖に対して、平衡状態の全原子モデルの配置を効率よく作成するための方法であって、前記ホモポリマー高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルをコンピュータに設定するステップ、前記コンピュータが、分子動力学計算に基づいて、前記粗視化モデルの平衡状態を計算するステップ、及び、前記コンピュータが、前記計算された平衡状態の粗視化モデルの前記各ビーズに、前記モノマー又はモノマーの一部分を表すモノマーモデルを割り当てて平衡状態の全原子モデルを作成するリバースマッピングステップを含み、前記モノマーモデルは、複数の炭素原子と、該炭素原子を連結する結合鎖とを含み、前記ホモポリマー高分子鎖は、前記モノマーの繰り返しによって大部分が構成されるものであり、前記リバースマッピングステップは、最初の１個目のビーズである第１ビーズから第Ｎビーズまでのモノマー配置を、前記ビーズ及び前記モノマーモデルの前記最初からの順序を示す定数Ｍ（但し、１≦Ｍ≦Ｎ−１の整数）を用いて順次決定するものであり、前記リバースマッピングステップは、Ｍ個目のビーズである第Ｍビーズを構成する第Ｍモノマーモデル、及び、Ｍ＋１個目のビーズである第Ｍ＋１ビーズを構成する第Ｍ＋１モノマーモデルを、前記第Ｍビーズ及び前記第Ｍ＋１ビーズに配置する第１ステップ、前記第Ｍビーズと前記第Ｍモノマーモデルの前記炭素原子との間、及び、前記第Ｍ＋１ビーズと前記第Ｍ＋１モノマーモデルの前記炭素原子との間に、引力が生じるポテンシャルを設定して分子動力学計算を行う第２ステップ、並びに、前記第２ステップ終了後に、前記第Ｍモノマーモデルを固定する第３ステップを含み、前記定数ＭがＮ−１に等しくなるまで、前記定数Ｍを１増加させて、前記第１ステップないし前記第３ステップが繰り返し実施されることを特徴とする。

本発明に係る前記全原子モデルの作成方法において、前記第２ステップは、前記第Ｍビーズと前記第Ｍモノマーモデルの前記炭素原子との間の距離、及び、前記第Ｍ＋１ビーズと前記第Ｍ＋１モノマーモデルの前記炭素原子との間の距離が、前記ホモポリマー高分子鎖のチューブ半径以下になるまで、前記分子動力学計算を行うのが望ましい。

本発明に係る前記全原子モデルの作成方法において、前記第３ステップ終了後に、前記第Ｍモノマーモデルを一時的に削除して、前記第Ｍモノマーモデルを前記分子動力学計算での対象外とするステップを含むのが望ましく、前記チューブ半径は、前記ホモポリマー高分子鎖に対して有効なレプテーション理論又は前記レプテーション理論から派生した理論で定義されるのが望ましい。

本発明に係る前記全原子モデルの作成方法において、前記リバースマッピングステップ終了後に、前記モノマーモデルの位置の固定を解除して、前記粗視化モデルを削除するステップをさらに含むのが望ましい。

本発明の全原子モデルの作成方法は、ホモポリマー高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルをコンピュータに設定するステップ、及び、コンピュータが、分子動力学計算に基づいて、粗視化モデルの平衡状態を計算するステップを含む。さらに、本発明の作成方法は、コンピュータが、計算された平衡状態の粗視化モデルの各ビーズに、モノマー又はモノマーの一部分を表すモノマーモデルを割り当てて平衡状態の全原子モデルを作成するリバースマッピングステップを含む。このような作成方法では、例えば、全原子モデルを用いて平衡状態を計算する方法に比べて、平衡状態の全原子モデルの配置を短時間に作成することができる。

さらに、本発明のリバースマッピングステップは、ホモポリマー高分子鎖のチューブ半径に基づいて、各ビーズにモノマーモデルを割り当てる。このホモポリマー高分子鎖のチューブ半径は、現実のホモポリマー高分子鎖の物性や運動を予測するのに用いられるものである。従って、本発明では、全原子モデルを、現実のホモポリマー高分子鎖に精度よく近づけることができる。

本実施形態の全原子モデルの作成方法を実行するコンピュータの斜視図である。ポリブタジエンの構造式である。本実施形態の作成方法のフローチャートである。粗視化モデルを示す概念図である。本実施形態の平衡状態計算ステップのフローチャートである。粗視化モデルのポテンシャルを説明する概念図である。仮想空間に配置された粗視化モデルを示す概念斜視図である。本実施形態のリバースマッピングステップのフローチャートである。粗視化モデル及びモノマーモデルを示す概念図である。本実施形態の初期マッピングステップのフローチャートである。初期計算ステップを説明する概念図である。本実施形態の主マッピングステップのフローチャートである。第２モノマーモデルに延長して配置された第３モノマーモデルを示す概念図である。主計算ステップを説明する概念図である。第１モノマーモデル乃至第Ｎモノマーモデルからなる全原子モデルの概念図である。全原子モデルのポテンシャルを説明する概念図である。炭素原子及び水素原子からなる全原子モデルの概念図である。

以下、本発明の実施の一形態が図面に基づき説明される。
本実施形態の全原子モデルの作成方法（以下、単に「作成方法」ということがある）は、コンピュータを用いて、高分子材料を構成する任意のホモポリマー高分子鎖に対して、モノマー単位での幾何異性体（シス・トランス異性体）の特徴を再現し、かつ、平衡状態に近似する全原子モデルを作成するための方法である。

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

高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。本実施形態の高分子材料は、図２に示されるように、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある）である場合が例示される。このポリブタジエンは、ホモポリマー高分子鎖（単一重合体）から構成される。このホモポリマー高分子鎖は、メチレン基（−ＣＨ_２−）とメチン基（−ＣＨ−）とからなるモノマー｛−［ＣＨ_２−ＣＨ＝ＣＨ−ＣＨ_２］−｝の繰り返し（重合度ｎで連結すること）によって、その大部分（本実施形態では、全体）が構成されている。このホモポリマー高分子鎖は、モノマー単位での幾何異性体に特徴を有している。ここで、「モノマー単位での幾何異性体に特徴を有する」とは、モノマー単位の幾何異性体（シス・トランス異性体）が存在する場合、幾何異性体の割合やホモポリマー高分子鎖中での配置等に何らかの規則性があることを意味している。本実施形態のホモポリマー高分子鎖は、モノマーが１００％シス型の構造からなるという特徴を有している。なお、ホモポリマー高分子鎖は、モノマーが１００％シス型ではない任意の幾何異性体を有するホモポリマー高分子鎖に対しても、本発明の作成方法は有効である。

また、ホモポリマー高分子鎖は、後述するレプテーション理論、又は、レプテーション理論から派生した理論（以下、これらの理論を総括して「レプテーション理論」という。）が用いられることにより、物性を予測しうることが知られている。なお、これらの理論に基づくシミュレーションは、例えば、ソフトマテリアル統合シミュレーターOCTAのPASTA（ Polymer rheology Analyzer with Slip-link model of enTAnglement )を用いて実施することができる。

なお、高分子材料としては、ポリブタジエン以外の高分子材料が採用されてもよい。例えば、ホモポリマー高分子鎖が、１００％シス等のモノマー単位での幾何異性体の特徴を有している天然ゴムや、モノマー単位での幾何異性体の特徴を人工的に付与した高分子材料が採用されてもよい。

図３には、本実施形態の作成方法の処理手順の一例が示されている。本実施形態では、先ず、コンピュータ１に、ホモポリマー高分子鎖をモデル化した粗視化モデルが設定される（ステップＳ１）。なお、粗視化モデルとしては、全原子モデルの分子動力学計算よりも計算コストが小さいものであれば、任意の粗視化ポテンシャル（例えば、angle、torsionのポテンシャル形式）を含むものでもよい。

ステップＳ１では、図２及び図４に示されるように、先ず、ホモポリマー高分子鎖のモノマー又はモノマーの一部分をなす構造単位２を、直径を持った球で表現されるビーズ（粒子）３で置換する。本実施形態では、ホモポリマー高分子鎖のチューブ半径Ｒ１に基づいて、１個のビーズ３に対応する構造単位２のモノマーの個数が決定される。このチューブ半径Ｒ１は、現実のホモポリマー高分子鎖の物性や運動を予測するのに用いられるものである。また、本実施形態のチューブ半径Ｒ１は、ホモポリマー高分子鎖に対して有効なレプテーション理論で定義される。

レプテーション理論は、現実のホモポリマー高分子鎖の物性や運動を予測するのに用いられるものである。レプテーション理論では、一本の高分子鎖をチューブ半径の細長い“ひも”で粗視化して、物性挙動を解析することに成功している。この“ひも”は、内部で高分子鎖の局所的相関を維持し、他の高分子鎖を排除するチューブ半径程度のブロッブと呼ばれる“糸まり”の連結で近似することができる。本発明では、ブロッブの内部に全原子模型を適用することで、粗視化した高分子鎖のモデルから、分子的な視点で解析を行うことができる全原子模型の配置を効率良く求めている。

レプテーション理論において、チューブ半径Ｒ１は、ホモポリマー高分子鎖がチューブ状にモデル化されるときの半径である。チューブ半径Ｒ１は、下記式（１）で定義される。

ここで、各変数は、次のとおりである。
＜Ｒ²＞₀：ホモポリマー高分子鎖の末端間距離の自乗平均
Ｍ：ホモポリマー高分子鎖の分子量
ρ：ホモポリマー高分子鎖の密度（ｋｇ／m³）
Ｔ：絶対温度（Ｋ）
Ｇ_ｅ：平衡弾性率（Ｐａ）
Ｒ：気体定数（ｍ^２ｋｇｓ^−２K^−１ｍｏｌ^−１）

上記式（１）において、ホモポリマー高分子鎖の密度ρ、絶対温度Ｔ、平衡弾性率Ｇ_ｅ及び＜Ｒ²＞₀／Ｍは、実験値である。特に、＜Ｒ²＞₀／Ｍは、高分子材料（例えば、ポリブタジエン）に中性子ビームを照射する中性子小角散乱（ＳＡＮＳ）を実施することによって求められる実験値である。なお、各変数は、例えば、論文（ L,J.Fetters , D.J.Lohse and R.H.Colby 著、「Chain Dimension and Entanglement Spacings」、Physical Properties of Polymers Handbook Second Edition P445-P452 ）に基づいて設定することもできる。

例えば、ホモポリマー高分子鎖がポリブタジエンである場合、チューブ半径Ｒ１は、４．７１nmである。このチューブ半径Ｒ１に基づいて、１．５５個分のモノマーを構造単位２として、対応する一個のビーズ３で当該構造単位２を置換している。なお、１．５５個分のモノマーを構造単位２としたのは、論文（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）のチューブ半径に関する記載に基づくものである。また、ホモポリマー高分子鎖がポリブタジエン以外の場合でも、例えば、上記論文に基づいて、構造単位２を設定することができる。

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

また、ビーズ３、３間には、平衡長が定義されたポテンシャルが設定される。これにより、ビーズ３、３間には、該ビーズ３、３を拘束する結合鎖４が定義される。このように、ステップＳ１では、ビーズ３及び結合鎖４が定義されることにより、粗視化モデル６を設定することができる。

ここで、「平衡長」とは、ビーズ３、３間の結合距離である。この結合距離が変化した場合は、結合鎖４によって、元の平衡長が保持される。これにより、粗視化モデル６は、直鎖状の三次元構造を維持することができる。なお、平衡長は、隣り合うビーズ３、３の中心間の距離として定義される。また、結合鎖４のポテンシャルには、例えば、論文（ 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 ）に基づいて設定されるのが望ましい。

次に、分子動力学計算に基づいて、粗視化モデルの平衡状態が計算される（平衡状態計算ステップＳ２）。図５には、本実施形態の平衡状態計算ステップＳ２の処理手順の一例が示されている。

本実施形態の平衡状態計算ステップＳ２では、先ず、図６に示されるように、同一の粗視化モデル６に含まれるビーズ３、３間、及び異なる粗視化モデル６、６のビーズ３、３間に、ビーズ間の距離の関数である相互ポテンシャルＰ１がそれぞれ定義される（ステップＳ２１）。本実施形態の相互ポテンシャルＰ１は、ＬＪ（Lennard-Jones）ポテンシャルであり、下記式（２）で定義される。

ここで、各定数及び変数は、次のとおりである。
ｒ_ij：各ビーズ間距離
ｒ_c：カットオフ距離
ε：各ビーズ間に定義される相互ポテンシャルの強度
σ：各ビーズの直径に相当
なお、ビーズ間の距離ｒ_ij及びカットオフ距離ｒ_cは、各ビーズ３、３の中心３ａ、３ａ間の距離として定義される。

本実施形態において、相互ポテンシャルＰ１の強度ε、相互ポテンシャルＰ１が作用する距離σ、及び、カットオフ距離ｒ_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 に基づいて、下記のように設定される。
強度ε：１．０
距離σ：１．０
カットオフ距離ｒ_c：2^1/6σ

これにより、相互ポテンシャルＰ１は、ビーズ間距離ｒ_ijが2^1/6σ未満になる場合にのみ、ビーズ３、３間に斥力を生じさせる。なお、ビーズ間距離ｒ_ijが2^1/6σ以上になった場合には、相互ポテンシャルＰ１がゼロになり、ビーズ３、３間に斥力が生じない。このような相互ポテンシャルＰ１は、分子動力学計算において、実際の高分子材料の分子運動に近似させることができる。このような相互ポテンシャルＰ１は、コンピュータ１に記憶される。

次に、図７に示されるように、予め定められた体積を持つ仮想空間８内に、粗視化モデル６が配置される（ステップＳ２２）。この仮想空間８は、解析対象の高分子材料の微小構造部分に相当する。

本実施形態の仮想空間８は、１辺の長さＬｓが、例えば２０〜１００［σ］の立方体として定義される。また、仮想空間８には、例えば、粗視化モデル６が１００〜８００本程度配置される。このような仮想空間８は、コンピュータ１に記憶される。

次に、コンピュータ１が、粗視化モデル６の分子動力学計算を行う（ステップＳ２３）。本実施形態の分子動力学計算では、例えば、仮想空間８について所定の時間、配置した全ての粗視化モデル６が古典力学に従うものとして、ランジュバンの運動方程式が適用される。そして、各時刻での全てのビーズ３（図６に示す）の動きが追跡される。また、分子動力学計算を行うに際しては、系内の粒子、体積及び温度は一定で行われる。

次に、粗視化モデル６の初期配置を十分に緩和できたか否かが判断される（ステップＳ２４）。このステップＳ２４では、粗視化モデル６の初期配置を十分に緩和できたと判断された場合、次のリバースマッピングステップＳ３が実施される。一方、粗視化モデル６の初期配置を十分に緩和できていないと判断された場合は、単位ステップを進めて（ステップＳ２５）、ステップＳ２３（分子動力学計算）が再度実施される。これにより、平衡状態計算ステップＳ２では、粗視化モデル６の平衡状態（構造が緩和した状態）を確実に計算することができる。

次に、コンピュータ１が、計算された平衡状態の粗視化モデル６の各ビーズ３に、モノマー又はモノマーの一部分をなすモノマーモデル１０を割り当てて、平衡状態の全原子モデルを作成する（リバースマッピングステップＳ３）。このようなリバースマッピングステップＳ３は、例えば、全原子モデルを仮想空間８に配置して、始めから平衡状態を計算する方法に比べて、モノマー単位での幾何異性体の特徴が再現された平衡状態の全原子モデルの配置を短時間に作成するのに役立つ。図８には、本実施形態のリバースマッピングステップＳ３の処理手順の一例が示されている。

本実施形態のリバースマッピングステップＳ３は、図９に示されるように、ホモポリマー高分子鎖に有効なレプテーション理論で定義されるホモポリマー高分子鎖のチューブ半径Ｒ１に基づいて、各ビーズ３にモノマーモデル１０を割り当てる。なお、上述したように、ホモポリマー高分子鎖がポリブタジエンである場合、チューブ半径Ｒ１は、４．７１nmである。また、チューブ半径Ｒ１の基準となる中心線１３は、粗視化モデル６の結合鎖４の中心線４ｃと一致している。

また、モノマーモデル１０は、粒子で表現された少なくとも一つ、本実施形態では１つの構造単位２に含まれる炭素原子の個数分の炭素原子１１を含んで構成される。即ち、一つのモノマーモデル１０には、１．５５個分のモノマーに含まれる炭素原子の個数（例えば、６個又は７個）分の炭素原子１１が含まれる。なお、本実施形態のモノマーモデル１０には、ポリブタジエンを構成する水素原子が省略されている。

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

また、炭素原子１１、１１間には、平衡長が定義されたポテンシャルが設定される。これにより、炭素原子１１、１１を連結（拘束）する結合鎖１２が定義される。なお、結合鎖１２に定義されるポテンシャルは、例えば、論文（ 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 が採用されるのが望ましい。なお、結合鎖１２に定義されるポテンシャルとしては、汎用パラメータである Dreiding のみならず、例えば、量子化学計算や第一原理計算などに基づいて得られた精度の高いパラメータが採用されてもよい。

また、結合鎖１２には、ホモポリマー高分子鎖の構造（本実施形態では、シス構造）に基づいて、結合鎖１２を介して連続する３つの炭素原子１１、１１がなす角度である結合角、及び、結合鎖１２を介して連続する４つの炭素原子１１において、隣り合う４つの粒子が作る二面角などがそれぞれ定義される。

本実施形態のリバースマッピングステップＳ３は、先ず、モノマー単位での幾何異性体の特徴が再現されるように、最初の１個目のビーズである第１ビーズ３ａ、及び、最初から２個目のビーズである第２ビーズ３ｂのモノマー配置が決定される（初期マッピングステップＳ３１）。ここで、「モノマー単位での幾何異性体の特徴を再現する」とは、例えば、１００％シス型のモノマーを対象とする場合、シス−トランスの選択肢のあるすべての結合部において、シス型の構造に再現することを意味している。本実施形態では、結合鎖１２に定義されている結合角及び二面角等を維持することにより、モノマー単位での幾何異性体の特徴を再現している。図１０には、本実施形態の初期マッピングステップＳ３１の処理手順の一例が示されている。

初期マッピングステップＳ３１では、先ず、図９に示されるように、第１ビーズ３ａを構成する第１モノマーモデル１０ａ、及び、第２ビーズ３ｂを構成する第２モノマーモデル１０ｂを、モノマー単位での幾何異性体の特徴を再現して、第１ビーズ３ａ及び第２ビーズ３ｂに配置する（ステップＳ３１１）。

第１モノマーモデル１０ａ及び第２モノマーモデル１０ｂは、モノマー単位での幾何異性体の構造（シス構造）を維持した状態で、結合鎖１２を介して一体に連結されている。本実施形態のステップＳ３１１では、第１モノマーモデル１０ａ及び第２モノマーモデル１０ｂを、第１ビーズ３ａ及び第２ビーズ３ｂの方向ベクトルＶ１に沿って配置している。これにより、ステップＳ３１１では、第１モノマーモデル１０ａ及び第２モノマーモデル１０ｂを、モノマー単位での幾何異性体の特徴を再現しつつ、第１ビーズ３ａと第２ビーズ３ｂとの間をのびる結合鎖４と略平行に配置することができる。なお、第１モノマーモデル１０ａの最初の一個目の炭素原子１１ｓは、第１ビーズ３ａの表面に固定されるのが望ましい。

次に、第１ビーズ３ａと第１モノマーモデル１０ａの炭素原子１１との間、及び、第２ビーズ３ｂと第２モノマーモデル１０ｂの炭素原子１１との間に、引力が生じるポテンシャル（以下、単に「引力ポテンシャル」ということがある）Ｐ２を設定して分子動力学計算を行う（初期計算ステップＳ３１２）。本実施形態の引力ポテンシャルＰ２は、例えば、下記式（３）で定義される。
Ｐ２＝−５００／ｒ … （３）
ここで、各定数及び変数は次のとおりである。
ｒ：ビーズと炭素原子との距離
なお、ビーズと炭素原子との距離ｒは、ビーズ３の中心と、炭素原子１１の中心との間の距離として定義される。

引力ポテンシャルＰ２によって生じる引力の大きさは、ビーズ３と炭素原子１１との距離ｒが小さいほど大きくなる。これにより、初期計算ステップＳ３１２では、図１１に示されるように、第１ビーズ３ａ及び第１モノマーモデル１０ａの炭素原子１１、並びに、第２ビーズ３ｂ及び第２モノマーモデル１０ｂの炭素原子１１を接近させることができる。

従って、初期計算ステップＳ３１２では、第１モノマーモデル１０ａ及び第２モノマーモデル１０ｂを、平衡状態の粗視化モデル６の第１ビーズ３ａ及び第２ビーズ３ｂに沿って配置することができる。また、初期計算ステップＳ３１２では、各モノマーモデル１０ａ、１０ｂの全ての炭素原子１１に、引力ポテンシャルＰ２を均等に作用させているため、例えば、人為的に配置される場合に比べて、各モノマーモデル１０ａ、１０ｂの構造が不自然になるのを防ぐことができる。従って、初期計算ステップＳ３１２では、各モノマーモデル１０ａ、１０ｂの炭素原子１１の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に近づけることができる。

上記のような作用を効果的に発揮させるために、分子動力学計算が実施される前に、モノマーモデル１０の結合鎖１２に定義されるポテンシャルのうち、炭素原子１１、１１間の二重結合に関する二面角のポテンシャルを大きくするのが望ましい。これにより、モノマーモデル１０は、引力ポテンシャルＰ２の影響を受けて、二面角（炭素原子１１、１１間の二重結合）が回転するのを防ぐことができる。従って、モノマーモデル１０は、ホモポリマー高分子鎖の構造（本実施形態では、シス構造）を維持することができ、炭素原子１１の配列を、ホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。なお、初期計算ステップＳ３１２において、結合鎖１２に定義される二面角のポテンシャルの大きさは、例えば、汎用ポテンシャル Dreiding の５０倍〜１５０倍（本実施形態では、１００倍）が望ましい。

次に、第１ビーズ３ａと第１モノマーモデル１０ａの結合鎖１２との間の距離Ｌ１、及び、第２ビーズ３ｂと第２モノマーモデル１０ｂの結合鎖１２との間の距離Ｌ２が、チューブ半径Ｒ１以下であるか否かが判断される（ステップＳ３１３）。なお、各距離Ｌ１、Ｌ２は、結合鎖１２の中心線１２ｃと、ビーズ３の中心ｃとの最短距離とする。

このステップＳ３１３では、全ての第１モノマーモデル１０ａ及び全ての第２モノマーモデル１０ｂを構成する結合鎖１２において、各距離Ｌ１、Ｌ２がチューブ半径Ｒ１以下であると判断された場合、次のステップＳ３１４が実施される。

一方、各距離Ｌ１、Ｌ２がチューブ半径Ｒ１以下でないと判断された場合は、単位ステップを進めて（ステップＳ３１５）、ステップＳ３１２〜ステップＳ３１３が再度実施される。これにより、初期マッピングステップＳ３１では、第１モノマーモデル１０ａ及び第２モノマーモデル１０ｂを、各距離Ｌ１、Ｌ２がチューブ半径Ｒ１以下になるように確実に配置することができる。従って、初期マッピングステップＳ３１では、レプテーション理論に基づいて、各モノマーモデル１０ａ、１０ｂの炭素原子１１を配置することができるため、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。

次に、初期計算ステップＳ３１２が終了した後に、第１モノマーモデル１０ａの位置を固定する（ステップＳ３１４）。これにより、第１モノマーモデル１０ａは、その位置の固定が解除されるまでの間、第１ビーズ３ａと第１モノマーモデル１０ａの結合鎖１２との間の距離Ｌ１を、チューブ半径Ｒ１以下に保持することができる。なお、第１モノマーモデル１０ａの位置を固定は、該第１モノマーモデル１０ａの位置情報を、変更不能に固定することにより設定することができる。

また、第１モノマーモデル１０ａは、その位置が固定されても、以後の分子動力学計算の計算対象となり、計算時間が増大するおそれがある。このため、固定された第１モノマーモデル１０ａは、その位置情報が記憶された後に、仮想空間８（図７に示す）から一時的に削除されてもよい。これにより、第１モノマーモデル１０ａは、計算対象外となるため、計算時間の増大を効果的に防ぐことができる。

次に、第２ビーズ３ｂから最後のビーズである第Ｎビーズ３ｎまでのモノマー配置を順次決定する（主マッピングステップＳ３２）。図１２には、本実施形態の主マッピングステップＳ３２の処理手順の一例が示されている。

主マッピングステップＳ３２では、先ず、初期条件として、定数Ｍに２が設定される（ステップＳ３２１）。この定数Ｍは、コンピュータ１に記憶される。

次に、図１３に示されるように、最初からＭ個目のビーズである第Ｍビーズ（例えば、第２ビーズ３ｂ）を構成する第Ｍモノマーモデル（例えば、第２モノマーモデル１０ｂ）に１モノマー分延長して配される第Ｍ＋１モノマーモデル（例えば、第３モノマーモデル１０ｃ）を設定する（ステップＳ３２２）。

ステップＳ３２２では、第Ｍモノマーモデル（例えば、第２モノマーモデル１０ｂ）に、モノマー単位での幾何異性体の特徴を維持しつつ、その構造が安定する方向に、第Ｍ＋１モノマーモデル（例えば、第３モノマーモデル１０ｃ）が延長して配される。ここで、構造が安定する方向とは、例えば、第２モノマーモデル１０ｂ及び第３モノマーモデル１０ｃが、予め定められた平衡長、結合角及び二面角を維持できる方向を意味している。また、第Ｍモノマーモデル及び第Ｍ＋１モノマーモデルは、結合鎖１２を介して一体に連結されている。また、「モノマー単位での幾何異性体の特徴を維持する」とは、例えば、１００％シス型のモノマーを対象とする場合、シス−トランスの選択肢のあるすべての結合部において、シス型の構造を維持することを意味している。

次に、第Ｍビーズ（例えば、第２ビーズ３ｂ）と第Ｍモノマーモデル（例えば、第２モノマーモデル１０ｂ）の炭素原子１１との間、及び、第Ｍ＋１ビーズ（例えば、第３ビーズ３ｃ）と第Ｍ＋１モノマーモデル（例えば、第３モノマーモデル１０ｃ）の炭素原子１１との間に、引力ポテンシャルＰ２を設定して分子動力学計算を行う（主計算ステップＳ３２３）。引力ポテンシャルＰ２は、初期計算ステップＳ３１２と同様に、上記式（３）で定義されるのが望ましい。

主計算ステップＳ３２３では、図１４に示されるように、引力ポテンシャルＰ２により、第２ビーズ３ｂ及び第２モノマーモデル１０ｂの炭素原子１１、並びに、第３ビーズ３ｃ及び第３モノマーモデル１０ｃの炭素原子１１を接近させることができる。従って、主計算ステップＳ３２３では、第２モノマーモデル１０ｂ及び第３モノマーモデル１０ｃを、第２ビーズ３ｂ及び第３ビーズ３ｃに沿って配置することができる。

さらに、主計算ステップＳ３２３では、初期マッピングステップＳ３１で分子動力学計算が行われた第２モノマーモデル１０ｂの位置が固定されていないため、引力ポテンシャルＰ２による第３モノマーモデル１０ｃの動きが阻害されない。このため、第３モノマーモデル１０ｃは、シス構造を維持しつつ、第３ビーズ３ｃに円滑に接近することができる。従って、主計算ステップＳ３２３では、各モノマーモデル１０ｂ、１０ｃの炭素原子１１の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に近づけることができる。

また、主計算ステップＳ３２３では、初期計算ステップＳ３１２と同様に、モノマーモデル１０の結合鎖１２に定義されるポテンシャルのうち、炭素原子１１、１１間の二重結合に関する二面角のポテンシャルを大きくするのが望ましい。これにより、モノマーモデル１０は、ホモポリマー高分子鎖の構造（シス構造）を確実に維持することができる。

次に、第Ｍビーズ（例えば、第２ビーズ３ｂ）と第Ｍモノマーモデル（例えば、第２モノマーモデル１０ｂ）の結合鎖１２との間の距離Ｌ２、及び、第Ｍ＋１ビーズ（例えば、第３ビーズ３ｃ）と第Ｍ＋１モノマーモデル（例えば、第３モノマーモデル１０ｃ）の結合鎖１２との間の距離Ｌ３が、チューブ半径Ｒ１以下であるか否かが判断される（ステップ３２４）。

このステップＳ３２４では、例えば、第２モノマーモデル１０ｂ及び第３モノマーモデル１０ｃを構成する全ての炭素原子１１において、各距離Ｌ２、Ｌ３がチューブ半径Ｒ１以下であると判断された場合、次のステップＳ３２５が実施される。

一方、各距離Ｌ２、Ｌ３がチューブ半径Ｒ１以下でないと判断された場合は、単位ステップを進めて（ステップＳ３２６）、ステップＳ３２３〜ステップＳ３２４が実施される。これにより、主マッピングステップＳ３２では、各距離Ｌ２、Ｌ３をチューブ半径Ｒ１以下に確実にすることができる。従って、主マッピングステップＳ３２では、レプテーション理論に基づいて、各モノマーモデル１０ｂ、１０ｃの炭素原子１１を配列することができるため、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。

次に、主計算ステップＳ３２３が終了した後に、第Ｍモノマーモデル（例えば、第２モノマーモデル）を固定する（ステップＳ３２５）。これにより、第２モノマーモデル１０ｂは、その位置の固定が解除されるまでの間、第２ビーズ３ｂと第２モノマーモデル１０ｂの結合鎖１２との間の距離Ｌ２を、チューブ半径Ｒ１以下に保持することができる。また、固定された第２モノマーモデル１０ｂは、第１モノマーモデル１０ａと同様に、その位置情報が記憶された後に、一時的に削除してもよい。

次に、第２ビーズ３ｂから第Ｎビーズ（図示省略）までのモノマー配置が決定されたか否か（定数ＭがＮ−１であるか否か）が判断される（ステップＳ３２７）。このステップＳ３２７では、第２ビーズ３ｂから第Ｎビーズまでのモノマー配置が決定されたと判断された場合、次のステップＳ３３が実施される。

一方、第２ビーズ３ｂから第Ｎビーズまでのモノマー配置が決定されていない（定数Ｍ＜Ｎ−１）と判断された場合は、定数Ｍを１増加させて（ステップＳ３２８）、ステップＳ３２２〜Ｓ３２７が再度実施される。これにより、主マッピングステップＳ３２では、第２モノマーモデル１０ｂ〜第Ｎモノマーモデル（図示省略）を、ビーズ３とモノマーモデル１０の結合鎖１２との間の距離Ｌがチューブ半径Ｒ１以下になるように配置することができる。従って、主マッピングステップＳ３２では、ホモポリマー高分子鎖に対して有効なレプテーション理論に基づいて、第２モノマーモデル１０ｂ〜第Ｎモノマーモデルの炭素原子１１の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。

次に、主マッピングステップＳ３２が終了した後に、第１モノマーモデル１０ａ〜第Ｎモノマーモデル１０ｎの位置の固定を解除して、粗視化モデル６を削除する（ステップＳ３３）。これにより、図１５に示されるように、ステップＳ３３では、第１モノマーモデル１０ａ〜第Ｎモノマーモデル１０ｎの炭素原子１１のみからなる原子モデル１６を設定することができる。なお、各モノマーモデル１０が一時的に削除されている場合は、ステップＳ３３に先立って、モノマーモデル１０を表示させるのが望ましい。

次に、原子モデル１６の炭素原子１１に、相互ポテンシャルを設定して、原子モデル１６の分子動力学計算を行う（ステップＳ３４）。図１６に示されるように、相互ポテンシャルＰ３は、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（２）で定義される。また、相互ポテンシャルＰ３は、同一の原子モデル１６に含まれる炭素原子１１、１１間、及び異なる原子モデル１６、１６の炭素原子１１、１１間にそれぞれ定義される。

このステップＳ３４では、相互ポテンシャルＰ３を設定して分子動力学計算が実施されることにより、原子モデル１６を効果的に緩和することができる。なお、原子モデル１６は、平衡状態の粗視化モデル６に基づいて設定されているため、非常に少ないステップ数（例えば、１００００〜１０００００ステップ）で緩和することができる。

なお、主マッピングステップＳ３２が終了した直後の原子モデル１６は、その炭素原子１１が、隣接する原子モデル１６の炭素原子１１と重なって配置される場合がある。上記式（２）で定義される相互ポテンシャルＰ３は、粒子間距離が小さくなるほど大きくなる。このため、互いに重なる炭素原子１１、１１間では、相互ポテンシャルＰ３が非常に大きくなり、計算が強制的に中断されるおそれがある。

従って、ステップＳ３４では、炭素原子１１の直径に相当する相互ポテンシャルＰ３のパラメータσを小さくした状態で、分子動力学計算が開始されるのが望ましい。このパラメータσには、例えば、最も隣接する炭素原子１１、１１の粒子間距離が設定されるのが望ましい。これにより、ステップＳ３４では、相互ポテンシャルＰ３が大きくなることに起因する計算の中断を防ぐことができる。また、ステップＳ３４では、炭素原子１１の直径に相当する相互ポテンシャルＰ３のパラメータσを徐々に大きくして、分子動力学計算が実施されるのが望ましい。これにより、ステップＳ３４では、平衡状態の原子モデル１６を効率的に計算することができる。

次に、図１７に示されるように、ホモポリマー高分子鎖の構造に基づいて、原子モデル１６の炭素原子１１に、水素原子１４が付加される（ステップＳ３５）。これにより、炭素原子１１及び水素原子１４を含む全原子モデル２０を設定することができる。

水素原子１４は、炭素原子１１と同様に、分子動力学計算において、運動方程式の質点として取り扱われる。このため、水素原子１４にも、質量、体積、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ１に記憶される。

また、水素原子１４と炭素原子１１との間には、平衡長が定義されたポテンシャルが設定される。これにより、水素原子１４と炭素原子１１との間を拘束する結合鎖１８が定義される。なお、結合鎖１８に定義されるポテンシャルは、炭素原子１１、１１間の結合鎖１２に定義されるポテンシャルと同様のものが定義されるのが望ましい。また、結合鎖１８には、ホモポリマー高分子鎖の構造に基づいて、結合角及び二面角が設定される。また、結合鎖１８には、ホモポリマー高分子鎖の構造に基づいて、結合鎖１８を介して連続する３つの炭素原子１１又は水素原子１４がなす角度である結合角、及び、結合鎖１８を介して連続する４つの炭素原子１１又は水素原子１４が作る二面角などがそれぞれ定義される。

なお、結合鎖１８に当初の平衡長が定義された状態で水素原子１４が付加されると、隣接する全原子モデル２０、２０の水素原子１４、１４同士、又は、水素原子１４と炭素原子１１とが衝突するおそれがある。このため、本実施形態のステップＳ３５では、先ず、結合鎖１８の平衡長を小さくした状態で、水素原子１４を炭素原子１１に付加している。なお、結合鎖１８の平衡長には、例えば、当初の平衡長の１／３０〜１／１０程度に設定されるのが望ましい。

次に、隣接する全原子モデル２０、２０の水素原子１４、１４間、及び、水素原子１４と炭素原子との間に、相互ポテンシャル（図示省略）を設定する。この相互ポテンシャルは、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（２）で定義される。そして、分子動力学計算を行いながら、結合鎖１８の平衡長を徐々に大きくして、当初の平衡長を定義する。これにより、水素原子１４は、隣接する全原子モデル２０の水素原子１４及び炭素原子１１を離間させながら、結合鎖１８に当初の平衡長を定義することができるため、水素原子１４、１４同士の衝突、又は水素原子１４と炭素原子１１との衝突を防ぐことができる。

このように、リバースマッピングステップＳ３では、初期マッピングステップＳ３１及び主マッピングステップＳ３２において、ホモポリマー高分子鎖に対して有効なレプテーション理論に基づいて、各モノマーモデル１０の炭素原子１１の配列を、現実のホモポリマー高分子鎖の炭素原子の配列に精度よく近づけることができる。さらに、ステップ３５では、炭素原子１１の配列に基づいて、水素原子１４を精度よく配置することができる。従って、本発明の作成方法では、現実のホモポリマー高分子鎖に精度よく近づけた全原子モデル２０を確実に作成することができる。

図３に示される手順に従って、ホモポリマー高分子鎖をビーズでモデル化した粗視化モデルを設定し、分子動力学計算に基づいて、粗視化モデルの平衡状態が計算された。また、平衡状態の粗視化モデルの各ビーズに、モノマーモデルを割り当てて平衡状態の全原子モデルが作成された。なお、各モノマーモデルは、各ビーズと各モノマーモデルの結合鎖との間の距離が、レプテーション理論で定義されるホモポリマー高分子鎖のチューブ半径以下になるように割り当てられた（実施例）。

また、比較のために、ホモポリマー高分子鎖のチューブ半径を考慮することなく、平衡状態の粗視化モデルの各ビーズに、モノマーモデルを割り当てて平衡状態の全原子モデルが作成された（比較例１）。さらに、ホモポリマー高分子鎖をモデル化した全原子モデルを設定し、分子動力学計算に基づいて、平衡状態の全原子モデルが計算された（比較例２）。

そして、実施例、比較例１及び比較例２の作成方法が夫々実施され、全原子モデルの緩和（平衡化）に要する時間が測定された。さらに、実施例１及び比較例１では、下記式（４）に示されるリバースマッピング率が計算された。
リバースマッピング率（％）＝Ｒｎ／Ｒａ×１００…（４）
ここで、各変数は次のとおりである。
Ｒｎ：リバースマッピングに成功した粗視化モデルの個数
Ｒａ：粗視化モデルの全個数

なお、各ポテンシャルパラメータは、明細書中の記載通りであり、共通仕様は次のとおりである。
コンピュータ：ＳＧＩ社製のＡｌｔｉｘ−３４０
ＣＰＵのコア数：８コア（８プロセス並列計算）
ホモポリマー高分子鎖
構造：cis-１，４ポリブタジエン（重合度：２０７０）
個数：１２４本
仮想空間：
１辺の長さＬｓ：４８．２σ

テストの結果、実施例では、全原子モデルの緩和に要した時間が３週間であった。一方、比較例２では、最長緩和時間（全原子モデルが緩和する時間）まで計算するのに必要な時間を見積もった結果、約７００万年であり、全原子モデルの平衡状態を計算できなかった。従って、実施例では、平衡状態の全原子モデルを短時間かつ確実に計算しうることが確認できた。

また、実施例のリバースマッピング率は１００％であったのに対し、比較例１のリバースマッピング率は、実施例の緩和計算が終了した３週間を経過した時点で、１２％であった。従って、実施例は、比較例１に比べて、短時間かつ確実に計算しうることが確認できた。

３ビーズ
６粗視化モデル
１０モノマーモデル
２０全原子モデル

Claims

任意のホモポリマー高分子鎖に対して、平衡状態の全原子モデルの配置を効率よく作成するための方法であって、
前記ホモポリマー高分子鎖を、モノマー又はモノマーの一部分をなす構造単位をビーズで置換した粗視化モデルをコンピュータに設定するステップ、
前記コンピュータが、分子動力学計算に基づいて、前記粗視化モデルの平衡状態を計算するステップ、及び、
前記コンピュータが、前記計算された平衡状態の粗視化モデルの前記各ビーズに、前記モノマー又はモノマーの一部分を表すモノマーモデルを割り当てて平衡状態の全原子モデルを作成するリバースマッピングステップを含み、
前記モノマーモデルは、複数の炭素原子と、該炭素原子を連結する結合鎖とを含み、
前記ホモポリマー高分子鎖は、前記モノマーの繰り返しによって大部分が構成されるものであり、
前記リバースマッピングステップは、最初の１個目のビーズである第１ビーズから第Ｎビーズまでのモノマー配置を、前記ビーズ及び前記モノマーモデルの前記最初からの順序を示す定数Ｍ（但し、１≦Ｍ≦Ｎ−１の整数）を用いて順次決定するものであり、
前記リバースマッピングステップは、
Ｍ個目のビーズである第Ｍビーズを構成する第Ｍモノマーモデル、及び、Ｍ＋１個目のビーズである第Ｍ＋１ビーズを構成する第Ｍ＋１モノマーモデルを、前記第Ｍビーズ及び前記第Ｍ＋１ビーズに配置する第１ステップ、
前記第Ｍビーズと前記第Ｍモノマーモデルの前記炭素原子との間、及び、前記第Ｍ＋１ビーズと前記第Ｍ＋１モノマーモデルの前記炭素原子との間に、引力が生じるポテンシャルを設定して分子動力学計算を行う第２ステップ、並びに、
前記第２ステップ終了後に、前記第Ｍモノマーモデルを固定する第３ステップを含み、
前記定数ＭがＮ−１に等しくなるまで、前記定数Ｍを１増加させて、前記第１ステップないし前記第３ステップが繰り返し実施されることを特徴とする全原子モデルの作成方法。
前記第２ステップは、前記第Ｍビーズと前記第Ｍモノマーモデルの前記炭素原子との間の距離、及び、前記第Ｍ＋１ビーズと前記第Ｍ＋１モノマーモデルの前記炭素原子との間の距離が、前記ホモポリマー高分子鎖のチューブ半径以下になるまで、前記分子動力学計算を行う請求項１記載の全原子モデルの作成方法。
前記第３ステップ終了後に、前記第Ｍモノマーモデルを一時的に削除して、前記第Ｍモノマーモデルを前記分子動力学計算での対象外とするステップを含む請求項１又は２記載の全原子モデルの作成方法。
前記チューブ半径は、前記ホモポリマー高分子鎖に対して有効なレプテーション理論又は前記レプテーション理論から派生した理論で定義される請求項１乃至３のいずれかに記載の全原子モデルの作成方法。
前記リバースマッピングステップ終了後に、前記モノマーモデルの位置の固定を解除して、前記粗視化モデルを削除するステップをさらに含む請求項１乃至４のいずれかに記載の全原子モデルの作成方法。