(57)【要約】      （修正有）【課題】分子動力学を用いた新規なシミュレーション方法を提供する。【解決手段】Ｎ個の粒子を含み各質量がｍの粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表せる粒子系Ｓを考え、粒子系Ｓが配置される空間の次元数ｄを用い、１より大きい第１の繰り込み因子αと０以上ｄ以下の第２の繰り込み因子γと０以上の第３の繰り込み因子δとを決定し、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／αｄで求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍαδ／αγで求め、繰り込まれた相互作用係数ε´を、変換式ε´＝εαγで求める。繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、無次元化関数ｆと、繰り込まれた相互作用係数ε´との積ε´ｆで表される粒子系Ｓ´について、分子動力学計算を実行する。【選択図】図２
2010146368.tif 000002
【特許請求の範囲】
【請求項１】
（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓを考えたとき、該粒子系Ｓが配置されている空間の次元数ｄを用い、１より大きい第１の繰り込み因子αと、０以上ｄ以下の第２の繰り込み因子γと、０以上の第３の繰り込み因子δとを決定して、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／αｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍαδ／αγにより求め、繰り込まれた相互作用係数ε´を、変換式ε´＝εαγにより求める工程と、
（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、繰り込まれた相互作用係数ε´との積ε´ｆで表される粒子系Ｓ´について、分子動力学計算を実行する工程と
を有するシミュレーション方法。
【請求項２】
さらに、
（ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算におけるある時間間隔をｔ´と表すとき、前記工程（ａ）で決定された、前記第１の繰り込み因子α、第２の繰り込み因子γ、及び第３の繰り込み因子δのうちの少なくとも１つを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´αにより求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／αδ／２により求める計算、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´α（−γ＋δ／２）により求める計算、及び、前記粒子系Ｓについて実行すると仮定した分子動力学計算におけるある時間間隔ｔを、変換式ｔ＝ｔ´α（１＋γ−δ／２）により求める計算のうち、少なくとも１つの計算を実行する工程
を有する請求項１に記載のシミュレーション方法。
【請求項３】
前記第２の繰り込み因子γが０である請求項１または２に記載のシミュレーション方法。
【請求項４】
前記第１の繰り込み因子が、ｎを正整数として、２ｎである請求項１〜３のいずれかに記載のシミュレーション方法。
【請求項５】
前記次元数ｄが３である請求項１〜４のいずれかに記載のシミュレーション方法。
【請求項６】
（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓを考えたとき、該粒子系Ｓが配置されている空間の次元数ｄを用い、１より大きい第１の繰り込み因子αと、０以上の第３の繰り込み因子δとを決定して、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／αｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍαδにより求める工程と、
（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、分子動力学計算を実行する工程とを有するシミュレーション方法。
【請求項７】
さらに、
（ｃ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｂ）における分子動力学計算で求められた速度ベクトルをｖ´と表すとき、前記工程（ａ）で決定された前記第１の繰り込み因子αと前記第３の繰り込み因子δの少なくとも一方を用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑを、変換式ｑ＝ｑ´α  により求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐを、変換式ｐ＝ｐ´／αδ／２により求める計算、及び前記粒子系Ｓに含まれるある粒子の速度ベクトルｖを、変換式ｖ＝ｖ´αδ／２により求める計算のうち、少なくとも１つの計算を実行する工程
を有する請求項６に記載のシミュレーション方法。
【請求項８】
（ｄ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓと、前記粒子系Ｓが配置されている空間にＸＹＺ直交座標系を考えたとき、それぞれが１より大きいＸ方向の繰り込み因子αＸ、Ｙ方向の繰り込み因子αＹ、及びＺ方向の繰り込み因子αＺを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／（αＸαＹαＺ）により求め、更に０以上の繰り込み因子δを決定して、Ｘ方向に関する繰り込まれた粒子の質量ｍＸ´を、変換式ｍＸ´＝ｍαＸδにより求め、Ｙ方向に関する繰り込まれた粒子の質量ｍＹ´を、変換式ｍＹ´＝ｍαＹδにより求め、Ｚ方向に関する繰り込まれた粒子の質量ｍＺ´を、変換式ｍＺ´＝ｍαＺδにより求める工程と、
（ｅ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、Ｘ方向の運動方程式にはＸ方向に関する繰り込まれた粒子の質量ｍＸ´を用い、Ｙ方向の運動方程式にはＹ方向に関する繰り込まれた粒子の質量ｍＹ´を用い、Ｚ方向の運動方程式にはＺ方向に関する繰り込まれた粒子の質量ｍＺ´を用いて、分子動力学計算を実行する工程と
を有するシミュレーション方法。
【請求項９】
さらに、
（ｆ）前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた位置ベクトルをｑ´と表し、該位置ベクトルｑ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｑＸ´、ｑＹ´、及びｑＺ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた運動量ベクトルをｐ´と表し、該運動量ベクトルｐ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｐＸ´、ｐＹ´、及びｐＺ´と表し、前記粒子系Ｓ´に含まれるある粒子の、前記工程（ｅ）における分子動力学計算で求められた速度ベクトルをｖ´と表し、該速度ベクトルｖ´のＸ、Ｙ及びＺ方向それぞれの成分を、ｖＸ´、ｖＹ´、及びｖＺ´と表すとき、前記Ｘ方向の繰り込み因子αＸ、Ｙ方向の繰り込み因子αＹ、Ｚ方向の繰り込み因子αＺ、及び前記工程（ｄ）で決定された繰り込み因子δを用いて、前記粒子系Ｓに含まれるある粒子の位置ベクトルｑのＸ、Ｙ及びＺ方向それぞれの成分ｑＸ、ｑＹ、及びｑＺを、それぞれ、変換式ｑＸ＝ｑＸ´αＸ、ｑＹ＝ｑＹ´αＹ、及びｑＺ＝ｑＺ´αＺにより求める計算、前記粒子系Ｓに含まれるある粒子の運動量ベクトルｐのＸ、Ｙ及びＺ方向それぞれの成分ｐＸ、ｐＹ、及びｐＺを，それぞれ、変換式ｐＸ＝ｐＸ´／αＸδ／２、ｐＹ＝ｐＹ´／αＹδ／２、及びｐＺ＝ｐＺ´／αＺδ／２により求める計算、及び、前記粒子系Ｓに含まれるある粒子の速度ベクトルｖのＸ、Ｙ及びＺ方向それぞれの成分ｖＸ、ｖＹ、及びｖＺを，それぞれ、変換式ｖＸ＝ｖＸ´αＸδ／２、ｖＹ＝ｖＹ´αＹδ／２、及びｖＺ＝ｖＺ´αＺδ／２により求める計算のうち、少なくとも１つの計算を実行する工程
を有する請求項８に記載のシミュレーション方法。
【請求項１０】
前記Ｘ方向の繰り込み因子αＸ、Ｙ方向の繰り込み因子αＹ、及びＺ方向の繰り込み因子αＺのうち、少なくとも２つが互いに異なる請求項８または９に記載のシミュレーション方法。
【請求項１１】
請求項１〜１０のいずれかに記載のシミュレーション方法を、コンピュータで実行するためのプログラム。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、シミュレーション方法及びプログラムに関し、特に、分子動力学を用いたシミュレーション方法、及び、それをコンピュータに実行させるためのプログラムに関する。
【背景技術】
【０００２】
分子動力学を用いたコンピュータシミュレーションが行われている。分子動力学では、シミュレーションの対象となる系を構成する粒子の運動方程式が、数値的に解かれる。シミュレーション対象の系の粒子数が増大すれば、必要な計算量が増大する。現行のコンピュータの演算能力では、例えば、粒子数１０万程度の系のシミュレーションまでしか行うことができない。
【０００３】
シミュレーションに必要な計算量を減らすために、例えば非特許文献１において、分子動力学に繰り込みの手法を応用しようとする試みがなされている。
【０００４】
なお、本願に関連する発明が本願発明者らによってなされ、特許文献１に開示されている。
【０００５】
【特許文献１】特開２００６−２８５８６６号公報
【非特許文献１】稲村，武澤，社本，「繰り込みの手法に基づく可変スケール分子動力学シミュレーションについて」，日本機械学会論文集（Ａ編），１９９７年，第６３巻，６０８号，ｐ．２０２−２０７
【発明の開示】
【発明が解決しようとする課題】
【０００６】
非特許文献１の「３．数値計算法」の欄に記載されているように、非特許文献１のシミュレーション方法では、温度を陽に評価できない。
【０００７】
本発明の一目的は、分子動力学を用いた新規なシミュレーション方法、及びその方法をコンピュータで実行するためのプログラムを提供することである。
【０００８】
本発明の他の目的は、繰り込み群の手法に基づき、温度を陽に評価することが可能なシミュレーション方法、及びその方法をコンピュータで実行するためのプログラムを提供することである。
【課題を解決するための手段】
【０００９】
本発明の一観点によれば、（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓを考えたとき、該粒子系Ｓが配置されている空間の次元数ｄを用い、１より大きい第１の繰り込み因子αと、０以上ｄ以下の第２の繰り込み因子γと、０以上の第３の繰り込み因子δとを決定して、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／αｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍαδ／αγにより求め、繰り込まれた相互作用係数ε´を、変換式ε´＝εαγにより求める工程と、（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、繰り込まれた相互作用係数ε´との積ε´ｆで表される粒子系Ｓ´について、分子動力学計算を実行する工程とを有するシミュレーション方法が提供される。
【００１０】
また、本発明の他の観点によれば、（ａ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓを考えたとき、該粒子系Ｓが配置されている空間の次元数ｄを用い、１より大きい第１の繰り込み因子αと、０以上の第３の繰り込み因子δとを決定して、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／αｄにより求め、繰り込まれた粒子の質量ｍ´を、変換式ｍ´＝ｍαδにより求める工程と、（ｂ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、分子動力学計算を実行する工程とを有するシミュレーション方法が提供される。
【００１１】
更に、本発明の他の観点によれば、（ｄ）Ｎ個の粒子を含み、各粒子の質量がｍであり、粒子間の相互作用ポテンシャルエネルギを、粒子間距離に対する依存性を表す無次元化関数ｆ及び相互作用係数εの積εｆで表すことができる粒子系Ｓと、前記粒子系Ｓが配置されている空間にＸＹＺ直交座標系を考えたとき、それぞれが１より大きいＸ方向の繰り込み因子αＸ、Ｙ方向の繰り込み因子αＹ、及びＺ方向の繰り込み因子αＺを用いて、繰り込まれた粒子数Ｎ´を、変換式Ｎ´＝Ｎ／（αＸαＹαＺ）により求め、更に０以上の繰り込み因子δを決定して、Ｘ方向に関する繰り込まれた粒子の質量ｍＸ´を、変換式ｍＸ´＝ｍαＸδにより求め、Ｙ方向に関する繰り込まれた粒子の質量ｍＹ´を、変換式ｍＹ´＝ｍαＹδにより求め、Ｚ方向に関する繰り込まれた粒子の質量ｍＺ´を、変換式ｍＺ´＝ｍαＺδにより求める工程と、（ｅ）繰り込まれた粒子数Ｎ´個の粒子を含み、各粒子が、繰り込まれた粒子の質量ｍ´を有し、粒子間の相互作用ポテンシャルエネルギが、前記無次元化関数ｆと、前記相互作用係数εとの積εｆで表される粒子系Ｓ´について、Ｘ方向の運動方程式にはＸ方向に関する繰り込まれた粒子の質量ｍＸ´を用い、Ｙ方向の運動方程式にはＹ方向に関する繰り込まれた粒子の質量ｍＹ´を用い、Ｚ方向の運動方程式にはＺ方向に関する繰り込まれた粒子の質量ｍＺ´を用いて、分子動力学計算を実行する工程とを有するシミュレーション方法が提供される。
【発明の効果】
【００１２】
繰り込まれた粒子数Ｎ´は、粒子数Ｎよりも少ない。このため、粒子系Ｓに対して分子動力学計算を実行するよりも少ない計算量で、粒子系Ｓ´に対する分子動力学計算を実行することができる。
【００１３】
また、０以上の繰り込み因子δを任意に決定して計算を行うため、シミュレーションを柔軟に実行することができる。
【発明を実施するための最良の形態】
【００１４】
まず、分子動力学（Ｍｏｌｅｃｕｌａｒ  Ｄｙｎａｍｉｃｓ；ＭＤ）について簡潔に説明する。Ｎ個の粒子（例えば原子）からなり、ハミルトニアンＨが、
2010146368.tif 000003
【００１５】
・・・（１）
と表される粒子系について考える。ここで、各粒子の質量がｍであり、粒子ｊの運動量ベクトルがｐｊであり、粒子ｊの位置ベクトル（位置座標）がｑｊであり、粒子間の相互作用ポテンシャルエネルギがφである。
【００１６】
ハミルトニアンＨを、ハミルトンの正準方程式に代入することにより、各粒子に対する運動方程式
2010146368.tif 000004
【００１７】
・・・（２）
及び
2010146368.tif 000005
【００１８】
・・・（３）
が得られる。ここで式（３）において、ｖｊは粒子ｊの速度ベクトルである。分子動力学では、粒子系を構成する各粒子に対し、式（２）及び（３）で表される運動方程式を数値的に積分して解くことにより、各時刻における各粒子の運動量ベクトル（または速度ベクトル）及び位置ベクトルが求められる。なお、数値積分には、Ｖｅｒｌｅｔ法が用いられることが多い。Ｖｅｒｌｅｔ法については、例えば、Ｊ．Ｍ．Ｔｈｉｊｓｓｅｎ，“Ｃｏｍｐｕｔａｔｉｏｎａｌ  Ｐｈｙｓｉｃｓ”，Ｃａｍｂｒｉｄｇｅ Ｕｎｉｖｅｒｓｉｔｙ Ｐｒｅｓｓ（１９９９）のｐ．１７５に説明されている。分子動力学計算で得られた各粒子の位置、速度に基づき、系の種々の物理量を算出することができる。
【００１９】
次に、繰り込み群の手法を応用した分子動力学（これを、繰り込み群分子動力学と呼ぶこととする）について、概念的に説明する。
【００２０】
繰り込み群分子動力学では、物理量を求めたい所望の系Ｓを、系Ｓよりも少ない粒子数で構成される系（これを、繰り込まれた系Ｓ´と呼ぶこととする）に対応付ける。次に、繰り込まれた系Ｓ´に対する分子動力学計算を実行する。そして、繰り込まれた系Ｓ´に対する計算結果を、所望の系Ｓに対応付ける。これにより、所望の系Ｓに対する分子動力学計算を直接的に実行するよりも少ない計算量で、所望の系Ｓに対する物理量を算出することが可能となる。所望の系Ｓにおける物理量（例えば、粒子数、各粒子の質量等）と、繰り込まれた系Ｓ´における物理量とを互いに対応付けるための変換則を、スケール変換則と呼ぶ。
【００２１】
次に、本願発明者が見出したスケール変換則について説明する。物理量を算出したい所望の粒子系Ｓが、Ｎ個の粒子から構成されるとする。粒子系Ｓの全ハミルトニアンＨが、式（１）に示したように、
2010146368.tif 000006
【００２２】
と表されるとする。ある粒子ｊのハミルトニアンＨｊは、
2010146368.tif 000007
【００２３】
・・・（４）
と表される。
【００２４】
繰り込みの議論のために、相互作用ポテンシャルエネルギφを、粒子間距離に対する依存性を表し、無次元化された関数ｆと、相互作用の強さを表し、エネルギの次元を持つ係数ε（これを相互作用係数εと呼ぶこととする）との積により、以下の様に表す。
2010146368.tif 000008
【００２５】
・・・（５）
例えば、不活性原子に対しては、
2010146368.tif 000009
【００２６】
・・・（６）
という相互作用ポテンシャルが用いられる。ここでｒは原子間距離である。相互作用係数ε及び定数σが、原子の種類に対応する値を取る。
【００２７】
次に、ハミルトニアンＨから、繰り込まれたハミルトニアンＨ´を導出する方法について説明する。以下詳細に示すように、繰り込まれたハミルトニアンＨ´は、粒子系Ｓに対する分配関数Ｚ（β）における積分の一部を実行し、ハミルトニアンを粗視化し、続いて積分変数を再スケールすることで得られる。
【００２８】
粒子数一定における正準集団（Ｃａｎｏｎｉｃａｌ  ｅｎｓｅｍｂｌｅ）に対する分配関数Ｚ(β)は、粒子系Ｓに対して、
2010146368.tif 000010
【００２９】
・・・（７）
と表される。ここで、係数βは、系の温度Ｔとボルツマン定数ｋＢとによりβ＝１／（ｋＢＴ）と定義される。ｄΓＮは位相空間内の体積要素であり、より詳細には、
2010146368.tif 000011
【００３０】
・・・（８）
と表される。ここで、ＷＮは、Ｎ！ｈ３Ｎ（ｈはプランク定数）である。また、ＤｐＮ及びＤｑＮをそれぞれ、
2010146368.tif 000012
2010146368.tif 000013
【００３１】
と定義している。
【００３２】
調和振動子以外のポテンシャルの繰り込みを議論することは困難である。そこで、鞍点を持つ相互作用ポテンシャルを３つの領域に分けて考える。粒子間距離ｒ→∞とｒ→０は、固定点であるから繰り込み変換に際しポテンシャルは不変である。そこで、鞍点近傍の繰り込み変換がすべての粒子間距離ｒについても成立すると仮定する。鞍点の周りをテーラー展開し、２次の揺らぎまで残す。鞍点の位置をｒ０とし、相対変位をδｕとすると、δｕの１次の項は消えて、
2010146368.tif 000014
【００３３】
・・・（９）
を得る。従って、粒子ｊについて、鞍点近傍を記述するハミルトニアンは、
2010146368.tif 000015
【００３４】
・・・（１０）
となる。ここで、ｕｉ及びｕｊは、それぞれ、粒子ｉ及びｊの鞍点からの変位を表す。φ″は、粒子間距離ｒについての２階微分である。
【００３５】
次に、分配関数Ｚ（β）の粒子間相互作用部分の粗視化について説明する。まず、一次元鎖上に配置された粒子系における粗視化について考察した後、単純立方格子上に配置された粒子系へ拡張する。
【００３６】
図１（Ａ）に示すように、１次元鎖上の格子点（ラティスポイント）に、粒子ｉとｊとが互いに最隣接して配置され、粒子ｊとｋとが互いに最隣接して配置されている。粒子ｉとｋとの中間に、粒子ｊが配置されている。最隣接する粒子同士の相互作用（二アレストネイバーカップリング）を実線で示す。粒子ｊを消去し、粗視化することを考える。粒子ｊが関与する相互作用を書き出し、粒子ｉとｋとの中間にある粒子ｊの変位ｕｊについて積分を実行すれば、粒子ｊの影響が繰り込まれた後の、粒子ｉとｋとの相互作用ポテンシャルφ´（ｑｉ−ｑｋ）について、
2010146368.tif 000016
【００３７】
・・・（１１）
という式を得る。
【００３８】
この粗視化により、格子定数がα（＝２）倍になるから、ｕ´＝ｕ／αという式でスケール変換された変数ｕ´により、元のポテンシャル関数と同形に表すことができる。つまり、以下の相似則が得られる。
2010146368.tif 000017
【００３９】
・・・（１２）
ただし、ここで、式（１２）の２行目へ移るに際して、鞍点近傍以外においても同じ繰り込み変換が成立すると仮定した。αを、第１の繰り込み因子と呼ぶこととする。
【００４０】
ｄ次元単純立方格子における粗視化は、ポテンシャルムービング（Ｐｏｔｅｎｔｉａｌ   Ｍｏｖｉｎｇ）の方法により実現できる。Ｐｏｔｅｎｔｉａｌ  Ｍｏｖｉｎｇについては、例えば、Ｌｅｏ Ｐ．Ｋａｄａｎｏｆｆ，“ＳＴＡＴＩＳＴＩＣＡＬ ＰＨＹＳＩＣＳ，Ｓｔａｔｉｃｓ，Ｄｙｎａｍｉｃｓ ａｎｄ Ｒｅｎｏｒｍａｌｉｚａｔｉｏｎ”，Ｗｏｒｌｄ  Ｓｃｉｅｎｔｉｆｆｉｃ（１９９９） の１４章に解説されている。
【００４１】
図１（Ｂ）を参照して、Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇの考え方について説明する。Ｐｏｔｅｎｔｉａｌ  Ｍｏｖｉｎｇでは、多次元格子を１次元鎖へ帰着させる。図１（Ｂ）には、２次元格子の場合を示す。
【００４２】
図１（Ｂ）の上段の図に示すように、２次元正方格子の格子点上に粒子が配置されている。最隣接する粒子同士の相互作用（最隣接間相互作用）を実線で示す。格子の有る辺に平行な方向をＸ方向とし、Ｘ方向に直交する方向をＹ方向とする。
【００４３】
図１（Ｂ）の中段の図に示すように、Ｘ方向に並んだ粒子について、１つ置きに粒子を積分することを考える（白抜きの丸で示した格子点（ラティスポイントフォーサメイション）上の粒子を積分することを考える）。積分したい（消去したい）粒子を、被積分粒子と呼ぶこととする。被積分粒子同士の最隣接間相互作用を波線で示す。もし被積分粒子同士の最隣接間相互作用が存在しなければ、被積分粒子については、Ｘ方向に延在する１次元鎖と見なせる。
【００４４】
Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇでは、図１（Ｂ）の下段の図に示すように、被積分粒子同士の最隣接間相互作用を、Ｘ方向について被積分粒子の両隣に配置された粒子に振り分ける。振り分けられた相互作用が加算された粒子間相互作用（ダブルカップリング）を、二重線で示す。二重線で示す相互作用は、元の最隣接間相互作用の２倍の強さになる。この様にして、被積分粒子に対して、２次元格子から１次元鎖に変換できる。そこで、被積分粒子に対して、上述した１次元鎖に対する粗視化の方法が実行できる。例えば３次元格子の場合、以上説明したのと同様な手順を、Ｙ、Ｚ方向についても繰り返す。このようにして、多次元格子における粗視化が実行される。
【００４５】
Ｐｏｔｅｎｔｉａｌ Ｍｏｖｉｎｇの方法を用いることにより、ｄ次元単純立方格子における分配関数Ｚ（β）の相互作用部分の粗視化について、
2010146368.tif 000018
【００４６】
・・・（１３）
という結果が導かれる。ここで、繰り込まれた粒子数Ｎ´及び第１の繰り込み因子αが、
2010146368.tif 000019
2010146368.tif 000020
【００４７】
という式で表される。ｎは正整数であり、粗視化の回数を表す。
【００４８】
次に、分配関数Ｚ（β）の運動エネルギ部分の粗視化について説明する。まずは、１次元鎖において考える。互いに最隣接する粒子対である粒子ｊ−１と粒子ｊとを、１つの粒子へ粗視化する。Ｋａｄａｎｏｆｆのブロックスピン法の考え方を運動量に適用し、粒子ｊと粒子ｊ−１とからなる粗視化された粒子の運動量ｇｊを以下の様に定義する。
2010146368.tif 000021
【００４９】
また、重み関数Ｔ｛Ｐ，ｐ｝を以下の様に定義する。
2010146368.tif 000022
【００５０】
ここで、δはＤｉｒａｃのδ-関数である。明らかに、
2010146368.tif 000023
【００５１】
を満たす。この条件から、分配関数の運動エネルギ部分Ｚｐへ重み関数Ｔ｛Ｐ，ｐ｝を挿入でき、
2010146368.tif 000024
【００５２】
・・・（１４）
という式を得る。δ-関数を積分表示すれば、分配関数の運動エネルギ部分Ｚｐは、以下の様に表せる。
2010146368.tif 000025
【００５３】
・・・（１５）
この積分は、以下の公式を使えば容易に実行できる。
2010146368.tif 000026
【００５４】
ここで、＜Ａ，Ｂ＞は、ベクトルＡとＢとの内積を表す。この公式を用いると、分配関数の運動エネルギ部分Ｚｐは、
2010146368.tif 000027
【００５５】
・・・（１６）
と表すことができる。
【００５６】
この粗視化により、運動エネルギがα（＝２）倍になる。なお、２つの粒子の全運動エネルギを重心の運動エネルギと相対運動の運動エネルギへ分け、相対運動量について積分を実行した場合も同様な結果が得られる。
【００５７】
次に、ｄ次元単純立方格子における粗視化について考察する。まず、結晶軸ａ方向について粗視化すると、粒子の持つ運動エネルギＥｋｉｎはα倍になる。引き続いてｂ軸方向において粗視化すると、粒子の持つ運動エネルギＥｋｉｎαは、（Ｅｋｉｎα）のα倍、即ちＥｋｉｎα２になる。従って、ｄ次元格子ではＥｋｉｎαｄになるから、
2010146368.tif 000028
【００５８】
・・・（１７）
という関係式が得られる。ｎ回の粗視化では、α＝２ｎとなることは自明である。
【００５９】
以上の議論から、粗視化された全ハミルトニアンに対する分配関数Ｚ（β）は、結果に影響しない係数を除いて、
2010146368.tif 000029
【００６０】
・・・（１８）
と表すことができる。
【００６１】
次に、粗視化された系においてもハミルトニアンを同形とするために、“繰り込まれた”変数を導入する。以下のように、繰り込まれた位置ベクトルｑ´、繰り込まれた運動量ベクトルｐ´、繰り込まれた質量ｍ´、繰り込まれた相互作用係数ε´、及び、繰り込まれた係数β´を定義する。ハミルトニアンの形が不変となるようにスケール変換する。
2010146368.tif 000030
2010146368.tif 000031
2010146368.tif 000032
2010146368.tif 000033
2010146368.tif 000034
【００６２】
ここで、γは、０以上ｄ以下である。γを、第２の繰り込み因子と呼ぶこととする。また、δは０以上の値をとる。δを第３の繰り込み因子と呼ぶ。
【００６３】
繰り込まれた変数をこのように定めることにより、ハミルトニアンＨと同形な、繰り込まれたハミルトニアンＨ´が得られる。繰り込まれたハミルトニアンＨ´は、
2010146368.tif 000035
【００６４】
・・・（１９）
と表される。
【００６５】
なお、式（１８）から、分配関数は以下の様になる。
2010146368.tif 000036
【００６６】
・・・（２０）
ここで、位相空間の体積素は、結果に影響しない定数を除いて、
2010146368.tif 000037
【００６７】
・・・（２１）
である。
【００６８】
繰り込まれたハミルトニアンＨ´で記述される系（繰り込まれた粒子系Ｓ´）における運動方程式は、正準方程式
2010146368.tif 000038
【００６９】
・・・（２２）
及び
2010146368.tif 000039
【００７０】
・・・（２３）
から与えられ、
2010146368.tif 000040
【００７１】
・・・（２４）
及び
2010146368.tif 000041
【００７２】
・・・（２５）
となる。ここで、繰り込まれた時間ｔ´が、繰り込まれた変数ｑ´、ｐ´、ｍ´及びε´に整合するように、
2010146368.tif 000042
【００７３】
という式で定められる。運動方程式（２４）及び（２５）に基づいて、繰り込まれた粒子系Ｓ´に対する分子動力学計算を実行することができる。
【００７４】
以上の考察をまとめる。粒子系Ｓと、繰り込まれた粒子系Ｓ´とは、以下の変換式で表されるスケール変換則により対応付けられる。
【００７５】
粒子系Ｓの粒子数Ｎと、繰り込まれた粒子系Ｓ´の繰り込まれた粒子数Ｎ´とは、変換式
2010146368.tif 000043
【００７６】
・・・（変換式Ｎ）
により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑと、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´とは、変換式
2010146368.tif 000044
【００７７】
・・・（変換式ｑ）
により対応付けられる。粒子系Ｓのある粒子の運動量ベクトルｐと、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´とは、変換式
2010146368.tif 000045
【００７８】
・・・（変換式ｐ）
により対応付けられる。粒子系Ｓの各粒子の質量ｍと、繰り込まれた粒子系Ｓ´の各粒子の質量ｍ´とは、変換式
2010146368.tif 000046
【００７９】
・・・（変換式ｍ）
により対応付けられる。粒子系Ｓにおける相互作用係数εと、繰り込まれた粒子系Ｓ´における相互作用係数ε´とは、変換式
2010146368.tif 000047
【００８０】
・・・（変換式ε）
により対応付けられる。
【００８１】
粒子系Ｓに対して運動方程式を解く場合の時間スケールｔ（粒子系Ｓに対して実行すると仮定した分子動力学計算におけるある時間間隔ｔ）と、繰り込まれた粒子系Ｓ´に対して運動方程式を解く場合の時間スケールｔ´（繰り込まれた粒子系Ｓ´に対して実行する分子動力学計算におけるある時間間隔ｔ´）とは、変換式
2010146368.tif 000048
【００８２】
・・・（変換式ｔ）
により対応付けられる。
【００８３】
なお、粒子系Ｓのある粒子の速度ベクトルｖと、繰り込まれた粒子系Ｓ´のある粒子の速度ベクトルｖ´とは、速度ベクトルｖがｖ＝ｐ／ｍと表されることから、変換式
2010146368.tif 000049
【００８４】
・・・（変換式ｖ）
により対応付けられる。
【００８５】
ここで、粒子系Ｓの配置されている空間の次元数がｄである。第１の繰り込み因子αは２ｎ（ｎは正整数）であり、第２の繰り込み因子γは０以上ｄ以下である。また、第３の繰り込み因子δは０以上である。
【００８６】
第２の繰り込み因子γの値が大きく、例えばｄに等しいとき、極めて大きな、繰り込まれた相互作用係数ε´と、小さな、繰り込まれた質量ｍ´とが数値計算に現れる。数値計算上、第２の繰り込み因子γは０に設定することが好ましい。このとき更に、第３の繰り込み因子δを０とすれば、時刻ｔのみが変換され、従来のＭＤが利用できる。
【００８７】
なお、たとえばδ＝２のとき、質量ｍとバネ定数ｓを持つ１次元鎖のバネ・マス系を考えたとき、繰り込まれたハミルトニアンから分散関係が以下の様に得られる。
2010146368.tif 000050
【００８８】
・・・（２６）
ここで、角周波数がωであり、波数がｋである。平衡状態における原子間距離がａである。なお、第２の繰り込み因子γは０とした。式（２６）は、長波長の極限において、
2010146368.tif 000051
【００８９】
・・・（２７）
となり、第１の繰り込み因子αには無関係になる。
【００９０】
次に、繰り込み群分子動力学によるコンピュータシミュレーションの実施例について説明する。
【００９１】
まず、アルミニウムに熱を流入させて昇温させ、融点及び潜熱を算出したシミュレーションについて説明する。第１のシミュレーションでは、比較のため、１６０００個のアルミニウム原子からなる系に対して、従来の（繰り込みを応用しない）分子動力学計算を行った。アルミニウム原子は３次元空間内に配置されており、次元数ｄは３である。
【００９２】
第２及び第３のシミュレーションでは、第１のシミュレーションに対応する系を繰り込んだ系に対する分子動力学計算を行った。第１の繰り込み因子αは、第２のシミュレーションでは２と決定し、第３のシミュレーションでは２２（＝４）と決定した。（変換式Ｎ）より、第２のシミュレーションでは、系の粒子数が２０００個となり、第３のシミュレーションでは、系の粒子数が２５０個となる。第２の繰り込み因子γは０と決定した。また、第３の繰り込み因子δは２と決定した。なお、従来の手法（第１のシミュレーション）は、第１の繰り込み因子αを１と設定した場合に対応する。
【００９３】
第１のシミュレーションにおける原子間ポテンシャル（つまり、元の系における粒子間ポテンシャル）として、モース（Ｍｏｒｓｅ）ポテンシャルを用いた。このポテンシャルは、
2010146368.tif 000052
【００９４】
・・・（２８）
と表される。ここで、相互作用係数εは１．９２×１０−２０[Ｊ]であり、定数Ａは２．３５×１０１０[ｍ−１]であり、定数ｒ０は２．８６×１０−１０[ｍ]である。粒子間距離がｒである。
【００９５】
第２及び第３のシミュレーションにおける粒子間ポテンシャル（つまり、繰り込まれた系における粒子間ポテンシャル）としても、式（２８）と同形のものが用いられる。第２及び第３のシミュレーションにおける繰り込まれた相互作用係数は、それぞれ、（変換式ε）より求められる。第２の繰り込み因子γが０であるので、第２及び第３のシミュレーションでの相互作用係数は、双方とも、第１のシミュレーションの相互作用係数εと等しい。定数ｒ０及び定数Ａは第１のシミュレーションのそれと等しい。
【００９６】
第１のシミュレーションにおけるアルミニウム原子の質量（つまり、元の系における粒子の質量）は、４．４８×１０−２６[ｋｇ]である。第２及び第３のシミュレーションにおける粒子の繰り込まれた質量は、それぞれ、（変換式ｍ）より求められる。第２のシミュレーションにおける質量は、第１のシミュレーションでの質量に２２（＝４）を乗じた１７．９×１０−２６[ｋｇ]となり、第３のシミュレーションにおける質量は、第１のシミュレーションでの質量に４２（＝１６）を乗じた７１．７×１０−２６[ｋｇ]となる。
【００９７】
次に、粒子の初期配置について説明する。初期において、粒子は面心立方格子の格子点上に配置される。第１のシミュレーションでは、幅方向に２０個、奥行き方向に２０個の粒子が並んだ層が、４０層重なる配置とした。元の系の、幅方向について１／α、奥行き方向について１／α、及び高さ方向について１／αの部分（元の系の１／α３の部分）が、繰り込まれた系に対応する。第２のシミュレーションでは、幅方向に１０個、奥行き方向に１０個の粒子が並んだ層が、２０層重なった配置となり、第３のシミュレーションでは、幅方向に５個、奥行き方向に５個の粒子が並んだ層が、１０層重なった配置となる。
【００９８】
次に、数値積分の時間刻み幅について説明する。元の系における時間刻み幅をΔｔとし、粒子の速さをｖとする。時間刻み幅Δｔは、例えば  
2010146368.tif 000053
【００９９】
を満たすように定められる（「＜＜」は、式の左辺が右辺より充分小さいという意）。ここでｒ０は粒子間ポテンシャルにおける定数であり、粒子間距離の目安を与える。繰り込まれた系における時間刻み幅をΔｔ´とし、粒子の速さをｖ´とする。（変換式ｖ）よりｖ´＝ｖ／αであるので、繰り込まれた系における時間刻み幅Δｔ´は、例えば、αΔｔに設定される。
【０１００】
数値積分の時間刻み幅は、第１のシミュレーションにおいては５．０[ｆｓ]とし、第２のシミュレーションにおいては１０[ｆｓ]とし、第３のシミュレーションにおいては２０[ｆｓ]とした。
【０１０１】
次に、粒子に与える初期速度について説明する。まず、元の系の温度と、繰り込まれた系の温度との関係について説明する。元の系の温度Ｔは、
2010146368.tif 000054
【０１０２】
・・・（２９）
という式で定義される。ここで、ｋＢはボルツマン定数であり、ｖｉは粒子ｉの速度ベクトルであり、ｖｇは粒子系の重心の速度ベクトルである。温度は、つまり、速度の“揺らぎ”である。繰り込まれた系の温度Ｔ´は、
2010146368.tif 000055
【０１０３】
・・・（３０）
という式から定まる。
【０１０４】
元の系において、質量がｍで速さがｖの粒子の運動エネルギは、１／２ｍｖ２である。繰り込まれた系において、質量がｍ´で速さがｖ´の粒子の運動エネルギは、１／２ｍ´ｖ´２である。（変換式ｍ）及び（変換式ｖ）から、
2010146368.tif 000056
【０１０５】
・・・（３１）
という関係が得られる。よって、元の系の温度Ｔと、繰り込まれた系の温度Ｔ´とは等しくなる。第１〜第３のシミュレーションにおいて、粒子系の温度が３００[Ｋ]となるように、初期速度が与えられる。
【０１０６】
次に、粒子系に熱を流入させる方法について説明する。第１〜第３のシミュレーションでは、粒子系の底面上に配置された粒子に運動エネルギを与えて、系への熱の流入を模擬した。熱の流入する面を加熱面と呼ぶこととする。
【０１０７】
まず、第１のシミュレーションにおける熱の流入方法（つまり、元の系における熱の流入方法）ついて説明する。加熱面上の粒子１個当りに、時間刻み幅Δｔの間に、運動エネルギΔＥが与えられるとする。粒子の運動量ベクトルの変化は、下式（３２）より求められる。ここで、ある時刻ｔにおける運動量ベクトルがｐｔであり、時刻ｔ＋Δｔにおける運動量ベクトルがｐｔ＋Δｔである。
2010146368.tif 000057
【０１０８】
・・・（３２）
運動エネルギΔＥによって、運動量ベクトルがλ倍されるとする。つまり、
2010146368.tif 000058
【０１０９】
・・・（３３）
と表されるとする。式（３３）を式（３２）に代入することにより、λは、
2010146368.tif 000059
【０１１０】
・・・（３４）
と表される。時間刻み幅Δｔごとに、式（３３）により運動量ベクトルを更新することにより、時間刻み幅Δｔ当たりに運動エネルギΔＥが流入するような熱の流入を模擬できる。
【０１１１】
加熱面上の粒子数をＮｈとし、加熱面の面積をＳｈとすると、単位面積当たりに流入するパワーｆｈは、
2010146368.tif 000060
【０１１２】
・・・（３５）
と表される。第１のシミュレーションでは、加熱面上の粒子数Ｎｈが４００個であり、加熱面の面積Ｓｈが２．８３×１０−１７[ｍ２]であり、運動エネルギΔＥが４．１４×１０−２４[Ｊ]であり、時間刻み幅Δｔが５．０[ｆｓ]であり、パワーｆｈが１．１７×１０１０[Ｗ／ｍ２]である。
【０１１３】
第２及び第３のシミュレーションにおいて（繰り込まれた系において）、単位面積当たりに流入するパワーｆｈ´は、
2010146368.tif 000061
【０１１４】
・・・（３６）
という式で与えられる。運動量ベクトルを更新するための係数λ´は、
2010146368.tif 000062
【０１１５】
・・・（３７）
という式で与えられる。なお式（３１）から、与えられる運動エネルギΔＥは元の系と繰り込まれた系とで等しくなる。
【０１１６】
次に、図２を参照して、シミュレーション結果について説明する。図２に示すグラフのうち、下段のグラフが第１のシミュレーション（従来の分子動力学計算）の結果を表し、中段のグラフが第２のシミュレーション（繰り込み因子αを２とした繰り込み群分子動力学計算）の結果を表し、上段のグラフが第３のシミュレーション（繰り込み因子αを４とした繰り込み群分子動力学計算）の結果を表す。３つのグラフにおいて、縦軸が温度（単位[Ｋ]）を示し、横軸が流入した熱量、つまり流入したエネルギ（単位１０−１６[Ｊ]）を示す。
【０１１７】
まず、下段のグラフを参照して、第１のシミュレーションの結果について説明する。系にエネルギが流入していないとき、系の温度は３００[Ｋ]である。流入したエネルギが増加するにつれ、系の温度が上昇する。流入エネルギが約５×１０−１６[Ｊ]に達すると、温度の上昇が、約８５０[Ｋ]で一旦止まり、流入エネルギが約９×１０−１６[Ｊ]に達するまで、温度がほぼ一定となる。この温度が、融点を示す。流入エネルギが約９×１０−１６[Ｊ]を超えると、再び温度が上昇する。
【０１１８】
温度一定の期間中に系に流入したエネルギが潜熱に対応する。このシミュレーションから求められる潜熱は、１４．７[ｋＪ／ｍｏｌ]となる。なお、実験から求められたアルミニウムの融点は９３３[Ｋ]であり、潜熱は１１[ｋＪ／ｍｏｌ]である。シミュレーションで得られた値と実験値との相違は、用いた原子間ポテンシャルの精度等に起因すると考えられる。
【０１１９】
次に、中段及び上段のグラフを参照して、第２及び第３のシミュレーションの結果について説明する。第２及び第３のシミュレーションの結果も、第１のシミュレーションの結果と同様な傾向を示す。すなわち、流入エネルギが約５×１０−１６[Ｊ]から約９×１０−１６[Ｊ]までの期間中に、温度が８５０[Ｋ]付近でほぼ一定となる傾向を示す。これにより、第２及び第３のシミュレーションからも、融点が約８５０[Ｋ]であり、潜熱が１４．７[ｋＪ／ｍｏｌ]であると求めることができる。
【０１２０】
粒子系に流入したエネルギは、シミュレーションした期間分だけ流入パワーを積分すれば求められる。繰り込まれた系に流入したエネルギをＥ´とすると、元の系へ流入したエネルギＥは、
2010146368.tif 000063
【０１２１】
・・・（３８）
という式から求められる。なお、粒子１つ当たりの流入エネルギは、元の系と繰り込まれた系とで等しい。
【０１２２】
繰り込まれた系の温度は、式（３０）から求められる。上述の議論より、元の系の温度は、繰り込まれた系の温度と等しい。
【０１２３】
なお、上述のシミュレーションでは、元の系における流入エネルギ及び温度を求めたが、元の系における他の物理量、例えば粒子の位置、運動量、速度等も、必要に応じて、上述のスケール変換則に基づいて求めることができる。
【０１２４】
以上説明したように、繰り込み群分子動力学では、繰り込む前の系（上記実施例の第１のシミュレーションに対応する系）に対する分子動力学計算で得られると期待される物理量（上記実施例では融点及び潜熱）を、繰り込む前の系よりも粒子数の少ない系に対する分子動力学計算を行うことにより、算出することが可能となる。繰り込み群分子動力学によるシミュレーションに要する種々の計算は、プログラムにより、コンピュータで実行することができる。
【０１２５】
なお、上記実施例のシミュレーションにおいて、アルミニウムは面心立方構造を取る。スケール変換則は、単純立方格子上に配置された粒子系に対する考察から導出したが、単純立方構造に限らない他の結晶構造に対しても有効である。
【０１２６】
次に、シリコンに熱を流入させて昇温させ、融点及び潜熱を算出したシミュレーションについて説明する。シリコンは、ダイヤモンド構造を取る。アルミニウムに対するシミュレーションと同様に、３つのシミュレーションを行った。第１のシミュレーションは、比較のために行った従来の分子動力学手法によるものである。第２及び第３のシミュレーションにおいて、第１のシミュレーションに対応する系を繰り込んだ系に対する分子動力学計算を行った。繰り込み因子αは、第２のシミュレーションでは２とし、第３のシミュレーションでは２２（＝４）とした。系の粒子数は、第１のシミュレーションでは８１９２個であり、第２のシミュレーションでは１０２４個であり、第３のシミュレーションでは１２８個である。
【０１２７】
第１のシミュレーションにおける原子間ポテンシャル（つまり、元の系における粒子間ポテンシャル）として、Ｓ−Ｗ（Ｓｔｉｌｌｉｎｇｅｒ−Ｗｅｂｅｒ）ポテンシャルを用いた。このポテンシャルは、
2010146368.tif 000064
【０１２８】
・・・（３９）
と表される。ここで、
2010146368.tif 000065
【０１２９】
であり、相互作用係数εは、２．１６７[ｅＶ]であり、σは２．０９５[Å]であり、定数Ａ、Ｂ、ｐ、ｑ、ａ、λ、及びγは、それぞれ、７．０５、０．６０２、４．０、０．０、１．８、２１．０、及び１．２である。
【０１３０】
また、図３（Ａ）に示すように、粒子ｉ及び粒子ｋの結合と、粒子ｉ及び粒子ｊの結合とのなす角がθｉｊｋであり、粒子ｉから粒子ｊまでの粒子間距離がｒｉｊである。第２及び第３のシミュレーションの粒子間ポテンシャルにおける繰り込まれた相互作用係数は、上述のアルミニウムに対するシミュレーションと同様な方法で求められる。
【０１３１】
第１のシミュレーションにおけるシリコン原子の質量（つまり、元の系における粒子の質量）は、４．６７×１０−２６[ｋｇ]とした。第２及び第３のシミュレーションにおける粒子の繰り込まれた質量は、上述のアルミニウムに対するシミュレーションと同様な方法で求められる。
【０１３２】
次に、図３（Ｂ）を参照して、シミュレーション結果について説明する。図３（Ｂ）に示すグラフのうち、下段のグラフが第１のシミュレーション（従来の分子動力学計算）の結果を表し、中段のグラフが第２のシミュレーション（繰り込み因子αを２とした繰り込み群分子動力学計算）の結果を表し、上段のグラフが第３のシミュレーション（繰り込み因子αを４とした繰り込み群分子動力学計算）の結果を表す。
【０１３３】
まず、下段のグラフを参照して、第１のシミュレーションの結果について説明する。流入エネルギが約１．３×１０−１６[Ｊ]から約４．６×１０−１６[Ｊ]に達するまでの期間中、温度が１７５０[Ｋ]程度でほぼ一定となる。潜熱は、２５．６[ｋＪ／ｍｏｌ]と求められる。なお、実験から求められたシリコンの融点は１６８７[Ｋ]であり、潜熱は５０．６[ｋＪ／ｍｏｌ]である。シミュレーションで得られた値と実験値との相違は、原子間ポテンシャルの精度等に起因すると考えられる。
【０１３４】
次に、中段及び上段のグラフを参照して、第２及び第３のシミュレーションの結果について説明する。第２及び第３のシミュレーションの結果も、第１のシミュレーションの結果とほぼ同様な傾向を示す。これにより、第２及び第３のシミュレーションに基づいて、第１のシミュレーションで得られた値に近い融点及び潜熱を得ることができる。
【０１３５】
次に、図４を参照して、分散関係を繰り込み群分子動力学で求めたシミュレーションについて説明する。アルミニウムの[１１１]方向に平行に振動を加えて、[１１１]方向の波数を算出した。原子間はバネで接続した。バネ定数ｓは、４０．３[Ｎ／ｍ]とした。なお、バネ定数ｓは、アルミニウムの原子間ポテンシャルの最小近傍を２次関数で近似することにより求められる。アルミニウム原子の質量ｍは、４．４８×１０−２６[ｋｇ]とした。原子間距離ａは、２．８６×１０−１０[ｍ]とした。分散関係の厳密解（解析的に求められた解）は、加振の角周波数をωとし、[１１１]方向の波数をｋとして、
2010146368.tif 000066
【０１３６】
という式で与えられる。なお、この厳密解は、式（２６）においてα＝１として求めることができる。
グラフ中の実線Ｌが、厳密解を示す。点Ｐ１が、従来の分子動力学計算（繰り込み因子α＝１の場合に対応）の結果を示す。点Ｐ２〜点Ｐ５が、それぞれ、繰り込み因子αを、２７、２１５、２２１、及び２２８とした場合の結果を示す。繰り込み因子αを大きくすることは、元の系（繰り込まれる前の系）のサイズを大きくすることに対応する。繰り込み因子αを大きくすることにより、小さい波数を得ることができる。繰り込み群分子動力学で得られた結果は、厳密解と良く一致している。このように、繰り込み群分子動力学を用いれば、従来の分子動力学で実行するのが困難な長波長域の計算（連続体とみなせるようなマクロな系に対する計算）を、精度良く行うことができる。
【０１３７】
粒子系が配置されている空間に、ＸＹＺ直交座標系を考える。上述の実施例では、元の系の位置ベクトルｑから、繰り込み因子αを用いて、繰り込まれた位置ベクトルｑ´が、
2010146368.tif 000067
【０１３８】
・・・（変換式ｑ）
と求められた。これは、Ｘ、Ｙ及びＺのすべての方向に関する空間的なスケール変換が、共通の繰り込み因子αにより行われることを示す。繰り込み因子αを、Ｘ、Ｙ及びＺの各方向について独立に設定することも可能である。つまり、異方性を有するスケール変換を行うことも可能である。
【０１３９】
第２の繰り込み因子γが０、及び第３の繰り込み因子δが２と決定された場合について、異方性を有するスケール変換則は、以下に示すような変換式で表される。Ｘ方向の繰り込み因子をαＸとし、Ｙ方向の繰り込み因子をαＹとし、Ｚ方向の繰り込み因子をαＺとする。各繰り込み因子αＸ〜αＺは、それぞれ、例えば２ｎ（ｎは正整数）である。
【０１４０】
粒子系Ｓの粒子数Ｎと、繰り込まれた粒子系Ｓ´の繰り込まれた粒子数Ｎ´とは、変換式
2010146368.tif 000068
【０１４１】
・・・（変換式Ｎａ）
により対応付けられる。粒子系Ｓのある粒子の位置ベクトルｑの各成分と、繰り込まれた粒子系Ｓ´のある粒子の位置ベクトルｑ´の各成分とは、変換式
2010146368.tif 000069
【０１４２】
・・・（変換式ｑａ）
により対応付けられる。ここで、ηは、ベクトルのＸ，Ｙ、Ｚのいずれかの成分を示す（以下、（変換式ｐａ）、（変換式ｖａ）、及び（変換式ｍａ）においても同様である。）。粒子系Ｓのある粒子の運動量ベクトルｐの各成分と、繰り込まれた粒子系Ｓ´のある粒子の運動量ベクトルｐ´の各成分とは、変換式
2010146368.tif 000070
【０１４３】
・・・（変換式ｐａ）
により対応付けられる。
【０１４４】
なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。
2010146368.tif 000071
【０１４５】
粒子系Ｓのある粒子の速度ベクトルｖの各成分と、繰り込まれた粒子系Ｓ´のある粒子の速度ベクトルｖ´の各成分とは、変換式
2010146368.tif 000072
【０１４６】
・・・（変換式ｖａ）
により対応付けられる。
【０１４７】
なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。
2010146368.tif 000073
【０１４８】
方向ごとに繰り込み因子が設定されているので、粒子の質量も方向ごとにスケール変換される。粒子系Ｓの各粒子の質量ｍと、繰り込まれた粒子系Ｓ´の各粒子の、Ｘ方向の質量ｍ´Ｘ、Ｙ方向の質量ｍ´Ｙ、及びＺ方向の質量ｍ´Ｚとは、変換式
2010146368.tif 000074
【０１４９】
・・・（変換式ｍａ）
により対応付けられる。繰り込まれた系における運動方程式（２５）で、方向ごとに、対応する質量が用いられる。
【０１５０】
なお、γ＝０、δが未決定の場合は、δを決定し、以下の式に代入して対応付けを行う。
2010146368.tif 000075
【０１５１】
粒子系Ｓにおける相互作用係数εと、繰り込まれた粒子系Ｓ´における相互作用係数ε´とは、変換式
2010146368.tif 000076
【０１５２】
・・・（変換式εａ）
により対応付けられ、粒子系Ｓにおける時間スケールｔと、繰り込まれた粒子系Ｓ´におけるスケールｔ´とは、変換式
2010146368.tif 000077
【０１５３】
・・・（変換式ｔａ）
により対応付けられる。
【０１５４】
例えば、薄膜に対するシミュレーションについて考える。薄膜は、膜厚方向の寸法に比べて、膜面に平行な方向の寸法が大きい。このため。膜厚に平行な方向について、膜厚方向よりも繰り込み因子を大きくするとよい。例えば、膜厚方向をＺ方向とするならば、膜厚方向の繰り込み因子αＺに比べて、膜面に平行な方向の繰り込み因子αＸ及びαＹを大きく設定するとよい。
【０１５５】
なお、第２のスケーリングパラメータγは、（変換式ｍ）、（変換式ε）、及び（変換式ｔ）に、αγという形で含まれている。γ＝０と設定した場合、αγは単に「１」となる。よって、たとえばδ＝２と決定された場合、粒子の質量、相互作用係数、及び時間に関して、第２のスケーリングパラメータγを用いずに（γという第２のスケーリングパラメータを導入せずに）、
2010146368.tif 000078
2010146368.tif 000079
2010146368.tif 000080
【０１５６】
というスケール変換を施したとしても、第２のスケーリングパラメータγを０と設定したのと同様なスケール変換が行われる。なお、γ＝０、δ＝２のとき、速度ベクトルについては、
2010146368.tif 000081
【０１５７】
というスケール変換が行われる。
【０１５８】
なお、繰り込み因子α、αＸ、αＹ、及びαＺが、それぞれ１より大きければ、繰り込まれた系における粒子数が元の系の粒子数よりも少なくなる。
【０１５９】
なお、必要に応じ、シミュレーション対象の系のある部分と他の部分とで、繰り込み因子（例えばα）を異ならせることも可能である。
【０１６０】
なお、上述の実施例では、繰り込み群分子動力学の例として、融点、潜熱及び分散関係を求めた。繰り込み群分子動力学は、これらに限らない種々の物理現象のシミュレーションに用いることができる。繰り込み群分子動力学は、分子動力学を用いるので、分子動力学で扱える様々な問題に適用可能である。例えば、流れ、構造、熱、反応等の問題、さらには、機構、磨耗や潤滑の問題等を扱うことができる。繰り込み群分子動力学は、実際に分子動力学計算を実行する系より粒子数の多い系に対する物理量を求めることを可能にする。よって、本願発明者が見出した繰り込み群分子動力学を用いれば、従来シミュレーションが困難であったスケールの問題に対して、精度良いシミュレーションを行うことが可能となる。なお、繰り込み群分子動力学において実行される分子動力学計算の手法としては、公知の種々のものを用いることができる。
【０１６１】
なお、本願明細書の「背景技術」の欄に記載した非特許文献１では、繰り込みの手法を分子動力学に応用する試みがなされている。しかし、非特許文献１では、ハミルトニアンの相互作用部分の繰り込みがなされているのみであり、運動エネルギ部分の繰り込みはなされていない。非特許文献１では、運動量（または速度）に対して、繰り込みではなく、単純平均を取っている。このため、静的な解析に留まっており、温度を陽に考慮できず、強制的にエネルギの散逸を起こさせるための人為的な手法（非特許文献１の「３．数値計算法」の欄に記載の式（２３）、（２４）、（２５））が必要となっている。そのため、有限温度の計算や振動・騒音、相転移が伴う動的な問題へ適用する際、信頼性に欠ける等の課題があった。非特許文献１が開示するスケール変換則は、本願発明のγ＝０、δ＝３の場合に相当する。
【０１６２】
本願発明者のスケール変換則の導出過程においては、ハミルトニアンの運動エネルギ部分の繰り込みもなされている。これにより、運動量ベクトル（または速度ベクトル）の適正なスケール変換則が導かれる。よって、シミュレーションにおいて、温度を陽に考慮することが可能となる。
【０１６３】
以上実施例に沿って本発明を説明したが、本発明はこれらに制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【図面の簡単な説明】
【０１６４】
【図１】図１（Ａ）及び図１（Ｂ）は、スケール変換則の導出方法を説明するための図である。
【図２】繰り込み群分子動力学を用いたアルミニウムに対するシミュレーションの結果を示すグラフである。
【図３】図３（Ａ）は、Ｓ−Ｗポテンシャルについて説明するための図であり、図３（Ｂ）は、繰り込み群分子動力学を用いたシリコンに対するシミュレーションの結果を示すグラフである。
【図４】繰り込み群分子動力学を用いて分散関係を求めたシミュレーションの結果を示すグラフである。
【符号の説明】
【０１６５】
α  第１の繰り込み因子
Ｎ  粒子数
Ｎ´  繰り込まれた粒子数
【図１】
2010146368.tif 000082
【図２】
2010146368.tif 000083
【図３】
2010146368.tif 000084
【図４】
2010146368.tif 000085