JP6651254B2 - シミュレーション方法、シミュレーションプログラム、及びシミュレーション装置 - Google Patents

Description

本発明は、くりこみ群を応用した分子動力学によるシミュレーション方法、シミュレーションプログラム、及びシミュレーション装置に関する。

分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算能力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。

下記の特許文献１に、必要な計算量を減らすためにくりこみ変換の手法を用いたシミュレーション方法が開示されている。以下、特許文献１に開示されたくりこみ変換の手法について説明する。

シミュレーション対象の粒子系Ｓの粒子数をＮ、各粒子の質量をｍ、粒子間の相互作用ポテンシャルをφ（ｒ）で表す。ここで、ｒは、粒子間の距離である。相互作用ポテンシャルφ（ｒ）が、相互作用係数εと関数ｆ（ｒ）との積で表される。相互作用係数εは、相互作用の強さを表し、エネルギの次元を持つ。関数ｆ（ｒ）は、粒子間距離に対する依存性を表し、無次元化されている。

第１のくりこみ因子α、第２のくりこみ因子γ、及び第３のくりこみ因子δが決定される。第１のくりこみ因子αは１より大きい。第２のくりこみ因子γは０以上であり、かつ空間の次元数ｄ以下である。第３のくりこみ因子δは０以上である。くりこみの回数をｎで表したとき、第１のくりこみ因子αは、
α＝２^ｎ
と表される。

くりこみ手法を用いてくりこみ変換された粒子系Ｓ’の粒子数をＮ’、各粒子の質量をｍ’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系Ｓ’の粒子数Ｎ’、質量ｍ’、及び相互作用係数ε’を求める。
Ｎ’＝Ｎ／α^ｄ
ｍ’＝ｍα^δ−γ
ε’＝εα^γ

くりこみ変換された粒子系Ｓ’について、分子動力学計算を行なう。分子動力学計算で求められた各粒子の位置ベクトルをｑ’、運動量ベクトルをｐ’で表す。元の粒子系Ｓの各粒子の位置ベクトルｑ、運動量ベクトルｐは、以下の式を用いて算出することができる。
ｑ＝ｑ’α
ｐ＝ｐ’／α^δ／２

特許第５２４１４６８号公報

粒子系Ｓが単原子の集合体である場合、特許文献１に開示されたくりこみ変換により、くりこまれた粒子系Ｓ’に対して分子動力学計算を行うことができる。多数のポリマー分子からなる粒子系Ｓに対して分子動力学計算を行う手法として、クレマーグレストモデル（Kremer Grest Model）が提案されている。クレマーグレストモデルにくりこみ変換処理を適用してシミュレーションを行う手法は確立されていない。

本発明の目的は、ポリマー分子からなる粒子系に対してくりこみ群の考え方を適用し、ポリマー分子の挙動をシミュレーションする方法を提供することである。本発明の他の目的は、このシミュレーション方法を実行するコンピュータプログラムを提供することである。本発明のさらに他の目的は、このシミュレーション方法を実行するシミュレーション装置を提供することである。

本発明の一観点によると、
くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
粒子間の相互作用ポテンシャルφ（ｒ）として、
任意のモノマー粒子間においては、

が適用され、ポテンシャルＵ_０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
次元を持つパラメータσと、無次元関数ｆを用いて

と表すことができ、
同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
Ｕ_０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵ_ｃｈ
（ｒ：ε，σ）がポテンシャルＵ_０（ｒ）に追加されて

が適用され、
粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
表したとき、下記の変換則

を適用し、
粒子系ＳにおけるポテンシャルＵ_０（ｒ）、Ｕ_ｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ
’におけるポテンシャルＵ_０’（ｒ）、Ｕ_ｃｈ’（ｒ：ε，σ）として、

を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法が提
供される。

本発明の他の観点によると、上記シミュレーションを実行するコンピュータプログラム、コンピュータプログラムを記録した記録媒体、及びシミュレーション装置が提供される。

複数のモノマー粒子が鎖状に繋がったポリマーからなる粒子系にくりこみ変換を適用して、分子動力学計算を行うことができる。これにより、計算時間を短縮することができる。

