JP5241468B2 - シミュレーション方法及びプログラム - Google Patents

Description

本発明は、シミュレーション方法及びプログラムに関し、特に、分子動力学を用いたシミュレーション方法、及び、それをコンピュータに実行させるためのプログラムに関する。

分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が、数値的に解かれる。シミュレーション対象の系の粒子数が増大すれば、必要な計算量が増大する。現行のコンピュータの演算能力では、例えば、粒子数１０万程度の系のシミュレーションまでしか行うことができない。

シミュレーションに必要な計算量を減らすために、例えば非特許文献１において、分子動力学に繰り込みの手法を応用しようとする試みがなされている。

なお、本願に関連する発明が本願発明者らによってなされ、特許文献１に開示されている。

特開２００６−２８５８６６号公報 稲村，武澤，社本，「繰り込みの手法に基づく可変スケール分子動力学シミュレーションについて」，日本機械学会論文集（Ａ編），１９９７年，第６３巻，６０８号，ｐ．２０２−２０７

非特許文献１の「３．数値計算法」の欄に記載されているように、非特許文献１のシミュレーション方法では、温度を陽に評価できない。

本発明の一目的は、分子動力学を用いた新規なシミュレーション方法、及びその方法をコンピュータで実行するためのプログラムを提供することである。

本発明の他の目的は、繰り込み群の手法に基づき、温度を陽に評価することが可能なシミュレーション方法、及びその方法をコンピュータで実行するためのプログラムを提供することである。

本発明の一観点によれば、
（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、１より大きい第１の繰り込み因子αと、０以上ｄ以下の第２の繰り込み因子γと、０以上の第３の繰り込み因子δとを決定する工程と、
（ａ１）前記粒子系Ｓが配置されている空間の次元数ｄ、工程（ａ）で決定された第１の繰り込み因子α、第２の繰り込み因子γ、及び第３の繰り込み因子δを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／α^ｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍα^δ／α^γにより求め、繰り込まれた相互作用係数ε´を、変換式ε´＝εα^γにより求める工程と、
（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、繰り込まれた相互作用係数ε´との積ε´ｆで表される粒子系Ｓ´について、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
（ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系
Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、前記工程（ｂ）における分子動力学計算におけるある時間間隔をｔ´と表すとき、前記工程（ａ）で決定された、前記第１の繰り込み因子α、第２の繰り込み因子γ、及び第３の繰り込み因子δのうちの少なくとも１つを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´αにより求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／α ^δ／２ により求める計算、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´α ^{（−γ＋δ／２）} により求める計算、及び、前記粒子系Ｓについて実行すると仮定した分子動力学計算におけるある時間間隔ｔを、変換式ｔ＝ｔ´α ^{（１＋γ−δ／２）} により求める計算のうち、少なくとも１つの計算を実行する工程と
を有するシミュレーション方法が提供される。

また、本発明の他の観点によれば、
（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、１より大きい第１の繰り込み因子αと、０以上の第３の繰り込み因子δとを決定する工程と、
（ａ１）前記粒子系Ｓが配置されている空間の次元数ｄ、前記工程（ａ）で決定された第１の繰り込み因子α、及び第３の繰り込み因子δを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／α^ｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍα^δにより求める工程と、
（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
（ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表すとき、前記工程（ａ）で決定された前記第１の繰り込み因子αと前記第３の繰り込み因子δの少なくとも一方を用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´α により求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／α ^δ／２ により求める計算、及び前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´α ^δ／２ により求める計算のうち、少なくとも１つの計算を実行する工程と
を有するシミュレーション方法が提供される。

更に、本発明の他の観点によれば、
（ｄ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、前記粒子系Ｓが配置されている空間にＸＹＺ直交座標系を考えたとき、それぞれが１より大きいＸ方向の繰り込み因子α_Ｘ、Ｙ方向の繰り込み因子α_Ｙ、及びＺ方向の繰り込み因子α_Ｚ 、及び０以上の繰り込み因子δを決定する工程と、
（ｄ１）前記工程（ｄ）で決定された繰り込み因子α _Ｘ 、α _Ｙ 、α _Ｚ 、及びδを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／（α_Ｘα_Ｙα_Ｚ）により求め、Ｘ方向に関する繰り込まれた粒子の質量ｍ_Ｘ´を、変換式ｍ_Ｘ´＝ｍα_Ｘ ^δにより求め、Ｙ方向に関する繰り込まれた粒子の質量ｍ_Ｙ´を、変換式ｍ_Ｙ´＝ｍα_Ｙ ^δにより求め、Ｚ方向に関する繰り込まれた粒子の質量ｍ_Ｚ´を、変換式ｍ_Ｚ´＝ｍα_Ｚ ^δにより求める工程と、
（ｅ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、Ｘ方向の運動方程式にはＸ方向に関する繰り込まれた粒子の質量ｍ_Ｘ´を用い、Ｙ方向の運動方程式にはＹ方向に関する繰り込まれた粒子の質量ｍ_Ｙ´を用い、Ｚ方向の運動方程式にはＺ方向に関する繰り込まれた粒子の質量ｍ_Ｚ´を用いて、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
（ｆ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、該位置ベクトルｑ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｑ _Ｘ ´、ｑ _Ｙ ´、及びｑ _Ｚ ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、該運動量ベクトルｐ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｐ _Ｘ ´、ｐ _Ｙ ´、及びｐ _Ｚ ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、該速度ベクトルｖ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｖ _Ｘ ´、ｖ _Ｙ ´、及びｖ _Ｚ ´と表すとき、前記Ｘ方向の繰り込み因子α _Ｘ 、Ｙ方向の繰り込み因子α _Ｙ 、Ｚ方向の繰り込み因子α _Ｚ 、及び前記工程（ｄ）で決定された繰り込み因子δを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑのＸ、Ｙ及びＺ方向それぞれの成分ｑ _Ｘ 、ｑ _Ｙ 、及びｑ _Ｚ を、それぞれ、変換式ｑ _Ｘ ＝ｑ _Ｘ ´α _Ｘ 、ｑ _Ｙ ＝ｑ _Ｙ ´α _Ｙ 、及びｑ _Ｚ ＝ｑ _Ｚ ´α _Ｚ により求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐのＸ、Ｙ及びＺ方向それぞれの成分ｐ _Ｘ 、ｐ _Ｙ 、及びｐ _Ｚ を，それぞれ、変換式ｐ _Ｘ ＝ｐ _Ｘ ´／α _Ｘ ^δ／２ 、ｐ _Ｙ ＝ｐ _Ｙ ´／α _Ｙ ^δ／２ 、及びｐ _Ｚ ＝ｐ _Ｚ ´／α _Ｚ ^δ／２ により求める計算、及び、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖのＸ、Ｙ及びＺ方向それぞれの成分ｖ _Ｘ 、ｖ _Ｙ 、及びｖ _Ｚ を，それぞれ、変換式ｖ _Ｘ ＝ｖ _Ｘ ´α _Ｘ ^δ／２ 、ｖ _Ｙ ＝ｖ _Ｙ ´α _Ｙ ^δ／２ 、及びｖ _Ｚ ＝ｖ _Ｚ ´α _Ｚ ^δ／２ により求める計算のうち、少なくとも１つの計算を実行する工程と
を有するシミュレーション方法が提供される。

