(57)【要約】【課題】  物質の物理的、化学的性質を計算するためのサンプリングをシステマティックに効率よく行うことができるサンプリングシミュレーション装置を提供する。【解決手段】  決定論的微分方程式を数値積分することによって、解の軌跡を求め、この軌跡に沿ってサンプリングを行う際に、該微分方程式の解がBoltzmann-Gibbs(ＢＧ)分布より広いエネルギー範囲をカバーする複数のTsallis分布を合成した分布を再現するように設定する。Tsallis分布のパラメータには、いくつかの値を設定し、それぞれ異なったパラメータ値に対応するTsallis分布を合成したものを作る。複数のTsallis分布を合成した分布から得られた軌跡からサンプリングしたサンプリング点を使って、ＢＧ分布に従う物理系の物理的及び化学的性質を、統計力学の方法によって、計算する。【特許請求の範囲】
【請求項１】
決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション装置において、
Boltzmann-Gibbs分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成する複数分布生成手段と、
数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分する数値積分手段と、
該数値積分の結果得られた軌跡に沿ってサンプリングするサンプリング手段と、
を備えることを特徴とするサンプリングシミュレーション装置。
【請求項２】
前記エネルギー分布がBoltzmann-Gibbs分布より広い分布は、Tsallis分布であることを特徴とする請求項１に記載のサンプリングシミュレーション装置。
【請求項３】
複数の異なるパラメータについて生成される分布は、温度に対応するパラメータのみを異なる値に設定することによって得られることを特徴とする請求項１に記載のサンプリングシミュレーション装置。
【請求項４】
前記合成された分布は、Boltzmann-Gibbs分布と、Boltzmann-Gibbs分布よりエネルギー分布が広い分布を合成して得られることを特徴とする請求項１に記載のサンプリングシミュレーション装置。
【請求項５】
決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション方法において、
Boltzmann-Gibbs分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成し、
数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分し、
該数値積分の結果得られた軌跡に沿ってサンプリングする、
ことを特徴とするサンプリングシミュレーション方法。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、物質の物理的、化学的な性質を計算するための決定論的方式によるサンプリングシミュレーション装置に関する。
【背景技術】
【０００２】
物質の物理的・化学的な性質は、その系を構成する系の自由度と同じ数のｎ個の座標ｘ＝（ｘ１、・・・、ｘｎ）およびｎ個の運動量ｐ＝（ｐ１、・・・、ｐｎ）で記述することができる。多くの系は、構成要素であるＮ個の原子の座標で決まるポテンシャルエネルギー関数と運動量から計算される運動エネルギー関数の和であるエネルギー関数で定義できる。従ってこれらの関数形を与えてやれば、計算機シミュレーションにて物質の性質を調べることが可能である。これにより、通常行われる実験に対する補完、さらに進んで、実験困難な状況に対する対応、実験に先立つ予測や、大量・分散処理等による実験コスト削減、といったことが期待され、多くの研究開発がなされてきている。ここで、特に注目するのは、薬物の対象となる分子と、それが働く環境を構成する分子とで表される古典力学系である。このような系におけるポテンシャルエネルギー関数の典型的な例は、次式のようなものである。
【０００３】
【数１】
2006221578.tif 000002
【０００４】
第１項が、各原子間の距離で決まるようなポテンシャルエネルギーを表し、第２項は、分子内原子の結合長、結合角などの平衡値を表現するようなポテンシャルエネルギーの総和である。
【０００５】
物理的・化学的な性質を知るためには、Boltzmann-Gibbs (BG)分布などで代表される分布における、物理量の期待値で表わされるような熱力学量を計算する必要がある。例えば、固体の典型的な熱力学量としては、比熱、弾性係数などがある。一方、生体高分子である蛋白質がとる安定な立体構造は、各蛋白質を規定する固有のアミノ酸配列をもつポリペプチド鎖が構築できる様々な立体構造の中から、ＢＧ分布に基づく自由エネルギーの最も低い状態が選択されていると考えられる。また、薬物と受容体蛋白質との親和力は、両者が結合した複合体として存在する系と、両者が離散した場合の系との自由エネルギーの差によって定量的に説明される。このように、物理量の期待値を算出して熱力学量を求めたり、自由エネルギーを算出して蛋白質の立体構造や薬物と受容体蛋白質の親和力を求めたりする際に重要なのが、（ｘ、ｐ）で表わされる状態を効率的にサンプリングすることである。
【０００６】
状態サンプリングによく用いられる方法として、モンテカルロ(MC)法がある。アルゴリズムが簡単で、一般化も進んでいる。しかし、ここで主に対象とするような複数分子で表わされる系を、この方法で取り扱おうとすると、局所変位によるエネルギー増大が起きやすく、従って効率的な状態更新の受理がされにくいという困難を生じる。一方、分子動力学(MD)法では、複数分子で構成されるような現実的な系に対しても対応できる。従って、MD法により、状態サンプリングすることが望ましい。所謂、定温の BG分布を実現するMDの手法としては、Nose-Hoover (NH)の手法（非特許文献１、２）やその発展形がある（非特許文献３、４）。
【０００７】
しかしながら、自由度が高く、上記(0)式で表される系などは、一般的にエネルギー面が複雑であり、その状態空間を効率的にサンプリングすることは従来のMD法技術では困難であった。例えば、任意に選んだ初期構造に近い構造に非常にしばしばtrap されてしまい、正しい構造予測ができない、といういわゆる局所トラップの問題がある。
【０００８】
これを解決する手段として、Tsallis 分布密度(規格化因子を除く)（非特許文献５、６）、
【０００９】
【数２】
2006221578.tif 000003
【００１０】
がエネルギーEに関し冪則的であり、ｑ＞１の場合は、指数関数的なBG分布密度
【００１１】
【数３】
2006221578.tif 000004
【００１２】
と比較してなだらかであるので局所トラップ回避に有効であることに着目したのが、特許文献１であり、(1)式の分布を決定論的方程式で実現する（非特許文献７）。ここで、Ｅ（ｘ，ｐ）＝Ｕ（ｘ）＋Ｋ（ｐ）は全エネルギーであり、Ｕはポテンシャルエネルギー, Ｋは運動エネルギー
【００１３】
【数４】
2006221578.tif 000005
【００１４】
であり(質量は全て1とする),β＝１／ｋＴ、Ｔは温度に相当するパラメータ(以下、温度)、ｋはBoltzmann定数である。ｑは実数パラメータであり、q,βの値は、Ｅの値の変域に応じて(1)式が実数として定義されるような範囲で与えられる。ｑ→１の極限が、カノニカル分布を特徴付けるものとして従来から用いられてきた BG分布の密度(2)式に等しくなる。特許文献１の技術を用いて物理系のシミュレーションが可能になることが非特許文献７に示されている。さらに、従来のシミュレーション手法では、その状態サンプリング性能に限界があり、BG分布が生成できない、といった取り扱いが困難な系（非特許文献８、９）に対しても、効率良く正確な結果を与えること、またpeptide 系においてエネルギー空間の幅広いサンプリングが可能なことが示された（非特許文献１０、１１）。
【００１５】
他のMDサンプリング手法としては、Replica exchange molecular dynamics (REMD)（非特許文献１２） やthe multicanonical molecular dynamics (MCMD)（非特許文献１３）がある。
【００１６】
他方、MDではないが、GA法などのいわゆるヒューリスックな最適化手法では、上記のような局所trapを回避することに対しては有効である。
【非特許文献１】S. Nose, J. Chem. Phys. , 511 (1984).
【非特許文献２】W. G. Hoover, Phys. Rev. A , 1695 (1985).
【非特許文献３】S.Nose, Prog. Theor. Phys. Suppl. , 1 (1991).
【非特許文献４】W. G. Hoover and B. L. Holian, Phys. Lett. A , 253 (1996), and the references therein.
【非特許文献５】C. Tsallis, J. Stat. Phys. , 479 (1988).
【非特許文献６】C. Tsallis, Braz. J. Phys. , 1 (1999); for an updated bibliography, cf. http://tsallis.cat.cbpf.br/biblio.htm.
【非特許文献７】I. Fukuda and H. Nakamura, Phys. Rev. E 65 (2002) 026105.
【非特許文献８】D. Kusnezov, A. Bulgac, and W. Bauer, Ann. of Phys. , 155 (1990).
【非特許文献９】I. L’Heureux and I. Hamilton, Phys. Rev. E , 1411 (1993).
【非特許文献１０】I. Fukuda and H. Nakamura, Chem. Phys. Lett. , 367 (2003).
【非特許文献１１】I. Fukuda and H. Nakamura, J. Phys. Chem. B, , 4162 (2004).
【非特許文献１２】Y. Sugita and Y. Okamoto, Chem. Phys. Lett. , 141 (1999).
【非特許文献１３】N.Nakajima, H. Nakamura, and A.Kidera, J. Phys. Chem. B 101 (1997) 817.
【特許文献１】特開２００３−４４５２４号公報
【発明の開示】
【発明が解決しようとする課題】
【００１７】
・特許文献１では、Tsallis分布を生成する決定論的方程式を導き、それを用いて物理系のシミュレーションを可能にする技術が示された。しかしながら、本方程式を状態サンプリング手法として用いる場合、最も有効な結果を与えるTsallis分布のパラメータ値（文献１）をどのようにして求めるかという問題が残っている。
【００１８】
このパラメータの設定・選択は、通常の系では非常に注意深く行われる必要があり、例えばパラメータを少し動かすだけで、実効的にカバーするエネルギー領域が大きく変わってしまう、という場合があった。これに対し、ヒューリスティックな方法や試行錯誤的方法も可能であるが、対象とする系が大きくなると、計算コストの増大が避けられない。したがって、よりスムースな応用を目指すためには、システマティックな方法が必要になってくる。
【００１９】
また、たとえ、最良と思われるパラメータが見つかったとしても、さらに効率的に広いサンプリングが可能となるような分布が求めよ、といった要請には容易には対処できなかった。
・REMDでは、いくつかのBG分布を、予め定めた適当なタイミングで、適当な確率をもって交換、といった作業が必要である。したがって、タイミングを決めるパラメータ値の設定の違いによるシミュレーション結果への依存性といった問題がある。
・GA法などの最適化手法では、上記のような局所trapを回避することに対しては有効かもしれないが、統計力学的にはっきり定まった分布を生成しているわけではないので、直接には、熱力学量の計算などができない。そのため、現実の問題、例えば、新規薬物の受容体蛋白質への親和力を推定するin silico 薬物設計などには、容易には対応できない。
【００２０】
本発明の課題は、物質の物理的、化学的性質を計算するためのサンプリングをシステマティックに効率よく行うことができるサンプリングシミュレーション装置を提供することである。
【課題を解決するための手段】
【００２１】
本発明のサンプリングシミュレーション装置は、決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション装置において、BG分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成する複数分布生成手段と、数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分する数値積分手段と、該数値積分の結果得られた軌跡に沿ってサンプリングするサンプリング手段とを備えることを特徴とする。
【発明の効果】
【００２２】
本発明によれば、薬物やその他の蛋白質のような複雑な物理系の物理的、化学的性質を精度良く計算することを可能とするサンプリングをシスタマティックに実行することができるサンプリングシミュレーション装置を提供することができる。
【発明を実施するための最良の形態】
【００２３】
本発明の実施形態では、決定論的方程式を数値的に解くことによって得られる軌跡に沿ってサンプリングを行う手法を適用するが、特に、本発明の実施形態のシミュレーション技術により、この軌跡が局所trap されることを回避し、かつ、熱力学量の計算を可能にする。本発明の実施形態では、適切なパラメータ値をひとつ選ぶ代わりに、各パラメータ値を伴った多数の分布を同時に生成するという方式を採用し、それを実現する手法を提供する。
【００２４】
対象となる分布がＭ個あるとし、その密度関数を、ρａ：ＲＮ→Ｒ、ａ＝１、・・・、Ｍとする。物理系を運動エネルギーＫ（ｐ）≡（１／２）‖ｐ‖２とポテンシャルエネルギーＵ（ｘ）であらわすことにしているので、ρａを次のような形で与えられる分布となるようにする。
【００２５】
【数５】
2006221578.tif 000006
【００２６】
ここでρＺは新たに導入された実変数ζの適当な実数値関数である。従って、目的は、
【００２７】
【数６】
2006221578.tif 000007
【００２８】
の形の分布密度を実現することである。このために、次のような常微分方程式を解く：
【００２９】
【数７】
2006221578.tif 000008
【００３０】
ここで、
【００３１】
【数８】
2006221578.tif 000009
【００３２】
である。Ｔは温度の次元をもつ定数であり、各物理変数の次元を自然なものにし、時間発展のスピードも定める。以下簡単のため、ＲＮ≡Ｒ２ｎ＋１≡Γとする。
次に、各ρａが
【００３３】
【数９】
2006221578.tif 000010
【００３４】
の形で規格化されている場合を考える。ここで主に対象にしている分子系にて、このＺａを求めるのは非常に困難である場合がしばしばであるので、次のような対策をとる：
(a)
【００３５】
【数１０】
2006221578.tif 000011
【００３６】
の代わりに
【００３７】
【数１１】
2006221578.tif 000012
【００３８】
を用いる。
(b)ＹａはＺａを計算するより容易であり、一般的な方法として、free energy perturbation methodや thermodynamic integration methodなどが知られている他、本発明の実施形態で改めて提案するような方法（後述）も用いることができる。 
(c)ρから
2006221578.tif 000013
への変更により、式(4)-(8)で変更を受けるのは、(7)のみで、それは以下のようになる：
【００３９】
【数１２】
2006221578.tif 000014
【００４０】
<(b)の詳細>
free energy perturbation methodや thermodynamic integration methodにしても、以下に述べる簡便な方法においても、Ｙａ≡Ｚａ／Ｚｃを求める場合は、分布
2006221578.tif 000015
が十分に重なっていることが必要である。このためにまず、各ａ＝１、・・・Ｍ−１について
2006221578.tif 000016
が十分に重なるように並べられたと仮定する。するとＹ１＝１(i.e.,Ｚｃ＝Ｚ１)とすることにより、
【００４１】
【数１３】
2006221578.tif 000017
【００４２】
として求めることができる。このＺｊ／Ｚｊ−１を求めるための簡便な方法として、
【００４３】
【数１４】
2006221578.tif 000018
【００４４】
を用いることができる。ここで、Ｐａはρａから得られたエネルギー分布密度である。
さて、前述したようにＹａを求める際には、分布間の重なりが重要になってくる。Tsallis 分布[(1)式]では、通常対象となるBG分布に比べ、分布の広がりが大きいので、これを複数個利用すれば分布間の重なりをより簡便に実現することができるので有利である。そこで、分布ρａの具体形としてTsallis 分布を用いたものを考察する。この場合、(9)式は次のようになる。
【００４５】
【数１５】
2006221578.tif 000019
【００４６】
ここで、
【００４７】
【数１６】
2006221578.tif 000020
【００４８】
であり、αａ≡（ｑａ−１）／Ｔａ、及び、ｂａ≡１／（１−ｑａ）［もしくは、ｑａ／（１−ｑａ）］とおき、εａは１＋αａ（ｅ＋εａ）＞０を保障するために導入したパラメータである（文献２、３）。
【００４９】
実際にどのようにεａ、ｑａ、Ｔａを設定すればよいかは、後に、具体例をもって説明する。
本発明の、式(4)-(6), (8), (13)で、Ｍ≡１とおくことは、単一のTsallis 分布を生成することになるが、そのときこれらの式[でα１≡（ｑ−１）／Ｔ、ｂ１≡ｑ／（１−ｑ）、ε１≡０としたもの]は、特許文献１で開示した方程式
【００５０】
【数１７】
2006221578.tif 000021
【００５１】
に本質的に同じ（上式右辺にTをかけたものに同じ）になることがわかる。
なお、ここでは分布の広がりが十分確保できることが重要であり、必ずしもTsallis分布を使用しなくてはならないということではない。したがって、式(4)-(9)は、複数のTsallis分布の適用に限定したものではない。また、Tsallis分布であっても、（１）式の他にいくつかの種類があり、それらを利用することも勿論可能である。
【００５２】
<本発明の実施形態の動作>
Liouvilleの定理 により、測度
2006221578.tif 000022
が式(4)-(6), (8), (9)[もしくは(13)]の常微分方程式の定める位相空間内の解の流れの不変測度になっていることが証明できる。従って、適当な数学的条件の仮定の下で、Birkhoff のergodic Theorem により、
2006221578.tif 000023
に関して殆ど全ての初期値に関し、
【００５３】
【数１８】
2006221578.tif 000024
【００５４】
な関数ｆの長時間平均が存在する。もし系がこの測度に関してergodic ならば、
【００５５】
【数１９】
2006221578.tif 000025
【００５６】
となる。つまりｆの長時間平均は、‘multi’分布
【００５７】
【数２０】
2006221578.tif 000026
【００５８】
における期待値に等しくなる。また、任意領域Ａの実現確率は、上式と
【００５９】
【数２１】
2006221578.tif 000027
【００６０】
を用いて、
【００６１】
【数２２】
2006221578.tif 000028
【００６２】
とあらわすことができるので、状態ω∈Ωの実現確率密度は、与えた‘multi’分布ρに比例することになる。これらが、式(4)-(6), (8), (9)が‘multi’分布
【００６３】
【数２３】
2006221578.tif 000029
【００６４】
を実現することを表す。
同様に、ポテンシャルエネルギー分布密度ＰＵ（ｕ）については、各ａのρａに対応するポテンシャルエネルギー分布密度をＰａＵ（ｕ）とすると、以下のようになることが証明される：
【００６５】
【数２４】
2006221578.tif 000030
【００６６】
これは、本手法で得るポテンシャルエネルギー分布密度が、各ａのポテンシャルエネルギー分布密度の平均になることを示している。
（ｘ、ｙ）の関数である物理量Ｏ（ｘ、ｙ）に関しては、
【００６７】
【数２５】
2006221578.tif 000031
【００６８】
ただし、ここで、
【００６９】
【数２６】
2006221578.tif 000032
【００７０】
とした。従ってまた、BG分布を生成するためには、
【００７１】
【数２７】
2006221578.tif 000033
【００７２】
となることが示せるので、最左辺の長時間平均を見積もればよい。ただし、
【００７３】
【数２８】
2006221578.tif 000034
【００７４】
とした（Ｙａ＝Ｚａ／Ｚｃ＝Ｚ’ａ／Ｚ’ｃを使った）。ここで、ρＢＧは（２）式のものであるが、逆温度パラメータβは任意の正値にとれる。また、異なったβについても並列に計算できるので、ひとつのシミュレーションで任意個数のBG分布の再現が可能となる。
【００７５】
前述の手法（特許文献１）では、有効な分布ρａ／Ｚａをひとつ選択する（つまりρａを定義しているパラメータをひとつ設定する）ものであった。一方、本手法では、有効な分布ρａ／Ｚａをひとつ選択するのではなく、多数の分布
【００７６】
【数２９】
2006221578.tif 000035
【００７７】
を（同時に）生成する。そのため、最良のサンプリング結果をもたらす分布を狙って、パラメータを注意深くひとつ選択するという労力が削減される。さらに、実効的にカバーするエネルギー領域が異なる分布を任意個数加え合わせることができるので、容易に広い領域のサンプリングが可能になるという効果がある。
【００７８】
また、本手法では「同時に」生成するため、上述のREMDのように、いくつかの分布を、適当なタイミングで切り替える、といった作業が不要である。従って、タイミングを決めるパラメータ値の設定の違いによるシミュレーション結果への依存性といった問題はなくなる。
【００７９】
以下に、アラニンぺプチド、Ac-Ala-Ala-NMe (AcとNMe は、それぞれAcetyl基と N-Methyl基)からなる系に本手法を適用した例を示す。ポテンシャル関数は、AMBER force-field（文献４、５）を用いた。
【００８０】
【数３０】
2006221578.tif 000036
【００８１】
とし、Ｔ＝２００Kとし5ns のシミュレーションを行った。
この系のシミュレーション例で、εａ、ｑａ、Ｔａの具体的設定例を述べる。まず、εａと ｑａは、ａによらず以下のような一定値を与え、異なるＴａに関しての分布を加え合わせる。  
【００８２】
【数３１】
2006221578.tif 000037
【００８３】
とする。下限は厳密には求められないが、これにより、通常（マージンがあるので）は、
【００８４】
【数３２】
2006221578.tif 000038
【００８５】
が保障され、各ａに対し密度（(1)式）が良く定義される。この候補値は、低温度のBG分布を生成する（通常の）MDのシミュレーション結果から得たポテンシャルエネルギー値の最小値とする。本例では、
2006221578.tif 000039
【００８６】
Kのシミュレーション結果から
2006221578.tif 000040
を得た。
【００８７】
【数３３】
2006221578.tif 000041
【００８８】
とする。単一のTsallis 分布に対しｑ＝１＋１／ｎと設定することは、Hansmann-Okamoto（文献６）が最初に提案したものであり、BG分布に比べ一般に広い分布を生成することが可能である。従って、本実施形態では、この設定方法を使用する。
(iii) どのようなＴａを足し合わせるかは、最終的におのおの分布の重なりを確認する必要があるので、予め、具体的に単一のTsallis 分布を計算するほうが良い。
【００８９】
図１は、複数のＴの値に対するTsallisエネルギー分布を示す図である。
図１に、(i),(ii)の設定の下、いくつかのＴａでの単一のTsallis分布の全エネルギー分布密度を、特許文献１のシミュレーション手法を用いて算出したものを示した。Ｔａ＝１５Ｋの分布とＴａ＝６０Ｋの分布で十分な重なりが確認され、それらの２つの和で十分広いエネルギー領域がカバーできるので、本実施形態の例では、Ｍ＝２；Ｔ１＝１５Ｋ、Ｔ２＝６０Ｋとした。
【００９０】
このＴａの設定をスムースに行うために、以下に述べるようなＴａの基準値を定めることができる。ある基準値が求まれば、Tsallisエネルギー分布密度Ｐの値が大きい部分は、Ｔａの増加とともに右側に移動する、という経験則に基づいて、ある基準値の近傍でTsallisエネルギー分布密度がどのように変化するかを考察することが可能になる。前記は経験則ではあるが、エネルギー状態密度関数に適当な仮定、例えばΩ（Ｅ）＝（Ｅ−Ｅｍｉｎ）ν、νは適当な冪数、とすれば、エネルギー期待値がＴの増加関数となることが知られている。
【００９１】
Ｔａの基準値を定めるためにいくつか仮定をする。温度１／（ｋＢβ）のBGエネルギー分布密度ＰＢＧについてｌｎＰＢＧが上に凸な関数で、エネルギーＥ１＜Ｅ２でＰＢＧ（Ｅ１）＝ＰＢＧ（Ｅ２）だと仮定する。単一のTsallisエネルギー分布密度ＰについてもｌｎＰが上に凸と仮定する。このとき、もし、Ｐ（Ｅ１）＝Ｐ（Ｅ２）が成立すれば、ＰＢＧの最大値もＰの最大値もＥ１とＥ２の間に存在することが言える。これは、Ｐの有効な台がＰＢＧの有効な台の主領域を十分含んでいることを示唆する。従って、Ｐ（Ｅ１）＝Ｐ（Ｅ２）が成立するようなＴを基準値とすることができる。このようなＴは以下となることが証明される：
【００９２】
【数３４】
2006221578.tif 000042
【００９３】
すなわち、(i)と(ii)により設定したεとｑを用いた(22)式のＴをＴａの基準値とする。ペプチド系の例では、NH法（非特許文献１、２）を用いて、１／（ｋＢβ）＝３００ＫのBG分布を発生させ,
【００９４】
【数３５】
2006221578.tif 000043
【００９５】
(エネルギー平均値＋その標準偏差値) とし、ＰＢＧ（Ｅ１）＝ＰＢＧ（Ｅ２）となるようなＥ１、により、Ｔ＝１９．９Ｋを得た。
図２に、３００ＫのBGエネルギー分布と、Ｔ＝１９．９Ｋの単一の Tsallisエネルギー分布Ｐを示す。
【００９６】
非常によい近似でＰ（Ｅ１）＝Ｐ（Ｅ２）が成立し、Ｐの有効な台がＰＢＧの有効な台を十分含んでいることがわかる。これは、前記で述べたTaの基準値の定め方の妥当性を示すものである。
【００９７】
図３に、２つのTsallisポテンシャルエネルギー分布と、これらの平均の分布、及び２つのTsallis分布を(i)-(iii)に従って合成させた場合のポテンシャルエネルギー分布(‘multi’と記した)を示す。図４に、２つのTsallis分布とそれを合成した分布、及び各温度でのＢＧ分布を示す。
【００９８】
図３より、合成した分布は、単一のTsallis 分布の平均値とよく一致しており、式(12)が成立していることを示す。また、非常に広い温度領域をカバーしている(図４参照)。これは、単一のTsallis 分布を考えていただけでは、容易にはカバーできない領域であり、これが精度の良い物理的、化学的性質の計算に必要なサンプリングに有効に働く。
【００９９】
また、本発明の実施形態では、異なるエネルギー領域間の遷移が、自動的に起こる。遷移のタイミングを人工的に設定しないのでこれによる artifact がない。これは、REMD法 と異なるものである。
【０１００】
図５に、本実施形態で得た、ポテンシャルエネルギー値の軌跡を示す。これは、単一のTsallis 分布を生成させたときのポテンシャルエネルギー値の軌跡と本質的に異なるものであり、各々の主だった領域間で自然に起こっているジャンプが観察される。すなわち、図５（ａ）は、温度パラメータが小さく設定された単一Tsallis分布のポテンシャルエネルギー値の軌跡を示したものであり、軌跡がポテンシャルエネルギーの低い部分に局在しているのが分かる。また、図５（ｂ）は、温度パラメータが大きく設定された単一Tsallis分布のポテンシャルエネルギー値の軌跡を示したものであり、軌跡がポテンシャルエネルギーの低い部分には現れず、ポテンシャルエネルギーの高い部分に広がっているのが分かる。図５（ｃ）は、本発明の実施形態にしたがって温度パラメータが大きいTsallis分布と小さいTsallis分布を合成した分布の軌跡を示す。図５（ｃ）から分かるように、ポテンシャルエネルギーの大きい部分から小さい部分まで、非常に広い範囲にわたって軌跡が広がっているのがわかる。
【０１０１】
次の計算例は、上述の２つのTsallis 分布にさらに、３００K のBG分布を足したものである。実際に計算したいのが３００Kの状態だとすれば、このような合成の仕方は、特有の効果がある。すなわち、３００KのBG分布が支配する領域を効率的にサンプリングするためには、通常複雑なエネルギー面形状から規定される位相空間内の非連結領域のサンプリングが必要であり、そのためには、その領域より低エネルギーの領域、および高エネルギーの領域の両方への到達を達成する必要があるのである。
【０１０２】
上述の（ｑ、ε、Ｔ１）のTsallis 分布、（ｑ、ε、Ｔ２）のTsallis 分布、３００K のBG分布、の3つの分布を本手法により合成した結果を図６に示す。
図６によれば、本手法の結果が、3つの分布の平均を与えていることがわかる。そして、３００K のBG分布の領域の達成は勿論、（ｑ、ε、Ｔ１）のTsallis 分布がカバーする低エネルギー領域への到達と、（ｑ、ε、Ｔ２）に対応する高エネルギー領域への到達が共に達成されていることがわかる。
【０１０３】
図７は、本発明の実施形態に従ったサンプリングシミュレーション装置への入力データを説明する図である。
自由度数ｎは、シミュレーション対象の物理系の自由度の数である。density数は、いくつのTsallis分布を合成するかを示す数である。multi-Tsallis分布パラメータは、各Tsallis分布を規定するパラメータであり、それぞれのTsallis分布にそれぞれ（ｑa、εa、Ｔa)の組を設定する。なお、これらのパラメータをそれぞれ別個に設定する代わりに、上述した実施形態のように、ｑとεを全てに共通して設定し、ＴをそれぞれのTsallis分布に個別に設定するようにしても良い。
【０１０４】
次に、シミュレーション条件として、前述の常微分方程式の解を数値的に求める場合の時間のステップ数、時間刻み幅を設定し、同時に境界条件も設定する。また、ρｚ関数のパラメータとして、関数番号、設定係数を設定する。この場合、ｚ−密度値計算部４１に複数の異なる関数を登録しておき、関数番号により、使用する関数を選択する。また、ＢＧ分布の温度リストとして、計算するＢＧ分布の温度の数と、温度の値を設定する。度数計算パラメータは、求めるエネルギー分布のエネルギーの値の離散幅と変域のパラメータを設定する。更に、モニター変量として、特定の２原子間の距離や、特定の３原子間で定義される角度などを設定する。
【０１０５】
シミュレーション制御変数としては、出力タイミング等の出力制御やリスタート制御のための変数を設定する。また、運動エネルギーデータとしては、各粒子の質量等を設定する。ポテンシャル関数Ｕのデータとしては、関数形及びパラメータ、長距離力のカットオフの方法やカットオフ長などを設定する。また、各自由度の初期値として、座標、運動量、ζの初期値を設定する。
【０１０６】
図８は、ＢＧ分布を生成するための入力データを示す図である。
自由度の数ｎ、ポテンシャル関数Ｕのデータとして、関数形及びパラメータ、長距離力カットオフの方法及びカットオフ長、シミュレーション条件であるステップ数、時間刻み幅、境界条件、度数計算パラメータとしての離散幅、変域パラメータ、温度Ｔまたは逆温度β、各自由度の初期値としての座標、運動量、運動エネルギーデータとしての各粒子の質量等を設定する。
【０１０７】
図９は、パラメータ列生成部入力データを示す図である。
自由度の数ｎ、後述のＵｉｎｆ計算部で発生させる、通常低温に設定するＢＧ分布の逆温度β０、Ｅ１、Ｅ２計算部で発生させる、通常常温に設定するＢＧ分布の逆温度β、発生させるTsallis分布の数Ｍ、分布重なり計算部で用いる閾値Ｗａ、並び替え部で用いる閾値ｕ、Ｔａ計算部で用いるパラメータＣ、ポテンシャルエネルギーの最小値の候補Ｕｉｎｆを設定する。
【０１０８】
図１０は、本発明の実施形態に従った出力データを示す図である。
ポテンシャルエネルギー等の物理量やζ等の変数値の時間変化、変数値の期待値（有限時間平均）の時間変化、また、各温度におけるＢＧ分布におけるエネルギーや比熱等の物理量の期待値（有限時間平均）の時間変化及び最終値、ポテンシャルエネルギー、運動エネルギー、全エネルギー、モニタする物理量である、特定の２原子間の距離、特定の３原子間で定義される角度等の変数値の実現確率（頻度）が出力となる。
【０１０９】
図１１は、本発明の実施形態のシミュレーション装置のブロック図である。
入力部１０へは、入力データとTsallis-parameter列生成部からデータが入力され、計算部１１に送られる。計算部１１では、微分方程式を逐次数値解法により積分し、物理量の値及び期待値を計算する。計算結果は、出力部１２から出力データとして、出力される。
【０１１０】
図１２は、入力部のブロック図である。
入力部１０では、入力データをデータ入力部１３で受け、ファイル出力準備部１４で、必要なファイルの準備を行って、データの整合性をエラーチェック部１５でチェックして、計算部１１に渡す。
【０１１１】
図１３は、Tsallis-parameter列生成部の構成及びフローを示す図である。
Ｕｉｎｆ計算部では、前述のβ０をＢＧＤ入力部２０で受け、ＢＧＤ計算部２１でNose-Hoover法等に従い、シミュレーションを行って、ＢＧ分布を発生させる。Ｕｉｎｆ計算部２２でポテンシャルエネルギーの最小値を計算しこれをＵｉｎｆとする。
【０１１２】
Ｕｉｎｆと自由度数ｎと、上述の（ｉ）と（ｉｉ）で求められたｑとεは、ＴＤ入力部２５とＴｒｅｆ計算部２３に入力される。また、前述のβは、Ｅ１、Ｅ２計算部のＢＧＤ入力部２８とＴｒｅｆ計算部２３に入力される。Ｅ１、Ｅ２計算部では、βから、ＢＧＤ計算部２９で、Nose-Hoover法等によって、シミュレーションを行い、ＢＧ分布を発生させ、ＢＧＤ出力部３０がＢＧエネルギー分布を出力する。そして、ブロック３１において、エネルギー分布の平均、分散等から、Ｅ１、Ｅ２を計算し、結果をＴｒｅｆ計算部２３に入力する。
【０１１３】
Ｔｒｅｆ計算部２３では、温度変数の基準値Ｔｒｅｆを式（２２）にしたがって算出し、Ｔａ計算部２４において、Ｔｒｅｆ、Ｃ、Ｍより、複数の温度変数Ｔ１〜ＴＭを求める。これらの複数の温度変数は、ＴＤ入力部２５に送られる。
【０１１４】
ＴＤ入力部２５で入力された(ｑ1、ε1、Ｔ1),…, (ｑM、εM、ＴM)の各々のパラメータ値に対し、ＴＤ計算部２６で、Tsallis dynamics（ＴＤ）法等に従い、シミュレーションを行い、Tsallis分布を発生させ、ＴＤ出力部２７において、Tsallisエネルギー分布を出力する。出力された複数のTsallis分布は、前述のｕに基づいて並び替え部３２において並び替えられ、分布重なり計算部３３において、前述のｗａに基づいて重なり部分が計算され、表示部３４が生成された複数のTsallis分布を表示する。計算された重なり部分の結果を基にさらにユーザが目でTsallis分布の重なりを確認して、重なりが十分か否かを判断する。重なりが十分でない場合には、Ｍを大きくして、Ｔａ計算部２４で再度複数の温度パラメータを計算して、Tsallis分布を計算しなおす。重なりが十分である場合には、各Tsallis分布の（ｑ、ε、Ｔ）の変数組を決定し、Ｚｉの比で定義されるＹａをＹａ計算部３５で計算して、入力部１０へ入力する。
【０１１５】
Ｔａ計算部２４では、Ｔ＝Ｔｒｅｆを基準値として、Ｔ１〜ＴＭを生成する。例えば、０＜Ｃ×Ｔｒｅｆを（Ｍ−１）等分して、各々をＴ１、Ｔ２、・・・ＴＭとする。Ｃは入力値である。例えば、Ｃ＝２の場合は、Ｔｒｅｆを中心とした等分割されたＭ個のＴａを得る。
【０１１６】
分布重なり計算部３３では、例えば、隣り合う２つのエネルギー分布値があるＥで、Ｐａ（Ｅ）＝Ｐａ−１（Ｅ）となるとき、Ｐａ（Ｅ）がＰａに対し相対的に大きいかどうかを調べる。同様にａ−１についても調べる。これをａ＝１〜Ｍについて行う。例えば、ｗａ＝Ｐａ（Ｅ）／Ｐａ（Ｅｍａｘ）と所定の閾値ｗとの比較を行う。ｗとしては、例えば、０．５と設定する。あるいは、全てのエネルギー分布のペアに対し、エネルギー分布値があるＥでＰａ（Ｅ）＝Ｐｂ（Ｅ）となるとき、Ｐａ（Ｅ）がＰａに対し、相対的に大きいかどうかを上記と同様に調べ、大きなＭの場合、いらない分布を除去する。
【０１１７】
Ｙａ計算部３５は、（Ｐｊ−１（Ｅ）／Ｐｊ（Ｅ））（ρｊ（Ｅ）／ρｊ−１（Ｅ））＝Ｚｊ−１／Ｚｊにしたがって、分配関数比を求め、Ｙａを計算する。
並び替え部３２は、以下の処理を行う。
（ｉ）ｌｎＰａ（Ｅ）＝ｕ（閾値）なるＥ（通常２つある）の最小値Ｅａが一番小さいａを１番とする。
（ｉｉ）｛２、・・・、Ｍ｝内のａで、ＰａとＰ１の交点Ｅａ（複数あれば最小値）が最小を取るＰａをＰ２とする。
（ｉｉｉ）｛３、・・・、Ｍ｝内のａでＰａとＰ２の交点Ｅａ（複数あれば最小値）が最小を取るＰａをＰ３とする。
（ｉｖ）以下同様に、Ｐ４、・・・、ＰＭを決定する。
【０１１８】
図１４は、計算部の構成及び処理フローである。
入力部１０からのデータを得て、数値積分が開始される。ステップＳ１０において、ｔをｔ＋Δｔとする。ステップＳ１１において、数値積分を実行する。この段階において、積分手法を選択・追加可能とする。数値積分においては、ポテンシャル関数計算部４０がポテンシャル関数及びその偏導関数の値を計算し、方程式右辺計算部４３が、数値積分すべき方程式のポテンシャル関数を含む右辺を計算する。このとき、ｚ−密度値計算部４１が、ρｚを計算する。このとき、関数の選択・追加を可能とする。この計算結果は、方程式右辺計算部４３に入力される。ポテンシャル関数計算部４０の計算結果は、エネルギー計算部４２に入力され、エネルギーの計算に使用される。数値積分を実行すると、ステップＳ１２において、境界条件の処理を行い、ステップＳ１３において、エネルギー計算を行い、ステップＳ１４において、密度比計算を行う。ステップＳ１４においては、Tsallis分布密度値計算部４４とＢＧ分布密度値計算部４５のそれぞれが計算した、Tsallis分布密度値とＢＧ分布密度値を使用する。ステップＳ１５において、数値積分正常に終了したか否かをチェックする。ステップＳ１６において、数値積分の結果の数値にエラーが含まれていないかを判断する。ステップＳ１６の判断の結果、エラーがあるとステップＳ２０に進んで、数値積分を終了し、最終出力部４６がエラーを出力データとして出力して終了する。ステップＳ１６の判断の結果、エラーが無いと判断された場合には、ステップＳ１７において、度数計算を行う。すなわち、数値計算の結果からサンプリングを行い、ポテンシャルエネルギーＵ、運動エネルギーＫ、モニターする物理量等の頻度を計算する。逐次出力部４７が計算結果を逐次、出力データとして出力すると共に、ステップＳ１９において、終了条件を判断する。ステップＳ１９の判断の結果、終了でないと判断された場合には、ステップＳ１０に戻って、計算を続ける。ステップＳ１９の判断の結果、計算が終了する条件が満たされたと判断された場合には、ステップＳ２０において、数値積分を終了し、最終出力部４６が出力データを出力して、数値計算を終了する。
【０１１９】
図１５は、逐次出力部の構成を示す図である。
計算部１１から計算結果を取得した逐次出力部４７は、変数値を算出し、変数値の期待値（有限時間平均）を算出し、各温度におけるＢＧ分布における物理量の期待値（有限時間平均）を算出し、出力データの可視化を行って、出力データとして出力する。
【０１２０】
また、ポテンシャル関数値計算部４０は、更新されたｘ（ｔ）に基づき、ポテンシャル関数値、ポテンシャル関数の偏導関数値を求め、また、必要に応じて解を拘束するためのポテンシャル関数の追加等を行う。このとき、関数の選択・追加が可能となるように、ポテンシャル関数値計算部４０は構成される。
【０１２１】
図１６は、最終出力部の構成を示す図である。
計算部１１から計算結果を受け取った最終出力部４６は、変数値の実現確率（頻度）を算出し、エラー処理を行い、出力データを可視化して、出力データを出力する。
参考文献：
文献１： I. Fukuda and H. Nakamura, in Proceedings of the Third International Symposium on Slow Dynamics in Complex Systems, Sendai, 2003, edited by M. Tokuyama and I. Oppenheim (AIP, 2004), pp. 356-357. 
文献２：  A. R. Plastino and A. Plastino, Phys. Lett. A , 140 (1994). 
文献３：  C. Tsallis, R. S. Mendes and A. R. Plastino, Physica A , 534 (1998). 
文献４：  W. D. Cornell, P. Cieplak, C. I. Bayly, I.R.Gould, K.M.Merz,Jr., D.M.Ferguson, D.C.Spellmeyer, T.Fox, J.W.Caldwell, and P.A.Kollman, J. Am. Chem. Soc. 117 (1995) 5179. 
文献５：  P. Kollman, R. Dixon, W. Cornell, T. Fox, C. Chipot, and A. Pohorille, in: W. F. van Gunsteren, P. K. Weiner, and A. J. Wilkinson (Eds.), Computer Simulation of Biomolecular Systems, Vol. 3, Kluwer, Netherlands, 1997, pp. 83-96. 
文献６：  U. H. E. Hansmann and Y. Okamoto, Phys. Rev. E 56 (1997) 2228. 
（付記１）
決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション装置において、
Boltzmann-Gibbs分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成する複数分布生成手段と、
数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分する数値積分手段と、
該数値積分の結果得られた軌跡に沿ってサンプリングするサンプリング手段と、
を備えることを特徴とするサンプリングシミュレーション装置。
【０１２２】
（付記２）
前記エネルギー分布がBoltzmann-Gibbs分布より広い分布は、Tsallis分布であることを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１２３】
（付記３）
複数の異なるパラメータについて生成される分布は、温度に対応するパラメータのみを異なる値に設定することによって得られることを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１２４】
（付記４）
前記合成された分布は、Boltzmann-Gibbs分布と、Boltzmann-Gibbs分布よりエネルギー分布が広い分布を合成して得られることを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１２５】
（付記５）
前記決定論的微分方程式は、ｘ、ｐ、ζを変数とし、M、Ｔ、ｎをパラメータとし、Uをポテンシャルエネルギー関数、Kを運動エネルギー関数とし、ρａＰを該広い分布とし、ρｚを任意の分布とし、Ｄｉをｘのｉ成分についての偏微分としたときに、
【０１２６】
【数３６】
2006221578.tif 000044
【０１２７】
ここで、
【０１２８】
【数３７】
2006221578.tif 000045
【０１２９】
で与えられることを特徴とする請求項１に記載のサンプリングシミュレーション装置。
（付記６）
前記Boltzmann-Gibbs分布よりエネルギー分布の広い分布の温度に対応するパラメータの基準値は、所望の温度の該Boltzmann-Gibbs分布でのエネルギー分布の値が同じ値を示す２つのエネルギー値において、該広い分布のエネルギー分布の２つの値が同じ値となるように設定されることを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１３０】
（付記７）
Boltzmann-Gibbs分布を計算することによってポテンシャルエネルギーの最小値を算出し、
Boltzmann-Gibbs分布から、よりエネルギー分布の広い分布の温度に対応するパラメータの基準値を算出し、
上記で求められた最小値と基準値を用いて、複数の温度に対応するパラメータについて、よりエネルギー分布の広い分布を計算し、
複数の温度に対応するパラメータに対する、よりエネルギー分布の広い分布を表示し、
表示された分布から、異なる温度に対応するパラメータに対応するエネルギー分布の重なり方を判断した結果により、規格化係数を計算するための分布を決定し、
規格化係数を計算する
ことを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１３１】
（付記８）
前記決定論的微分方程式を数値積分して、サンプリングを行い、
境界条件を処理し、
エネルギーの計算を行い、
サンプリング結果を出力する
ことを特徴とする付記１に記載のサンプリングシミュレーション装置。
【０１３２】
（付記９）
決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション方法において、
Boltzmann-Gibbs分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成し、
数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分し、
該数値積分の結果得られた軌跡に沿ってサンプリングする、
ことを特徴とするサンプリングシミュレーション方法。
【０１３３】
（付記１０）
決定論的微分方程式を数値積分した結果を用いてサンプリングを行い、物質の物理的、化学的性質を計算するサンプリングシミュレーション方法をコンピュータに実現させるプログラムを記録した記録媒体において、
Boltzmann-Gibbs分布より、カバーするエネルギー範囲が広い分布を異なる複数のパラメータについて複数生成し、
数値積分した結果得られた軌跡に沿ってサンプリングした結果が該複数の分布を合成した分布を再現する決定論的微分方程式を数値積分し、
該数値積分の結果得られた軌跡に沿ってサンプリングする、
ことを特徴とするサンプリングシミュレーション方法をコンピュータに実現させるプログラムを記録した記録媒体。
【図面の簡単な説明】
【０１３４】
【図１】複数のＴの値に対するTsallisエネルギー分布を示す図である。
【図２】３００ＫのBGエネルギー分布と、Ｔ＝１９．９Ｋの単一の Tsallisエネルギー分布Ｐを示す図である。
【図３】２つのTsallisポテンシャルエネルギー分布と、これらの平均の分布と、２つのTsallis分布を合成させた場合のポテンシャルエネルギー分布を示す図である。
【図４】２つのTsallis分布とそれを合成した分布、及び各温度でのＢＧ分布を示す図である。
【図５】本実施形態で得た、ポテンシャルエネルギー値の軌跡を示す図である。
【図６】（ｑ、ε、Ｔ１）のTsallis 分布、（ｑ、ε、Ｔ２）のTsallis 分布、３００K のBG分布、の3つの分布を本手法により合成した結果を示す図である。
【図７】本発明の実施形態に従ったサンプリングシミュレーション装置への入力データを説明する図である。
【図８】ＢＧ分布を生成するための入力データを示す図である。
【図９】パラメータ列生成部入力データを示す図である。
【図１０】本発明の実施形態に従った出力データを示す図である。
【図１１】本発明の実施形態のシミュレーション装置のブロック図である。
【図１２】入力部のブロック図である。
【図１３】Tsallis-parameter列生成部の構成及びフローを示す図である。
【図１４】計算部の構成及び処理フローである。
【図１５】逐次出力部の構成を示す図である。
【図１６】最終出力部の構成を示す図である。
【符号の説明】
【０１３５】
１０      入力部
１１      計算部
１２      出力部
１３      データ入力部
２０      ＢＧＤ入力部
２１      ＢＧＤ計算部
２２      Ｕｉｎｆ計算部
２３      Ｔｒｅｆ計算部
２４      Ｔａ計算部
２５      ＴＤ入力部
２６      ＴＤ計算部
２７      ＴＤ出力部
２８      ＢＧＤ入力部
２９      ＢＧＤ計算部
３０      ＢＧＤ出力部
３２      並び替え部
３３      分布重なり計算部
３４      表示部
３５      Ｙａ計算部
４０      ポテンシャル関数計算部
４１      ｚ−密度値計算部
４２      エネルギー計算部
４３      方程式右辺計算部
４４      Tsallis分布密度計算部
４５      ＢＧ分布密度値計算部
４６      最終出力部
４７      逐次出力部
【図１】
2006221578.tif 000046
【図２】
2006221578.tif 000047
【図３】
2006221578.tif 000048
【図４】
2006221578.tif 000049
【図６】
2006221578.tif 000050
【図７】
2006221578.tif 000051
【図８】
2006221578.tif 000052
【図９】
2006221578.tif 000053
【図１０】
2006221578.tif 000054
【図１１】
2006221578.tif 000055
【図１２】
2006221578.tif 000056
【図１３】
2006221578.tif 000057
【図１４】
2006221578.tif 000058
【図１５】
2006221578.tif 000059
【図１６】
2006221578.tif 000060
【図５】
2006221578.tif 000061