図１は、実施例によるシミュレーション方法のフローチャートである。 図２は、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている状態を示す模式図である。 図３は、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示すグラフである。 図４Ａ〜図４Ｃは、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明するための図である。 図５は、２φ（（ｑ_ｉ−ｑ_ｋ）／２）の計算結果、及びδφの数値積分結果を示すグラフである。 図６は、直線状のポリマーが３次元方向に規則正しく配置されている例を示す模式図である。 図７は、高分子材料の粘度のひずみ速度依存性をシミュレーションによって求めた結果を示すグラフである。 図８は、実施例によるシミュレーション装置のブロック図である。

本願発明の実施例で適用される分子動力学（Molecular Dynamics）について簡単に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが下記の式で表される粒子系について考える。

ここで、
ｍは粒子の質量を表し、
φは粒子間の相互作用ポテンシャルを表し、
ベクトルｐ_ｊは粒子の運動量ベクトルを表し、
ベクトルｑ_ｊは粒子の位置ベクトル（位置座標）を表す。

ハミルトニアンＨを、ハミルトニアンの正準方程式に代入することにより、粒子ｊに対する下記の運動方程式が得られる。

分子動力学では、粒子系を構成する各粒子に対し、式（７）及び式（８）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトルｐ_ｊ及び位置ベクトルｑ_ｊが求められる。数値積分には、Verlet法が用いられることが多い。Verlet法については、例えば、J.M.Thijssen,「Computational Physics」,Cambridge University Press (1999)の１７５頁に説明されている。分子動力学計算で得られた各粒子の運動量ベクトル及び位置ベクトルに基づいて、粒子系の種々の物理量を算出することができる。

次に、くりこみ群の手法を応用した分子動力学（以下、くりこみ群分子動力学という。）について概念的に説明する。

くりこみ群分子動力学では、シミュレーション対象の粒子系Ｓを、粒子系Ｓの粒子数よりも少ない数の粒子で構成された粒子系Ｓ’（以下、くりこまれた粒子系Ｓ’という。）に対応づける。次に、くりこまれた粒子系Ｓ’に対して分子動力学計算を実行する。くりこまれた粒子系Ｓ’に対する計算結果を、シミュレーション対象の粒子系Ｓに対応付ける。これにより、シミュレーション対象の粒子系Ｓに対して分子動力学計算を直接的に実行する場合よりも計算量を少なくすることができる。シミュレーション対象の粒子系Ｓにおける物理量（例えば、粒子数、粒子の質量等）と、くりこまれた粒子系Ｓ’における物理量とを相互に対応付けるための変換則を、くりこみ変換測という。

図１に、実施例によるシミュレーション方法のフローチャートを示す。ステップＳ１において、シミュレーション対象の粒子系Ｓに対して、くりこみ変換測に基づいてくりこみ変換処理を実行することにより、くりこまれた粒子系Ｓ’を定義する。粒子系Ｓは、複数のポリマーで構成され、複数のポリマーの各々は、一次元状に繋がった複数のモノマー粒子で構成される。

粒子系Ｓのモノマー粒子間の相互作用ポテンシャルφ（ｒ）として、任意のモノマー粒子間においては、

が適用される。ここで、ｒは、粒子系Ｓを構成する各粒子からの距離を表す。

ポテンシャルＵ_０（ｒ）は、基本的に下記の式で表される。

ここで、
ｆは無次元化された関数を表し、ε及びσは、モノマー粒子を特徴づけるパラメータである。パラメータεはエネルギの次元を持ち、相互作用係数と呼ばれる。パラメータσは距離の次元を持ち、粒子の大きさに依存する。

同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ_０（ｒ）に有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε，σ）が追加されて

が適用される。有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε，σ）は、ポテンシャルＵ_０（ｒ）を定義するパラメータε及びσに依存するパラメータを含む。