繰り込まれた粒子数Ｎ´は、粒子数Ｎよりも少ない。このため、粒子系Ｓに対して分子動力学計算を実行するよりも少ない計算量で、粒子系Ｓ´に対する分子動力学計算を実行することができる。

また、０以上の繰り込み因子δを任意に決定して計算を行うため、シミュレーションを柔軟に実行することができる。

まず、分子動力学（Ｍｏｌｅｃｕｌａｒ Ｄｙｎａｍｉｃｓ；ＭＤ）について簡潔に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが、

・・・（１）

と表される粒子系について考える。ここで、各粒子の質量がｍであり、粒子ｊの運動量ベクトルがｐ_ｊであり、粒子ｊの位置ベクトル（位置座標）がｑ_ｊであり、粒子間の相互作用ポテンシャルエネルギがφである。

ハミルトニアンＨを、ハミルトンの正準方程式に代入することにより、各粒子に対する運動方程式

・・・（２）

及び

・・・（３）

が得られる。ここで式（３）において、ｖ_ｊは粒子ｊの速度ベクトルである。分子動力学では、粒子系を構成する各粒子に対し、式（２）及び（３）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトル（または速度ベクトル）及び位置ベクトルが求められる。なお、数値積分には、Ｖｅｒｌｅｔ法が用いられることが多い。Ｖｅｒｌｅｔ法については、例えば、Ｊ．Ｍ．Ｔｈｉｊｓｓｅｎ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ Ｐｈｙｓｉｃｓ”，Ｃａｍｂｒｉｄｇｅ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ（１９９９）のｐ．１７５に説明されている。分子動力学計算で得られた各粒子の位置、速度に基づき、系の種々の物理量を算出することができる。

次に、繰り込み群の手法を応用した分子動力学（これを、繰り込み群分子動力学と呼ぶこととする）について、概念的に説明する。

繰り込み群分子動力学では、物理量を求めたい所望の系Ｓを、系Ｓよりも少ない粒子数で構成される系（これを、繰り込まれた系Ｓ´と呼ぶこととする）に対応付ける。次に、繰り込まれた系Ｓ´に対する分子動力学計算を実行する。そして、繰り込まれた系Ｓ´に対する計算結果を、所望の系Ｓに対応付ける。これにより、所望の系Ｓに対する分子動力学計算を直接的に実行するよりも少ない計算量で、所望の系Ｓに対する物理量を算出することが可能となる。所望の系Ｓにおける物理量（例えば、粒子数、各粒子の質量等）と、繰り込まれた系Ｓ´における物理量とを互いに対応付けるための変換則を、スケール変換則と呼ぶ。

次に、本願発明者が見出したスケール変換則について説明する。物理量を算出したい所望の粒子系Ｓが、Ｎ個の粒子から構成されるとする。粒子系Ｓの全ハミルトニアンＨが、式（１）に示したように、

と表されるとする。ある粒子ｊのハミルトニアンＨ_ｊは、

・・・（４）

と表される。

繰り込みの議論のために、相互作用ポテンシャルエネルギφを、粒子間距離に対する依存性を表し、無次元化された関数ｆと、相互作用の強さを表し、エネルギの次元を持つ係数ε（これを相互作用係数εと呼ぶこととする）との積により、以下の様に表す。

・・・（５）

例えば、不活性原子に対しては、

・・・（６）

という相互作用ポテンシャルが用いられる。ここでｒは原子間距離である。相互作用係数ε及び定数σが、原子の種類に対応する値を取る。

次に、ハミルトニアンＨから、繰り込まれたハミルトニアンＨ´を導出する方法について説明する。以下詳細に示すように、繰り込まれたハミルトニアンＨ´は、粒子系Ｓに対する分配関数Ｚ（β）における積分の一部を実行し、ハミルトニアンを粗視化し、続いて積分変数を再スケールすることで得られる。

粒子数一定における正準集団（Ｃａｎｏｎｉｃａｌ ｅｎｓｅｍｂｌｅ）に対する分配関数Ｚ(β)は、粒子系Ｓに対して、

・・・（７）

と表される。ここで、係数βは、系の温度Ｔとボルツマン定数ｋ_Ｂとによりβ＝１／（ｋ_ＢＴ）と定義される。ｄΓ_Ｎは位相空間内の体積要素であり、より詳細には、

・・・（８）

と表される。ここで、Ｗ_Ｎは、Ｎ！ｈ^３Ｎ（ｈはプランク定数）である。また、Ｄ_ｐ ^Ｎ及びＤ_ｑ ^Ｎをそれぞれ、

と定義している。

調和振動子以外のポテンシャルの繰り込みを議論することは困難である。そこで、鞍点を持つ相互作用ポテンシャルを３つの領域に分けて考える。粒子間距離ｒ→∞とｒ→０は、固定点であるから繰り込み変換に際しポテンシャルは不変である。そこで、鞍点近傍の繰り込み変換がすべての粒子間距離ｒについても成立すると仮定する。鞍点の周りをテーラー展開し、２次の揺らぎまで残す。鞍点の位置をｒ_０とし、相対変位をδｕとすると、δｕの１次の項は消えて、

・・・（９）

を得る。従って、粒子ｊについて、鞍点近傍を記述するハミルトニアンは、

・・・（１０）

となる。ここで、ｕ_ｉ及びｕ_ｊは、それぞれ、粒子ｉ及びｊの鞍点からの変位を表す。φ″は、粒子間距離ｒについての２階微分である。

次に、分配関数Ｚ（β）の粒子間相互作用部分の粗視化について説明する。まず、一次元鎖上に配置された粒子系における粗視化について考察した後、単純立方格子上に配置された粒子系へ拡張する。

図１（Ａ）に示すように、１次元鎖上の格子点（ラティスポイント）に、粒子ｉとｊとが互いに最隣接して配置され、粒子ｊとｋとが互いに最隣接して配置されている。粒子ｉとｋとの中間に、粒子ｊが配置されている。最隣接する粒子同士の相互作用（二アレストネイバーカップリング）を実線で示す。粒子ｊを消去し、粗視化することを考える。粒子ｊが関与する相互作用を書き出し、粒子ｉとｋとの中間にある粒子ｊの変位ｕ_ｊについて積分を実行すれば、粒子ｊの影響が繰り込まれた後の、粒子ｉとｋとの相互作用ポテンシャルφ´（ｑ_ｉ−ｑ_ｋ）について、

