(57)【要約】【課題】ＰＩＭＤ法を用いて、時間によって変化する物理量の解析を行うことができるシミュレーション方法を提供する。【解決手段】本シミュレーション方法は、複数の粒子からなる系の粒子の各々に対応して複数のビーズを定義し、複数のビーズの挙動を解析する経路積分分子動力学法を用いたものである。１つの粒子に対応する複数のビーズの質量が同一であり、１つの粒子に対応する複数のビーズの質量の合計が、対応する粒子の質量と等しくなるように複数のビーズの質量を設定する。複数の粒子に対応する複数のビーズの挙動を、分子動力学法を用いてシミュレーションする。【選択図】図１
2019046313.tif 000002
【特許請求の範囲】
【請求項１】
複数の粒子からなる系の前記粒子の各々に対応して複数のビーズを定義し、複数の前記ビーズの挙動を解析する経路積分分子動力学法を用いたシミュレーション方法であって、
１つの前記粒子に対応する複数の前記ビーズの質量が同一であり、１つの前記粒子に対応する複数の前記ビーズの質量の合計が、対応する前記粒子の質量と等しくなるように複数の前記ビーズの質量を設定し、
複数の前記粒子に対応する複数の前記ビーズの挙動を、分子動力学法を用いてシミュレーションするシミュレーション方法。
【請求項２】
分子動力学法を用いてシミュレーションする際に、１つの前記粒子に対応する複数の前記ビーズの間に、前記粒子の間に作用するポテンシャルに基づいて得られるポテンシャルと、前記ビーズの間に引力を生じさせるポテンシャルとを重畳したポテンシャルを導入する請求項１に記載のシミュレーション方法。
【請求項３】
入出力装置と、処理装置とを有し、
前記処理装置は、
前記入出力装置から入力されたシミュレーション条件に基づき、複数の粒子からなる系の前記粒子の各々に対して複数のビーズを定義し、
１つの前記粒子に対応する複数の前記ビーズの質量が同一であり、１つの前記粒子に対応する複数の前記ビーズの質量の合計が、対応する前記粒子の質量と等しくなるように複数の前記ビーズの質量を設定し、
複数の前記粒子に対応する複数の前記ビーズの挙動を、分子動力学法を用いてシミュレーションし、
シミュレーション結果を前記入出力装置に出力するシミュレーション装置。
【請求項４】
前記処理装置は、前記粒子の各々に対応する複数の前記ビーズの重心位置及び平均速度が、対応する前記粒子の位置及び速度の初期条件に一致するように、複数の前記ビーズを定義する請求項３に記載のシミュレーション装置。
【請求項５】
前記処理装置は、分子動力学法を用いてシミュレーションする際に、１つの前記粒子に対応する複数の前記ビーズの間に、前記粒子の間に作用するポテンシャルに基づいて得られるポテンシャルと、前記ビーズの間に引力を生じさせるポテンシャルとを重畳したポテンシャルを導入する請求項３または４に記載のシミュレーション装置。
【請求項６】
前記入出力装置に出力される前記シミュレーション結果は、複数の前記粒子の各々に対応する複数の前記ビーズの重心位置の時間変化を含む請求項３乃至５のいずれか１項に記載のシミュレーション装置。
【発明の詳細な説明】
【技術分野】
【０００１】
本発明は、分子動力学法を利用したシミュレーション方法及びシミュレーション装置に関する。
【背景技術】
【０００２】
量子力学に従う系を古典力学で表現して解析する経路積分分子動力学法（ＰＩＭＤ法）が知られている。ＰＩＭＤ法においては、量子的な１個の粒子が古典系ではＰ個のビーズで表現される。これらのビーズは、レプリカと呼ばれる場合もある。Ｐ個のビーズは、例えばバネによって環状に接続されている。複数のビーズの間には、バネ力の他に、元の粒子に作用するポテンシャルと同形のポテンシャルによる力が作用する。このような条件の下で、分子動力学法を用いて複数のビーズの挙動を解析することにより、物理量の統計平均を求める。このＰＩＭＤ法によって求められた古典的統計平均が、元の量子系における統計平均に対応する（非特許文献１）。
【先行技術文献】
【非特許文献】
【０００３】
【非特許文献１】David Chandler and Peter G. Wolynes, “Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids”, The Jounal of Chemical Physics 74.7, 4078-4095 (1981)
【発明の概要】
【発明が解決しようとする課題】
【０００４】
元の粒子の各々を表す複数のビーズに付される通し番号は量子系の虚時間発展と関係しており、ＰＩＭＤ法における時間発展の物理的な意味は明らかではなかった。このため、従来のＰＩＭＤ法では、時間に依存しない物理量（静的または定常的な量、統計平均）を計算することは可能であったが、動的な物理量を計算することはできなかった。
【０００５】
本発明の目的は、ＰＩＭＤ法を用いて、時間によって変化する物理量の解析を行うことができるシミュレーション方法及びシミュレーション装置を提供することである。
【課題を解決するための手段】
【０００６】
本発明の一観点によると、
複数の粒子からなる系の前記粒子の各々に対応して複数のビーズを定義し、複数の前記ビーズの挙動を解析する経路積分分子動力学法を用いたシミュレーション方法であって、
１つの前記粒子に対応する複数の前記ビーズの質量が同一であり、１つの前記粒子に対応する複数の前記ビーズの質量の合計が、対応する前記粒子の質量と等しくなるように複数の前記ビーズの質量を設定し、
複数の前記粒子に対応する複数の前記ビーズの挙動を、分子動力学法を用いてシミュレーションするシミュレーション方法が提供される。
【０００７】
本発明の他の観点によると、
入出力装置と、処理装置とを有し、
前記処理装置は、
前記入出力装置から入力されたシミュレーション条件に基づき、複数の粒子からなる系の前記粒子の各々に対して複数のビーズを定義し、
１つの前記粒子に対応する複数の前記ビーズの質量が同一であり、１つの前記粒子に対応する複数の前記ビーズの質量の合計が、対応する前記粒子の質量と等しくなるように複数の前記ビーズの質量を設定し、
複数の前記粒子に対応する複数の前記ビーズの挙動を、分子動力学法を用いてシミュレーションし、
シミュレーション結果を前記入出力装置に出力するシミュレーション装置が提供される。
【発明の効果】
【０００８】
１つの粒子に対応する複数のビーズの質量の合計が、対応する粒子の質量と等しくなるように複数のビーズの質量を設定することにより、ＰＩＭＤ法を用いて、時間によって変化する物理量の解析を行うことができる。
【図面の簡単な説明】
【０００９】
【図１】図１は、本実施例によるシミュレーション方法の手順のフローチャートである。
【図２】図２は、複数の粒子の各々に対して複数（Ｐ個）のビーズを定義することを示す図である。
【図３】図３Ａ、図３Ｂ、及び図３Ｃは、それぞれ第１のシミュレーション結果、第２のシミュレーション結果、及び第３のシミュレーション結果を表すグラフである。
【図４】図４は、実施例によるシミュレーション装置のブロック図である。
【発明を実施するための形態】
【００１０】
まず、本願発明の実施例について説明する前に、量子力学で取り扱われる量子系の分配関数ＺＱＭと、古典力学で取り扱われる古典的分配関数ＺＰとの関係について説明する。
【００１１】
量子系の分配関数ＺＱＭは、以下の式で表される。
【数１】
2019046313.tif 000003
ここで、Ｈハットはハミルトニアンであり、Ｔは温度であり、ｋＢはボルツマン定数である。記号Ｔｒは、行列の対角和を意味する。
【００１２】
量子的な１個の粒子は、古典系ではＰ個のビーズ（レプリカ）によって表現される。ｓ番目のビーズの座標をｒ（ｓ）（ｓ＝１、２、・・・Ｐ）と表記する。便宜上、ｒ（Ｐ＋１）＝ｒ（１）とする。ポテンシャルＶ（ｒ）中の１粒子系では、古典的分配関数ＺＰは以下の式で表される。
【数２】
2019046313.tif 000004
ここで、Ｐは整数、ｍは量子的な粒子の質量を表す。１つの粒子に対応するＰ個のビーズは、有効ポテンシャルＶｅｆｆを表す式の右辺の第二項に示すように、バネによって環状につながっている。
【００１３】
式（１）に示した量子系の分配関数ＺＱＭは、式（２）に示した古典的分配関数ＺＰにおいてパラメータＰが無限大の極限を取ったものと一致する。すなわち、以下の関係が成立する。
【数３】
2019046313.tif 000005
ＰＩＭＤ法では、古典的分配関数ＺＰのパラメータＰを有限の値に設定して古典的運動方程式を解くことにより、物理量の統計平均を求める。古典的手法で求められた統計平均は、元の量子系における統計平均に対応する。ところが、ＰＩＭＤ法における時間発展の物理的意味は明らかではなかった。
【００１４】
以下に説明する本実施例では、ＰＩＭＤ法における時間発展が、実時間の時間発展に対応すると解釈される。
【００１５】
上述の式（２）の右辺には、系のポテンシャルのみが明示的に表れている。ｓ番目のビーズの仮想的運動量ｐ（ｓ）と各ビーズの質量Ｍを導入して、式（２）は以下のように書き換えることができる。
【数４】
2019046313.tif 000006
ここで、Ｈｅｆｆは、以下の式（５）で表される。
【数５】
2019046313.tif 000007
【００１６】
式（５）で表されたＨｅｆｆは、元の量子系に対応した古典系のハミルトニアンとみなすことができる。式（５）の右辺の第１項は、ビーズの運動エネルギを表し、第２項は元の粒子間に作用するポテンシャルエネルギをビーズ間に作用するように置き換えたポテンシャルエネルギを表し、第３項は、１つの粒子に対応するＰ個のビーズの間に作用するバネ力（引力）によるポテンシャルエネルギを表している。このように、１つの粒子に対応するＰ個のビーズの間に、元の粒子の間に作用するポテンシャルＶ（ｒ）に基づいて得られるポテンシャル（式（５）の右辺第二項）と、ビーズの間に引力を生じさせるポテンシャル（式（５）の右辺第三項）とを重畳した相互作用ポテンシャルが作用する。以下、式（５）のＨｅｆｆを古典系のハミルトニアンと解釈して、量子系と古典系とを比較する。
【００１７】
時間発展に対応するパラメータをτと表記すると、ＰＩＭＤ法で使用する運動方程式は以下のように表すことができる。
【数６】
2019046313.tif 000008
【００１８】
量子力学における時間に依存せず粒子の位置のみに依存する演算子Ａハットに対応して、ＰＩＭＤ法で計算できる以下の量〈Ａ〉ＰＩＭＤを、演算子Ａハットの期待値として定義する。
【数７】
2019046313.tif 000009
【００１９】
演算子Ａハットを粒子の位置と定義すると、式（６）及び式（７）から、以下の式が導出される。
【数８】
2019046313.tif 000010
式（８）は、元の粒子に対応するＰ個のビーズの重心が従う運動方程式である。
【００２０】
量子力学において、波動関数ψの下での時間に依存しない演算子Ａハットの期待値は、以下の式で定義される。
【数９】
2019046313.tif 000011
【００２１】
粒子の位置の期待値の従う方程式は、式（９）とシュレディンガー方程式とから、以下の式で表される（エーレンフェストの定理）。
【数１０】
2019046313.tif 000012
【００２２】
Ｐ個のビーズの重心が従う運動方程式（８）と、量子系の粒子の位置の期待値が従う方程式（１０）とを比較すると、Ｍ＝ｍ／Ｐのとき、両者は同一の形をしていることがわかる。従って、量子系における時刻ｔ＝０での各粒子の初期条件を与え、対応してＰＩＭＤ法における時刻τ＝０でのビーズの初期条件が与えられれば、量子系での粒子の位置の期待値と、ＰＩＭＤ法でのビーズの重心の位置とはＰが無限大の極限において同一の値になる。すなわち、ＰＩＭＤ法における時間τが、量子系における実時間ｔと同一の物理的な意味を持つと解釈することができる。
【００２３】
このように、ＰＩＭＤ法におけるビーズの質量Ｍを、Ｍ＝ｍ／Ｐを満足するように設定してビーズの挙動をシミュレーションすることにより、量子系の粒子の動力学的性質を調べることができる。等式Ｍ＝ｍ／Ｐは、１つの粒子に対応する複数のビーズの質量が同一であり、１つの粒子に対応する複数のビーズの質量の合計が、対応する粒子の質量と等しくなるように複数のビーズの各々の質量を設定することを意味する。
【００２４】
図１は、本実施例によるシミュレーション方法の手順を示すフローチャートである。以下の各処理は、シミュレーション装置により実行される。まず、複数の粒子からなるシミュレーション対象の系の初期条件を設定する（ステップＳ１）。初期条件には、粒子の位置及び速度が含まれる。さらに、粒子間に作用する相互作用ポテンシャル、境界条件等を設定する。初期条件、相互作用ポテンシャル、境界条件等を設定する基礎となる情報は、ユーザがシミュレーション装置に入力することにより与えられる。
【００２５】
シミュレーション装置は、図２に示すように、複数の粒子１０の各々に対して複数（Ｐ個）のビーズ１１を定義する（ステップＳ２）。１つの粒子に対応するビーズの個数Ｐは２以上の整数である。さらに、複数のビーズの各々の質量Ｍを設定する（ステップＳ３）。実施例においては、ＰＩＭＤ法における時間τに、量子系における実時間ｔと同一の物理的な意味を持たせるために、Ｍ＝ｍ／Ｐを満たすように質量Ｍを設定する。ここで、ｍはＰ個のビーズに対応する元の粒子の質量を表す。
【００２６】
その後、シミュレーション装置は、元の粒子の初期条件に基づいて、複数のビーズの各々の初期条件を設定する（ステップＳ４）。例えば、Ｐ個のビーズの重心位置が、元の粒子の位置に一致し、Ｐ個のビーズの平均速度が元の粒子の速度に一致するように、各ビーズの初期条件を設定する。
【００２７】
次に、複数のビーズの挙動を、分子動力学法を用いてシミュレーションする（ステップＳ５）。このとき、式（６）の運動方程式を用いる。シミュレーションが終了すると、シミュレーション結果を出力する（ステップＳ６）。例えば、粒子の各々に対応するＰ個のビーズの重心の位置の時間変動を図形として表示する。
【００２８】
次に、図３を参照して、陽子に電子を衝突させるシミュレーションを例にとって、上記実施例の持つ優れた効果について説明する。
【００２９】
まず、シミュレーションの条件について説明する。座標ｒｐ、質量ｍｐの陽子と、座標ｒｅ、質量ｍｅの電子との運動を、その重心座標Ｒと、相対座標ｒに座標変換する。重心座標Ｒ及び相対座標ｒは、以下の式で表される。
【数１１】
2019046313.tif 000013
【００３０】
陽子と電子との重心は、等速直線運動をする。相対座標ｒは、ポテンシャルＶ（ｒ）中における換算質量μの粒子の運動方程式に従う。ポテンシャルＶ（ｒ）及び換算質量μは、以下の式で表される。
【数１２】
2019046313.tif 000014
ここで、ｅは電気素量、ε０は真空の誘電率を表す。
【００３１】
本シミュレーションでは、相対座標ｒの時間変化を求めた。ただし、ｒ＝０でポテンシャルＶ（ｒ）及び粒子に作用する力が発散してしまうため、シミュレーションではｒ＝０の近傍ではポテンシャル及び力が発散しないように、誤差関数ｅｒｆとボーア半径ａ０を用いて、式（１２）のポテンシャルＶ（ｒ）に代えて、以下のように丸めたポテンシャルＶティルダを用いた。
【数１３】
2019046313.tif 000015
【００３２】
以下の３種類のシミュレーションを行った。
第１のシミュレーションでは、換算質量μの古典的粒子の位置の時間変化を、分子動力学法により求めた。古典的粒子の位置をｒ１、速度をｖ１で表す。
【００３３】
第２のシミュレーションでは、量子系をＰＩＭＤ法で記述し、Ｐ個のビーズの位置の時間変化を求めた。Ｐ＝３２とし、ビーズの各々の質量をμ／Ｐに設定した。ビーズの位置をｒ２（ｓ）、速度をｖ２（ｓ）で表す。ここで、ｓはビーズに付された通し番号であり、１以上Ｐ以下の整数である。
【００３４】
第３のシミュレーションでは、ビーズの各々の質量をμに設定した。その他の条件は第２のシミュレーションの条件と同一である。
【００３５】
第１のシミュレーションにおける古典的粒子の位置の初期値、第２のシミュレーションにおけるＰ個のビーズの重心の位置の初期値、及び第３のシミュレーションにおけるＰ個のビーズの重心の位置の初期値を同一にした。さらに、第１のシミュレーションにおける古典的粒子の速度の初期値、第２のシミュレーションにおけるＰ個のビーズの平均速度の初期値、及び第３のシミュレーションにおけるＰ個のビーズの平均速度の初期値を同一にした。第２のシミュレーションにおけるＰ個のビーズの位置の初期値と第３のシミュレーションにおけるＰ個のビーズの位置の初期値とを同一とした。第２のシミュレーションにおけるＰ個のビーズの各々の速度の初期値と第３のシミュレーションにおけるＰ個のビーズの各々の速度の初期値とを同一とした。
【００３６】
すなわち、初期条件は以下の式を満たす。
【数１４】
2019046313.tif 000016
【００３７】
図３Ａ、図３Ｂ、及び図３Ｃは、それぞれ第１のシミュレーション結果、第２のシミュレーション結果、及び第３のシミュレーション結果を表すグラフである。図３Ａ、図３Ｂ、及び図３Ｃにおいて、重心座標Ｒを原点に固定した。重心座標Ｒは、ほぼ陽子の位置に相当する。図３Ａ、図３Ｂ、及び図３Ｃのグラフは、ポテンシャルＶティルダ中における相対座標ｒの時間変化を示している。相対座標ｒは、ほぼ電子の位置を表す。図３Ａ、図３Ｂ、及び図３Ｃにおいて横方向に並ぶ５個のグラフは、左側のグラフから右側のグラフに向かって、時間が経過していることを示している。
【００３８】
図３Ｂに示したＰ個のビーズ１１の重心の軌跡は、図３Ａに示した粒子１０の古典的解析による軌跡とほぼ一致する。ところが、図３Ｃに示したＰ個のビーズ１１の重心は、図３Ａに示した粒子１０の古典的解析による挙動とは全く異なる挙動を示している。このことから、量子系の粒子に対応する複数のビーズの各々の質量を、元の粒子の質量の１／Ｐ倍に設定してＰＩＭＤ法を用いてシミュレーションした結果は、古典系に対応した正しい時間発展を与えることがわかる。これに対し、量子系の粒子に対応する複数のビーズの各々の質量をその他の値に設定してＰＩＭＤ法によりシミュレーションした場合には、ＰＩＭＤ法における時間発展が、元の粒子系における時間発展に対応しているということはできない。
【００３９】
次に、図４を参照して実施例によるシミュレーション装置について説明する。
図４は、実施例によるシミュレーション装置のブロック図である。シミュレーション装置は一般的なコンピュータで構成され、中央処理ユニット（ＣＰＵ、処理装置）５０、メモリ５１、入出力装置５２、及び補助記憶装置５３を含む。これらは、バス５４で相互に接続されている。
【００４０】
入出力装置５２は、キーボード、マウス等のポインティングデバイス、ディスプレイ、リムーバブルメディアのリーダライタ、通信装置等を含む。ディスプレイは、ユーザの操作に必要な各種ウィンドウやデータ等を表示する。通信装置は、外部機器とデータ通信を行う。ユーザは、入出力装置を操作することにより、コンピュータに対する指示、各処理に必要なデータの入力等を行う。また、処理結果が、入出力装置に出力される。
【００４１】
補助記憶装置５３は、ハードディスク等で構成され、コンピュータのＯＳ、シミュレーションプログラム、シミュレーションに必要なデータや処理結果等を蓄積する。
【００４２】
メモリ５１は、ＲＯＭ（Ｒｅａｄ  Ｏｎｌｙ  Ｍｅｍｏｒｙ）やＲＡＭ（Ｒａｎｄｏｍ  Ａｃｃｅｓｓ  Ｍｅｍｏｒｙ）等で構成される。メモリ５１は、ＣＰＵ５０により補助記憶装置５３から読み出されたシミュレーションプログラム等を格納する。
【００４３】
ＣＰＵ５０は、コンピュータのＯＳ、メモリ５１に格納されているシミュレーションプログラム等に基づいて、各種演算や入出力装置５２及び補助記憶装置５３の制御を行う。
【００４４】
次に、シミュレーション装置の動作について説明する。ＣＰＵ５０は、ステップＳ１（図１）の処理を実行するために、初期条件等のシミュレーション条件を入出力装置５２から取得する。例えば、ユーザが入出力装置５２を操作することにより、初期条件等のシミュレーション条件を入力する。または、ＣＰＵ５０は、初期条件等のシミュレーション条件をリムーバブルメディアから読み込む。その後、ＣＰＵ５０は、ステップＳ２からステップＳ５までの処理を実行する。ステップＳ６において、ＣＰＵ５０は、例えば解析結果を図形としてディスプレイに表示させる。
【００４５】
上記実施例は例示であり、本発明は上記実施例に制限されるものではない。例えば、種々の変更、改良、組み合わせ等が可能なことは当業者に自明であろう。
【符号の説明】
【００４６】
１０  粒子
１１  ビーズ
５０  中央処理ユニット（ＣＰＵ、処理装置）
５１  メモリ
５２  入出力装置
５３  補助記憶装置
５４  バス
【図１】
2019046313.tif 000017
【図２】
2019046313.tif 000018
【図３】
2019046313.tif 000019
【図４】
2019046313.tif 000020