ポテンシャルＵ_０（ｒ）として、例えばレナード・ジョーンズ型のポテンシャルを適用することができる。ポテンシャルＵ_０（ｒ）に、シフトさせたレナード・ジョーンズ型のポテンシャルを適用した場合、ポテンシャルＵ_０（ｒ）は、例えば下記の式で定義することができる。

ポテンシャルＵ_０（ｒ）として、モース型ポテンシャルを適用することも可能である。さらに、ポテンシャルＵ_０（ｒ）として、これらの関数をテーラー展開した近似関数を用いることも可能である。本実施例においては、ポテンシャルＵ_０（ｒ）としてレナード・ジョーンズ型のポテンシャルを適用することとして説明を続ける。

くりこみ変換処理により、シミュレーション対象の粒子系Ｓのモノマー粒子の位置ベクトルｑ、モノマー粒子の質量ｍ、ポリマー内のモノマー粒子の個数Ｎｍ、ポリマーの数Ｎｐ、ポテンシャルＵ_０（ｒ）を定義するパラメータε、及びσが、それぞれ繰り込まれた粒子系Ｓ’のモノマー粒子の位置ベクトルｑ’、モノマー粒子の質量ｍ’、ポリマー内のモノマー粒子の個数Ｎｍ’、ポリマーの数Ｎｐ’、ポテンシャルＵ_０’（ｒ）を定義するパラメータε’、及びσ’に変換される。ステップＳ１のくりこみ変換処理では、下記のくりこみ変換測が適用される。

ここで、ｄはシミュレーション対象の粒子系Ｓが配置された空間の次元数を表す。λは繰り込み回数に依存するくりこみ因子である。くりこみ回数がｎ回のとき、くりこみ因子λは以下の式で表される。

粒子系ＳにおけるポテンシャルＵ_０（ｒ）及びＵ_ｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ’におけるポテンシャルＵ_０’（ｒ）及びＵ_ｃｈ’（ｒ：ε，σ）として、下記のポテンシャルを用いる。

粒子系Ｓ’の任意のモノマー粒子間に適用される相互作用ポテンシャルφ’（ｒ）は、下記の式で定義される。

同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ_０’（ｒ）に有限伸長非線形弾性ポテンシャルＵ_ｃｈ’（ｒ：ε，σ）が追加されて、相互作用ポテンシャルφ’として、

が適用される。

上述のくりこみ変換測を適用することにより、粒子の個数が１／λ^ｄ倍になり、粒子の質量がλ^ｄ倍になる。このため、粒子系Ｓの全体の質量と、粒子系Ｓ’の全体の質量とは同一である。さらに、粒子系Ｓ’におけるモノマー粒子の位置ベクトルｑ’が、粒子系Ｓにおけるモノマー粒子の位置ベクトルｑと同一であるため、粒子系Ｓとくりこまれた粒子系Ｓ’との巨視的な寸法は同一である。くりこみ変換処理の前後で、粒子系の全質量及び寸法が不変であるため、粒子系の密度も不変である。

次に、ステップＳ２において、シミュレーションの初期条件を設定する。初期条件には、各粒子の位置ベクトルｑ_ｊ及び運動量ベクトルｐ_ｊの初期値が含まれる。運動量ベクトルｐ_ｊは、くりこまれた粒子系Ｓ’の温度Ｔ’に基づいて設定される。シミュレーション対象の粒子系Ｓの温度をＴで表したとき、温度に関して下記のくりこみ変換測が適用される。

次に、ステップＳ３において、くりこまれた粒子系Ｓ’に関して分子動力学計算を実行する。運動方程式は、下記の式で表される。

ここで、Ｎ’は粒子系Ｓ’のモノマー粒子の全個数を表しており、Ｎ’＝Ｎｍ’×Ｎｐ’である。

この運動方程式を、数値的に積分して解く。これにより、くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’及び運動量ベクトルｐ’の時刻歴が求まる。

くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’と、シミュレーション対象の粒子
系Ｓの各粒子の位置ベクトルｑとの対応関係は、式（１３）に示したとおりである。くりこまれた粒子系Ｓ’における運動量ベクトルｐ’と、シミュレーション対象の粒子系Ｓにおける運動量ベクトルｐとは、下記の関係を有する。