・・・（１１）

という式を得る。

この粗視化により、格子定数がα（＝２）倍になるから、ｕ´＝ｕ／αという式でスケール変換された変数ｕ´により、元のポテンシャル関数と同形に表すことができる。つまり、以下の相似則が得られる。

・・・（１２）

ただし、ここで、式（１２）の２行目へ移るに際して、鞍点近傍以外においても同じ繰り込み変換が成立すると仮定した。αを、第１の繰り込み因子と呼ぶこととする。

ｄ次元単純立方格子における粗視化は、ポテンシャルムービング（Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇ）の方法により実現できる。Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇについては、例えば、Ｌｅｏ Ｐ．Ｋａｄａｎｏｆｆ，“ＳＴＡＴＩＳＴＩＣＡＬ ＰＨＹＳＩＣＳ，Ｓｔａｔｉｃｓ，Ｄｙｎａｍｉｃｓ ａｎｄ Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ”，Ｗｏｒｌｄ Ｓｃｉｅｎｔｉｆｆｉｃ（１９９９） の１４章に解説されている。

図１（Ｂ）を参照して、Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇの考え方について説明する。Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇでは、多次元格子を１次元鎖へ帰着させる。図１（Ｂ）には、２次元格子の場合を示す。

図１（Ｂ）の上段の図に示すように、２次元正方格子の格子点上に粒子が配置されている。最隣接する粒子同士の相互作用（最隣接間相互作用）を実線で示す。格子の有る辺に平行な方向をＸ方向とし、Ｘ方向に直交する方向をＹ方向とする。

図１（Ｂ）の中段の図に示すように、Ｘ方向に並んだ粒子について、１つ置きに粒子を積分することを考える（白抜きの丸で示した格子点（ラティスポイントフォーサメイション）上の粒子を積分することを考える）。積分したい（消去したい）粒子を、被積分粒子と呼ぶこととする。被積分粒子同士の最隣接間相互作用を波線で示す。もし被積分粒子同士の最隣接間相互作用が存在しなければ、被積分粒子については、Ｘ方向に延在する１次元鎖と見なせる。

Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇでは、図１（Ｂ）の下段の図に示すように、被積分粒子同士の最隣接間相互作用を、Ｘ方向について被積分粒子の両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を、二重線で示す。二重線で示す相互作用は、元の最隣接間相互作用の２倍の強さになる。この様にして、被積分粒子に対して、２次元格子から１次元鎖に変換できる。そこで、被積分粒子に対して、上述した１次元鎖に対する粗視化の方法が実行できる。例えば３次元格子の場合、以上説明したのと同様な手順を、Ｙ、Ｚ方向についても繰り返す。このようにして、多次元格子における粗視化が実行される。

Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇの方法を用いることにより、ｄ次元単純立方格子における分配関数Ｚ（β）の相互作用部分の粗視化について、

・・・（１３）

という結果が導かれる。ここで、繰り込まれた粒子数Ｎ´及び第１の繰り込み因子αが、

という式で表される。ｎは正整数であり、粗視化の回数を表す。

次に、分配関数Ｚ（β）の運動エネルギ部分の粗視化について説明する。まずは、１次元鎖において考える。互いに最隣接する粒子対である粒子ｊ−１と粒子ｊとを、１つの粒子へ粗視化する。Ｋａｄａｎｏｆｆのブロックスピン法の考え方を運動量に適用し、粒子ｊと粒子ｊ−１とからなる粗視化された粒子の運動量ｇ_ｊを以下の様に定義する。

また、重み関数Ｔ｛Ｐ，ｐ｝を以下の様に定義する。

ここで、δはＤｉｒａｃのδ-関数である。明らかに、

を満たす。この条件から、分配関数の運動エネルギ部分Ｚ_ｐへ重み関数Ｔ｛Ｐ，ｐ｝を挿入でき、

・・・（１４）

という式を得る。δ-関数を積分表示すれば、分配関数の運動エネルギ部分Ｚ_ｐは、以下の様に表せる。

・・・（１５）

この積分は、以下の公式を使えば容易に実行できる。

ここで、＜Ａ，Ｂ＞は、ベクトルＡとＢとの内積を表す。この公式を用いると、分配関数の運動エネルギ部分Ｚ_ｐは、

・・・（１６）

と表すことができる。

この粗視化により、運動エネルギがα（＝２）倍になる。なお、２つの粒子の全運動エネルギを重心の運動エネルギと相対運動の運動エネルギへ分け、相対運動量について積分を実行した場合も同様な結果が得られる。

次に、ｄ次元単純立方格子における粗視化について考察する。まず、結晶軸ａ方向について粗視化すると、粒子の持つ運動エネルギＥ_ｋｉｎはα倍になる。引き続いてｂ軸方向において粗視化すると、粒子の持つ運動エネルギＥ_ｋｉｎαは、（Ｅ_ｋｉｎα）のα倍、即ちＥ_ｋｉｎα^２になる。従って、ｄ次元格子ではＥ_ｋｉｎα^ｄになるから、

・・・（１７）

という関係式が得られる。ｎ回の粗視化では、α＝２^ｎとなることは自明である。

以上の議論から、粗視化された全ハミルトニアンに対する分配関数Ｚ（β）は、結果に影響しない係数を除いて、

・・・（１８）

と表すことができる。

次に、粗視化された系においてもハミルトニアンを同形とするために、“繰り込まれた”変数を導入する。以下のように、繰り込まれた位置ベクトルｑ´、繰り込まれた運動量ベクトルｐ´、繰り込まれた質量ｍ´、繰り込まれた相互作用係数ε´、及び、繰り込まれた係数β´を定義する。ハミルトニアンの形が不変となるようにスケール変換する。

ここで、γは、０以上ｄ以下である。γを、第２の繰り込み因子と呼ぶこととする。また、δは０以上の値をとる。δを第３の繰り込み因子と呼ぶ。

繰り込まれた変数をこのように定めることにより、ハミルトニアンＨと同形な、繰り込まれたハミルトニアンＨ´が得られる。繰り込まれたハミルトニアンＨ´は、

・・・（１９）

と表される。

なお、式（１８）から、分配関数は以下の様になる。

・・・（２０）

ここで、位相空間の体積素は、結果に影響しない定数を除いて、

・・・（２１）

である。

繰り込まれたハミルトニアンＨ´で記述される系（繰り込まれた粒子系Ｓ´）における運動方程式は、正準方程式

