(57)【特許請求の範囲】
【請求項１】
コンピュータに、
入力部より、ユーザから物理系あるいは化学系の熱力学量を算出するためのシミュレーションに必要なパラメータ、及び、少なくともポテンシャル関数を含む微分方程式を設定するためのデータの入力を受け付け、入力された微分方程式を設定するためのデータに基づき、統計力学によって定義される統計分布をエネルギーの関数としてみたとき、エネルギーが大きくなるにつれての分布の変化が、指数関数より遅く０に漸近する一般化分布を再現する微分方程式を設定する処理と、
入力されたパラメータに基づいて、設定された微分方程式を積分し、得られた軌道に沿って、サンプリングを行い、該サンプリングに基づいて物理量の長時間平均を算出する処理と、
を実行させることを特徴とするシミュレーションプログラム。
【請求項２】
前記一般化分布は、Tsallis分布であることを特徴とする請求項１に記載のシミュレーションプログラム。
【請求項３】
前記微分方程式は、Tsallis分布を再現する連立常微分方程式であることを特徴とする請求項２に記載のシミュレーションプログラム。
【請求項４】
前記連立微分方程式は、ｘiを一般化座標、ｐiを一般化座標に共役な運動量、Ｕをポテンシャル関数、ｎを整数、Ｔを温度、ｋをボルツマン定数、β＝１／ｋＴ、Ｅをエネルギー関数、Ｄを微分演算子、ｑを１以上の定数、ζをTsallis分布を適切に作るための制御変数、ρzをユーザが与える任意関数とした場合に、
【数１】
0004733307.tif 000032
で与えられることを特徴とする請求項３に記載のシミュレーションプログラム。
【請求項５】
物理系あるいは化学系の熱力学量を算出するためのシミュレーション装置であって、
ユーザから、物理系あるいは化学系の熱力学量を算出するためのシミュレーションに必要なパラメータ、及び、少なくともポテンシャル関数を含む微分方程式を設定するためのデータの入力を受け付け、入力された微分方程式を設定するためのデータに基づき、統計力学によって定義される統計分布をエネルギーの関数としてみたとき、エネルギーが大きくなるにつれての分布の変化が、指数関数より遅く０に漸近する一般化分布を再現する微分方程式を設定する入力手段と、
入力されたパラメータに基づいて、設定された微分方程式を積分し、得られた軌道に沿って、サンプリングを行い、該サンプリングに基づいて物理量の長時間平均を算出する物理量算出手段と、
を備えることを特徴とするシミュレーション装置。
【請求項６】
シミュレーション装置に、
入力部により、ユーザから物理系あるいは化学系の熱力学量を算出するためのシミュレーションに必要なパラメータ、及び、少なくともポテンシャル関数を含む微分方程式を設定するためのデータの入力を受け付け、入力された微分方程式を設定するためのデータに基づき、統計力学によって定義される統計分布をエネルギーの関数としてみたとき、エネルギーが大きくなるにつれての分布の変化が、指数関数より遅く０に漸近する一般化分布を再現する微分方程式を設定する処理と、
入力されたパラメータに基づいて、設定された微分方程式を積分し、得られた軌道に沿って、サンプリングを行い、該サンプリングに基づいて物理量の長時間平均を算出する処理と、
を実行させることを特徴とするシミュレーション方法。
【請求項７】
コンピュータに、
入力部により、ユーザから物理系あるいは化学系の熱力学量を算出するためのシミュレーションに必要なパラメータ、及び、少なくともポテンシャル関数を含む微分方程式を設定するためのデータの入力を受け付け、入力された微分方程式を設定するためのデータに基づき、統計力学によって定義される統計分布をエネルギーの関数としてみたとき、エネルギーが大きくなるにつれての分布の変化が、指数関数より遅く０に漸近する一般化分布を再現する微分方程式を設定する処理と、
入力されたパラメータに基づいて、設定された微分方程式を積分し、得られた軌道に沿って、サンプリングを行い、該サンプリングに基づいて物理量の長時間平均を算出する処理と、
を実行させることを特徴とするシミュレーションプログラムを記録した、コンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【０００１】
【発明の属する技術分野】
本発明は、薬剤などの開発などにおいて有用なシミュレーション装置に関する。
【０００２】
【従来の技術】
Tsallis分布［Ｔ］は、非加法的に拡張した一般化されたカノニカル分布であり、長距離相互作用や長時間の記憶効果を持つ系などの特徴である非加法性を記述するのに適した統計力学の理論体系として、近年注目されている。以下の説明では、変数ｘ、ｐあるいはこれらの関数であるＥ（ｘ、ｐ）が連続な値をとる連続型のTsallis分布を対象とする。このような連続型Tsallis分布は、次の確率密度により定義される（なお、この定義は一義てきなものではなく、様々な変形がしられている）。
【０００３】
【数２】
0004733307.tif 000002
【０００４】
ここで、ｘ＝（ｘ1、・・・、ｘn）は（一般化）座標の組であり、ｐ＝（ｐ1、・・・、ｐn）は、（それに共役な）運動量の組を表す。また、Ｅ（ｘ、ｐ）＝Ｕ（ｘ）＋Ｋ（ｐ）は全エネルギーであり、Ｕはポテンシャルエネルギー、Ｋは運動エネルギー
【０００５】
【数３】
0004733307.tif 000003
【０００６】
である。（質量は全て１とした）。、ｃは規格化因子、β＝１／ｋＴ、Ｔは温度に相当するパラメータ（以下、温度）、ｋはボルツマン定数である。ｑは実数パラメータであり、Ｅの値の変域に応じてρTsallisが実数として定義されるような範囲で与えられる。ｑ→１の極限において、カノニカル分布を特徴付けるものとして従来から用いられてきたボルツマン−ギプス（ＢＧ）の密度
【０００７】
【数４】
0004733307.tif 000004
【０００８】
に等しくなる。これが、（１．１）で定義される分布が一般化されたカノニカル分布と呼ばれる理由の一つである。
物質の物理的・化学的な性質に関して、実験室で通常得られる情報は、一定の温度の下の熱浴内にある系の情報と言うことができる。以下の説明の適用対象として注目する物理系は、原子数で数個から数百万個の原子・分子の集合で表される系である。特に注目するのは、薬物の対象となる（生体系の）高分子と、それが働く環境を構成する分子とで表される系である。このような系におけるポテンシャルエネルギー関数の典型的な例は、次式のようなものである。
【０００９】
【数５】
0004733307.tif 000005
【００１０】
第１項が、各原子間の距離で決まるようなポテンシャルエネルギーを表し、第２項は、分子内原子の結合長、結合角などの平衡値を表現するようなポテンシャルエネルギーの総和である。
【００１１】
ＢＧ分布は、熱浴内の物理系の状態を表す分布としてよく研究されてきた。特に、ＢＧ分布に基づいた自由エネルギー等の熱力学諸量によって、上記（１．３）で記述される高分子とその環境を構成する系の動的性質を良く説明できることが、従来から多くの研究で明らかにされている。例えば、生体高分子であるタンパク質が取る安定な立体構造は、各タンパク質を規定する固有のアミノ酸配列をもつポリペプチド鎖が構築できる様々な立体構造の中から、ＢＧ分布に基づく自由エネルギーの最も低い状態が選択されていると考えられる。また、薬物と受容体タンパク質との親和力は、両者が結合した複合体として存在する系と、両者が離散した場合の系との自由エネルギーの差によって定量的に説明される。
【００１２】
しかしながら、上記（１．３）で表される系は、一般的にエネルギー面が複雑であり、そのＢＧ分布を正確に得ることは従来の技術では困難であった。更に、Tsallis分布（１．１）はべき則的であり、指数関数的なＢＧ分布（１．２）と比較して明らかなように、これらの分布をエネルギーＥの関数と見たときのエネルギーＥへの依存性は、ｑ＞１のTsallis分布の場合、エネルギーＥが大きくなるにつれてＢＧ分布よりもなだらかに減少しており、広い範囲のエネルギー値を有する状態の集合となっている。すなわち、Tsallis分布に基づいて相空間内でサンプリングを行い、（ｘ、ｐ）の組の集合を作った場合、ｑ＞１のTsallis分布においては、ＢＧ分布よりも遙かに広い相空間がサンプリングされていることになる。従って、数値計算によって、熱力学量を計算する場合に、このサンプリングデータをもとに演算を行うと、Tsallis分布を用いたときの方がより広い相空間をサンプリングするので、より良い結果が得られる。
実際、効果的サンプリングのためにTsallis分布の考え方が利用された例もある［ＡＳ］。効果的サンプリングがTsallis分布の実現によって可能となり、更にＢＧ分布をTsallis分布から正確に再構築することができれば、ゲノム配列からタンパク質の安定な立体構造を予測したり、新規薬物の受容体タンパク質への親和力を推定する、コンピュータ上の薬物設計への応用が期待される。
【００１３】
更に、Tsallis分布そのものが、従来のＢＧ分布と同様に、熱浴内の現実の物理系の状態を表す分布の候補と考えられ始めている［ＰＰ，ＡＲ１，ＡＲ２］。従って、Tsallis分布を計算機上で実現することによっても、熱浴内の物質の情報を計算機にて得ることが可能であり、更に、ＢＧ加法的統計力学ではうまく取り扱えない現実の物理現象に対し、説明能力を持ちうると言う報告もなされている［Ａ］。従って、Tsallis分布の計算機上での実現により、熱浴内の物理系の振る舞いに関する全く新しい知見を得ることが期待される。対象の系としては、様々なレベルのものが考えられる：
・電子・陽電子散乱により生成されるハドロン・ジェットの横運動量分布といった素粒子レベルのミクロなもの［ＢＣＭ］から、
・銀河クラスターの速度分布といった大きなスケールの系のもの［Ｌ  et  al  ］にまで、
それらの特徴がTsallis分布で記述できる例が報告されている。基本的には、エネルギー関数が定義されている系がその対象であり、特に、長距離相互作用や長時間の記憶効果を持つ系では、Tsallis分布の効果が期待される。また、式（１．３）の第１項は一般に長距離相互作用であり、これらの系においてTsallis分布を実現することにより、複雑な（生体系）高分子の動的挙動や構造形成に関する全く新しい知見を得て、新たな視点による薬物設計なども期待できる。
【００１４】
Tsallis分布を実現する従来手法として以下の種類のものが知られている：
（１）（ｘ、ｐ）両変数を用いたｑ＝１のTsallis分布（すなわちＢＧ分布）
（２）ｘ変数を用いたｑ≧１のTsallis分布
（３）（ｘ、ｐ）両変数を用いたｑ≦１のTsallis分布
各々を実現する方法の内、決定論的手法としては、以下が代表的なものである。
【００１５】
（方法１）Nos&eacute;とHooverはエネルギー関数が定義された物理系に対してあるＯＤＥ（Ordinary  Differential  Equation：常微分方程式）を提唱し、単に一自由度を追加するのみで温度ＴのＢＧ分布が実現可能であることを示した［Ｎ，Ｈ］。これにより熱浴内の物理系のシミュレーションを現実的なものにし、シミュレーションにより得られるカノニカル（ボルツマン−ギプス）集団に関する知識を飛躍的に増大させた。Nos&eacute;−Hoover（ＮＨ）の方法は、分子動力学（ＭＤ：Molecular  Dynamics）の分野で広く用いられている。この方程式はNos&eacute;−Hoover方程式と呼ばれ以下で表される。
【００１６】
【数６】
0004733307.tif 000006
【００１７】
ここで、Ｕはポテンシャル関数であり、Ｑは正のパラメータである。
（方法２）AndricioaeiとStraubはTsallis分布を実現する方法を提案し、化学系に適用し、それらの有効性を示した［ＡＳ］。効果的サンプリングを目的にしたこのアプローチでは、Tsallis分布としては座標変数ｘのみ考慮し、ニュートン方程式を基礎にした方法が提案された。
【００１８】
（方法３）PlastinoとAnteneodoの手法［ＰＡ］は、Tsallis分布を生成する方法であるが、方法２と異なり、（ｘ、ｐ）両変数を考慮したものである。この手法がサポートするのは、主に、０＜ｑ≦１の場合である。（ｑ＜１なる分布を用いた解析が有効であるような現象には、前記銀河クラスターの速度分布などが報告されている。）
一方、ＢＧ分布におけるサンプリングを行う従来手法としては、ＢＧ分布及びその拡張であるTsallis分布を実現して行う手法がある。この内、決定論的方法には、やはり上記方法１−３を代表的手法として挙げることができる。
【００１９】
【発明が解決しようとする課題】
前記（１）−（３）に対するTsallis分布を実現する従来方法としては方法１−３があるが、（ｘ、ｐ）両変数を用いたｑ≧１のTsallis分布を実現する方法は未だ成功していない。本発明はこの分布を実現する。
【００２０】
次に、ＢＧ分布におけるサンプリングを行う際の従来手法の課題点を述べる。（方法１）ＢＧ分布を生成する方法として良く用いられる方法であるが、パラメータＱの設定が簡単ではない［Ｎ２，ＰＨＶ］。Ｑ→∞で定摩擦付きの純粋なニュートン方程式に近づくので、これが大きすぎると温度Ｔへの温度制御の効果が落ちる。逆に、Ｑが小さすぎると、摩擦項の非線形性が大きくなりすぎて、数値積分アルゴリズムが破綻することが知られている。また、扱う系により、十分に点（ｘ、ｐ）のサンプリングができず、熱力学量を正しく計算できないことがある。
【００２１】
なお、ＢＧ分布を生成する方法としては、ＮＨの方法以外にもいろいろ知られている［Ｎ２］。しかし、（ｘ、ｐ）両変数を用いた決定論的な方程式を用いたＢＧ分布を実現する方法では、対象としている温度に対し高いポテンシャル障壁を超えることは、本質的に困難である。すなわち、障壁Ｕ（ｘw）が高くなるにつれ、確率Ｐ（ｘw、ｐ）∝ｅｘｐ［−Ｅ（ｘw、ｐ）／ｋＴ］≦ｅｘｐ［−Ｕ（ｘw）／ｋＴ］は指数的に小さくなるので、障壁点ｘwに到達するには、極めてまれに起こる現象を統計的に有意な量サンプルする必要がある。従って、そのような極めてまれに起こる現象に対応するサンプリングを相当数行うためには、本質的に長い計算時間が必要である。
【００２２】
（方法２）Tsallis分布をｘ変数のみを考慮して実現しようとする決定論的な方法である。これは、基本的に、系を表すポテンシャルＵに対し、ｘを変数とするTsallis分布用に変形されたポテンシャルＵ’を用いたニュートン方程式を用いる方法である。Ｕと比較してＵ’の障壁で働く力が小さくなることによる効果が期待でき、実際に、従来手法に比べサンプリング効果の向上が示された［ＡＳ，ＰＷ１，ＰＷ２，ＨＯ］。化学系に応用する際には、適当な温度コントロールを用いることが試みられた［ＰＷ１，ＰＷ２］。これらは、Tsallis分布をｘ変数のみで、すなわち、
【００２３】
【数７】
0004733307.tif 000007
【００２４】
のような形で、実現しようとする試みである。なお、前記Ｕ’は、
【００２５】
【数８】
0004733307.tif 000008
【００２６】
にて定義されるものである。
しかしながら、それらの試みに用いられている基本的方程式においては、エネルギー保存則が成立するので、Ｕ’（ｘ）＞Ｅ（ｘ（０）、ｐ（０））なるＵ’の障壁に対してはこれを超えることができず、局所的なエネルギー極小値状態へのトラップが起こり得る。従って、このような障壁の内部に閉じこめられてしまい、十分なサンプリングに至らない。これは、初期値（ｘ（０）、ｐ（０））に依存する問題であり、改良の余地がある。応用に際しては、温度コントロールの手法が明示されている報告と明示されていない報告とがある。しかしながら、いずれの報告にしても、（３．１）が真に実現されているか否かは、障壁のない単純な場合（１次元調和振動子）を除いては明らかにされていない。
【００２７】
（方法３）この手法が主にサポートする０＜ｑ≦１の場合、Ｅが有限の値Ｅ0≡１／（１−ｑ）β以上では、
【００２８】
【数９】
0004733307.tif 000009
【００２９】
と定義される。しかし、方法１にて述べたのと同様に、これはサンプリングに対しては不利となる（なお［ＰＡ］はサンプリングを目的にはしていない）。高いエネルギーの実現確率が小さく、もしくは零（Ｕ（ｘ）≧Ｅcのとき）になることによって、ポテンシャルＵの障壁を越すことが困難になるからである。これは、（
ｘ、ｐ）−空間の十分なサンプリングができず、正しい熱力学関数の計算ができないことを意味する。一方、ｑ＞１の場合は、数値解法におけるオーバフローの問題が指摘されている。その詳細は明らかにされていない。
【００３０】
本発明は、従来の方法に比べて、サンプリング能力の向上を図ることを目的とするが、これは上述した（ｘ、ｐ）両変数を用いたｑ≧１のTsallis分布を実現する方法の実現によって、自ずと達成される。この達成によりＢＧ分布での熱力学量の正確な値を求めることが可能になる。
【００３１】
これにより、計算機シミュレーションによる薬物設計の効率化・精密化を図ることができる。
本発明の課題は、薬物設計などにおけるシミュレーションに有益なシミュレーション装置を提供することである。
【００３２】
【課題を解決するための手段】
本発明のシミュレーション装置は、物理系あるいは化学系の熱力学量を精度良く算出するためのシミュレーション装置であって、シミュレーションに必要なパラメータを設定するパラメータ設定手段と、ボルツマン−ギプス分布に比べ、エネルギーの関数としてみたとき、よりゆっくり０に漸近する一般化分布を設定する一般化分布設定手段と、該一般化分布を再現する微分方程式を設定する微分方程式設定手段と、該パラメータに基づいて、該微分方程式を積分し、得られた軌道に沿って、サンプリングを行い、該サンプリングに基づいて物理量の長時間平均を算出する物理量算出手段とを備えることを特徴とする。
【００３３】
本発明によれば、ボルツマン−ギプス分布を用いて熱力学量を算出する場合、ボルツマン−ギプス分布ではなく、これよりもより、０への漸近の仕方が弱い一般化分布を用いてサンプリングを行うので、エネルギー関数の形が、極小値を有する場合でも、この極小値を囲む障壁を乗り越えてサンプリングを有効にすることができる。従って、サンプリングがエネルギー関数の極小値にトラップされてしまう可能性を低くし、より精度の良い熱力学量を算出を可能とする。
【００３４】
【発明の実施の形態】
以下、本発明の実施形態の説明においては、Tsallis分布を仮定して説明するが、本発明は、Tsallis分布に限定されるものではなく、ＢＧ分布をエネルギーＥの関数と見たとき、エネルギーＥが大きくなるにつれての分布の変化が、指数関数よりも遅く０に近づくような関数を分布として持つものであっても同じく適用が可能である。
【００３５】
まず、本発明の実施形態においては、（ｘ、ｐ）両変数を用いたｑ≧１のTsallis分布を、決定論的アルゴリズム、具体的には常微分方程式にて実現する。常微分方程式の解の相空間内の微小領域Δ（ｘ、ｐ）への滞在時間率の長時間極限が、（１、１）式のρTsallis（ｘ、ｐ）に比例するように、常微分方程式を定義する。ここで、Δ（ｘ、ｐ）は、任意の点（ｘ、ｐ）に対し、それを含む微小領域を意味する。本発明の実施形態においては、次のような常微分方程式を用いる。
【００３６】
【数１０】
0004733307.tif 000010
【００３７】
ここで、Ｄの記号は、微分演算子を表す。また、ζは、Tsallis分布を適切に生成するための制御用変数であり、ρz（ζ）は、ユーザが与える任意関数である。
【００３８】
後述するように、この常微分方程式を解くと、解はTsallis分布密度に比例する出現頻度にて時間変化することになり、Tsallis分布が実現される（（５．４）式）。また、Tsallis分布における物理量の期待値も、この常微分方程式の軌道にそって相空間内の点をサンプリングし、これを該常微分方程式に従った時間発展の方向にサンプリングによる長時間平均を行うことによって計算できる。なお、長時間平均は、通常初期値に依存しないことが示される。
【００３９】
このように、Tsallis分布を再現したサンプリングにより、物理量ＡのＢＧ分布での期待値を求めるには、
【００４０】
【数１１】
0004733307.tif 000011
【００４１】
のように２つの長時間平均の比を求めることにより、計算できる。ここでは、長時間平均は、上付線で表している（以下同じ）。ここで、ρTsallis及びρBGは各々、（１．１）式及び次式
【００４２】
【数１２】
0004733307.tif 000012
【００４３】
で与えられるものである。ここで、β’は、１／ｋＴ’であるが、（１．１）式のβと同じ値である必要はない。すなわち、ＢＧ分布の温度Ｔ’は、Tsallis分布の温度Ｔと異なって良い。なお、Tsallis分布を求めたのに、ＢＧ分布によって熱力学量（物理量の長時間平均）を求めるのは、ＢＧ分布を用いて計算した熱力学量から対象とする系の物理的、あるいは、化学的性質を推測する方法が確立されているからである。
【００４４】
方程式（４．１）が拡張可能であることを示すために、これがある方程式の特別な場合であることを説明する。該方程式は、Ｒ2n+1（２ｎ＋１次元空間（相空間に対応する））の領域Γにおいて任意に与えた密度関数ρに対し定義される：
【００４５】
【数１３】
0004733307.tif 000013
【００４６】
ここで、
【００４７】
【数１４】
0004733307.tif 000014
【００４８】
にて関数を定義した。また、外１  は各々、
【００４９】
【外１】
0004733307.tif 000015
【００５０】
関数Θのｘi成分、ｐi成分、ζ成分の点ω≡（ｘ、ｐ、ζ）での偏微分を表す。（４．２）式において、リユビル方程式、
【００５１】
【数１５】
0004733307.tif 000016
【００５２】
が満足されていることが確認できる（Ｘは（４．２）式右辺を表すものとする）。ρが適当な数学的条件を満たすと仮定すれば、リユビルの定理とバーコフのエルゴード定理［Ｍ］により、領域Ａ⊂Γに対し、軌道の滞在時間率の長時間極限が、測度ρｄωに関してほとんど全ての初期値に対し存在する（相空間内にある点を初期値とする軌道に沿った長時間平均の長時間極限が有限の値に収束するということが、ほとんどの相空間内の点について成立する）。もし、系がこの測度に関してエルゴード的ならば、
【００５３】
【数１６】
0004733307.tif 000017
【００５４】
となる。すなわち、状態ω∈Γの実現確率密度は、与えたρに比例する。（４．１）式は、（４．２）式と（４．３）式において、
【００５５】
【数１７】
0004733307.tif 000018
【００５６】
なる代入を行って得られるものである。これにより、（ｘ、ｐ）の実現確率密度を、定義したρTsallis（ｘ、ｐ）に比例させようとするものである。
なお、Tsallis分布を定義するものとして（１．１）式以外の密度が用いられることがある［ＴＭＰ］が、本発明でも対応可能（（４．５）式内のρTsallis（ｘ、ｐ）の定義を変更すれば良い）であり、（１．１）式の形に限定されない。ところで、通常は（１．１）式でなく、
【００５７】
【数１８】
0004733307.tif 000019
【００５８】
を用いて定義されることが多い。しかし、物理量の期待値を求める際には、［ＴＭＰ］で提案されているように、上記をｑ乗した（１．１）式の形で用いることも多い。本実施形態は物理量の期待値を求めることを目的としているので、便宜的に、（１．１）をその定義として採用したものである。また、Ｅ（ｘ、ｐ）も、更に一般的な形の関数として取り扱うことが可能である（（４．５）式のρTsallis（ｘ、ｐ）の定義内容が変更を受けることに相当する）。また、（４．４）式を満たすＸには自明な任意性（Ｘを定数倍できる）があるので、これを利用して、例えば、（４．３）式のΘが物理次元を持つ時に対応できる（Θを、これと同じ次元を持つ定数εによりスケーリングして、
【００５９】
【数１９】
0004733307.tif 000020
【００６０】
に対する式（４．２）を考え、その右辺をε倍すると、次元を持つΘで方程式を表せる）。これは同時に（４．１）式にも当てはまる。
数学的付加的条件の下で、方程式（４．１）に対し、バーコフのエルゴード定理［Ｍ］により、関数ｆの長時間極限が、測度ρｄωに関してほとんど全ての初期値に対し存在する。ここで、
【００６１】
【数２０】
0004733307.tif 000021
【００６２】
である。もし系がこの測度に関してエルゴード的ならば、
【００６３】
【数２１】
0004733307.tif 000022
【００６４】
特にまた、ｘ、ｐの関数としての物理量Ａの長時間平均値は、Tsallis分布におけるＡの期待値（重みρTsallisの空間平均）に等しくなる。すなわち、
【００６５】
【数２２】
0004733307.tif 000023
【００６６】
同様に、微小領域Δ（ｘ、ｐ）に対しては、次が成立する：
【００６７】
【数２３】
0004733307.tif 000024
【００６８】
左辺はΔ（ｘ、ｐ）への軌道訪問頻度の長時間極限である。従って、点（ｘ、ｐ）の実現確率密度はρTsallisに比例する。
なお、計算機シミュレーションでは、（５．３）の左辺の値を、十分大きい時間ステップＮに対し、
【００６９】
【数２４】
0004733307.tif 000025
【００７０】
などの近似値で評価する。
また、（５．３）式より、次式が成立する：
【００７１】
【数２５】
0004733307.tif 000026
【００７２】
これにより、ＢＧ分布における物理量Ａの期待値（右辺）は、この左辺の２つの長時間平均により計算できる。
図１〜図２０は、数値計算と理論値の比較による本実施形態の手法の有効性を示す図である。
【００７３】
最初に、（ｘ、ｐ）両変数を用いたｑ≧１のTsallis分布を実現した例を示す。次に、ＢＧ分布におけるサンプリング能力が従来法より向上した例を示す。そして、異なる温度でのＢＧ分布におけるサンプリング能力の検証例を示す。
【００７４】
Tsallis分布を実現した第一の例として、ポテンシャル関数Ｕを一次元の調和振動子
【００７５】
【数２６】
0004733307.tif 000027
【００７６】
としたものを図１〜図４に示す。温度Ｔは、１．０とした。簡単のため、全ての量は無次元で表す。ＯＤＥの初期値はｘ（０）＝０．１、ｐ（０）＝０．０、ζ（０）＝０．０とし、数値積分は４次のルンゲ−クッタ法で行い、いくつかのｑの値に対して計算した。シミュレーションに使ったパラメータ値を表１に示す。
【００７７】
【表１】
0004733307.tif 000028
【００７８】
Δｔは、時間刻み幅、Ｎは、シミュレーションのタイムステップを表す。ｑの値を３．０としたとき（case  １）のρTsallisをヒストグラムにて求めたシミュレーション値と理論値（（１．１）式）を図１及び図２に示した。これらは、良く一致している。同様に、ｑの値を２．０（case  ２）としたときの結果を図３及び図４に示す。シミュレーション値と理論値は良く一致していることが分かる。
【００７９】
Tsallis分布を実現した第二の例は、ポテンシャル関数Ｕを二重井戸振動子型（図３）としたものである（図５〜図９に示す）：
【００８０】
【数２７】
0004733307.tif 000029
【００８１】
ポテンシャル関数は、上記式で表され、その図を図５に示す。
温度Ｔ＝１．０に対し、障壁の高さは１０．０である。ｑの値を３．０（case  ４）、２．０（case  ５）としたのときの結果を、各々図６〜図９で示す。調和振動子の場合と同様に、シミュレーション値と理論値は良く一致している。
【００８２】
ＢＧ分布におけるサンプリング能力を検証するため、上記case  ４のシミュレーションにおいて、（５．６）式を用いて求めた結果を示し、従来法で得られる結果と比較した。まず、Ｔ＝１．０のＢＧ分布の密度ρBGをヒストグラムにて求めたシミュレーション値と理論値（（１．２）式）を図１０、及び図１１に示した。これらは良く一致している。次に、物理量としてエネルギーＥ（ｘ、ｐ）のＢＧ分布における期待値のシミュレーション結果（（５．５）式）を示す（図１２）。エネルギー期待値の理論値は、定積分
【００８３】
【数２８】
0004733307.tif 000030
【００８４】
の値を見積もると、１．０２４８となる。従って、シミュレーション結果は正しい値に収束しているのが分かる。比較として、方法１の従来法にて、同様に熱力学関数を計算した結果を示す。温度Ｔ、初期値、数値積分法、時間刻み幅はcase４と全て同じものである。方程式（２．１）のパラメータ１／２Ｑの値を、１０-2、１０-1、・・・、１０4、１０5、と２．０×１０5の９通り行った。全ての場合において、系の温度
【００８５】
【数２９】
0004733307.tif 000031
【００８６】
のコントロールは成功したが、解の座標値ｘは初期座標０．１を含むポテンシャルの左側の谷領域に局在していて、右側に移動が起こらなかった。シミュレーション結果の幾つかを図１３〜図１８に示す。ここで、図１３は、パラメータ１／２Ｑ＝１０-2、初期値ｘ＝０．１、図１４は、初期値は同じで、パラメータ１／２Ｑ＝１０、図１５は、初期値は同じで、パラメータ１／２Ｑ＝１０3、図１６は、初期値が同じで、パラメータ１／２Ｑ＝１０5、図１７は、初期値は同じで、パラメータ１／２Ｑ＝２．０×１０5、図１８は、初期値がｘ＝０．９で、パラメータ１／２Ｑ＝２．０×１０5としている。
従来法である方法１においては、サンプリングが二重井戸の一方の井戸へ局在した結果、ＢＧ分布は正しく生成されていない。従って、また、熱力学関数値も正しく計算されないことになる。パラメータ１／２Ｑ＝２．０×１０5に対するエネルギー期待値のシミュレーション結果を図１７に示す。サンプリングが不十分であるため、正しい値に収束していない。前述のように、パラメータ１／Ｑの値を更に小さくすると、ニュートン方程式に近づくのでＢＧ分布の生成能力が落ち、これらの結果は改善されない。逆に、パラメータ１／Ｑの値を大きくしすぎると、数値積分アルゴリズムが破綻することが知られている。実際、今回の計算でも、１／２Ｑ＝１０6とすると座標ｘが数値オーバフローを起こし、数値積分アルゴリズムが破綻した。また、初期座標値を右側の谷領域ｘ＝０．９としても、同様に局在化し、サンプリングが不十分であった（１／２Ｑ＝２．０×１０5の場合を図１８に示す）。
【００８７】
異なる温度でのＢＧ分布におけるサンプリング能力を検証するため、上記case４と同条件を用い、ただし２倍の温度２．０のＢＧ分布に対し本実施形態に従って、計算を行った。密度ρBGをヒストグラムにて求めたシミュレーション値と理論値を図１９と図２０に示した。これらの一致を見ることができる。
【００８８】
図２１は、本発明の実施形態に従ったシミュレーション装置のブロック構成図である。
シミュレーション装置１は、入力部１０、計算部１１、及び出力部１２によって構成される。入力部１０には、演算に必要な入力データが入力される。入力部１０に入力された入力データは、計算部１１に送られる。計算部１１では、入力データに含まれているパラメータ値などを使って、実施形態の説明で示したようなTsallis分布を再現する（４．１）式、あるいは、ＢＧ分布よりエネルギーについて、よりゆっくり（エネルギーＥの関数として指数関数よりゆっくりであって、例えば、Ｅの有限べき乗に従って漸近的振る舞いが支配されるような関数の漸近的振る舞いを「ゆっくり」としている）０に漸近する分布を再現する微分方程式をルンゲークッタ法などの逐次数値解法により積分する（解を求める）。そして、微分方程式を積分した結果えられた軌道上の点からサンプリングを行い、サンプリングされた点を用いて、物理量の値及び期待値を計算する。
【００８９】
計算された物理量の値あるいは期待値は、出力データとして出力部１２に渡され、不図示の表示装置などによって、ユーザに提示される。
図２２は、入力データの例を示した図である。
【００９０】
同図に示されるように、入力データの例としては、自由度数ｎ、ポテンシャル関数Ｎｏ（この場合、計算部に複数の異なるポテンシャル関数を登録しておき、これらの中から選択するようにすることを念頭に置いている。なお、関数の形を直接入力し、シミュレーション装置に解釈させる様にしても良い。）、シミュレーション条件、境界条件、Tsallis分布パラメータ、ρz関数パラメータ、ＢＧ分布温度リスト、及び度数計算パラメータなどがある。
【００９１】
シミュレーション条件は、シミュレーションのステップ数、時間刻み幅、出力指定パラメータなどがある。また、Tsallis分布パラメータとしては、ｑや温度Ｔあるいはβがある。ρz関数パラメータとしては、ポテンシャル関数Ｎｏと同様の意味での関数Ｎｏ、及び設定係数などがある。ＢＧ分布温度リストとしては、計算するＢＧ分布の温度の数、及び、計算するＢＧ分布の温度の値がある。度数計算パラメータとしては、分布のサンプリングをヒストグラムとして形成する場合の離散幅、変域を制御する変域パラメータなどである。
【００９２】
図２３は、出力データの例を示す図である。
出力データとしては、微分方程式を積分して解いた状態変数ｘ、ｐの値の時間変化や、その他の変数値の時間変化（ζ、ポテンシャルエネルギーなど）、物理量の期待値（時間平均）の時間変化、各温度におけるＢＧ分布における物理量の期待値（時間平均）の時間変化（及び最終値）（エネルギー、比熱等）、ポテンシャルエネルギー値の度数（実現確率）、運動エネルギー値の度数（実現確率）、Tsallis分布の密度（密度、周辺分布）、ＢＧ分布の密度（密度、周辺分布）等が考えられる。上記データの内、熱力学量となるのは、物理量の期待値、各温度におけるＢＧ分布における物理量の期待値である。その他のデータは、シミュレーションの精度を評価するなどのために用いることができる。
【００９３】
図２４、及び図２５は、本発明の実施形態の全体の処理の流れを示すフローチャートである。
まず、ステップＳ１において、データ入力を行う。次に、ステップＳ２において、定数定義を行う。そして、ステップＳ３において、計算結果を出力するファイルへの出力準備を行う。ステップＳ４では、変数の初期設定を定義する。この場合、どのように設定するかは、ステップＳ５において自動生成（ステップＳ６）と外部ファイルからの読み込み（ステップＳ７）との間で選択する。変数の初期設定が完了すると、ステップＳ８において、度数計算の準備を行う。そして、ステップＳ９で数値積分を開始する。
【００９４】
ステップＳ１０においては、ｔをｔ＋Δｔとして、ステップＳ１１において、数値積分を実行する。数値計算を実行するに当たっては、ｚ−密度値計算部（ステップＳ１６）やポテンシャル関数計算部（ステップＳ１７）などの結果を解くべき方程式の（本発明の実施形態の場合、（４．１）式の）右辺を計算し（ステップＳ１８）、この結果を使って、数値積分を行う。この場合、積分手法の選択や追加を可能とすることが可能である。また、ポテンシャル関数計算部の演算結果は、ステップＳ１５のエネルギー計算部に渡され、エネルギー計算に使用される。ステップＳ１２において、数値積分するための境界条件を処理する。ステップＳ１３で、エネルギー計算を行う。ステップＳ１４において、分布の密度計算を行う。このとき、分布としては、Tsallis分布やＢＧ分布がある。密度計算は、それぞれの分布の計算部からのデータを使用する。
【００９５】
ステップＳ１９においては、ｚ−密度値計算部からのデータに基づいて、数値積分チェックを行う。そして、ステップＳ２０において、数値エラーチェックを行い、エラーがある場合には、数値積分を終了する。エラーがない場合には、ステップＳ２１において、度数計算を行い、ポテンシャル、運動量エネルギーの度数を求めたり、Tsallis分布の密度を求める。ステップＳ２２においては、物理量の期待値を計算する。物理量の例としては、運動エネルギー、全エネルギー、ＢＧ分布におけるエネルギー、ＢＧ分布における比熱などがある。
【００９６】
そして、ステップＳ２３において、上記計算結果を出力データとして出力し、ステップＳ２４において、ユーザが予め定めた計算終了条件を満たしているか否かを判断する。満たしていない場合には、ステップＳ１０に戻って計算を繰り返し、満たしている場合には、数値積分を終了する。そして、最後に、数値積分終了後、出力データの確認として、Tsallis分布密度（周辺分布）やＢＧ分布密度（周辺分布）の理論値を可能な場合に計算し、出力データと比較して出力データの精度を吟味し、適切な分布が生成されたかを確認する。これにより、演算された物理量が正しい値を与えるか否かを評価する。
【００９７】
図２６は、入力部の処理の流れを示すフローチャートである。
まず、ステップＳ３０において、データの入力を行う。ステップＳ３１では、数値積分のための定数の定義を行い、ステップＳ３２では、ファイル出力準備を行う。そして、ステップＳ３３において、変数の初期値を定義する。このとき、変数の初期値を自動生成するか外部ファイルから読み込むか選択するようにする。あるいは、ユーザが手で入力しても良い。そして、ステップＳ３４において、密度・度数計算の準備を行い、計算部へと進む。
【００９８】
図２７は、計算部の処理の流れを示すフローチャートである。
ステップＳ４０において、数値積分を開始する。ステップＳ４１では、時間パラメータｔをｔ＋ΔｔとΔｔだけ進める。ステップＳ４２では、数値積分を実行する。このとき、ｚ−密度値計算部において、関数の選択、追加を行って適切な関数を使って計算を行った結果を、ポテンシャル関数計算部の計算結果と共に、当該方程式の右辺に代入し、計算する。その結果をステップＳ４２において、数値積分計算に使用する。数値積分計算においては、計算手法を選択、追加できるようにしておくのが便利である。ステップＳ４３においては、境界条件を処理し、ステップＳ４４においては、ポテンシャル関数計算部からのデータをつかったエネルギー計算がエネルギー計算部において行われものをつかって、エネルギー値の計算を行う。ステップＳ４５では、密度比を計算し、ステップＳ４６において数値積分をチェックし、ステップＳ４７において、数値エラーがあるか否かを判断する。エラーがあった場合には、数値積分を終了する。エラーがなかった場合には、ステップＳ４８において、度数計算を行い、ステップＳ４９において、物理量期待値計算を行う。そして、ステップＳ５０において、データを出力し、ステップＳ５１において、計算終了条件を判断し、終了条件が満たされていない場合には、ステップＳ４１にもどって計算を繰り返すが、終了条件が満たされた場合には、数値積分を終了する。
【００９９】
図２８は、出力部の処理の流れを示すフローチャートである。
ステップＳ６０においては、計算部で計算された出力データの確認を行う。すなわち、可能な場合には、Tsallis分布の密度やＢＧ分布の密度の周辺分布の理論値を計算し、ステップＳ６１において出力データを可視化してユーザが正しいシミュレーションが行われたか否かを判断し、出力データとする。
【０１００】
図２９は、ポテンシャル関数値計算部の演算内容を説明する図である。
当該微分方程式の積分結果としてｘ（ｔ）が更新された場合、以下の計算を行う。すなわち、ポテンシャル関数値、ポテンシャル関数の偏導関数（Force  Function）値（ポテンシャルの偏導関数は系に働く力を示すことが知られている）、エネルギー原点の調整、及び必要ならば、解を拘束させるためのポテンシャル関数の追加を行う。ここで、ポテンシャル関数としてどのようなものを使うかは、複数の候補の中から選択できるようにし、追加も可能としておくと便利である。
【０１０１】
図３０は、本発明の実施形態のシミュレーション装置をプログラムで実現する場合のコンピュータのハードウェア環境を示す図である。
ＣＰＵ２１は、バス２０によって接続されたＲＯＭ２２、ＲＡＭ２３に格納された当該プログラムを実行することによって、本発明の実施形態の処理を実行する。ＲＡＭ２３にコピーされる当該プログラムは、ハードディスクなどの記憶装置２７に格納されるか、フロッピーディスク、ＣＤ−ＲＯＭ、ＤＶＤなどの可搬記録媒体２９に格納される。
【０１０２】
可搬記録媒体２９に格納された当該プログラムは、読み取り装置２８によって読み取られ、バス２０を介して、ＲＡＭ２３にコピーされるか、記憶装置２７に格納されてコンピュータにインストールされる。
【０１０３】
入出力装置３０は、ディスプレイ、キーボード、マウス、テンプレート、などからなり、ユーザからの入力をＣＰＵ２１に伝えたり、ＣＰＵ２１の演算結果をユーザに提示するために使用される。
【０１０４】
通信インターフェース２４は、ネットワーク２５を介して情報提供者２６と接続し、当該プログラムをダウンロードして、ＣＰＵ２１が実行可能とすることができる。あるいは、情報提供者２６が有する当該プログラムをネットワーク環境下で実行するようにしても良い。
＜参照文献＞
・[A]阿部純義、数理科学、No.439、No.440、No.441；No.442（2000）
・[AR1] S.Abe  and  A.K.Rajagopal,  Nonuniqueness  of  canonical  ensemble  theory  arising  from  microcanonical  basis,  Phys.  Lett.  A272,341(2000)
・[AR2] S.Abe  and  A,K,Rajagopal,  Microcanonical  foundation  for  systems  with  power-law  distributions,  J.Phys.  A33,8733(2000)
・[AS]I.Andricioaei  and  J.E.Straub,  On  Monte  Carlo  and  molecular  dynamics  method  inspired  by  Tsallis  statistics:  Methodology,  optimization  and  applications  to  atomic  clusters,  J.Chem.Phys.  107,9117(1997)
・[BCM]I.Bediaga,  E.M.F.Curado  and  J.Miranda,  A  nonextensive  thermodynamical  equilibrium  approach  in  e+e_!  hadrons,  Physica  A  286,156(2000)
・[H]W.G.Hoover, Phys. Rev. A 31 1695 (1985)
・[HO] U.H.E.Hansmann and Y. Okamoto, Tackling the protein folding problem by ageneralized-ensemble approach with Tsallis statistics, in NonextensiveStatistical Mechanics and Thermodynamics, eds. S.R.A. Salinas and C. Tsallis,Braz. J.Phys. 29, 187(1999)
・[L et al] A.Lavagno, G. Kaniadakis, M. Rego-Monteiro, P. Quarati and C. Tsallis,Nonextensive thermostatistical approach of the peculiar velocity function ofgalaxy clusters, Astrophysical Letters and Communications 35, 449(1989)
・[M] R. Mane,"Ergodic Theory and Differentiable Dynamics," p.35, Springer-Verlag(1987)
・[N] S. Nos&eacute;, J. Chem.Phys. 81, 511 (1984)
・[N2] S. Nos&eacute;, Prog. Theor.Phys. Suppl. 103, 1 (1991)
・[PA] A.R.Plastino and C. Anteneodo, A dynamical thermostatting approach to non-extensivecanonical ensembles, Ann. Phys. 255, 250 (1997)
・[PHV] H.A.Posh, W.G. Hoover, and F.j.Vesely, Phys. Rev. A 33, 4253 (1986)
・[PP] A.R.Plastino and A. Plastino, From Gibbs microcanonical ensemble to Tsallisgeneralized canonical distribution, Phys. Lett. A 193, 140 (1994)
・[PW1] Y. Pakand S. Wang, Folding of a 16_residue helical peptide using molecular dynamicssimulation with Tsallis effective potential, J. Chem. Phys. 111,4356 (1999)
・[PW2] Y. Pakand S. Wang, J. Phys. Chem. B 2000, 104, 354 (2000)
・[T]C.Tsallis, Possble generalization of Boltzmann-Gibbs statistics, J. Stat. Phys.52, 479 (1988)
・[TMP] C.Tsallis, R.S.Mendes and A.R. Plastino, The  role of constraints within generalized nonextensivestatistics, Physica A 261,534 (1998)
【０１０５】
【発明の効果】
本発明によれば、以下のような効果を有している。
・（ｘ、ｐ）両変数を用いたｑ≧１のTsallis分布を生成する事が可能である。これは、従来手法では成功していない。
・任意の温度におけるＢＧ分布の生成が可能である。
・ＢＧ分布の生成により、熱力学量の計算が可能である。
・サンプリング能力が従来の方法よりも向上したことにより、従来法ではサンプリングが不十分で正しい分布や熱力学量値が計算できない場合でも、正しい値を計算することができる。
【図面の簡単な説明】
【図１】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１）である。
【図２】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その２）である。
【図３】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その３）である。
【図４】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その４）である。
【図５】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その５）である。
【図６】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その６）である。
【図７】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その７）である。
【図８】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その８）である。
【図９】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その９）である。
【図１０】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１０）である。
【図１１】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１１）である。
【図１２】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１２）である。
【図１３】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１３）である。
【図１４】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１４）である。
【図１５】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１５）である。
【図１６】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１６）である。
【図１７】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１７）である。
【図１８】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１８）である。
【図１９】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その１９）である。
【図２０】数値計算と理論値の比較による本実施形態の手法の有効性を示す図（その２０）である。
【図２１】本発明の実施形態に従ったシミュレーション装置のブロック構成図である。
【図２２】入力データの例を示した図である。
【図２３】出力データの例を示す図である。
【図２４】本発明の実施形態の全体の処理の流れを示すフローチャート（その１）である。
【図２５】本発明の実施形態の全体の処理の流れを示すフローチャート（その２）である。
【図２６】入力部の処理の流れを示すフローチャートである。
【図２７】計算部の処理の流れを示すフローチャートである。
【図２８】出力部の処理の流れを示すフローチャートである。
【図２９】ポテンシャル関数値計算部の演算内容を説明する図である。
【図３０】本発明の実施形態のシミュレーション装置をプログラムで実現する場合のコンピュータのハードウェア環境を示す図である。
【符号の説明】
１０      入力部
１１      計算部
１２      出力部
【図１】
0004733307.tif 000033
【図２】
0004733307.tif 000034
【図３】
0004733307.tif 000035
【図４】
0004733307.tif 000036
【図５】
0004733307.tif 000037
【図６】
0004733307.tif 000038
【図７】
0004733307.tif 000039
【図８】
0004733307.tif 000040
【図９】
0004733307.tif 000041
【図１０】
0004733307.tif 000042
【図１１】
0004733307.tif 000043
【図１２】
0004733307.tif 000044
【図１３】
0004733307.tif 000045
【図１４】
0004733307.tif 000046
【図１５】
0004733307.tif 000047
【図１６】
0004733307.tif 000048
【図１７】
0004733307.tif 000049
【図１８】
0004733307.tif 000050
【図１９】
0004733307.tif 000051
【図２０】
0004733307.tif 000052
【図２１】
0004733307.tif 000053
【図２２】
0004733307.tif 000054
【図２３】
0004733307.tif 000055
【図２４】
0004733307.tif 000056
【図２５】
0004733307.tif 000057
【図２６】
0004733307.tif 000058
【図２７】
0004733307.tif 000059
【図２８】
0004733307.tif 000060
【図２９】
0004733307.tif 000061
【図３０】
0004733307.tif 000062