ステップＳ４において、シミュレーション結果を出力する。例えば、位置ベクトルｑ’及び運動量ベクトルｐ’を、そのまま数値で出力してもよいし、位置ベクトルｑ’に基づいて、粒子系Ｓ’の複数の粒子の空間内の分布を画像として表示してもよい。

［φ（ｒ）＝Ｕ_０（ｒ）の場合について考察］
次に、相互作用ポテンシャルφ（ｒ）＝Ｕ_０（ｒ）の場合について、くりこみ変換則を導き出す。まず、くりこまれた粒子系Ｓ’のハミルトニアンＨ’が

となることを示す。くりこみ前の粒子系Ｓの粒子数をＮで表し、くりこまれた粒子系Ｓ’の粒子数をＮ’で表す。このハミルトニアンＨ’を正準方程式に代入して運動方程式を得ることができる。ハミルトニアンＨ’のパラメータから、式（１３）のくりこみ変換則が導出される。

くりこまれた粒子系Ｓ’のハミルトニアンＨ’を得るためには、粒子系Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することによりハミルトニアンＨ’が得られる。

粒子数が一定の正準集団（Canonical ensemble）に対する分配関数Ｚ（β）は、以下の式で表される。

ここで、ｋ_Ｂはボルツマン定数であり、ｄΓ_Ｎは位相空間内の体積要素である。ｄΓ_Ｎは下記の式で表される。

ここで、ｈはプランク定数である。Ｗ_Ｎは、すべての状態の量子論的な和と、位相空間に亘る積分とが一致するように決められる。

まず、粒子間の相互作用ポテンシャルの粗視化について説明し、次に運動エネルギの粗視化について説明する。続いて、相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化を踏まえて、くりこみ変換測を定義する。

［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ_０の粗視化］
まず、一次元鎖状に配置された粒子系における相互作用ポテンシャルの粗視化について説明する。その後、単純立方格子状に配置された粒子における相互作用ポテンシャルについて説明する。

図２に示すように、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている。粒子ｊが関与する相互作用を書き出し、粒子ｉと粒子ｋとの中間にある粒子ｊの位置座標について積分を実行することにより、相互作用ポテンシャルの粗視化を行うことができる。まず、第二近接以上の粒子からの寄与を取り込むために、ポテンシャルムービング（Potential Moving）の方法を用いることができる。ポテンシャルムービングの方法については、Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999)に解説されている。

第二近接以上の粒子からの寄与を取り込んだ相互作用ポテンシャルφティルダは、下記の式で表すことができる。

ここで、ａは平衡状態における粒子間距離を表す。平衡状態における粒子間距離ａは、相互作用ポテンシャルφ（ｒ）が最小となる距離に等しい。

複数の粒子が一次元状に配置されているため、粒子ｊの位置ベクトルｑ_ｊは、一次元座標ｑ_ｊで表すことができる。粒子ｊの位置をｑ_ｊで表すと、粒子ｊに対して、最近接の粒子によって作られるケージポテンシャルは、下記の式で表される。

積分変数ｑ_ｊについて積分を実行すれば、式（２５）を用いて下記の式が得られる。

ここで、ｒ_ａは粒子の直径を表し、ｚ（ｑ_ｉ−ｑ_ｋ）及びＰ（ｑ_ｉ−ｑ_ｋ）は、下記の式で表される。

積分領域はケージポテンシャルの内部領域に制限される。

次に、ｚ（ｑ_ｉ−ｑ_ｋ）を具体的に計算する。相互作用ポテンシャルφがレナード・ジョーンズ型ポテンシャルの場合、φ^（２ｎ）は、下記の式で表される。

式（２９）を式（２８）に代入して数値積分を行う。式（２８）に式（２９）を代入する際には、式（２４）を用いる。