・・・（２２）

及び

・・・（２３）

から与えられ、

・・・（２４）

及び

・・・（２５）

となる。ここで、繰り込まれた時間ｔ´が、繰り込まれた変数ｑ´、ｐ´、ｍ´及びε´に整合するように、

という式で定められる。運動方程式（２４）及び（２５）に基づいて、繰り込まれた粒子系Ｓ´に対する分子動力学計算を実行することができる。

以上の考察をまとめる。粒子系Ｓと、繰り込まれた粒子系Ｓ´とは、以下の変換式で表されるスケール変換則により対応付けられる。

粒子系Ｓの粒子数Ｎと、繰り込まれた粒子系Ｓ´の繰り込まれた粒子数Ｎ´とは、変換式

・・・（変換式Ｎ）

により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑと、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´とは、変換式

・・・（変換式ｑ）

により対応付けられる。粒子系Ｓのある粒子の運動量ベクトルｐと、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´とは、変換式

・・・（変換式ｐ）

により対応付けられる。粒子系Ｓの各粒子の質量ｍと、繰り込まれた粒子系Ｓ´の各粒子の質量ｍ´とは、変換式

・・・（変換式ｍ）

により対応付けられる。粒子系Ｓにおける相互作用係数εと、繰り込まれた粒子系Ｓ´における相互作用係数ε´とは、変換式

・・・（変換式ε）

により対応付けられる。

粒子系Ｓに対して運動方程式を解く場合の時間スケールｔ（粒子系Ｓに対して実行すると仮定した分子動力学計算におけるある時間間隔ｔ）と、繰り込まれた粒子系Ｓ´に対して運動方程式を解く場合の時間スケールｔ´（繰り込まれた粒子系Ｓ´に対して実行する分子動力学計算におけるある時間間隔ｔ´）とは、変換式

・・・（変換式ｔ）

により対応付けられる。

なお、粒子系Ｓのある粒子の速度ベクトルｖと、繰り込まれた粒子系Ｓ´のある粒子の速度ベクトルｖ´とは、速度ベクトルｖがｖ＝ｐ／ｍと表されることから、変換式

・・・（変換式ｖ）

により対応付けられる。

ここで、粒子系Ｓの配置されている空間の次元数がｄである。第１の繰り込み因子αは２^ｎ（ｎは正整数）であり、第２の繰り込み因子γは０以上ｄ以下である。また、第３の繰り込み因子δは０以上である。

第２の繰り込み因子γの値が大きく、例えばｄに等しいとき、極めて大きな、繰り込まれた相互作用係数ε´と、小さな、繰り込まれた質量ｍ´とが数値計算に現れる。数値計算上、第２の繰り込み因子γは０に設定することが好ましい。このとき更に、第３の繰り込み因子δを０とすれば、時刻ｔのみが変換され、従来のＭＤが利用できる。

なお、たとえばδ＝２のとき、質量ｍとバネ定数ｓを持つ１次元鎖のバネ・マス系を考えたとき、繰り込まれたハミルトニアンから分散関係が以下の様に得られる。

・・・（２６）

ここで、角周波数がωであり、波数がｋである。平衡状態における原子間距離がａである。なお、第２の繰り込み因子γは０とした。式（２６）は、長波長の極限において、

・・・（２７）

となり、第１の繰り込み因子αには無関係になる。

次に、繰り込み群分子動力学によるコンピュータシミュレーションの実施例について説明する。

まず、アルミニウムに熱を流入させて昇温させ、融点及び潜熱を算出したシミュレーションについて説明する。第１のシミュレーションでは、比較のため、１６０００個のアルミニウム原子からなる系に対して、従来の（繰り込みを応用しない）分子動力学計算を行った。アルミニウム原子は３次元空間内に配置されており、次元数ｄは３である。

第２及び第３のシミュレーションでは、第１のシミュレーションに対応する系を繰り込んだ系に対する分子動力学計算を行った。第１の繰り込み因子αは、第２のシミュレーションでは２と決定し、第３のシミュレーションでは２^２（＝４）と決定した。（変換式Ｎ）より、第２のシミュレーションでは、系の粒子数が２０００個となり、第３のシミュレーションでは、系の粒子数が２５０個となる。第２の繰り込み因子γは０と決定した。また、第３の繰り込み因子δは２と決定した。なお、従来の手法（第１のシミュレーション）は、第１の繰り込み因子αを１と設定した場合に対応する。

第１のシミュレーションにおける原子間ポテンシャル（つまり、元の系における粒子間ポテンシャル）として、モース（Ｍｏｒｓｅ）ポテンシャルを用いた。このポテンシャルは、

・・・（２８）

と表される。ここで、相互作用係数εは１．９２×１０^−２０[Ｊ]であり、定数Ａは２．３５×１０^１０[ｍ^−１]であり、定数ｒ_０は２．８６×１０^−１０[ｍ]である。粒子間距離がｒである。

第２及び第３のシミュレーションにおける粒子間ポテンシャル（つまり、繰り込まれた系における粒子間ポテンシャル）としても、式（２８）と同形のものが用いられる。第２及び第３のシミュレーションにおける繰り込まれた相互作用係数は、それぞれ、（変換式ε）より求められる。第２の繰り込み因子γが０であるので、第２及び第３のシミュレーションでの相互作用係数は、双方とも、第１のシミュレーションの相互作用係数εと等しい。定数ｒ_０及び定数Ａは第１のシミュレーションのそれと等しい。

第１のシミュレーションにおけるアルミニウム原子の質量（つまり、元の系における粒子の質量）は、４．４８×１０^−２６[ｋｇ]である。第２及び第３のシミュレーションにおける粒子の繰り込まれた質量は、それぞれ、（変換式ｍ）より求められる。第２のシミュレーションにおける質量は、第１のシミュレーションでの質量に２^２（＝４）を乗じた１７．９×１０^−２６[ｋｇ]となり、第３のシミュレーションにおける質量は、第１のシミュレーションでの質量に４^２（＝１６）を乗じた７１．７×１０^−２６[ｋｇ]となる。

次に、粒子の初期配置について説明する。初期において、粒子は面心立方格子の格子点上に配置される。第１のシミュレーションでは、幅方向に２０個、奥行き方向に２０個の粒子が並んだ層が、４０層重なる配置とした。元の系の、幅方向について１／α、奥行き方向について１／α、及び高さ方向について１／αの部分（元の系の１／α^３の部分）が、繰り込まれた系に対応する。第２のシミュレーションでは、幅方向に１０個、奥行き方向に１０個の粒子が並んだ層が、２０層重なった配置となり、第３のシミュレーションでは、幅方向に５個、奥行き方向に５個の粒子が並んだ層が、１０層重なった配置となる。

