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