図３に、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示す。粒子ｉ、粒子ｊ、粒子ｋ、粒子ｌ、及び粒子ｍがこの順番に一次元状に配置されている場合の粒子ｋの位置座標がｑ_ｋで表される。図３の横軸は、ｑ_ｋ／２を表し、縦軸は、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）、及びｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）を、対数目盛で表す。図３に至る数値計算において、ε／ｋ_ＢＴ＝２．０とした。

図３から、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）の変化に比べて、ｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）の変化が緩やかであることがわかる。このため、Ｐ（ｑ_ｉ−ｑ_ｋ）Ｐ（ｑ_ｋ−ｑ_ｍ）に対してｚ（ｑ_ｉ−ｑ_ｋ）ｚ（ｑ_ｋ−ｑ_ｍ）は、ほぼ定数であると近似できる。

粒子ｋが位置座標ｑ_ｋに存在する確率ｐ（ｑ_ｋ）は、以下のように近似することができる。

従って、下記の式が導出される。

ここまでは、一次元状に複数の粒子が配置された粒子系の相互作用ポテンシャルの粗視化について説明した。多次元の粒子系の相互作用ポテンシャルは、ポテンシャルムービングの方法により実現できる。

図４Ａ〜図４Ｃを参照して、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する。

図４Ａに示すように、二次元正方格子の格子点の位置に粒子が配置されている。最近接する粒子間の相互作用（ニアレストネイバーカップリング）を実線で示す。粒子が配列する１つの方向をｘ方向とし、それに直交する方向をｙ方向とする。

図４Ｂに示すように、ｘ方向に並んだ粒子のうち、１つ置きに配置された粒子（中空の円で表された粒子）の変位について積分を実行することを考える。積分の対象となる粒子（消去したい粒子）を被積分粒子と呼ぶこととする。

図４Ｃに示すように、ポテンシャルムービングの方法では、被積分粒子同士のニアレストネイバーカップリングを、ｘ方向に関して両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を二重線で表す。二重線で表されたダブルカップリングは、元のニアレストネイバーカップリングの２倍の強さになる。この手法を用いて、二次元格子を一次元鎖に変換することができる。一次元鎖の粒子系においては、図２を参照して説明した手法により相互作用ポテンシャルの粗視化を行うことができる。シミュレーション対象の粒子系が三次元格子を構成する場合には、二次元格子を一次元鎖に変換する手順をｘ方向、ｙ方向、及びｚ方向の３方向に関して繰り返せばよい。このようにして、多次元格子を構成する粒子系の粗視化を行うことができる。

次元数ｄの多次元格子を構成する粒子系の相互作用ポテンシャルの粗視化は、下記の式で表される。

ここで、＜ｉ，ｊ＞は、最近接格子間で和を取ることを意味する。

最近接格子間の和を、すべての相互作用の和に変更すると、下記の式が得られる。

［運動エネルギの粗視化］
次に、運動エネルギの粗視化について説明する。運動エネルギについては、容易に積分を実行することができ、下記の式が導出される。

式（３４）の導出に際し、下記の式が利用される。ここで、運動量ベクトルｐ_ｊ ^２は、ベクトルの内積を意味する。

［くりこみ変換測の導出］
次に、上述の相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化から導出されるくりこみ変換測について説明する。

式（２２）に、式（３３）及び式（３４）を代入することにより、結果に影響しない係数を除いて下記の式が得られる。

式（３６）から、粗視化されたハミルトニアン（くりこまれた粒子系Ｓ’のハミルトニアン）Ｈ’は、以下の式で表される。

ハミルトニアンを粗視化する際の結合定数のリストをＫで表す。結合定数のリストＫは下記のように表される。

くりこみ変換Ｒを以下のように定義する。

ｎ回のくりこみ変換実行後の結合係数のリストＫ_ｎは、下記の式で表される。

従って、ｎ回のくりこみ変換が行われたハミルトニアンＨｎは、下記の式で表される。

ここで、くりこまれた粒子系Ｓ’における運動量ベクトルｐ_ｊ’は、下記の式で表される。

式（２１）と式（４１）とから、式（１３）に示したくりこみ変換測が導き出される。