次に、数値積分の時間刻み幅について説明する。元の系における時間刻み幅をΔｔとし、粒子の速さをｖとする。時間刻み幅Δｔは、例えば

を満たすように定められる（「＜＜」は、式の左辺が右辺より充分小さいという意）。ここでｒ_０は粒子間ポテンシャルにおける定数であり、粒子間距離の目安を与える。繰り込まれた系における時間刻み幅をΔｔ´とし、粒子の速さをｖ´とする。（変換式ｖ）よりｖ´＝ｖ／αであるので、繰り込まれた系における時間刻み幅Δｔ´は、例えば、αΔｔに設定される。

数値積分の時間刻み幅は、第１のシミュレーションにおいては５．０[ｆｓ]とし、第２のシミュレーションにおいては１０[ｆｓ]とし、第３のシミュレーションにおいては２０[ｆｓ]とした。

次に、粒子に与える初期速度について説明する。まず、元の系の温度と、繰り込まれた系の温度との関係について説明する。元の系の温度Ｔは、

・・・（２９）

という式で定義される。ここで、ｋ_Ｂはボルツマン定数であり、ｖ_ｉは粒子ｉの速度ベクトルであり、ｖ_ｇは粒子系の重心の速度ベクトルである。温度は、つまり、速度の“揺らぎ”である。繰り込まれた系の温度Ｔ´は、

・・・（３０）

という式から定まる。

元の系において、質量がｍで速さがｖの粒子の運動エネルギは、１／２ｍｖ^２である。繰り込まれた系において、質量がｍ´で速さがｖ´の粒子の運動エネルギは、１／２ｍ´ｖ´^２である。（変換式ｍ）及び（変換式ｖ）から、

・・・（３１）

という関係が得られる。よって、元の系の温度Ｔと、繰り込まれた系の温度Ｔ´とは等しくなる。第１〜第３のシミュレーションにおいて、粒子系の温度が３００[Ｋ]となるように、初期速度が与えられる。

次に、粒子系に熱を流入させる方法について説明する。第１〜第３のシミュレーションでは、粒子系の底面上に配置された粒子に運動エネルギを与えて、系への熱の流入を模擬した。熱の流入する面を加熱面と呼ぶこととする。

まず、第１のシミュレーションにおける熱の流入方法（つまり、元の系における熱の流入方法）ついて説明する。加熱面上の粒子１個当りに、時間刻み幅Δｔの間に、運動エネルギΔＥが与えられるとする。粒子の運動量ベクトルの変化は、下式（３２）より求められる。ここで、ある時刻ｔにおける運動量ベクトルがｐ_ｔであり、時刻ｔ＋Δｔにおける運動量ベクトルがｐ_ｔ＋Δｔである。

・・・（３２）

運動エネルギΔＥによって、運動量ベクトルがλ倍されるとする。つまり、

・・・（３３）

と表されるとする。式（３３）を式（３２）に代入することにより、λは、

・・・（３４）

と表される。時間刻み幅Δｔごとに、式（３３）により運動量ベクトルを更新することにより、時間刻み幅Δｔ当たりに運動エネルギΔＥが流入するような熱の流入を模擬できる。

加熱面上の粒子数をＮ_ｈとし、加熱面の面積をＳ_ｈとすると、単位面積当たりに流入するパワーｆ_ｈは、

・・・（３５）

と表される。第１のシミュレーションでは、加熱面上の粒子数Ｎ_ｈが４００個であり、加熱面の面積Ｓ_ｈが２．８３×１０^−１７[ｍ^２]であり、運動エネルギΔＥが４．１４×１０^−２４[Ｊ]であり、時間刻み幅Δｔが５．０[ｆｓ]であり、パワーｆ_ｈが１．１７×１０^１０[Ｗ／ｍ^２]である。

第２及び第３のシミュレーションにおいて（繰り込まれた系において）、単位面積当たりに流入するパワーｆ_ｈ´は、

・・・（３６）

という式で与えられる。運動量ベクトルを更新するための係数λ´は、

・・・（３７）

という式で与えられる。なお式（３１）から、与えられる運動エネルギΔＥは元の系と繰り込まれた系とで等しくなる。

次に、図２を参照して、シミュレーション結果について説明する。図２に示すグラフのうち、下段のグラフが第１のシミュレーション（従来の分子動力学計算）の結果を表し、中段のグラフが第２のシミュレーション（繰り込み因子αを２とした繰り込み群分子動力学計算）の結果を表し、上段のグラフが第３のシミュレーション（繰り込み因子αを４とした繰り込み群分子動力学計算）の結果を表す。３つのグラフにおいて、縦軸が温度（単位[Ｋ]）を示し、横軸が流入した熱量、つまり流入したエネルギ（単位１０^−１６[Ｊ]）を示す。

まず、下段のグラフを参照して、第１のシミュレーションの結果について説明する。系にエネルギが流入していないとき、系の温度は３００[Ｋ]である。流入したエネルギが増加するにつれ、系の温度が上昇する。流入エネルギが約５×１０^−１６[Ｊ]に達すると、温度の上昇が、約８５０[Ｋ]で一旦止まり、流入エネルギが約９×１０^−１６[Ｊ]に達するまで、温度がほぼ一定となる。この温度が、融点を示す。流入エネルギが約９×１０^−１６[Ｊ]を超えると、再び温度が上昇する。

温度一定の期間中に系に流入したエネルギが潜熱に対応する。このシミュレーションから求められる潜熱は、１４．７[ｋＪ／ｍｏｌ]となる。なお、実験から求められたアルミニウムの融点は９３３[Ｋ]であり、潜熱は１１[ｋＪ／ｍｏｌ]である。シミュレーションで得られた値と実験値との相違は、用いた原子間ポテンシャルの精度等に起因すると考えられる。

次に、中段及び上段のグラフを参照して、第２及び第３のシミュレーションの結果について説明する。第２及び第３のシミュレーションの結果も、第１のシミュレーションの結果と同様な傾向を示す。すなわち、流入エネルギが約５×１０^−１６[Ｊ]から約９×１０^−１６[Ｊ]までの期間中に、温度が８５０[Ｋ]付近でほぼ一定となる傾向を示す。これにより、第２及び第３のシミュレーションからも、融点が約８５０[Ｋ]であり、潜熱が１４．７[ｋＪ／ｍｏｌ]であると求めることができる。

粒子系に流入したエネルギは、シミュレーションした期間分だけ流入パワーを積分すれば求められる。繰り込まれた系に流入したエネルギをＥ´とすると、元の系へ流入したエネルギＥは、

・・・（３８）

