(57)【特許請求の範囲】
【請求項１】
くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
粒子間の相互作用ポテンシャルφ（ｒ）として、
任意のモノマー粒子間においては、
【数１】
0006651254.tif 000053
が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
次元を持つパラメータσと、無次元関数ｆを用いて
【数２】
0006651254.tif 000054
と表すことができ、
同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
Ｕ０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ
（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
【数３】
0006651254.tif 000055
が適用され、
粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
表したとき、下記の変換則
【数４】
0006651254.tif 000056
を適用し、
粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ
’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ｒ：ε，σ）として、
【数５】
0006651254.tif 000057
を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。
【請求項２】
粒子系ＳにおけるポテンシャルＵ０（ｒ）は、レナード・ジョーンズ型ポテンシャルの
形で表される請求項１に記載のシミュレーション方法。
【請求項３】
有限伸長非線形弾性ポテンシャルＵｃｈ（ε，σ，ｒ）は、パラメータε、σに依存す
るパラメータｋ及びＲ０を用いて、
【数６】
0006651254.tif 000058
と表される請求項１または２に記載のシミュレーション方法。
【請求項４】
請求項１乃至３のいずれか１項に記載のシミュレーション方法をコンピュータで実行す
るためのコンピュータプログラム。
【請求項５】
請求項４に記載のコンピュータプログラムが、コンピュータで読み取り可能に記録され
た記録媒体。
【請求項６】
各々が一次元状に繋がった複数のモノマー粒子からなる複数のポリマーで構成されたシ
ミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー
粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐ
の値、及び粒子系Ｓの初期条件が入力される入力装置と、
前記入力装置から入力されたデータに基づいてシミュレーション処理を実行する処理ユ
ニットと
を有し、
前記処理ユニットは、
粒子系Ｓに対して、くりこみ回数に依存するくりこみ因子λに基づいてくりこみ変換を
行う処理と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
粒子系Ｓ’のモノマー粒子の位置及び運動量を求める処理と
を実行し、
前記処理ユニットは、
粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、
粒子系Ｓのモノマー粒子の座標をｑで表し、モノマー粒子間の相互作用ポテンシャルφ
（ｒ）として、
任意のモノマー粒子間においては、
【数７】
0006651254.tif 000059
が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
次元を持つパラメータσと、無次元関数ｆを用いて
【数８】
0006651254.tif 000060
と表すことができ、
同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
Ｕ０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ
（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
【数９】
0006651254.tif 000061
が適用され、
粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
表したとき、下記の変換則
【数１０】
0006651254.tif 000062
を適用し、
粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ε，σ，ｒ）に対応する粒子系Ｓ
’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ε，σ，ｒ）として、
【数１１】
0006651254.tif 000063
を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション装置。
【請求項７】
くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
粒子間の相互作用ポテンシャルφ（ｒ）として、
任意のモノマー粒子間においては、
【数１２】
0006651254.tif 000064
が適用され、
同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
Ｕ０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ
（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
【数１３】
0006651254.tif 000065
が適用され、
ポテンシャルφ（ｒ）は、エネルギの次元を持つパラメータεと、長さの次元を持つパ
ラメータσと、無次元関数ｆを用いて
【数１４】
0006651254.tif 000066
と表すことができ、
粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
表したとき、下記の変換則
【数１５】
0006651254.tif 000067
を適用し、
粒子系Ｓにおけるポテンシャルφ（ｒ）に対応する粒子系Ｓ’におけるポテンシャルφ
’（ｒ）として、
【数１６】
0006651254.tif 000068
を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、くりこみ群を応用した分子動力学によるシミュレーション方法、シミュレーションプログラム、及びシミュレーション装置に関する。
【背景技術】
【０００２】
分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が数値的に解かれる。シミュレーション対象の系に含まれる粒子数が増加すれば、必要な計算量が増大する。現行のコンピュータの演算能力でシミュレーション可能な系の粒子数は、通常、数十万個程度である。
【０００３】
下記の特許文献１に、必要な計算量を減らすためにくりこみ変換の手法を用いたシミュレーション方法が開示されている。以下、特許文献１に開示されたくりこみ変換の手法について説明する。
【０００４】
シミュレーション対象の粒子系Ｓの粒子数をＮ、各粒子の質量をｍ、粒子間の相互作用ポテンシャルをφ（ｒ）で表す。ここで、ｒは、粒子間の距離である。相互作用ポテンシャルφ（ｒ）が、相互作用係数εと関数ｆ（ｒ）との積で表される。相互作用係数εは、相互作用の強さを表し、エネルギの次元を持つ。関数ｆ（ｒ）は、粒子間距離に対する依存性を表し、無次元化されている。
【０００５】
第１のくりこみ因子α、第２のくりこみ因子γ、及び第３のくりこみ因子δが決定される。第１のくりこみ因子αは１より大きい。第２のくりこみ因子γは０以上であり、かつ空間の次元数ｄ以下である。第３のくりこみ因子δは０以上である。くりこみの回数をｎで表したとき、第１のくりこみ因子αは、
α＝２ｎ
と表される。
【０００６】
くりこみ手法を用いてくりこみ変換された粒子系Ｓ’の粒子数をＮ’、各粒子の質量をｍ’、相互作用係数をε’で表す。以下の変換式を用いて、くりこみ変換された粒子系Ｓ’の粒子数Ｎ’、質量ｍ’、及び相互作用係数ε’を求める。
Ｎ’＝Ｎ／αｄ
ｍ’＝ｍαδ−γ
ε’＝εαγ
【０００７】
くりこみ変換された粒子系Ｓ’について、分子動力学計算を行なう。分子動力学計算で求められた各粒子の位置ベクトルをｑ’、運動量ベクトルをｐ’で表す。元の粒子系Ｓの各粒子の位置ベクトルｑ、運動量ベクトルｐは、以下の式を用いて算出することができる。
ｑ＝ｑ’α
ｐ＝ｐ’／αδ／２
【先行技術文献】
【特許文献】
【０００８】
【特許文献１】特許第５２４１４６８号公報
【発明の概要】
【発明が解決しようとする課題】
【０００９】
粒子系Ｓが単原子の集合体である場合、特許文献１に開示されたくりこみ変換により、くりこまれた粒子系Ｓ’に対して分子動力学計算を行うことができる。多数のポリマー分子からなる粒子系Ｓに対して分子動力学計算を行う手法として、クレマーグレストモデル（Kremer Grest Model）が提案されている。クレマーグレストモデルにくりこみ変換処理を適用してシミュレーションを行う手法は確立されていない。
【００１０】
本発明の目的は、ポリマー分子からなる粒子系に対してくりこみ群の考え方を適用し、ポリマー分子の挙動をシミュレーションする方法を提供することである。本発明の他の目的は、このシミュレーション方法を実行するコンピュータプログラムを提供することである。本発明のさらに他の目的は、このシミュレーション方法を実行するシミュレーション装置を提供することである。
【課題を解決するための手段】
【００１１】
本発明の一観点によると、
くりこみ回数に依存するくりこみ因子λに基づいて、各々が一次元状に繋がった複数の
モノマー粒子からなる複数のポリマーで構成されたシミュレーション対象の粒子系Ｓに対
してくりこみ変換処理を行う工程と、
くりこまれた粒子系Ｓ’について分子動力学計算を実行することにより、くりこまれた
粒子系Ｓ’のモノマー粒子の位置ベクトル及び運動量ベクトルを求める工程と
を有し、
粒子系Ｓ及び粒子系Ｓ’のモノマー間の距離をｒで表し、空間の次元数をｄで表し、
粒子系Ｓのモノマーの質量をｍで表し、１つのポリマーを構成しているモノマー粒子の
数をＮｍで表し、ポリマーの数をＮｐで表し、モノマー粒子の座標をｑで表し、モノマー
粒子間の相互作用ポテンシャルφ（ｒ）として、
任意のモノマー粒子間においては、
【数１】
0006651254.tif 000002
が適用され、ポテンシャルＵ０（ｒ）は、エネルギの次元を持つパラメータεと、長さの
次元を持つパラメータσと、無次元関数ｆを用いて
【数２】
0006651254.tif 000003
と表すことができ、
同一ポリマー内の隣接するモノマー粒子間においては、パラメータとしてポテンシャル
Ｕ０（ｒ）で用いられたパラメータε、σを持つ有限伸長非線形弾性ポテンシャルＵｃｈ
（ｒ：ε，σ）がポテンシャルＵ０（ｒ）に追加されて
【数３】
0006651254.tif 000004
が適用され、
粒子系Ｓ’のモノマー粒子の質量をｍ’で表し、１つのポリマーを構成しているモノマ
ー粒子の数をＮｍ’で表し、ポリマーの数をＮｐ’で表し、モノマー粒子の座標をｑ’で
表したとき、下記の変換則
【数４】
0006651254.tif 000005
を適用し、
粒子系ＳにおけるポテンシャルＵ０（ｒ）、Ｕｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ
’におけるポテンシャルＵ０’（ｒ）、Ｕｃｈ’（ｒ：ε，σ）として、
【数５】
0006651254.tif 000006
を用いて、粒子系Ｓ’について前記分子動力学計算を実行するシミュレーション方法が提
供される。
【００１２】
本発明の他の観点によると、上記シミュレーションを実行するコンピュータプログラム、コンピュータプログラムを記録した記録媒体、及びシミュレーション装置が提供される。
【発明の効果】
【００１３】
複数のモノマー粒子が鎖状に繋がったポリマーからなる粒子系にくりこみ変換を適用して、分子動力学計算を行うことができる。これにより、計算時間を短縮することができる。
【図面の簡単な説明】
【００１４】
【図１】図１は、実施例によるシミュレーション方法のフローチャートである。
【図２】図２は、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている状態を示す模式図である。
【図３】図３は、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示すグラフである。
【図４】図４Ａ〜図４Ｃは、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明するための図である。
【図５】図５は、２φ（（ｑｉ−ｑｋ）／２）の計算結果、及びδφの数値積分結果を示すグラフである。
【図６】図６は、直線状のポリマーが３次元方向に規則正しく配置されている例を示す模式図である。
【図７】図７は、高分子材料の粘度のひずみ速度依存性をシミュレーションによって求めた結果を示すグラフである。
【図８】図８は、実施例によるシミュレーション装置のブロック図である。
【発明を実施するための形態】
【００１５】
本願発明の実施例で適用される分子動力学（Molecular Dynamics）について簡単に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが下記の式で表される粒子系について考える。
【数６】
0006651254.tif 000007
ここで、
ｍは粒子の質量を表し、
φは粒子間の相互作用ポテンシャルを表し、
ベクトルｐｊは粒子の運動量ベクトルを表し、
ベクトルｑｊは粒子の位置ベクトル（位置座標）を表す。
【００１６】
ハミルトニアンＨを、ハミルトニアンの正準方程式に代入することにより、粒子ｊに対する下記の運動方程式が得られる。
【数７】
0006651254.tif 000008
【数８】
0006651254.tif 000009
【００１７】
分子動力学では、粒子系を構成する各粒子に対し、式（７）及び式（８）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトルｐｊ及び位置ベクトルｑｊが求められる。数値積分には、Verlet法が用いられることが多い。Verlet法については、例えば、J.M.Thijssen,「Computational Physics」,Cambridge University Press (1999)の１７５頁に説明されている。分子動力学計算で得られた各粒子の運動量ベクトル及び位置ベクトルに基づいて、粒子系の種々の物理量を算出することができる。
【００１８】
次に、くりこみ群の手法を応用した分子動力学（以下、くりこみ群分子動力学という。）について概念的に説明する。
【００１９】
くりこみ群分子動力学では、シミュレーション対象の粒子系Ｓを、粒子系Ｓの粒子数よりも少ない数の粒子で構成された粒子系Ｓ’（以下、くりこまれた粒子系Ｓ’という。）に対応づける。次に、くりこまれた粒子系Ｓ’に対して分子動力学計算を実行する。くりこまれた粒子系Ｓ’に対する計算結果を、シミュレーション対象の粒子系Ｓに対応付ける。これにより、シミュレーション対象の粒子系Ｓに対して分子動力学計算を直接的に実行する場合よりも計算量を少なくすることができる。シミュレーション対象の粒子系Ｓにおける物理量（例えば、粒子数、粒子の質量等）と、くりこまれた粒子系Ｓ’における物理量とを相互に対応付けるための変換則を、くりこみ変換測という。
【００２０】
図１に、実施例によるシミュレーション方法のフローチャートを示す。ステップＳ１において、シミュレーション対象の粒子系Ｓに対して、くりこみ変換測に基づいてくりこみ変換処理を実行することにより、くりこまれた粒子系Ｓ’を定義する。粒子系Ｓは、複数のポリマーで構成され、複数のポリマーの各々は、一次元状に繋がった複数のモノマー粒子で構成される。
【００２１】
粒子系Ｓのモノマー粒子間の相互作用ポテンシャルφ（ｒ）として、任意のモノマー粒子間においては、
【数９】
0006651254.tif 000010
が適用される。ここで、ｒは、粒子系Ｓを構成する各粒子からの距離を表す。
【００２２】
ポテンシャルＵ０（ｒ）は、基本的に下記の式で表される。
【数１０】
0006651254.tif 000011
ここで、
ｆは無次元化された関数を表し、ε及びσは、モノマー粒子を特徴づけるパラメータである。パラメータεはエネルギの次元を持ち、相互作用係数と呼ばれる。パラメータσは距離の次元を持ち、粒子の大きさに依存する。
【００２３】
同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ０（ｒ）に有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）が追加されて
【数１１】
0006651254.tif 000012
が適用される。有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε，σ）は、ポテンシャルＵ０（ｒ）を定義するパラメータε及びσに依存するパラメータを含む。
【００２４】
ポテンシャルＵ０（ｒ）として、例えばレナード・ジョーンズ型のポテンシャルを適用することができる。ポテンシャルＵ０（ｒ）に、シフトさせたレナード・ジョーンズ型のポテンシャルを適用した場合、ポテンシャルＵ０（ｒ）は、例えば下記の式で定義することができる。
【数１２】
0006651254.tif 000013
【００２５】
ポテンシャルＵ０（ｒ）として、モース型ポテンシャルを適用することも可能である。さらに、ポテンシャルＵ０（ｒ）として、これらの関数をテーラー展開した近似関数を用いることも可能である。本実施例においては、ポテンシャルＵ０（ｒ）としてレナード・ジョーンズ型のポテンシャルを適用することとして説明を続ける。
【００２６】
くりこみ変換処理により、シミュレーション対象の粒子系Ｓのモノマー粒子の位置ベクトルｑ、モノマー粒子の質量ｍ、ポリマー内のモノマー粒子の個数Ｎｍ、ポリマーの数Ｎｐ、ポテンシャルＵ０（ｒ）を定義するパラメータε、及びσが、それぞれ繰り込まれた粒子系Ｓ’のモノマー粒子の位置ベクトルｑ’、モノマー粒子の質量ｍ’、ポリマー内のモノマー粒子の個数Ｎｍ’、ポリマーの数Ｎｐ’、ポテンシャルＵ０’（ｒ）を定義するパラメータε’、及びσ’に変換される。ステップＳ１のくりこみ変換処理では、下記のくりこみ変換測が適用される。
【数１３】
0006651254.tif 000014
ここで、ｄはシミュレーション対象の粒子系Ｓが配置された空間の次元数を表す。λは繰り込み回数に依存するくりこみ因子である。くりこみ回数がｎ回のとき、くりこみ因子λは以下の式で表される。
【数１４】
0006651254.tif 000015
【００２７】
粒子系ＳにおけるポテンシャルＵ０（ｒ）及びＵｃｈ（ｒ：ε，σ）に対応する粒子系Ｓ’におけるポテンシャルＵ０’（ｒ）及びＵｃｈ’（ｒ：ε，σ）として、下記のポテンシャルを用いる。
【数１５】
0006651254.tif 000016
【００２８】
粒子系Ｓ’の任意のモノマー粒子間に適用される相互作用ポテンシャルφ’（ｒ）は、下記の式で定義される。
【数１６】
0006651254.tif 000017
【００２９】
同一ポリマー内の隣接するモノマー粒子間においては、ポテンシャルＵ０’（ｒ）に有限伸長非線形弾性ポテンシャルＵｃｈ’（ｒ：ε，σ）が追加されて、相互作用ポテンシャルφ’として、
【数１７】
0006651254.tif 000018
が適用される。
【００３０】
上述のくりこみ変換測を適用することにより、粒子の個数が１／λｄ倍になり、粒子の質量がλｄ倍になる。このため、粒子系Ｓの全体の質量と、粒子系Ｓ’の全体の質量とは同一である。さらに、粒子系Ｓ’におけるモノマー粒子の位置ベクトルｑ’が、粒子系Ｓにおけるモノマー粒子の位置ベクトルｑと同一であるため、粒子系Ｓとくりこまれた粒子系Ｓ’との巨視的な寸法は同一である。くりこみ変換処理の前後で、粒子系の全質量及び寸法が不変であるため、粒子系の密度も不変である。
【００３１】
次に、ステップＳ２において、シミュレーションの初期条件を設定する。初期条件には、各粒子の位置ベクトルｑｊ及び運動量ベクトルｐｊの初期値が含まれる。運動量ベクトルｐｊは、くりこまれた粒子系Ｓ’の温度Ｔ’に基づいて設定される。シミュレーション対象の粒子系Ｓの温度をＴで表したとき、温度に関して下記のくりこみ変換測が適用される。
【数１８】
0006651254.tif 000019
【００３２】
次に、ステップＳ３において、くりこまれた粒子系Ｓ’に関して分子動力学計算を実行する。運動方程式は、下記の式で表される。
【数１９】
0006651254.tif 000020
ここで、Ｎ’は粒子系Ｓ’のモノマー粒子の全個数を表しており、Ｎ’＝Ｎｍ’×Ｎｐ’である。
【００３３】
この運動方程式を、数値的に積分して解く。これにより、くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’及び運動量ベクトルｐ’の時刻歴が求まる。
【００３４】
くりこまれた粒子系Ｓ’の各粒子の位置ベクトルｑ’と、シミュレーション対象の粒子
系Ｓの各粒子の位置ベクトルｑとの対応関係は、式（１３）に示したとおりである。くりこまれた粒子系Ｓ’における運動量ベクトルｐ’と、シミュレーション対象の粒子系Ｓにおける運動量ベクトルｐとは、下記の関係を有する。
【数２０】
0006651254.tif 000021
【００３５】
ステップＳ４において、シミュレーション結果を出力する。例えば、位置ベクトルｑ’及び運動量ベクトルｐ’を、そのまま数値で出力してもよいし、位置ベクトルｑ’に基づいて、粒子系Ｓ’の複数の粒子の空間内の分布を画像として表示してもよい。
【００３６】
［φ（ｒ）＝Ｕ０（ｒ）の場合について考察］
次に、相互作用ポテンシャルφ（ｒ）＝Ｕ０（ｒ）の場合について、くりこみ変換則を導き出す。まず、くりこまれた粒子系Ｓ’のハミルトニアンＨ’が
【数２１】
0006651254.tif 000022
となることを示す。くりこみ前の粒子系Ｓの粒子数をＮで表し、くりこまれた粒子系Ｓ’の粒子数をＮ’で表す。このハミルトニアンＨ’を正準方程式に代入して運動方程式を得ることができる。ハミルトニアンＨ’のパラメータから、式（１３）のくりこみ変換則が導出される。
【００３７】
くりこまれた粒子系Ｓ’のハミルトニアンＨ’を得るためには、粒子系Ｓに対する分配関数Ｚ（β）の積分の一部を実行し、ハミルトニアンを粗視化することによりハミルトニアンＨ’が得られる。
【００３８】
粒子数が一定の正準集団（Canonical ensemble）に対する分配関数Ｚ（β）は、以下の式で表される。
【数２２】
0006651254.tif 000023
ここで、ｋＢはボルツマン定数であり、ｄΓＮは位相空間内の体積要素である。ｄΓＮは下記の式で表される。
【数２３】
0006651254.tif 000024
ここで、ｈはプランク定数である。ＷＮは、すべての状態の量子論的な和と、位相空間に亘る積分とが一致するように決められる。
【００３９】
まず、粒子間の相互作用ポテンシャルの粗視化について説明し、次に運動エネルギの粗視化について説明する。続いて、相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化を踏まえて、くりこみ変換測を定義する。
【００４０】
［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ０の粗視化］
まず、一次元鎖状に配置された粒子系における相互作用ポテンシャルの粗視化について説明する。その後、単純立方格子状に配置された粒子における相互作用ポテンシャルについて説明する。
【００４１】
図２に示すように、粒子ｉ、粒子ｊ、及び粒子ｋがこの順番に一次元状に並んで配置されている。粒子ｊが関与する相互作用を書き出し、粒子ｉと粒子ｋとの中間にある粒子ｊの位置座標について積分を実行することにより、相互作用ポテンシャルの粗視化を行うことができる。まず、第二近接以上の粒子からの寄与を取り込むために、ポテンシャルムービング（Potential Moving）の方法を用いることができる。ポテンシャルムービングの方法については、Leo P. Kadanoff, STATISTICAL PHYSICS Static, Dynamics and Renormalization, Chap. 14, World Scientific (1999)に解説されている。
【００４２】
第二近接以上の粒子からの寄与を取り込んだ相互作用ポテンシャルφティルダは、下記の式で表すことができる。
【数２４】
0006651254.tif 000025
ここで、ａは平衡状態における粒子間距離を表す。平衡状態における粒子間距離ａは、相互作用ポテンシャルφ（ｒ）が最小となる距離に等しい。
【００４３】
複数の粒子が一次元状に配置されているため、粒子ｊの位置ベクトルｑｊは、一次元座標ｑｊで表すことができる。粒子ｊの位置をｑｊで表すと、粒子ｊに対して、最近接の粒子によって作られるケージポテンシャルは、下記の式で表される。
【数２５】
0006651254.tif 000026
【００４４】
積分変数ｑｊについて積分を実行すれば、式（２５）を用いて下記の式が得られる。
【数２６】
0006651254.tif 000027
ここで、ｒａは粒子の直径を表し、ｚ（ｑｉ−ｑｋ）及びＰ（ｑｉ−ｑｋ）は、下記の式で表される。
【数２７】
0006651254.tif 000028
【数２８】
0006651254.tif 000029
積分領域はケージポテンシャルの内部領域に制限される。
【００４５】
次に、ｚ（ｑｉ−ｑｋ）を具体的に計算する。相互作用ポテンシャルφがレナード・ジョーンズ型ポテンシャルの場合、φ（２ｎ）は、下記の式で表される。
【数２９】
0006651254.tif 000030
【００４６】
式（２９）を式（２８）に代入して数値積分を行う。式（２８）に式（２９）を代入する際には、式（２４）を用いる。
【００４７】
図３に、相互作用ポテンシャルφ（ｒ）がレナード・ジョーンズ型ポテンシャルである場合の式（２７）の算出結果、及び式（２８）の数値積分結果を示す。粒子ｉ、粒子ｊ、粒子ｋ、粒子ｌ、及び粒子ｍがこの順番に一次元状に配置されている場合の粒子ｋの位置座標がｑｋで表される。図３の横軸は、ｑｋ／２を表し、縦軸は、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）、及びｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）を、対数目盛で表す。図３に至る数値計算において、ε／ｋＢＴ＝２．０とした。
【００４８】
図３から、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）の変化に比べて、ｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）の変化が緩やかであることがわかる。このため、Ｐ（ｑｉ−ｑｋ）Ｐ（ｑｋ−ｑｍ）に対してｚ（ｑｉ−ｑｋ）ｚ（ｑｋ−ｑｍ）は、ほぼ定数であると近似できる。
【００４９】
粒子ｋが位置座標ｑｋに存在する確率ｐ（ｑｋ）は、以下のように近似することができる。
【数３０】
0006651254.tif 000031
【００５０】
従って、下記の式が導出される。
【数３１】
0006651254.tif 000032
【００５１】
ここまでは、一次元状に複数の粒子が配置された粒子系の相互作用ポテンシャルの粗視化について説明した。多次元の粒子系の相互作用ポテンシャルは、ポテンシャルムービングの方法により実現できる。
【００５２】
図４Ａ〜図４Ｃを参照して、二次元格子を一次元格子に帰着させるポテンシャルムービングの方法について説明する。
【００５３】
図４Ａに示すように、二次元正方格子の格子点の位置に粒子が配置されている。最近接する粒子間の相互作用（ニアレストネイバーカップリング）を実線で示す。粒子が配列する１つの方向をｘ方向とし、それに直交する方向をｙ方向とする。
【００５４】
図４Ｂに示すように、ｘ方向に並んだ粒子のうち、１つ置きに配置された粒子（中空の円で表された粒子）の変位について積分を実行することを考える。積分の対象となる粒子（消去したい粒子）を被積分粒子と呼ぶこととする。
【００５５】
図４Ｃに示すように、ポテンシャルムービングの方法では、被積分粒子同士のニアレストネイバーカップリングを、ｘ方向に関して両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を二重線で表す。二重線で表されたダブルカップリングは、元のニアレストネイバーカップリングの２倍の強さになる。この手法を用いて、二次元格子を一次元鎖に変換することができる。一次元鎖の粒子系においては、図２を参照して説明した手法により相互作用ポテンシャルの粗視化を行うことができる。シミュレーション対象の粒子系が三次元格子を構成する場合には、二次元格子を一次元鎖に変換する手順をｘ方向、ｙ方向、及びｚ方向の３方向に関して繰り返せばよい。このようにして、多次元格子を構成する粒子系の粗視化を行うことができる。
【００５６】
次元数ｄの多次元格子を構成する粒子系の相互作用ポテンシャルの粗視化は、下記の式で表される。
【数３２】
0006651254.tif 000033
ここで、＜ｉ，ｊ＞は、最近接格子間で和を取ることを意味する。
【００５７】
最近接格子間の和を、すべての相互作用の和に変更すると、下記の式が得られる。
【数３３】
0006651254.tif 000034
【００５８】
［運動エネルギの粗視化］
次に、運動エネルギの粗視化について説明する。運動エネルギについては、容易に積分を実行することができ、下記の式が導出される。
【数３４】
0006651254.tif 000035
式（３４）の導出に際し、下記の式が利用される。ここで、運動量ベクトルｐｊ２は、ベクトルの内積を意味する。
【数３５】
0006651254.tif 000036
【００５９】
［くりこみ変換測の導出］
次に、上述の相互作用ポテンシャルの粗視化、及び運動エネルギの粗視化から導出されるくりこみ変換測について説明する。
【００６０】
式（２２）に、式（３３）及び式（３４）を代入することにより、結果に影響しない係数を除いて下記の式が得られる。
【数３６】
0006651254.tif 000037
【００６１】
式（３６）から、粗視化されたハミルトニアン（くりこまれた粒子系Ｓ’のハミルトニアン）Ｈ’は、以下の式で表される。
【数３７】
0006651254.tif 000038
【００６２】
ハミルトニアンを粗視化する際の結合定数のリストをＫで表す。結合定数のリストＫは下記のように表される。
【数３８】
0006651254.tif 000039
【００６３】
くりこみ変換Ｒを以下のように定義する。
【数３９】
0006651254.tif 000040
【００６４】
ｎ回のくりこみ変換実行後の結合係数のリストＫｎは、下記の式で表される。
【数４０】
0006651254.tif 000041
【００６５】
従って、ｎ回のくりこみ変換が行われたハミルトニアンＨｎは、下記の式で表される。
【数４１】
0006651254.tif 000042
ここで、くりこまれた粒子系Ｓ’における運動量ベクトルｐｊ’は、下記の式で表される。
【数４２】
0006651254.tif 000043
【００６６】
式（２１）と式（４１）とから、式（１３）に示したくりこみ変換測が導き出される。
【００６７】
［モノマー粒子間の相互作用ポテンシャルφ＝Ｕ０＋Ｕｃｈの場合について考察］
次に、ポリマー内の隣接するモノマー粒子間における相互作用ポテンシャルφ（ｒ）に関するくりこみ変換則の導出について説明する。
【００６８】
図２において、粒子ｉ、ｋが作る粒子ｊに対するケージポテンシャルは、式（２５）と同様に、
【数４３】
0006651254.tif 000044
と表すことができる。積分変数ｑｊについて、モノマー粒子ｊが一方のモノマー粒子ｉに接触した位置から、他方のモノマー粒子ｋに接触する位置まで積分する。モノマー粒子の直径をｒａで表すと、下記の式が得られる。
【数４４】
0006651254.tif 000045
【００６９】
上述の式のように、くりこみ後の相互作用ポテンシャルは、２φ（（ｑｉ−ｑｋ）／２）とδφとの関数の重ね合わせとなる。
【００７０】
有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε、σ）として、クレマーグレストモデルとして知られている以下のポテンシャルを採用することができる。
【数４５】
0006651254.tif 000046
ここで、パラメータｋ及びＲ０は、式（１０）に現れているパラメータε及びσに依存して決定される。一例として、パラメータｋ及びＲ０は下記の式で定義される。
【数４６】
0006651254.tif 000047
【００７１】
式（４６）の関係を適用した場合、有限伸長非線形弾性ポテンシャルＵｃｈ（ｒ：ε、σ）も、式（１０）に示したポテンシャルＵ０（ｒ）と同様に、下記の式で表すことができる。ここで、ｆｃｈは、無次元化された関数である。
【数４７】
0006651254.tif 000048
【００７２】
相互作用ポテンシャルφ（ｒ）は、
【数４８】
0006651254.tif 000049
と表される。ポテンシャルＵ０（ｒ）は、上述の式（１２）で示したとおりである。式（４４）の２φ（（ｑｉ−ｑｋ）／２）を実際に計算し、δφを数値積分により求めた。
【００７３】
図５に、２φ（（ｑｉ−ｑｋ）／２）の計算結果、及びδφの数値積分結果を示す。横軸は（ｑｉ−ｑｋ）／２を表し、縦軸はエネルギを単位「Ｊ」で表す。図５の実線は２φ（（ｑｉ−ｑｋ）／２）の計算結果を示し、破線はδφの数値積分結果を示す。図５の計算結果から、δφの変化量は、２φ（（ｑｉ−ｑｋ）／２）の変化量に比べて十分小さいことがわかる。従って、２φ（（ｑｉ−ｑｋ）／２）に対してδφは定数であると近似することができる。
【００７４】
上述の考察から、下記の式が導出される。
【数４９】
0006651254.tif 000050
【００７５】
この式は、相互作用ポテンシャルφ＝Ｕ０の粗視化の説明で示した式（３１）と等価である。従って、相互作用ポテンシャルφ＝Ｕ０＋Ｕｃｈの粗視化も、相互作用ポテンシャルφ＝Ｕ０の粗視化と同様に取り扱うことができる。
【００７６】
従って、くりこまれた粒子系Ｓ’における相互作用ポテンシャルφ’（ｒ）は、以下の式で表される。
【数５０】
0006651254.tif 000051
【００７７】
上式に、式（１３）のくりこみ変換則を適用すると、次の式が得られる。
【数５１】
0006651254.tif 000052
【００７８】
シミュレーションを行うときに、式（５１）を適用してもよい。この場合には、粒子系Ｓの相互作用ポテンシャルφ（ｒ）の変数ｒをｒ／λに置き換えてポテンシャルφ（ｒ／λ）を求め、その結果をλｄ倍することにより、粒子系Ｓ’について数値計算を行うことができる。
【００７９】
［ポリマー数、ポリマー内のモノマー粒子の個数のくりこみ変換則］
次に、図６を参照して、ポリマーの数Ｎｐ、及びポリマー内のモノマー粒子の個数Ｎｍのくりこみ変換則について説明する。
図６に、直線状のポリマーが３次元方向に規則正しく配置されている状態の模式図を示す。ポリマーの各々を構成するモノマー粒子がｙ方向に繋がっている。複数のポリマーが、ｘ、ｙ、ｚ方向に配置されている。
【００８０】
ｙ方向に関して１回のくりこみを行うと、ｙ方向に並ぶモノマー粒子が１つおきに消去される。例えば、矢印Ｐｙで示された面に配置されたモノマー粒子が消去される。これにより、ポリマー内のモノマー粒子の個数が１／２倍になる。次に、ｘ方向に関して１回のくりこみを行うと、ｘ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｘで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数が１／２倍になる。さらに、ｚ方向に関して１回のくりこみを行うと、ｚ方向に並ぶポリマーが１つおきに消去される。例えば、矢印Ｐｚで示された面に配置されたモノマー粒子が除去される。これにより、ポリマーの数がさらに１／２倍になる。
【００８１】
上述のように、３次元方向に１回のくりこみを行うと、モノマー粒子の個数が１／２倍になり、ポリマーの数が１／４倍になる。くりこみ回数をｎ回とすると、式（１３）に示したＮｍ及びＮｐのくりこみ変換則が得られる。
【００８２】
上記実施例によるシミュレーション方法、及びくりこみ変換を行わないシミュレーション方法を用いて、高分子に特有の現象である粘度のひずみ速度依存性を求めた。計算手法として、ＳＬＬＯＤ法を用いた。シミュレーション対象のパラメータは下記の通りである。
ｍ＝４２．３［ｇ／ｍｏｌ］
ε／ｋＢ＝４４３．０［Ｋ］
σ＝５．３×１０−１０［ｍ］
ｋ＝６．５３×１０−１［Ｎ／ｍ］
Ｒ０＝７．９５×１０−１０［ｍ］
Ｔ＝４４３．０［Ｋ］
Ｎｍ’＝６４
Ｎｐ’＝１０００
ここで、ｋＢはボルツマン定数である。ε、σ、ｋ、及びＲ０は、式（１２）及び式（４５）に現れているパラメータである。パラメータｋ及びＲ０は、レナード・ジョーンズ型ポテンシャル（式（１２））のε及びσによって決定される。くりこみ後の粒子系Ｓ’におけるｋ’及びＲ０’は、くりこみ変換後のε’及びＲ０’により決定される。
【００８３】
図７に、シミュレーション結果を示す。横軸はひずみ速度を表し、縦軸は粘度を表す。図７では、ひずみ速度及び粘度が無次元化されて表されている。図７中の丸記号はくりこみ変換を行わない粒子系Ｓについてシミュレーションを行った結果を示す、三角記号は、くりこみを行った粒子系Ｓ’についてシミュレーションを行った結果を示す。くりこみ後の粒子系Ｓ’について求めた粘度が、くりこみ前の粒子系Ｓについて求めた粘度を連ねる曲線に乗っていることがわかる。このように、くりこみを行った粒子系Ｓ’においても、粘度のひずみ速度依存性が保たれている。
【００８４】
図７に示したシミュレーション結果により、上記実施例によるくりこみ変換則を用いたシミュレーション方法を用いて、高分子材料のマクロな計算を行うことが可能であることが示された。
【００８５】
上記実施例によるシミュレーション方法は、コンピュータがコンピュータプログラムを実行することにより実現される。コンピュータプログラムは、例えばデータ記録媒体に記録された状態で提供される。または、コンピュータプログラムは、電気通信回線を介して提供することも可能である。
【００８６】
図８に、実施例によるシミュレーション装置２５のブロック図を示す。シミュレーション装置２５は、入力装置２１、処理ユニット１０、及び出力装置２２を含む。処理ユニット１０は、入力部１１、演算部１２、出力部１３、及び記憶装置１４を含む。
【００８７】
入力装置２１を通して入力部１１にシミュレーションの種々の条件が入力される。例えば、シミュレーション対象の粒子系Ｓを特徴づけるパラメータである空間の次元数ｄ、モノマー粒子の質量ｍ、１つのポリマーを構成しているモノマー粒子の数Ｎｍ、ポリマーの数Ｎｐの値、及び粒子系Ｓの初期条件が入力される、
【００８８】
記憶装置１４に、シミュレーションを実行するコンピュータプログラムが格納されている。演算部１２は、入力部１１に入力されたデータに基づいて、記憶装置１４に格納されているコンピュータプログラムを実行することにより、シミュレーション処理を実行する。シミュレーション結果が、出力部１３を介して出力装置２２に出力される。
【００８９】
演算部１２は、入力された粒子系Ｓのパラメータに対してくりこみ変換処理を行うことにより、くりこまれた粒子系Ｓ’のパラメータを算出する。算出されたパラメータを用いて、粒子系Ｓ’について分子動力学計算を行う。くりこまれた粒子系Ｓ’について分子動力学計算を行うため、計算時間を短縮することができる。
【００９０】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００９１】
１０  処理ユニット
１１  入力部
１２  演算部
１３  出力部
１４  記憶装置
２１  入力装置
２２  出力装置
２５  シミュレーション装置
【図１】
0006651254.tif 000069
【図２】
0006651254.tif 000070
【図３】
0006651254.tif 000071
【図４】
0006651254.tif 000072
【図５】
0006651254.tif 000073
【図６】
0006651254.tif 000074
【図７】
0006651254.tif 000075
【図８】
0006651254.tif 000076