［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ_０＋Ｕ_ｃｈの場合について考察］
次に、ポリマー内の隣接するモノマー粒子間における相互作用ポテンシャルφ（ｒ）に関するくりこみ変換則の導出について説明する。

図２において、粒子ｉ、ｋが作る粒子ｊに対するケージポテンシャルは、式（２５）と同様に、

と表すことができる。積分変数ｑ_ｊについて、モノマー粒子ｊが一方のモノマー粒子ｉに接触した位置から、他方のモノマー粒子ｋに接触する位置まで積分する。モノマー粒子の直径をｒ_ａで表すと、下記の式が得られる。

上述の式のように、くりこみ後の相互作用ポテンシャルは、２φ（（ｑｉ−ｑｋ）／２）とδφとの関数の重ね合わせとなる。

有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε、σ）として、クレマーグレストモデルとして知られている以下のポテンシャルを採用することができる。

ここで、パラメータｋ及びＲ_０は、式（１０）に現れているパラメータε及びσに依存して決定される。一例として、パラメータｋ及びＲ_０は下記の式で定義される。

式（４６）の関係を適用した場合、有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ｒ：ε、σ）も、式（１０）に示したポテンシャルＵ_０（ｒ）と同様に、下記の式で表すことができる。ここで、ｆ_ｃｈは、無次元化された関数である。

相互作用ポテンシャルφ（ｒ）は、

と表される。ポテンシャルＵ_０（ｒ）は、上述の式（１２）で示したとおりである。式（４４）の２φ（（ｑ_ｉ−ｑ_ｋ）／２）を実際に計算し、δφを数値積分により求めた。

図５に、２φ（（ｑ_ｉ−ｑ_ｋ）／２）の計算結果、及びδφの数値積分結果を示す。横軸は（ｑ_ｉ−ｑ_ｋ）／２を表し、縦軸はエネルギを単位「Ｊ」で表す。図５の実線は２φ（（ｑ_ｉ−ｑ_ｋ）／２）の計算結果を示し、破線はδφの数値積分結果を示す。図５の計算結果から、δφの変化量は、２φ（（ｑ_ｉ−ｑ_ｋ）／２）の変化量に比べて十分小さいことがわかる。従って、２φ（（ｑ_ｉ−ｑ_ｋ）／２）に対してδφは定数であると近似することができる。

上述の考察から、下記の式が導出される。

この式は、相互作用ポテンシャルφ＝Ｕ_０の粗視化の説明で示した式（３１）と等価である。従って、相互作用ポテンシャルφ＝Ｕ_０＋Ｕ_ｃｈの粗視化も、相互作用ポテンシャルφ＝Ｕ_０の粗視化と同様に取り扱うことができる。

従って、くりこまれた粒子系Ｓ’における相互作用ポテンシャルφ’（ｒ）は、以下の式で表される。

上式に、式（１３）のくりこみ変換則を適用すると、次の式が得られる。

シミュレーションを行うときに、式（５１）を適用してもよい。この場合には、粒子系Ｓの相互作用ポテンシャルφ（ｒ）の変数ｒをｒ／λに置き換えてポテンシャルφ（ｒ／λ）を求め、その結果をλ^ｄ倍することにより、粒子系Ｓ’について数値計算を行うことができる。

［ポリマー数、ポリマー内のモノマー粒子の個数のくりこみ変換則］
次に、図６を参照して、ポリマーの数Ｎｐ、及びポリマー内のモノマー粒子の個数Ｎｍのくりこみ変換則について説明する。
図６に、直線状のポリマーが３次元方向に規則正しく配置されている状態の模式図を示す。ポリマーの各々を構成するモノマー粒子がｙ方向に繋がっている。複数のポリマーが、ｘ、ｙ、ｚ方向に配置されている。

ｙ方向に関して１回のくりこみを行うと、ｙ方向に並ぶモノマー粒子が１つおきに消去される。例えば、矢印Ｐｙで示された面に配置されたモノマー粒子が消去される。これにより、ポリマー内のモノマー粒子の個数が１／２倍になる。次に、ｘ方向に関して１回のくりこみを行うと、ｘ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｘで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数が１／２倍になる。さらに、ｚ方向に関して１回のくりこみを行うと、ｚ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｚで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数がさらに１／２倍になる。