という式から求められる。なお、粒子１つ当たりの流入エネルギは、元の系と繰り込まれた系とで等しい。

繰り込まれた系の温度は、式（３０）から求められる。上述の議論より、元の系の温度は、繰り込まれた系の温度と等しい。

なお、上述のシミュレーションでは、元の系における流入エネルギ及び温度を求めたが、元の系における他の物理量、例えば粒子の位置、運動量、速度等も、必要に応じて、上述のスケール変換則に基づいて求めることができる。

以上説明したように、繰り込み群分子動力学では、繰り込む前の系（上記実施例の第１のシミュレーションに対応する系）に対する分子動力学計算で得られると期待される物理量（上記実施例では融点及び潜熱）を、繰り込む前の系よりも粒子数の少ない系に対する分子動力学計算を行うことにより、算出することが可能となる。繰り込み群分子動力学によるシミュレーションに要する種々の計算は、プログラムにより、コンピュータで実行することができる。

なお、上記実施例のシミュレーションにおいて、アルミニウムは面心立方構造を取る。スケール変換則は、単純立方格子上に配置された粒子系に対する考察から導出したが、単純立方構造に限らない他の結晶構造に対しても有効である。

次に、シリコンに熱を流入させて昇温させ、融点及び潜熱を算出したシミュレーションについて説明する。シリコンは、ダイヤモンド構造を取る。アルミニウムに対するシミュレーションと同様に、３つのシミュレーションを行った。第１のシミュレーションは、比較のために行った従来の分子動力学手法によるものである。第２及び第３のシミュレーションにおいて、第１のシミュレーションに対応する系を繰り込んだ系に対する分子動力学計算を行った。繰り込み因子αは、第２のシミュレーションでは２とし、第３のシミュレーションでは２^２（＝４）とした。系の粒子数は、第１のシミュレーションでは８１９２個であり、第２のシミュレーションでは１０２４個であり、第３のシミュレーションでは１２８個である。

第１のシミュレーションにおける原子間ポテンシャル（つまり、元の系における粒子間ポテンシャル）として、Ｓ−Ｗ（Ｓｔｉｌｌｉｎｇｅｒ−Ｗｅｂｅｒ）ポテンシャルを用いた。このポテンシャルは、

・・・（３９）

と表される。ここで、

であり、相互作用係数εは、２．１６７[ｅＶ]であり、σは２．０９５[Å]であり、定数Ａ、Ｂ、ｐ、ｑ、ａ、λ、及びγは、それぞれ、７．０５、０．６０２、４．０、０．０、１．８、２１．０、及び１．２である。

また、図３（Ａ）に示すように、粒子ｉ及び粒子ｋの結合と、粒子ｉ及び粒子ｊの結合とのなす角がθ_ｉｊｋであり、粒子ｉから粒子ｊまでの粒子間距離がｒ_ｉｊである。第２及び第３のシミュレーションの粒子間ポテンシャルにおける繰り込まれた相互作用係数は、上述のアルミニウムに対するシミュレーションと同様な方法で求められる。

第１のシミュレーションにおけるシリコン原子の質量（つまり、元の系における粒子の質量）は、４．６７×１０^−２６[ｋｇ]とした。第２及び第３のシミュレーションにおける粒子の繰り込まれた質量は、上述のアルミニウムに対するシミュレーションと同様な方法で求められる。

次に、図３（Ｂ）を参照して、シミュレーション結果について説明する。図３（Ｂ）に示すグラフのうち、下段のグラフが第１のシミュレーション（従来の分子動力学計算）の結果を表し、中段のグラフが第２のシミュレーション（繰り込み因子αを２とした繰り込み群分子動力学計算）の結果を表し、上段のグラフが第３のシミュレーション（繰り込み因子αを４とした繰り込み群分子動力学計算）の結果を表す。

まず、下段のグラフを参照して、第１のシミュレーションの結果について説明する。流入エネルギが約１．３×１０^−１６[Ｊ]から約４．６×１０^−１６[Ｊ]に達するまでの期間中、温度が１７５０[Ｋ]程度でほぼ一定となる。潜熱は、２５．６[ｋＪ／ｍｏｌ]と求められる。なお、実験から求められたシリコンの融点は１６８７[Ｋ]であり、潜熱は５０．６[ｋＪ／ｍｏｌ]である。シミュレーションで得られた値と実験値との相違は、原子間ポテンシャルの精度等に起因すると考えられる。

次に、中段及び上段のグラフを参照して、第２及び第３のシミュレーションの結果について説明する。第２及び第３のシミュレーションの結果も、第１のシミュレーションの結果とほぼ同様な傾向を示す。これにより、第２及び第３のシミュレーションに基づいて、第１のシミュレーションで得られた値に近い融点及び潜熱を得ることができる。

次に、図４を参照して、分散関係を繰り込み群分子動力学で求めたシミュレーションについて説明する。アルミニウムの[１１１]方向に平行に振動を加えて、[１１１]方向の波数を算出した。原子間はバネで接続した。バネ定数ｓは、４０．３[Ｎ／ｍ]とした。なお、バネ定数ｓは、アルミニウムの原子間ポテンシャルの最小近傍を２次関数で近似することにより求められる。アルミニウム原子の質量ｍは、４．４８×１０^−２６[ｋｇ]とした。原子間距離ａは、２．８６×１０^−１０[ｍ]とした。分散関係の厳密解（解析的に求められた解）は、加振の角周波数をωとし、[１１１]方向の波数をｋとして、

という式で与えられる。なお、この厳密解は、式（２６）においてα＝１として求めることができる。
グラフ中の実線Ｌが、厳密解を示す。点Ｐ１が、従来の分子動力学計算（繰り込み因子α＝１の場合に対応）の結果を示す。点Ｐ２〜点Ｐ５が、それぞれ、繰り込み因子αを、２^７、２^１５、２^２１、及び２^２８とした場合の結果を示す。繰り込み因子αを大きくすることは、元の系（繰り込まれる前の系）のサイズを大きくすることに対応する。繰り込み因子αを大きくすることにより、小さい波数を得ることができる。繰り込み群分子動力学で得られた結果は、厳密解と良く一致している。このように、繰り込み群分子動力学を用いれば、従来の分子動力学で実行するのが困難な長波長域の計算（連続体とみなせるようなマクロな系に対する計算）を、精度良く行うことができる。

粒子系が配置されている空間に、ＸＹＺ直交座標系を考える。上述の実施例では、元の系の位置ベクトルｑから、繰り込み因子αを用いて、繰り込まれた位置ベクトルｑ´が、

・・・（変換式ｑ）

