(57)【要約】      （修正有）【課題】粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算しうる。【解決手段】高分子材料のシミュレーション方法は、コンピュータが、粗視化モデルのRouseのパラメータを求める第１計算工程Ｓ３、コンピュータが、全原子モデルのRouseのパラメータを求める第２計算工程Ｓ６、及び、コンピュータが、粗視化モデルのRouseのパラメータと、全原子モデルのRouseのパラメータとに基づいて、粗視化モデルを用いた分子動力学計算での単位系を、任意の高分子鎖の単位系に換算する換算工程Ｓ７を含む。【選択図】図２
2014206913.tif 000002
【特許請求の範囲】
【請求項１】
コンピュータに、高分子鎖を、複数のビーズでモデル化した粗視化モデルを入力する工程、
前記コンピュータが、予め定めた仮想空間に配置された前記粗視化モデルを用いて分子動力学計算を行う工程、
前記コンピュータが、前記粗視化モデルのRouseのパラメータを求める第１計算工程、
前記コンピュータに、任意の高分子鎖を、複数の炭素を含む原子でモデル化した全原子モデルを入力する工程、
前記コンピュータが、予め定めた仮想空間に配置された前記全原子モデルを用いて分子動力学計算を行う工程、
前記コンピュータが、前記全原子モデルのRouseのパラメータを求める第２計算工程、及び、
前記コンピュータが、前記粗視化モデルのRouseのパラメータと、前記全原子モデルのRouseのパラメータとに基づいて、前記粗視化モデルを用いた分子動力学計算での単位系を、前記任意の高分子鎖の単位系に換算する換算工程を含み、
前記第２計算工程は、前記全原子モデルの全長Ｌfを求める全長計算工程を含み、
前記全長計算工程は、前記コンピュータが、前記全原子モデルの前記炭素原子の一端の固定する工程と、
前記コンピュータが、分子動力学計算を行いながら、前記全原子モデルの前記炭素原子の他端を、前記一端から離れる方向へ強制移動させる工程と、
前記コンピュータが、前記分子動力学計算によって強制的に引き伸ばされた前記全原子モデルを分子力学法による計算を行うことにより、前記全原子モデルの構造を安定化する工程と、
前記コンピュータが、前記全原子モデルの前記一端と前記他端との間の距離から前記全原子モデルの全長Ｌfを求める工程とを含むことを特徴とする高分子材料のシミュレーション方法。
【請求項２】
前記粗視化モデルのRouseのパラメータは、長さのパラメータｂcと、時間のパラメータτcと、摩擦のパラメータζcとを含み、
前記全原子モデルのRouseのパラメータは、長さのパラメータｂfと、時間のパラメータτfと、摩擦のパラメータζfとを含み、
前記換算工程は、前記粗視化モデルを用いた分子動力学計算での単位長さ、単位時間、又は、単位質量を、前記任意の高分子鎖の長さ、時間、又は、質量にそれぞれ換算する請求項１に記載の高分子材料のシミュレーション方法。
【請求項３】
前記粗視化モデルの前記長さのパラメータｂcは、下記式（１）で定義される請求項２に記載の高分子材料のシミュレーション方法。
【数１】
2014206913.tif 000012
ここで、
＜Ｒc2＞：粗視化モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌc：粗視化モデルの全長
Ａc：ビーズ間の平衡長
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
【請求項４】
前記粗視化モデルの前記時間のパラメータτcは、前記粗視化モデルの末端間ベクトルの時間相関関数が、１／ｅに減少するのに要する時間である請求項２又は３に記載の高分子材料のシミュレーション方法。
【請求項５】
前記粗視化モデルの前記摩擦のパラメータζcは、下記式（２）で定義される請求項２乃至４のいずれかに記載の高分子材料のシミュレーション方法。
【数２】
2014206913.tif 000013
ここで、
ｂc：粗視化モデルの長さのパラメータ
τc：粗視化モデルの時間のパラメータ
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
ｋB：ボルツマン定数
Ｔ：絶対温度
【請求項６】
前記全原子モデルの前記長さのパラメータｂfは、下記式（３）で定義される請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法。
【数３】
2014206913.tif 000014
ここで、
＜Ｒf2＞：全原子モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌf：全原子モデルの全長
【請求項７】
前記時間のパラメータτfは、前記全原子モデルの末端間ベクトルの時間相関関数が、１／ｅに減少するのに要する時間である請求項２乃至６のいずれかに記載の高分子材料のシミュレーション方法。
【請求項８】
前記摩擦のパラメータζfは、下記式（４）で定義される請求項２乃至７のいずれかに記載の高分子材料のシミュレーション方法。
【数４】
2014206913.tif 000015
ここで、
ｂf：全原子モデルの長さのパラメータ
τf：全原子モデルの時間のパラメータ
Ｎf：１つの全原子モデルあたりのRouseビーズの個数
Ｌf：全原子モデルの全長
ｋB：ボルツマン定数
Ｔ：絶対温度
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算しうる高分子材料のシミュレーション方法に関する。
【背景技術】
【０００２】
近年、ゴム配合の開発のために、高分子材料の反応を、コンピュータを用いて評価するためのシミュレーション方法（数値計算）が種々提案されている。この種のシミュレーション方法では、例えば、高分子鎖を、複数のビーズでモデル化した粗視化モデルが、コンピュータに入力される。そして、粗視化モデルが仮想空間に配置され、コンピュータによる分子動力学計算が行われる。関連する技術としては、次のものがある。
【先行技術文献】
【特許文献】
【０００３】
【非特許文献１】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
【発明の概要】
【発明が解決しようとする課題】
【０００４】
ところで、粗視化モデルを用いた分子動力学計算では、現実の高分子鎖とは異なる単位系（例えば、時間や長さ）が用いられている。このため、粗視化モデルを用いたシミュレーション結果から、現実の高分子鎖での時間や長さを求めるには、高分子材料を用いた実験により、高分子鎖の緩和時間等の物理量を求めて、分子動力学計算での単位系を、高分子鎖の単位系に換算していた。
【０００５】
しかしながら、上記のような方法では、高分子材料を用いた実験に要するコストや時間が増大するため、分子動力学計算での単位系を、高分子鎖の単位系に容易に換算できないという問題があった。また、従来の方法では、高分子材料を実際に準備する必要があるため、例えば、現実に存在しない高分子材料では換算できないという問題もあった。
【０００６】
また、全原子モデルを用いた分子動力学計算により、高分子鎖の単位系を求めることも考えられるが、全原子モデルは、複雑に屈曲して形成されるため、単位系の換算に必要となる全原子モデルの全長Ｌfを求めることが難しいという問題があった。
【０００７】
本発明は、以上のような実状に鑑み案出されたもので、粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算しうる高分子材料のシミュレーション方法を提供することを主たる目的としている。
【課題を解決するための手段】
【０００８】
本発明のうち請求項１記載の発明は、コンピュータに、高分子鎖を、複数のビーズでモデル化した粗視化モデルを入力する工程、前記コンピュータが、予め定めた仮想空間に配置された前記粗視化モデルを用いて分子動力学計算を行う工程、前記コンピュータが、前記粗視化モデルのRouseのパラメータを求める第１計算工程、前記コンピュータに、任意の高分子鎖を、複数の炭素を含む原子でモデル化した全原子モデルを入力する工程、前記コンピュータが、予め定めた仮想空間に配置された前記全原子モデルを用いて分子動力学計算を行う工程、前記コンピュータが、前記全原子モデルのRouseのパラメータを求める第２計算工程、及び、前記コンピュータが、前記粗視化モデルのRouseのパラメータと、前記全原子モデルのRouseのパラメータとに基づいて、前記粗視化モデルを用いた分子動力学計算での単位系を、前記任意の高分子鎖の単位系に換算する換算工程を含み、前記第２計算工程は、前記全原子モデルの全長Ｌfを求める全長計算工程を含み、前記全長計算工程は、前記コンピュータが、前記全原子モデルの前記炭素原子の一端の固定する工程と、  前記コンピュータが、分子動力学計算を行いながら、前記全原子モデルの前記炭素原子の他端を、前記一端から離れる方向へ強制移動させる工程と、前記コンピュータが、前記分子動力学計算によって強制的に引き伸ばされた前記全原子モデルを分子力学法による計算を行うことにより、前記全原子モデルの構造を安定化する工程と、前記コンピュータが、前記全原子モデルの前記一端と前記他端との間の距離から前記全原子モデルの全長Ｌfを求める工程とを含むことを特徴とする。
【０００９】
また、請求項２記載の発明は、前記粗視化モデルのRouseのパラメータは、長さのパラメータｂcと、時間のパラメータτcと、摩擦のパラメータζcとを含み、前記全原子モデルのRouseのパラメータは、長さのパラメータｂfと、時間のパラメータτfと、摩擦のパラメータζfとを含み、前記換算工程は、前記粗視化モデルを用いた分子動力学計算での単位長さ、単位時間、又は、単位質量を、前記任意の高分子鎖の長さ、時間、又は、質量にそれぞれ換算する請求項１に記載の高分子材料のシミュレーション方法である。
【００１０】
また、請求項３記載の発明は、前記粗視化モデルの前記長さのパラメータｂcは、下記式（１）で定義される請求項２に記載の高分子材料のシミュレーション方法である。
【数１】
2014206913.tif 000003
ここで、
＜Ｒc2＞：粗視化モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌc：粗視化モデルの全長
Ａc：ビーズ間の平衡長
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
【００１１】
また、請求項４記載の発明は、前記粗視化モデルの前記時間のパラメータτcは、前記粗視化モデルの末端間ベクトルの時間相関関数が、１／ｅに減少するのに要する時間である請求項２又は３に記載の高分子材料のシミュレーション方法である。
【００１２】
また、請求項５記載の発明は、前記粗視化モデルの前記摩擦のパラメータζcは、下記式（２）で定義される請求項２乃至４のいずれかに記載の高分子材料のシミュレーション方法である。
【数２】
2014206913.tif 000004
ここで、
ｂc：粗視化モデルの長さのパラメータ
τc：粗視化モデルの時間のパラメータ
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
ｋB：ボルツマン定数
Ｔ：絶対温度
【００１３】
また、請求項６記載の発明は、前記全原子モデルの前記長さのパラメータｂfは、下記式（３）で定義される請求項２乃至５のいずれかに記載の高分子材料のシミュレーション方法である。
【数３】
2014206913.tif 000005
ここで、
＜Ｒf2＞：全原子モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌf：全原子モデルの全長
【００１４】
また、請求項７記載の発明は、前記時間のパラメータτfは、前記全原子モデルの末端間ベクトルの時間相関関数が、１／ｅに減少するのに要する時間である請求項２乃至６のいずれかに記載の高分子材料のシミュレーション方法である。
【００１５】
また、請求項８記載の発明は、前記摩擦のパラメータζfは、下記式（４）で定義される請求項２乃至７のいずれかに記載の高分子材料のシミュレーション方法である。
【数４】
2014206913.tif 000006
ここで、
ｂf：全原子モデルの長さのパラメータ
τf：全原子モデルの時間のパラメータ
Ｎf：１つの全原子モデルあたりのRouseビーズの個数
Ｌf：全原子モデルの全長
ｋB：ボルツマン定数
Ｔ：絶対温度
【発明の効果】
【００１６】
本発明の高分子材料のシミュレーション方法は、コンピュータに、高分子鎖を、複数のビーズでモデル化した粗視化モデルを入力する工程、コンピュータが、予め定めた仮想空間に配置された粗視化モデルを用いて分子動力学計算を行う工程、及び、コンピュータが、粗視化モデルのRouseのパラメータを求める第１計算工程を含む。
【００１７】
また、本発明の高分子材料のシミュレーション方法では、コンピュータに、任意の高分子鎖を、複数の炭素原子でモデル化した全原子モデルを入力する工程、コンピュータが、予め定めた仮想空間に配置された全原子モデルを用いて分子動力学計算を行う工程、及び、コンピュータが、全原子モデルのRouseのパラメータを求める第２計算工程を含む。この全原子モデルのRouseのパラメータは、高分子鎖の単位系に一致する。
【００１８】
さらに、本発明の高分子材料のシミュレーション方法では、コンピュータが、粗視化モデルのRouseのパラメータと、全原子モデルのRouseのパラメータとに基づいて、粗視化モデルを用いた分子動力学計算での単位系を、任意の高分子鎖の単位系に換算する換算工程を含む。
【００１９】
このように、本発明の高分子材料のシミュレーション方法は、例えば、高分子鎖の物理量を実験によって求めることなく、粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算することができる。
【００２０】
また、全原子モデルは、例えば、現実に存在しない任意の高分子鎖をモデル化することができる。このため、第２計算工程では、現実に存在しない高分子鎖をモデル化した全原子モデルのRouseのパラメータを求めることができる。従って、本発明の高分子材料のシミュレーション方法では、粗視化モデルを用いた分子動力学計算での単位系を、現実に存在しない高分子鎖の単位系に換算することができる。
【００２１】
さらに、第２計算工程は、全原子モデルの全長Ｌfを求める全長計算工程を含んでいる。全長計算工程は、コンピュータが、全原子モデルの炭素原子の一端の固定する工程と、分子動力学計算を行いながら、全原子モデルの炭素原子の他端を、一端から離れる方向へ強制移動させる工程と、分子動力学計算によって強制的に引き伸ばされた全原子モデルを分子力学法による計算を行うことにより、全原子モデルの構造を安定化する工程と、全原子モデルの一端と他端との間の距離から全原子モデルの全長Ｌfを求める工程とを含んでいる。このような全長計算工程は、全原子モデルを鎖方向に伸ばしつつ、結合長、結合角、及び、二面角が安定した構造を有する全原子モデルを確実に得ることができるため、全原子モデルの全長Ｌfを効率よく計算することができる。従って、本発明の高分子材料のシミュレーション方法では、粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算しうる。
【図面の簡単な説明】
【００２２】
【図１】本実施形態の高分子材料のシミュレーション方法を実行するコンピュータの斜視図である。
【図２】本実施形態のシミュレーション方法のフローチャートである。
【図３】粗視化モデルの概念図である。
【図４】本実施形態の第１シミュレーション工程のフローチャートである。
【図５】粗視化モデルのポテンシャルを説明する概念図である。
【図６】粗視化モデルが配置された仮想空間の概念斜視図である。
【図７】本実施形態の第１計算工程のフローチャートである。
【図８】粗視化モデルの末端間ベクトルの時間相関関数と、時間との関係を示すグラフである。
【図９】全原子モデルの概念図である。
【図１０】本実施形態の第２シミュレーション工程のフローチャートである。
【図１１】全原子モデルのポテンシャルを説明する概念図である。
【図１２】全原子モデルが配置された仮想空間の概念斜視図である。
【図１３】本実施形態の第２計算工程のフローチャートである。
【図１４】本実施形態の全長計算工程のフローチャートである。
【図１５】一端が固定された全原子モデルの概念図である。
【図１６】他端が強制移動された全原子モデルの概念図である。
【図１７】分子力学法による計算後の全原子モデルの概念図である。
【図１８】全原子モデルの末端間ベクトルの時間相関関数と、時間との関係を示すグラフである。
【発明を実施するための形態】
【００２３】
以下、本発明の実施の一形態が図面に基づき説明される。
図１に示されるように、本実施形態の高分子材料のシミュレーション方法（以下、単に「シミュレーション方法」ということがある）は、高分子材料の反応を、コンピュータを用いて評価するためのものである。高分子材料としては、例えば、ゴム、樹脂又はエラストマー等が含まれる。
【００２４】
図１に示されるように、コンピュータ１は、本体１ａ、キーボード１ｂ、マウス１ｃ及びディスプレイ装置１ｄを含む。この本体１ａには、演算処理装置（ＣＰＵ）、ＲＯＭ、作業用メモリー、磁気ディスクなどの記憶装置及びディスクドライブ装置１ａ１、１ａ２などが設けられる。なお、記憶装置には、本実施形態のシミュレーション方法を実行するための処理手順（プログラム）が予め記憶される。
【００２５】
図２には、本実施形態のシミュレーション方法の具体的な処理手順が示されている。このシミュレーション方法では、先ず、図３に示されるように、高分子鎖を複数のビーズ３でモデル化した粗視化モデル６を、コンピュータ１に入力する（工程Ｓ１）。この粗視化モデル６は、高分子材料を分子動力学で取り扱うための数値データである。これらの数値データは、コンピュータ１に入力される。
【００２６】
ビーズ３は、分子動力学計算において、運動方程式の質点として取り扱われる。即ち、ビーズ３には、質量、直径、電荷又は初期座標などのパラメータが定義される。これらの各パラメータは、数値情報としてコンピュータ１に記憶される。
【００２７】
また、ビーズ３、３間には、平衡長が定義された結合ポテンシャルが設定される。これにより、ビーズ３、３間には、該ビーズ３、３を拘束する結合鎖４が定義される。
【００２８】
ここで、「平衡長」とは、ビーズ３、３間の結合距離である。この結合距離が変化した場合は、結合鎖４によって、元の平衡長が保持される。これにより、粗視化モデル６は、直鎖状の三次元構造を維持することができる。なお、平衡長は、隣り合うビーズ３、３の中心３ａ、３ａ間の距離として定義される。また、結合鎖４の結合ポテンシャルには、例えば、上記非特許文献１に基づいて設定されるのが望ましい。
【００２９】
次に、コンピュータ１が、予め定めた仮想空間に配置された粗視化モデル６を用いて分子動力学計算を行う（第１シミュレーション工程Ｓ２）。図４には、本実施形態の第１シミュレーション工程Ｓ２の処理手順の一例が示されている。
【００３０】
本実施形態の第１シミュレーション工程Ｓ２では、先ず、図５に示されるように、隣り合う粗視化モデル６、６のビーズ３、３間に、ビーズ３、３間の距離の関数である相互ポテンシャルＰ１をそれぞれ定義する（工程Ｓ２１）。本実施形態の相互ポテンシャルＰ１は、ＬＪ（Lennard-Jones）ポテンシャルであり、下記式（５）で定義される。
【００３１】
【数５】
2014206913.tif 000007
ここで、各定数及び変数は、次のとおりである。
ｒij：各ビーズ間距離
ｒc：カットオフ距離
ε：各ビーズ間に定義される相互ポテンシャルの強度
σ：各ビーズ間に定義される相互ポテンシャルの距離に関するパラメータ（分子動力学の分野では、ＬＪ球の直径と呼ばれる）
なお、ビーズ間の距離ｒij及びカットオフ距離ｒcは、各ビーズ３の中心３ａ（図３に示す）を基準として定義される。
【００３２】
相互ポテンシャルＰ１は、ビーズ３、３間距離ｒijが予め定められたカットオフ距離ｒc未満になる場合にのみ、ビーズ３、３間に引力又は斥力を生じさせる。なお、ビーズ間距離ｒijがカットオフ距離ｒc以上になった場合には、相互ポテンシャルＰ１がゼロになり、ビーズ３、３間に引力又は斥力が生じない。このような相互ポテンシャルＰ１は、分子動力学計算において、実際の高分子材料の分子運動に近似させることができる。なお、相互ポテンシャルＰ１の強度ε、相互ポテンシャルＰ１が作用する距離σ、及び、カットオフ距離ｒcは、例えば、論文（上記の非特許文献１）に基づいて、下記のように設定されるのが望ましい。このような相互ポテンシャルＰ１は、コンピュータ１に記憶される。
強度ε：１．０
距離σ：１．０
カットオフ距離ｒc：２1/6
【００３３】
次に、図６に示されるように、予め定められた体積を持つ仮想空間８内に、粗視化モデル６を配置する（工程Ｓ２２）。この仮想空間８は、解析対象の高分子材料の微小構造部分に相当する。
【００３４】
本実施形態の仮想空間８は、１辺の長さＬｓが、例えば非特許文献１に基づいて、系内のビーズの数密度０．８５と、仮想空間８内に配置するビーズ３の個数とに基づいて定義される。また、仮想空間８には、例えば、粗視化モデル６が１００〜４００本程度配置される。なお、粗視化モデル６は、例えば、モンテカルロ法に基づいて、仮想空間８内に配置されるのが望ましい。そして、粗視化モデル６のビーズ３の座標等が、コンピュータ１に記憶される。
【００３５】
次に、コンピュータ１が、粗視化モデル６の分子動力学計算を行う（工程Ｓ２３）。本実施形態の分子動力学計算では、例えば、仮想空間８について所定の時間、配置した全ての粗視化モデル６が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全てのビーズ３の動きが追跡され、コンピュータ１に記憶される。また、分子動力学計算の条件は、例えば、系内のビーズ３の個数、体積及び温度は一定で行われる。
【００３６】
次に、コンピュータ１が、粗視化モデル６の初期配置を十分に緩和できたか否かを判断する（工程Ｓ２４）。この工程Ｓ２４では、粗視化モデル６の初期配置を十分に緩和できたと判断された場合、次の第１計算工程Ｓ３が実施される。一方、粗視化モデル６の初期配置を十分に緩和できていないと判断された場合は、単位時間を進めて（工程Ｓ２５）、工程Ｓ２３及び工程Ｓ２４が再度実施される。これにより、第１シミュレーション工程Ｓ２では、粗視化モデル６の平衡状態（構造が緩和した状態）を確実に計算することができる。
【００３７】
次に、コンピュータ１が、粗視化モデル６のRouseのパラメータを求める（第１計算工程Ｓ３）。ここで、Rouseのパラメータとは、高分子物理で知られているRouseモデルのパラメータである。本実施形態では、長さのパラメータｂc（単位：σ）、時間のパラメータτc（単位：τ）、及び、摩擦のパラメータζc（単位：μ／τ）が求められ、コンピュータ１に記憶される。また、これらのパラメータｂc、τc、ζcは、工程Ｓ２３で実施された粗視化モデル６の分子動力学計算において、単位時間毎に記憶された数値データ等を用いて計算される。図７には、本実施形態の第１計算工程Ｓ３の処理手順の一例が示されている。
【００３８】
本実施形態の第１計算工程Ｓ３では、先ず、コンピュータ１が、粗視化モデル６の長さのパラメータｂc（単位：σ）を求める（工程Ｓ３１）。この長さのパラメータｂcは、下記式（１）で定義される。
【数１】
2014206913.tif 000008
ここで、
＜Ｒc2＞：粗視化モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌc：粗視化モデルの全長
Ａc：ビーズ間の平衡長
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
【００３９】
長さＲcは、図３に示すように、粗視化モデル６の末端６ｅ、６ｅ間を結ぶ末端間ベクトルＶ１の長さである。この長さＲcも、ビーズ３の中心３ａを基準に定義される。また、＜Ｒc2＞の計算方法としては、先ず、各粗視化モデル６において、単位時間毎の長さＲcの二乗について、分子動力学計算の全時間の平均をとる（時間平均）。そして、各粗視化モデル６の長さＲcの二乗の時間平均を、全ての粗視化モデル６について平均（アンサンブル平均）をとることにより、＜Ｒc2＞を求めることができる。
【００４０】
また、全長Ｌcは、１つの粗視化モデル６を構成する全ての結合鎖４の合計長さである。このような全長Ｌcは、ビーズ３、３間の平衡長Ａcに、結合鎖４（図３に示す）の個数（Ｎc−１）を乗じることにより求めることができる。なお、平衡長Ａcは、結合鎖４に定義される結合ポテンシャルの極小値が設定される。そして、上記式（１）で求められた粗視化モデル６の長さのパラメータｂcが、コンピュータ１に記憶される。
【００４１】
次に、コンピュータ１が、粗視化モデル６の時間のパラメータτc（単位：τ）を求める（工程Ｓ３２）。図８は、粗視化モデル６の末端間ベクトルＶ１の長さＲcの時間相関関数Ｃcと、時間ｔとの関係を示すグラフである。時間のパラメータτcは、末端間ベクトルＶ１の相関が無くなる時間、即ち、時間相関関数Ｃcが１／ｅに減少するのに要する時間（時定数）で定義することができる。これは、原子論の崩壊に基づくものである。なお、時間のパラメータτcは、時間相関関数Ｃcを、下記式（６）に近似させることにより、容易に求めることができる。そして、粗視化モデル６の時間のパラメータτcが、コンピュータ１に記憶される。
Ｃc＝Ａ×exp（−ｔ／τc）…（６）
ここで、
Ｃc：粗視化モデルの末端間ベクトルの時間相関関数
Ａ：フィッティングパラメータ
ｔ：時間
τc：粗視化モデルの時間のパラメータ
【００４２】
次に、コンピュータ１が、粗視化モデル６の摩擦のパラメータζc（単位：μ／τ）を求める（工程Ｓ３３）。この摩擦のパラメータζcは、下記式（２）で定義される。
【数２】
2014206913.tif 000009
ここで、
ｂc：粗視化モデルの長さのパラメータ
τc：粗視化モデルの時間のパラメータ
Ｎc：１つの粗視化モデルあたりのRouseビーズの個数
ｋB：ボルツマン定数
Ｔ：絶対温度
【００４３】
長さのパラメータｂc及び時間のパラメータτcは、工程Ｓ３１、及び、工程Ｓ３２で求められた数値である。ビーズの個数Ｎcは、粗視化モデル６のビーズ３の個数である。ただし、ビーズの個数Ｎcは、上記式（１）で求められる全長Ｌcを、長さのパラメータｂcで除することにより求めることができる。そして、上記式（２）で求められた粗視化モデル６の摩擦のパラメータζcが、コンピュータ１に記憶される。
【００４４】
このように、第１計算工程Ｓ３では、粗視化モデル６のRouseのパラメータ（長さのパラメータｂc、時間のパラメータτc、及び、摩擦のパラメータζc）を確実に求めることができる。
【００４５】
次に、コンピュータ１に、任意の高分子鎖をモデル化した全原子モデルを入力する（工程Ｓ４）。本実施形態の高分子鎖としては、例えば、cis-１，４ポリブタジエン（以下、単に「ポリブタジエン」ということがある）が例示される。図９に示されるように、本実施形態の全原子モデル１１は、炭素原子１２ｃ又は水素原子１２ｈを含む複数の粒子モデル１２でモデル化される。
【００４６】
粒子モデル１２は、分子動力学計算に基づいたシミュレーションにおいて、運動方程式の質点として取り扱われる。即ち、粒子モデル１２は、質量、直径、電荷、又は、初期座標などのパラメータが定義される。これらのパラメータは、コンピュータ１上で取り扱い可能な数値情報として、コンピュータ１に記憶される。
【００４７】
また、粒子モデル１２、１２間には、結合鎖１４が定義される。この結合鎖１４は、粒子モデル１２、１２間を拘束するものである。本実施形態の結合鎖１４は、炭素原子１２ｃ、１２ｃを連結する主鎖１４ａ、及び、炭素原子１２ｃと水素原子１２ｈとの間を連結する側鎖１４ｂとを含んでいる。これらの主鎖１４ａ及び側鎖１４ｂは、例えば、平衡長とバネ定数とが定義されたバネとして取り扱われる。これらの情報もコンピュータ１に記憶される。
【００４８】
全原子モデル１１は、各粒子モデル１２、１２間の結合長さである結合長、結合鎖１４を介して連続する３つの粒子モデル１２がなす角度である結合角、及び、結合鎖１４を介して連続する４つの粒子モデル１２において、隣り合う３つの粒子モデル１２が作る二面角などが定義される。これにより、全原子モデル１１は、３次元構造を有する。
【００４９】
また、全原子モデル１１は、慣例に従い、外力又は内力を受けることによって、結合長、結合角及び二面角が変化する。これにより、全原子モデル１１は、その三次元構造を変化させることができる。以上のようなモデル化は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡというソフトウエアを用いて処理することができる。
【００５０】
コンピュータ１が、予め定めた仮想空間に配置された全原子モデル１１を用いて分子動力学計算を行う（第２シミュレーション工程Ｓ５）。図１０には、本実施形態の第２シミュレーション工程Ｓ５の処理手順の一例が示されている。
【００５１】
本実施形態の第２シミュレーション工程Ｓ５では、先ず、図１１に示されるように、全原子モデル１１の各粒子モデル１２、１２間に相互ポテンシャルＰ２が定義される（工程Ｓ５１）。相互ポテンシャルＰ２は、ＬＪ（Lennard-Jones）ポテンシャルであり、上記式（５）で定義される。このような相互ポテンシャルＰ２は、コンピュータ１に入力される。
【００５２】
相互ポテンシャルＰ２は、炭素原子１２ｃ、１２ｃ間に設定される第１相互ポテンシャルＰ２ａ、水素原子１２ｈ、１２ｈ間に設定される第２相互ポテンシャルＰ２ｂ、及び、炭素原子１２ｃと水素原子１２ｈと間に設定される第３相互ポテンシャルＰ２ｃが含まれる。これらのポテンシャルＰ２ａ〜Ｐ２ｃにおいて、上記式（５）の距離σ、及び、カットオフ距離ｒcは、適宜設定することができる。
【００５３】
次に、図１２に示されるように、予め定められた体積を持つ仮想空間１６内に、全原子モデル１１が配置される（工程Ｓ５２）。この仮想空間１６は、解析対象の高分子材料の微小構造部分に相当し、図６に示した仮想空間８と同様に設定される。また、仮想空間１６には、例えば、全原子モデル１１が１０〜１００本程度配置される。なお、全原子モデル１１は、モンテカルロ法に基づいて、仮想空間１６内に配置されるのが望ましい。そして、全原子モデル１１の粒子モデル１２の座標等が、コンピュータ１に記憶される。
【００５４】
次に、コンピュータ１が、全原子モデル１１の分子動力学計算を行う（工程Ｓ５３）。本実施形態の分子動力学計算では、粗視化モデル６の分子動力学計算と同様に、例えば、仮想空間１６について所定の時間、配置した全ての全原子モデル１１が古典力学に従うものとして、ニュートンの運動方程式が適用される。そして、各時刻での全ての粒子モデル１２（図９に示す）の動きが追跡される。また、分子動力学計算の計算条件は、例えば系内の粒子、体積及び温度は一定で行われる。
【００５５】
次に、コンピュータ１が、全原子モデル１１の初期配置を十分に緩和できたか否かを判断する（工程Ｓ５４）。この工程Ｓ５４では、全原子モデル１１の初期配置を十分に緩和できたと判断された場合、次の第２計算工程Ｓ６が実施される。一方、全原子モデル１１の初期配置を十分に緩和できていないと判断された場合は、単位時間を進めて（工程Ｓ５５）、工程Ｓ５３（分子動力学計算）及び工程Ｓ５４が再度実施される。これにより、第２シミュレーション工程Ｓ５では、全原子モデル１１の平衡状態（構造が緩和した状態）を確実に計算することができる。
【００５６】
次に、コンピュータ１が、全原子モデルのRouseのパラメータを求める（第２計算工程Ｓ６）。本実施形態では、Rouseのパラメータのうち、長さのパラメータｂf（単位：ｍ）、時間のパラメータτf（単位：ｓ）、及び、摩擦のパラメータζf（単位：ｋｇ／ｓ）が求められる。これらのパラメータｂf、τf、ζfは、第２シミュレーション工程Ｓ５で実施された全原子モデル１１の分子動力学計算において、単位時間毎に記憶された数値データ等を用いて計算される。図１３には、本実施形態の第２計算工程Ｓ６の処理手順の一例が示されている。
【００５７】
本実施形態の第２計算工程Ｓ６では、先ず、コンピュータ１が、全原子モデル１１の全長Ｌfを求める（全長計算工程Ｓ６１）。図１４には、全長計算工程Ｓ６１の処理手順の一例が示されている。
【００５８】
本実施形態の全長計算工程Ｓ６１では、先ず、図１５に示されるように、１つの全原子モデル１１を、予め定められた体積を持つ仮想空間１８に配置する（工程Ｓ６１１）。この仮想空間１８としては、図１２に示した仮想空間１６と同一構成のものが用いられる。そして、全原子モデル１１の粒子モデル１２の座標等が、コンピュータ１に記憶される。
【００５９】
次に、コンピュータ１が、全原子モデル１１の炭素原子１２ｃの一端の固定する（工程Ｓ６１２）。この工程Ｓ６１２では、全原子モデル１１の末端１１ａ、１１ｂを構成する炭素原子１２ｃのうち、一端１１ａの炭素原子１２ｃの座標値が、変動不能に固定される。これにより、全原子モデル１１の炭素原子１２ｃの一端１１ａを、仮想空間１８において固定することができる。このような一端１１ａの炭素原子１２ｃの座標値は、コンピュータ１に記憶される。
【００６０】
次に、コンピュータ１が、分子動力学計算を行いながら、全原子モデル１１の他端１１ｂを、一端１１ａから離れる方向へ強制移動させる（工程６１３）。本実施形態の工程Ｓ６１３では、先ず、全原子モデル１１の末端１１ａ、１１ｂを構成する炭素原子１２ｃのうち、他端１１ｂの炭素原子１２ｃに、一端１１ａの炭素原子１２ｃから離れる方向へ移動する等速度運動Ｍが定義される。この等速度運動Ｍの速度としては、例えば、５００〜３５００ｍ／ｓが設定される。次に、全原子モデル１１の結合鎖１４に定義されるポテンシャル（図示省略）に基づいて、分子動力学計算を開始する。そして、全原子モデル１１の他端１１ｂに定義された等速度運動Ｍに基づいて、他端１１ｂの移動計算が実施される。
【００６１】
これにより、図１６に示されるように、全原子モデル１１は、各炭素原子１２ｃが略直線上に配置され、主鎖方向に強制的に引き伸ばされた構造が設定される。このように、本実施形態の工程Ｓ６１３では、分子動力学計算を行いながら、炭素原子１２ｃの移動計算が実施されるため、各炭素原子１２ｃの移動を円滑に計算することができる。なお、分子動力学計算のステップ数は、適宜設定することができるが、例えば、単位時間が１ｆｓ（フェムト秒）とすると、２５００〜７５００ステップが望ましい。
【００６２】
次に、コンピュータ１が、強制的に引き伸ばされた全原子モデル１１に対して、分子力学法による計算を行う（工程Ｓ６１４）。分子力学法（ＭＭ法）は、原子間に働く相互作用のポテンシャルエネルギーに基づいて、分子構造を最適化する方法である。このような分子力学法による計算は、例えば（株）ＪＳＯＬ社製のＪ−ＯＣＴＡに含まれるＣＯＧＮＡＣを用いて処理することができる。
【００６３】
これにより、図１７に示されるように、工程Ｓ６１４では、主鎖方向に伸ばしつつ、結合長、結合角、及び、二面角が安定した全原子モデル１１を得ることができる。なお、分子力学法による計算は、結合鎖１４の主鎖１４ａに定義されるポテンシャルが、例えば、６×１０−11〜８×１０−11Ｎに収束するまで行われるのが望ましい。これにより、上記のような全原子モデル１１を確実に得ることができる。
【００６４】
次に、コンピュータ１が、全原子モデル１１の一端１１ａと他端１１ｂとの間の距離を求める（工程Ｓ６１５）。これにより、全原子モデル１１の全長Ｌfを求めることができる。このように、本実施形態の全長計算工程Ｓ６１では、全原子モデル１１を主鎖方向に伸ばしつつ、結合長、結合角、及び、二面角が安定した構造を有する全原子モデル１１を確実に得ることができるため、全原子モデル１１の全長Ｌfを効率よく確実に計算することができる。そして、全原子モデル１１の全長Ｌfが、コンピュータ１に記憶される。
【００６５】
次に、コンピュータ１が、全原子モデル１１の長さのパラメータｂfを求める（工程Ｓ６２）。この長さのパラメータｂfは、下記式（３）で定義される。
【数３】
2014206913.tif 000010
ここで、
＜Ｒf2＞：全原子モデルの末端間ベクトルの長さの２乗のアンサンブル平均
Ｌf：全原子モデルの全長
【００６６】
図９に示されるように、長さＲfは、全原子モデル１１の末端１１ａ、１１ｂ間を結ぶ末端間ベクトルＶ２の長さである。この長さＲfは、炭素原子１２ｃの中心（図示省略）を基準に定義される。
【００６７】
また、＜Ｒf2＞の計算方法としては、先ず、各全原子モデル１１において、単位時間毎の長さＲfの二乗について、分子動力学計算の全時間の平均をとる（時間平均）。そして、各全原子モデル１１のＲf2の時間平均を、全ての全原子モデル１１について平均をとる（アンサンブル平均）ことにより、＜Ｒc2＞を求めることができる。なお、全長Ｌfは、上述した全長計算工程Ｓ６１で求められた数値である。そして、全原子モデル１１の長さのパラメータｂfが、コンピュータ１に記憶される。
【００６８】
次に、コンピュータ１が、全原子モデル１１の時間のパラメータτfを求める（工程Ｓ６３）。図１８は、全原子モデル１１の末端間ベクトルＶ２の長さＲfの時間相関関数Ｃfと、時間ｔとの関係を示すグラフである。時間のパラメータτfは、粗視化モデル６の時間のパラメータτcと同様に、時間相関関数Ｃfが１／ｅに減少するのに要する時間（時定数）で定義される。なお、時間のパラメータτfは、時間相関関数Ｃfを、下記式（７）に近似させることにより、容易に求めることができる。そして、全原子モデル１１の時間のパラメータτfが、コンピュータ１に記憶される。
Ｃf＝Ａ×exp（−ｔ／τf）…（６）
ここで、
Ｃf：全原子モデルの末端間ベクトルの長さの時間相関関数
Ａ：フィッティングパラメータ
ｔ：時間
τf：粗視化モデルの時間のパラメータ
【００６９】
次に、コンピュータ１が、全原子モデル１１の摩擦のパラメータζfを求める（工程Ｓ６４）。この摩擦のパラメータζfは、下記式（４）で定義される。
【数４】
2014206913.tif 000011
ここで、
ｂf：全原子モデルの長さのパラメータ
τf：全原子モデルの時間のパラメータ
Ｎf：１つの全原子モデルあたりのRouseビーズの個数
Ｌf：全原子モデルの全長
ｋB：ボルツマン定数
Ｔ：絶対温度
【００７０】
長さのパラメータｂf及び時間のパラメータτfは、上記の方法で求められた数値である。また、炭素原子の個数Ｎfには、全原子モデル１１の炭素原子１２ｃの個数が代入されるが、例えば、全長計算工程Ｓ６１で求められた全長Ｌfを、長さのパラメータｂfで除することにより求めることができる。そして、全原子モデル１１の摩擦のパラメータζfが、コンピュータ１に記憶される。
【００７１】
このように、第２計算工程Ｓ６では、全原子モデル１１のRouseのパラメータ（長さのパラメータｂf、時間のパラメータτf、及び、摩擦のパラメータζf）を確実に求めることができる。
【００７２】
次に、コンピュータ１が、粗視化モデル６を用いた分子動力学計算での単位系を、高分子鎖の単位系に換算する（換算工程Ｓ７）。
【００７３】
粗視化モデル６の分子動力学計算での単位系としては、σ（長さ）、τ（時間）、及び、μ（質量）が含まれる。一方、高分子鎖の単位系としては、ｍ（長さ）、ｓ（時間）、及び、kg（質量）が含まれる。本実施形態では、粗視化モデルを用いた分子動力学計算での単位長さ（１σ）、単位時間（１τ）、及び、単位質量（１μ）を、高分子鎖の長さ（ｍ）、時間（ｍ）、又は、質量（kg）にそれぞれ換算する。
【００７４】
換算工程Ｓ７では、粗視化モデル６のRouseのパラメータと、全原子モデル１１のRouseのパラメータとが同一であると仮定して、粗視化モデル６を用いた分子動力学計算での単位系を、任意の高分子鎖の単位系に換算する。
【００７５】
例えば、粗視化モデル６の長さのパラメータｂcが１σであり、かつ、全原子モデル１１の長さのパラメータｂfが１×１０−9ｍであったとする。本実施形態では、各長さのパラメータｂc、ｂfが同一であると仮定するため、分子動力学計算での単位長さ１σを、高分子鎖の長さ１×１０−9ｍに換算することができる。
【００７６】
また、粗視化モデル６の時間のパラメータτcが５τであり、かつ、全原子モデル１１の時間のパラメータτfが５×１０−11ｓであったとする。本実施形態では、各時間のパラメータτc、τfが同一であると仮定するため、分子動力学計算での単位時間１τを、高分子鎖での時間１×１０−11ｓに換算することができる。
【００７７】
さらに、粗視化モデル６の摩擦のパラメータζcが５０μ／τであり、かつ、全原子モデル１１の摩擦のパラメータζfが１×１０−11kg／ｓであったとする。この場合には、分子動力学計算での単位質量１μを、高分子鎖での質量２．０×１０−24kgに換算することができる。
【００７８】
このように、本実施形態のシミュレーション方法では、例えば、従来のように高分子鎖の物理量を実験によって求めることなく、粗視化モデル６を用いた分子動力学計算での単位系を、高分子鎖の単位系に換算することができる。従って、シミュレーション方法では、実験によって増大しがちなコストや時間を抑制でき、粗視化モデル６の単位系を、高分子鎖の単位系に容易に換算することができる。また、全長計算工程Ｓ６１では、全原子モデル１１の全長Ｌfを効率よく確実に計算することができるため、高分子鎖の単位系にさらに容易に換算することができる。
【００７９】
また、本実施形態では、ポリブタジエンをモデル化した全原子モデル１１のRouseのパラメータを求めるものが例示されたが、これに限定されるわけではない。例えば、ポリブタジエン以外の任意の高分子鎖をモデル化した全原子モデル１１に基づいて、Rouseのパラメータを求めてもよい。
【００８０】
換算工程Ｓ７では、このような全原子モデル１１においても、粗視化モデル６のRouseのパラメータと、全原子モデル１１のRouseのパラメータとが同一であると仮定される。これにより、本発明のシミュレーション方法では、粗視化モデル６を用いた分子動力学計算を別途行うことなく、粗視化モデル６（図３に示す）の単位系を、ポリブタジエンの単位系に限定されることなく、任意の高分子鎖の単位系に換算することができる。
【００８１】
さらに、全原子モデル１１は、例えば、現実に存在しない任意の高分子材料をモデル化することができる。このため、第２計算工程Ｓ６では、現実に存在しない高分子鎖をモデル化した全原子モデル１１のRouseのパラメータを求めることができる。これにより、本発明のシミュレーション方法では、上記の手順に従うことによって、粗視化モデル６を用いた分子動力学計算での単位系を、現実に存在しない高分子鎖の単位系に換算することができる。従って、本発明のシミュレーション方法では、未知の高分子材料の開発に役立つ。
【００８２】
次に、平衡状態の粗視化モデル６を含む仮想空間８（図６に示す）を用いて、高分子材料の変形シミュレーションを行う（工程Ｓ８）。この工程Ｓ８では、例えば、変形速度０．０５σ／τとして、変形量が１．０σになるまで変形計算する。そして、仮想空間８の変形による物理量（応力等）が計算される。
【００８３】
次に、変形シミュレーションによる物理量を、高分子鎖の単位系に換算する（工程Ｓ９）。この工程Ｓ９では、高分子鎖の長さ、時間、又は、質量に換算された粗視化モデル６の単位長さ、単位時間、又は、単位質量に基づいて、仮想空間８の変形による物理量（応力等）を、高分子鎖の単位系に換算する。
【００８４】
このように、本実施形態のシミュレーション方法では、例えば、全原子モデル１１を含む仮想空間１６（図１２に示す）を用いたり、高分子材料を用いた実験を行ったりすることなく、高分子鎖の単位系で表された物理量を容易に求めることができるため、計算時間を大幅に短縮しうる。
【００８５】
また、本実施形態のシミュレーション方法では、換算工程Ｓ７において、粗視化モデル６の単位系を、複数の高分子鎖の単位系に換算しておくことにより、粗視化モデル６を含む仮想空間８を用いた変形シミュレーションを１回行うだけで、複数の高分子鎖の単位系で表された物理量を同時に求めることができる。このため、本実施形態のシミュレーション方法では、変形シミュレーションの計算時間を大幅に短縮することができる。
【００８６】
次に、変形シミュレーションによる物理量が、許容範囲内であるかが判断される（工程Ｓ１０）。この工程Ｓ１０では、物理量が許容範囲内であると判断された場合、粗視化モデル６を含む仮想空間８の条件、及び、物理量が換算された高分子鎖の構造に基づいて、高分子材料が製造される（工程Ｓ１１）。一方、物理量が許容範囲内ではないと判断された場合には、粗視化モデル６を含む仮想空間８の諸条件を変更して（工程Ｓ１２）、工程Ｓ８〜Ｓ１１が再度行われる。
【００８７】
このように、本実施形態のシミュレーション方法では、変形シミュレーションによる物理量が許容範囲内になるまで、粗視化モデル６を含む仮想空間８の条件が変更されるため、所望の性能を有する高分子材料を、効率よく設計することができる。
【００８８】
以上、本発明の特に好ましい実施形態について詳述したが、本発明は図示の実施形態に限定されることなく、種々の態様に変形して実施しうる。
【実施例】
【００８９】
図２に示される手順に従って、粗視化モデルのRouseのパラメータ、及び、下記高分子鎖をモデル化した全原子モデルのRouseのパラメータが求められ、粗視化モデルを用いた分子動力学計算での単位系が、任意の高分子鎖の単位系に換算された（実施例）。この実施例では、全長計算工程により、全原子モデル１１の全長Ｌfが求められた。
【００９０】
また、比較のために、高分子材料を用いた実験により、下記高分子鎖の緩和時間等の物理量が求められ、分子動力学計算での単位系が、高分子鎖の単位系に換算された（比較例）。
【００９１】
そして、実施例及び比較例において、換算に要するコスト、及び時間が測定された。なお、各ポテンシャルのパラメータ等は、明細書中の記載通りであり、共通仕様は次のとおりである。
高分子鎖：
構造：cis-１，４ポリブタジエン
重合度：２０７０
仮想空間：
１辺の長さＬｓ：１８σ
【００９２】
テストの結果、比較例の換算に要するコストは、実施例の換算に要するコストを１００とする指数（小さいほど良い）で表すと、１００００であった。さらに、比較例の換算に要する時間は、実施例の換算に要する時間を１００とする指数（小さいほど良い）で表すと、５０００であった。このように、実施例は、比較例に比べて、換算に要するコスト及び時間を大幅に低減でき、粗視化モデルを用いた分子動力学計算での単位系を、高分子鎖の単位系に容易に換算しうることを確認できた。
【符号の説明】
【００９３】
１        コンピュータ
６        粗視化モデル
１１        全原子モデル
【図１】
2014206913.tif 000016
【図２】
2014206913.tif 000017
【図３】
2014206913.tif 000018
【図４】
2014206913.tif 000019
【図５】
2014206913.tif 000020
【図６】
2014206913.tif 000021
【図７】
2014206913.tif 000022
【図８】
2014206913.tif 000023
【図９】
2014206913.tif 000024
【図１０】
2014206913.tif 000025
【図１１】
2014206913.tif 000026
【図１２】
2014206913.tif 000027
【図１３】
2014206913.tif 000028
【図１４】
2014206913.tif 000029
【図１５】
2014206913.tif 000030
【図１６】
2014206913.tif 000031
【図１７】
2014206913.tif 000032
【図１８】
2014206913.tif 000033