上述のように、３次元方向に１回のくりこみを行うと、モノマー粒子の個数が１／２倍になり、ポリマーの数が１／４倍になる。くりこみ回数をｎ回とすると、式（１３）に示したＮｍ及びＮｐのくりこみ変換則が得られる。

上記実施例によるシミュレーション方法、及びくりこみ変換を行わないシミュレーション方法を用いて、高分子に特有の現象である粘度のひずみ速度依存性を求めた。計算手法として、ＳＬＬＯＤ法を用いた。シミュレーション対象のパラメータは下記の通りである。
ｍ＝４２．３［ｇ／ｍｏｌ］
ε／ｋ_Ｂ＝４４３．０［Ｋ］
σ＝５．３×１０^−１０［ｍ］
ｋ＝６．５３×１０^−１［Ｎ／ｍ］
Ｒ_０＝７．９５×１０^−１０［ｍ］
Ｔ＝４４３．０［Ｋ］
Ｎ_ｍ’＝６４
Ｎ_ｐ’＝１０００
ここで、ｋ_Ｂはボルツマン定数である。ε、σ、ｋ、及びＲ_０は、式（１２）及び式（４５）に現れているパラメータである。パラメータｋ及びＲ_０は、レナード・ジョーンズ型ポテンシャル（式（１２））のε及びσによって決定される。くりこみ後の粒子系Ｓ’におけるｋ’及びＲ_０’は、くりこみ変換後のε’及びＲ_０’により決定される。

図７に、シミュレーション結果を示す。横軸はひずみ速度を表し、縦軸は粘度を表す。図７では、ひずみ速度及び粘度が無次元化されて表されている。図７中の丸記号はくりこみ変換を行わない粒子系Ｓについてシミュレーションを行った結果を示す、三角記号は、くりこみを行った粒子系Ｓ’についてシミュレーションを行った結果を示す。くりこみ後の粒子系Ｓ’について求めた粘度が、くりこみ前の粒子系Ｓについて求めた粘度を連ねる曲線に乗っていることがわかる。このように、くりこみを行った粒子系Ｓ’においても、粘度のひずみ速度依存性が保たれている。

図７に示したシミュレーション結果により、上記実施例によるくりこみ変換則を用いたシミュレーション方法を用いて、高分子材料のマクロな計算を行うことが可能であることが示された。

上記実施例によるシミュレーション方法は、コンピュータがコンピュータプログラムを実行することにより実現される。コンピュータプログラムは、例えばデータ記録媒体に記録された状態で提供される。または、コンピュータプログラムは、電気通信回線を介して提供することも可能である。

図８に、実施例によるシミュレーション装置２５のブロック図を示す。シミュレーション装置２５は、入力装置２１、処理ユニット１０、及び出力装置２２を含む。処理ユニット１０は、入力部１１、演算部１２、出力部１３、及び記憶装置１４を含む。

入力装置２１を通して入力部１１にシミュレーションの種々の条件が入力される。例えば、シミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐの値、及び粒子系Ｓの初期条件が入力される、

記憶装置１４に、シミュレーションを実行するコンピュータプログラムが格納されている。演算部１２は、入力部１１に入力されたデータに基づいて、記憶装置１４に格納されているコンピュータプログラムを実行することにより、シミュレーション処理を実行する。シミュレーション結果が、出力部１３を介して出力装置２２に出力される。

演算部１２は、入力された粒子系Ｓのパラメータに対してくりこみ変換処理を行うことにより、くりこまれた粒子系Ｓ’のパラメータを算出する。算出されたパラメータを用いて、粒子系Ｓ’について分子動力学計算を行う。くりこまれた粒子系Ｓ’について分子動力学計算を行うため、計算時間を短縮することができる。

以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。

１０ 処理ユニット
１１ 入力部
１２ 演算部
１３ 出力部
１４ 記憶装置
２１ 入力装置
２２ 出力装置
２５ シミュレーション装置