と求められた。これは、Ｘ、Ｙ及びＺのすべての方向に関する空間的なスケール変換が、共通の繰り込み因子αにより行われることを示す。繰り込み因子αを、Ｘ、Ｙ及びＺの各方向について独立に設定することも可能である。つまり、異方性を有するスケール変換を行うことも可能である。

第２の繰り込み因子γが０、及び第３の繰り込み因子δが２と決定された場合について、異方性を有するスケール変換則は、以下に示すような変換式で表される。Ｘ方向の繰り込み因子をα_Ｘとし、Ｙ方向の繰り込み因子をα_Ｙとし、Ｚ方向の繰り込み因子をα_Ｚとする。各繰り込み因子α_Ｘ〜α_Ｚは、それぞれ、例えば２^ｎ（ｎは正整数）である。

粒子系Ｓの粒子数Ｎと、繰り込まれた粒子系Ｓ´の繰り込まれた粒子数Ｎ´とは、変換式

・・・（変換式Ｎａ）

により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑの各成分と、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´の各成分とは、変換式

・・・（変換式ｑａ）

により対応付けられる。ここで、ηは、ベクトルのＸ，Ｙ、Ｚのいずれかの成分を示す（以下、（変換式ｐａ）、（変換式ｖａ）、及び（変換式ｍａ）においても同様である。）。粒子系Ｓのある粒子の運動量ベクトルｐの各成分と、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´の各成分とは、変換式

・・・（変換式ｐａ）

により対応付けられる。

なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。

粒子系Ｓのある粒子の速度ベクトルｖの各成分と、繰り込まれた粒子系Ｓ´のある粒子の速度ベクトルｖ´の各成分とは、変換式

・・・（変換式ｖａ）

により対応付けられる。

なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。

方向ごとに繰り込み因子が設定されているので、粒子の質量も方向ごとにスケール変換される。粒子系Ｓの各粒子の質量ｍと、繰り込まれた粒子系Ｓ´の各粒子の、Ｘ方向の質量ｍ´_Ｘ、Ｙ方向の質量ｍ´_Ｙ、及びＺ方向の質量ｍ´_Ｚとは、変換式

・・・（変換式ｍａ）

により対応付けられる。繰り込まれた系における運動方程式（２５）で、方向ごとに、対応する質量が用いられる。

なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。

粒子系Ｓにおける相互作用係数εと、繰り込まれた粒子系Ｓ´における相互作用係数ε´とは、変換式

・・・（変換式εａ）

により対応付けられ、粒子系Ｓにおける時間スケールｔと、繰り込まれた粒子系Ｓ´におけるスケールｔ´とは、変換式

・・・（変換式ｔａ）

により対応付けられる。

例えば、薄膜に対するシミュレーションについて考える。薄膜は、膜厚方向の寸法に比べて、膜面に平行な方向の寸法が大きい。このため。膜厚に平行な方向について、膜厚方向よりも繰り込み因子を大きくするとよい。例えば、膜厚方向をＺ方向とするならば、膜厚方向の繰り込み因子α_Ｚに比べて、膜面に平行な方向の繰り込み因子α_Ｘ及びα_Ｙを大きく設定するとよい。

なお、第２のスケーリングパラメータγは、（変換式ｍ）、（変換式ε）、及び（変換式ｔ）に、α^γという形で含まれている。γ＝０と設定した場合、α^γは単に「１」となる。よって、たとえばδ＝２と決定された場合、粒子の質量、相互作用係数、及び時間に関して、第２のスケーリングパラメータγを用いずに（γという第２のスケーリングパラメータを導入せずに）、

というスケール変換を施したとしても、第２のスケーリングパラメータγを０と設定したのと同様なスケール変換が行われる。なお、γ＝０、δ＝２のとき、速度ベクトルについては、

というスケール変換が行われる。

なお、繰り込み因子α、α_Ｘ、α_Ｙ、及びα_Ｚが、それぞれ１より大きければ、繰り込まれた系における粒子数が元の系の粒子数よりも少なくなる。

なお、必要に応じ、シミュレーション対象の系のある部分と他の部分とで、繰り込み因子（例えばα）を異ならせることも可能である。

なお、上述の実施例では、繰り込み群分子動力学の例として、融点、潜熱及び分散関係を求めた。繰り込み群分子動力学は、これらに限らない種々の物理現象のシミュレーションに用いることができる。繰り込み群分子動力学は、分子動力学を用いるので、分子動力学で扱える様々な問題に適用可能である。例えば、流れ、構造、熱、反応等の問題、さらには、機構、磨耗や潤滑の問題等を扱うことができる。繰り込み群分子動力学は、実際に分子動力学計算を実行する系より粒子数の多い系に対する物理量を求めることを可能にする。よって、本願発明者が見出した繰り込み群分子動力学を用いれば、従来シミュレーションが困難であったスケールの問題に対して、精度良いシミュレーションを行うことが可能となる。なお、繰り込み群分子動力学において実行される分子動力学計算の手法としては、公知の種々のものを用いることができる。

なお、本願明細書の「背景技術」の欄に記載した非特許文献１では、繰り込みの手法を分子動力学に応用する試みがなされている。しかし、非特許文献１では、ハミルトニアンの相互作用部分の繰り込みがなされているのみであり、運動エネルギ部分の繰り込みはなされていない。非特許文献１では、運動量（または速度）に対して、繰り込みではなく、単純平均を取っている。このため、静的な解析に留まっており、温度を陽に考慮できず、強制的にエネルギの散逸を起こさせるための人為的な手法（非特許文献１の「３．数値計算法」の欄に記載の式（２３）、（２４）、（２５））が必要となっている。そのため、有限温度の計算や振動・騒音、相転移が伴う動的な問題へ適用する際、信頼性に欠ける等の課題があった。非特許文献１が開示するスケール変換則は、本願発明のγ＝０、δ＝３の場合に相当する。

本願発明者のスケール変換則の導出過程においては、ハミルトニアンの運動エネルギ部分の繰り込みもなされている。これにより、運動量ベクトル（または速度ベクトル）の適正なスケール変換則が導かれる。よって、シミュレーションにおいて、温度を陽に考慮することが可能となる。

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

図１（Ａ）及び図１（Ｂ）は、スケール変換則の導出方法を説明するための図である。 繰り込み群分子動力学を用いたアルミニウムに対するシミュレーションの結果を示すグラフである。 図３（Ａ）は、Ｓ−Ｗポテンシャルについて説明するための図であり、図３（Ｂ）は、繰り込み群分子動力学を用いたシリコンに対するシミュレーションの結果を示すグラフである。 繰り込み群分子動力学を用いて分散関係を求めたシミュレーションの結果を示すグラフである。

符号の説明

α 第１の繰り込み因子
Ｎ 粒子数
Ｎ´ 繰り込まれた粒子数

Claims (8)

1. （ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、１より大きい第１の繰り込み因子αと、０以上ｄ以下の第２の繰り込み因子γと、０以上の第３の繰り込み因子δとを決定する工程と、
   （ａ１）前記粒子系Ｓが配置されている空間の次元数ｄ、工程（ａ）で決定された第１の繰り込み因子α、第２の繰り込み因子γ、及び第３の繰り込み因子δを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／α^ｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍα^δ／α^γにより求め、繰り込まれた相互作用係数ε´を、変換式ε´＝εα^γにより求める工程と、
   （ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、繰り込まれた相互作用係数ε´との積ε´ｆで表される粒子系Ｓ´について、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
   （ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、前記工程（ｂ）における分子動力学計算におけるある時間間隔をｔ´と表すとき、前記工程（ａ）で決定された、前記第１の繰り込み因子α、第２の繰り込み因子γ、及び第３の繰り込み因子δのうちの少なくとも１つを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´αにより求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／α ^δ／２ により求める計算、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´α ^{（−γ＋δ／２）} により求める計算、及び、前記粒子系Ｓについて実行すると仮定した分子動力学計算におけるある時間間隔ｔを、変換式ｔ＝ｔ´α ^{（１＋γ−δ／２）} により求める計算のうち、少なくとも１つの計算を実行する工程と
   を有するシミュレーション方法。
2. 前記第２の繰り込み因子γが０である請求項１に記載のシミュレーション方法。
3. 前記第１の繰り込み因子が、ｎを正整数として、２^ｎである請求項１または２に記載のシミュレーション方法。
4. 前記次元数ｄが３である請求項１〜３のいずれかに記載のシミュレーション方法。
5. （ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、１より大きい第１の繰り込み因子αと、０以上の第３の繰り込み因子δとを決定する工程と、
   （ａ１）前記粒子系Ｓが配置されている空間の次元数ｄ、前記工程（ａ）で決定された第１の繰り込み因子α、及び第３の繰り込み因子δを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／α^ｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍα^δにより求める工程と、
   （ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
   （ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表すとき、前記工程（ａ）で決定された前記第１の繰り込み因子αと前記第３の繰り込み因子δの少なくとも一方を用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´α により求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／α ^δ／２ により求める計算、及び前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´α ^δ／２ により求める計算のうち、少なくとも１つの計算を実行する工程と
   を有するシミュレーション方法。
6. （ｄ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができるシミュレーション対象の粒子系Ｓに関し、前記粒子系Ｓが配置されている空間にＸＹＺ直交座標系を考えたとき、それぞれが１より大きいＸ方向の繰り込み因子α_Ｘ、Ｙ方向の繰り込み因子α_Ｙ、及びＺ方向の繰り込み因子α_Ｚ 、及び０以上の繰り込み因子δを決定する工程と、
   （ｄ１）前記工程（ｄ）で決定された繰り込み因子α _Ｘ 、α _Ｙ 、α _Ｚ 、及びδを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／（α_Ｘα_Ｙα_Ｚ）により求め、Ｘ方向に関する繰り込まれた粒子の質量ｍ_Ｘ´を、変換式ｍ_Ｘ´＝ｍα_Ｘ ^δにより求め、Ｙ方向に関する繰り込まれた粒子の質量ｍ_Ｙ´を、変換式ｍ_Ｙ´＝ｍα_Ｙ ^δにより求め、Ｚ方向に関する繰り込まれた粒子の質量ｍ_Ｚ´を、変換式ｍ_Ｚ´＝ｍα_Ｚ ^δにより求める工程と、
   （ｅ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、Ｘ方向の運動方程式にはＸ方向に関する繰り込まれた粒子の質量ｍ_Ｘ´を用い、Ｙ方向の運動方程式にはＹ方向に関する繰り込まれた粒子の質量ｍ_Ｙ´を用い、Ｚ方向の運動方程式にはＺ方向に関する繰り込まれた粒子の質量ｍ_Ｚ´を用いて、粒子の初期位置及び初期速度と、時間刻み幅とを与えて、分子動力学計算を実行する工程と、
   （ｆ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、該位置ベクトルｑ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｑ _Ｘ ´、ｑ _Ｙ ´、及びｑ _Ｚ ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、該運動量ベクトルｐ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｐ _Ｘ ´、ｐ _Ｙ ´、及びｐ _Ｚ ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、該速度ベクトルｖ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｖ _Ｘ ´、ｖ _Ｙ ´、及びｖ _Ｚ ´と表すとき、前記Ｘ方向の繰り込み因子α _Ｘ 、Ｙ方向の繰り込み因子α _Ｙ 、Ｚ方向の繰り込み因子α _Ｚ 、及び前記工程（ｄ）で決定された繰り込み因子δを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑのＸ、Ｙ及びＺ方向それぞれの成分ｑ _Ｘ 、ｑ _Ｙ 、及びｑ _Ｚ を、それぞれ、変換式ｑ _Ｘ ＝ｑ _Ｘ ´α _Ｘ 、ｑ _Ｙ ＝ｑ _Ｙ ´α _Ｙ 、及びｑ _Ｚ ＝ｑ _Ｚ ´α _Ｚ により求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐのＸ、Ｙ及びＺ方向それぞれの成分ｐ _Ｘ 、ｐ _Ｙ 、及びｐ _Ｚ を，それぞれ、変換式ｐ _Ｘ ＝ｐ _Ｘ ´／α _Ｘ ^δ／２ 、ｐ _Ｙ ＝ｐ _Ｙ ´／α _Ｙ ^δ／２ 、及びｐ _Ｚ ＝ｐ _Ｚ ´／α _Ｚ ^δ／２ により求める計算、及び、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖのＸ、Ｙ及びＺ方向それぞれの成分ｖ _Ｘ 、ｖ _Ｙ 、及びｖ _Ｚ を，それぞれ、変換式ｖ _Ｘ ＝ｖ _Ｘ ´α _Ｘ ^δ／２ 、ｖ _Ｙ ＝ｖ _Ｙ ´α _Ｙ ^δ／２ 、及びｖ _Ｚ ＝ｖ _Ｚ ´α _Ｚ ^δ／２ により求める計算のうち、少なくとも１つの計算を実行する工程と
   を有するシミュレーション方法。
7. 前記Ｘ方向の繰り込み因子α_Ｘ、Ｙ方向の繰り込み因子α_Ｙ、及びＺ方向の繰り込み因子α_Ｚのうち、少なくとも２つが互いに異なる請求項６に記載のシミュレーション方法。
8. 請求項１〜７のいずれかに記載のシミュレーション方法を、コンピュータで実行するためのプログラム。