Claims (7)

1. くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
   モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
   してくりこみ変換処理を行う工程と、
   くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
   粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
   を有し、
   粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
   粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
   数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
   粒子間の相互作用ポテンシャルφ（ｒ）として、
   任意のモノマー粒子間においては、
   

   

   
   が適用され、ポテンシャルＵ_０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
   次元を持つパラメータσと、無次元関数ｆを用いて
   

   

   
   と表すことができ、
   同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
   Ｕ_０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵ_ｃｈ
   （ｒ：ε，σ）がポテンシャルＵ_０（ｒ）に追加されて
   

   

   
   が適用され、
   粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
   ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
   表したとき、下記の変換則
   

   

   
   を適用し、
   粒子系ＳにおけるポテンシャルＵ_０（ｒ）、Ｕ_ｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ
   ’におけるポテンシャルＵ_０’（ｒ）、Ｕ_ｃｈ’（ｒ：ε，σ）として、
   

   

   
   を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。
2. 粒子系ＳにおけるポテンシャルＵ_０（ｒ）は、レナード・ジョーンズ型ポテンシャルの
   形で表される請求項１に記載のシミュレーション方法。
3. 有限伸長非線形弾性ポテンシャルＵ_ｃｈ（ε，σ，ｒ）は、パラメータε、_σに依存す
   るパラメータｋ及びＲ_０を用いて、
   

   

   
   と表される請求項１または２に記載のシミュレーション方法。
4. 請求項１乃至３のいずれか１項に記載のシミュレーション方法をコンピュータで実行す
   るためのコンピュータプログラム。
5. 請求項４に記載のコンピュータプログラムが、コンピュータで読み取り可能に記録され
   た記録媒体。
6. 各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシ
   ミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー
   粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐ
   の値、及び粒子系Ｓの初期条件が入力される入力装置と、
   前記入力装置から入力されたデータに基づいてシミュレーション処理を実行する処理ユ
   ニットと
   を有し、
   前記処理ユニットは、
   粒子系Ｓに対して、くりこみ回数に依存するくりこみ因子λに基づいてくりこみ変換を
   行う処理と、
   くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
   粒子系Ｓ’のモノマー粒子の位置及び運動量を求める処理と
   を実行し、
   前記処理ユニットは、
   粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、
   粒子系Ｓのモノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ
   （ｒ）として、
   任意のモノマー粒子間においては、
   

   

   
   が適用され、ポテンシャルＵ_０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
   次元を持つパラメータσと、無次元関数ｆを用いて
   

   

   
   と表すことができ、
   同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
   Ｕ_０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵ_ｃｈ
   （ｒ：ε，σ）がポテンシャルＵ_０（ｒ）に追加されて
   

   

   
   が適用され、
   粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
   ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
   表したとき、下記の変換則
   

   

   
   を適用し、
   粒子系ＳにおけるポテンシャルＵ_０（ｒ）、Ｕ_ｃｈ（ε，σ，ｒ）に対応する粒子系Ｓ
   ’におけるポテンシャルＵ_０’（ｒ）、Ｕ_ｃｈ’（ε，σ，ｒ）として、
   

   

   
   を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション装置。
7. くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
   モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
   してくりこみ変換処理を行う工程と、
   くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
   粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
   を有し、
   粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
   粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
   数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
   粒子間の相互作用ポテンシャルφ（ｒ）として、
   任意のモノマー粒子間においては、
   

   

   
   が適用され、
   同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
   Ｕ_０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵ_ｃｈ
   （ｒ：ε，σ）がポテンシャルＵ_０（ｒ）に追加されて
   

   

   
   が適用され、
   ポテンシャルφ（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパ
   ラメータσと、無次元関数ｆを用いて
   

   

   
   と表すことができ、
   粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
   ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
   表したとき、下記の変換則
   

   

   
   を適用し、
   粒子系Ｓにおけるポテンシャルφ（ｒ）に対応する粒子系Ｓ’におけるポテンシャルφ
   ’（ｒ）として、
   

   

   
   を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。

Images