Basis Quantum Monte Carlo 法による1次元及び3次元調和振動子内の2フェルミ粒子の計算

八木 徹, 長嶋 雲兵


Return

1 はじめに

多電子系の原子や分子の波動関数は、ab initio計算に代表される様々な手法により求めることができる。特にFull CI計算においては、Born-Oppenheimer近似と用いた基底関数をパラメータとした制限の中で、最も厳密な解を得ることができる。しかし、Full CI計算は電子数の増加とともに計算量が膨大なものとなる事がよく知られている[1]。
Full CI計算に代わり、多電子系の電子状態を求める手法として量子Monte Carloシミュレーションがあり、これまでに様々な改良、発展がなされている[2 - 5]。Monte Carloシミュレーションでは、電子数Nに対して計算量がN2程度の増加量となり、また、効率の良い並列化も可能である。このため、量子Monte CarloシミュレーションはFull CIとほぼ同等の効果を含んだ大規模系の電子状態の計算を行う可能性を持つ。
しかし、量子Monte Carloシミュレーションではフェルミ粒子の反対称性により負の確率密度が生じる問題が発生する。これに対しては、現在一般に用いられる拡散量子Monte Carlo法では、Fixed-node近似などの手法により、反対称性の問題を含む多電子系の計算を行うことができるようになっている[2 - 4]。
現在広く行われているFixed-node近似を用いた拡散量子Monte Carlo法では、事前に節面の情報が必要となる。このためガイド関数と呼ばれる節面を記述する関数を利用するが、ガイド関数を得るためにab initio分子軌道計算などを行うため、準備段階に高い計算コストが必要となる。また、ガイド関数の精度がMonte Carlo計算結果に大きく反映する。
ガイド関数を使うことなしに、多電子系の反対称性問題を取り扱うことが可能な別の量子Monte Carlo法として、OksuzによりBasis Quantum Monte Carlo(BQMC)法が提唱されている[6 - 8]。
BQMC法においては電子状態の符号を2粒子の座標の交換で表しており、事前の節面情報を与えるガイド関数を必要とせずに計算を行うことができる。このため、ガイド関数を求める手間やその精度に左右されずに、原子・分子の基底状態を求めることができる。また、BQMC法においては、Monte Carloシミュレーションの状態分布から、系の波動関数を求めることができ、エネルギー計算だけでなく、さまざまな物理量の計算や、他の問題に適用できる可能性を持っている。
本論文では、BQMC法を用いて1次元及び3次元の調和振動子内における2個のフェルミ粒子に対する計算を行った。また、Oksuzと異なるエネルギーの計算手法を用いて基底状態のエネルギーを求めて結果を比較した。

2 方法

BQMCの手法についてはOksuzの論文[6, 7]に記述されている。しかし、これらの文献は通常入手困難であるため、本節では通常の読者のために、概要を解説する。

2. 1 Basis Quantum Monte Carlo 法

対象とする系のハミルトニアン に対する固有状態を|i>とすると、任意の状態|x>は以下の式で展開される。

ここで、Ciは展開の係数である。この状態 |x> に、時間発展を表わす を作用させることにより、次式が得られる。

ここで、tは虚数時間、Eii番目の固有状態の固有値、及びETは試行エネルギーと呼ばれ、任意のエネルギーシフト量である。ここで、ETを基底状態の固有値E0と選ぶことにより、

となる。E0は基底状態のエネルギー固有値であり、E0 < Eiであるため、t®¥の時に |x>® |0 > が得られる。
Monte Carloシミュレーションでは、tは小さな値(t< 1)として、 を繰り返し作用させることにより基底状態を求める。
ある時間ステップtにおける状態を |x,t>とすると、以下の関係が得られる

このとき、

であり、基底状態のエネルギー固有値E0

として得られる。
BQMC法では、前述の繰り返し演算の関係を用い、Monte Carloシミュレーションにおける状態の時間発展を記述する[6, 7]。始めに、空間を等間隔の格子で分割し、各格子点を中心とする基底関数を導入する。この基底関数としては、中心となる格子点上での値が1、それ以外の場所では0〜1の間の正の値を取り、中心から離れると急速に減衰するような関数を用いる。計算上の取り扱いのしやすさから、以下に示すGauss型の関数を用いる。

ここで、bは格子点の間隔、添え字iはi番目の格子点、riはi番目の格子点の座標を表わす。(7)式は、i番目の格子点を中心としたGauss関数であり、等方的かつ中心から離れた位置では急速に減衰する。この基底関数を用い、2個のフェルミ粒子系の波動関数を以下のようにあらわす。

r1及びr2は、それぞれフェルミ粒子の空間座標を表わす。添え字ijに対する和の制限(i > j)は、同等の項を重複して加えることを避けるためのものである。(8)式の波動関数に対して、時間発展の演算子 を作用させる。その際に、tを十分小さな値として展開した

を用いる。反対称化された波動関数である(8)式に、(9)式を作用させた結果を整理すると、BQMCの時間発展は次式で表わされる[6 - 8]。

ここで、Uは粒子の拡散を表わす項であり、3次元と1次元でそれぞれ、

となる。ここで、中かっこ[ ] でくくられた項をフェルミ孔項と呼ぶことにする。これらは、2個の粒子が(rk, rl)から(rm, rn) (または(xk, xl)から(xm, xn))へ移動する確率を表わす。フェルミ孔項により、2個のフェルミ粒子が同じ座標を占める場合 (rk = rlrm = rnなど)の確率Uは0となり、2個の粒子が同時に同じ位置を占めることができない。これはパウリの排他原理を表している。また、フェルミ孔項は、2個の粒子間の距離が離れると急速に1に近づく。このとき並行スピン粒子間の相関はなくなり、通常の拡散と同じものになる。
フェルミ孔項は0から1の間の値を取るため、実際の計算では(11)、(12)式の評価にvon Neumannの採用-棄却法を用いる。まず、(11)および(12)式のフェルミ孔項を無視し、通常の拡散項を用いて配置の新しい座標を生成する。次に0と1の間の一様乱数rを生成してrがフェルミ孔項の値よりも小さい場合にその配置を採用し、それ以外の場合は棄却して新しい配置の生成に戻る。これにより、フェルミ孔項で重みづけされた拡散による配置を生成することができる。
Lは配置の生成・消滅を表わす分岐項であり、次式で与えられる。

ここで、Vはフェルミ粒子の座標に依存して決まるポテンシャルである。
dは時刻tにおける状態を表わし、

となる。ここで、添え字mnの和の制限については、様々な取り方が考えられるが、フェルミ粒子のx座標に注目し、粒子1と粒子2のx座標が常に、-¥ < x2 < x1 < ¥ となる状態のみを考慮する。
また、Uは粒子の拡散による遷移確率を与えるが、和の制限の存在のために規格化されていない。このため、以下の定数Sを導入してUを規格化する。

この規格化定数により式を

と表わす。U'、L'の個々の要素は

となる。
BQMC法でのMonte Carloシミュレーションは、(16)式の繰り返し計算を基に実行する。実際のシミュレーションにおいては、状態dはフェルミ粒子の座標(r)と重み(w)の組(r1a, r2a, wa)で表される配置となる。ここでaは個々の配置を区別する添え字であり、N個の配置がある場合にa=1,2,...,Nとなる。
系のエネルギーを求める際にOksuz は、各配置の重みを1とし、エネルギーを配置数の平均として求めている[6, 7]。本論文では、いくつかの重みつき平均エネルギーを比較するために、時刻ステップtにおける配置の重みを次式で定めた。

ここで、Dwa(t)は時間ステップtにおけるa番目の配置の重みを表わし、wa(t)は時間ステップtまでのすべての重みの積である。

2. 2 試行エネルギー

一般的な拡散量子モンテカルロシミュレーションにおいて、(2)式内の試行エネルギーETは、基底状態のエネルギーE0に対して ET - E00となるようにシミュレーション過程のエネルギー平均値などを用いて表す。また、配置の数Nをシミュレーションの間、ある一定の値付近に保つ役割も担っている。
Oksuzは、ETをポテンシャルの値をシフトさせる定数として、シミュレーション中の粒子数が安定となるように事前に決めるとしている[6, 7]。しかし、シミュレーションの全過程において定数を用いて粒子数を一定に保つことは困難であるため、一般的な量子モンテカルロの手法に従い、エネルギーのシフト量として以下を用いた[2, 5]。

ここで、Eestは、シミュレーションの各時間ステップにおけるエネルギー平均値で、Wtは時間ステップtにおける各配置の重みの総和、W0は目標とする重みの基準値である。Eestとしては、2.4節に示すエネルギーの平均値を用いた。

2. 3 反対称性

通常の拡散量子Monte Carlo計算においては、何らかの方法で事前に決定したガイド関数を用いて節面を表わす。また、フェルミ粒子の拡散において、粒子がガイド関数の節面を超えて移動した場合には、配置の符号が変化したとする。(実際には、符号が変化する場合には、配置を消滅させるか、粒子の移動をキャンセルする。)
文献[6]に従い、BQMC法における符号問題を以下に示す手順で取り扱う。
  1. 時間ステップtにおける2つのフェルミ粒子の座標を(rk, rl)、差分ベクトルをvi = rl - rkとする。
  2. 拡散過程による2粒子の移動後の座標を(rm, rn)、差分ベクトルをvf = rn - rmとする。
  3. 新しい粒子の座標は、(11)または(12)式に従い生成する。このとき、2つのベクトルの内積vi.vfが正となる(rm, rn)のみを生成する。
    1. vfx < 0の領域にある場合、基底関数の和の制限(2.1節参照)により、この配置を記述する基底関数が存在しない。このため粒子の座標を交換し、配置の符号を変える。
    2. vfx > 0の領域にある場合、その配置を採用してシミュレーションを継続する。
  4. 配置の符号が変化した場合、その配置は消滅させる。
上記の手順により、BQMC法においては、ガイド関数の仮定なしに、配置の符号を取り扱うことが可能となる。

2. 4 エネルギー計算

平衡化の過程を経たのちに、各配置の重みを用いてエネルギーの固有値を求める。本論文では以下に示す式を用いて基底状態のエネルギーを求め、結果を比較した。

ここで、シミュレーションは、タイムステップt = 0からTまで行うものとする。また、 は全ての配置に対する和を表わす。
エネルギーの統計誤差は以下の式で求めた。

ここで、Eは(22)又は(23)式で求めた平均エネルギー、Eiはあるステップでの平均エネルギー、Mは平均を計算したステップの数である。
(22)式のエネルギー(Egr1)に対する統計誤差を求める際には、(24)式のEiとして、各時刻ステップでのエネルギーの値を用いた。(23)式のエネルギー(Egr2)に対する統計誤差では、10個の時刻ステップ分のエネルギーの平均値をEiとして用いた。

3 シミュレーションの実際

本方法は無限回の行列積を行うことで全エネルギーを得る方法である。ここでのMonte Carloシミュレーションは、本研究で現れる性質を持つ行列の無限回の行列積をシミュレートする。本節では実際のMonte Carloシミュレーションについてアルゴリズムを説明する。

3. 1 Monte Carloアルゴリズム

Monte Carloシミュレーションは、以下の手順で実施した。
  1. 初期化
  2. 個々の配置に対して以下を実施
  3. 各配置を重みの値に応じて生成・消滅させる
  4. ステップ2.と3.を繰り返し、平衡化が達せられたのちに系のエネルギーを求める。
ステップ3における配置の生成・消滅にはUmrigarらの方法[5]を用いた。

3. 2 調和振動子ポテンシャル内のフェルミ粒子

相互作用をせず、並行スピンの関係にある2つのフェルミ粒子が、調和振動子ポテンシャル内に存在する系のシミュレーションを行った。ポテンシャルは1次元と3次元それぞれに対して以下を用いた。

3. 3 各種パラメータ

1次元と3次元の計算に対して、格子間隔としてそれぞれb = 0.1(a.u.) 及び b = 0.02 (a.u.)を用いた。格子間隔btt = 3b2の関係にあるとし[7]、それぞれの格子間隔に対応してt = 0.03 及びt = 0.0012とした。
粒子の初期配置はランダムに作成し、十分な平衡化の後に106ステップのシミュレーションを行った。平衡後の100ステップごとの配置に対して各種エネルギーの平均量を求めた。
試行エネルギーETは(21)式で計算した。ETの値の更新は、1次元調和振動子の系では100ステップに1回、3次元調和振動子では毎時間ステップ行った。

4 結果と考察

4. 1 1次元調和振動子

1次元調和振動子内の平衡スピンをもつ2個のフェルミ粒子の系に対して、BQMC法で求めた基底状態のエネルギーをTable 1に示す。Egr1及びEgr2はそれぞれ、(22)、(23)式のエネルギー計算式を用い、総ステップ数106のうちの104個の配置に対する計算結果である。これらは、ともに基底状態のエネルギーを精度よく再現している。
表のreference1及び2は、同じ系に対するOksuzの結果をまとめたものであり、それぞれエネルギーを今回の計算と同程度、及び10倍以上のステップ数で計算した結果である。Egr1による平均値では、同程度のステップ数での平均値(reference1)と比べて、同等以上の精度でエネルギーが得られている。また、Egr2による平均値では、Egr1と同じ長さの時間ステップの計算で、reference2に近い高精度の結果が得られている。

Table 1. Ground state energies of parallel spin two Fermi particles in a 1 dimensional harmonic potential (Hartree)
Energys
Egr12.0019470.00126
Egr22.0003530.000488
Reference1 1)2.005660.0077
Reference2 2)2.0001330.00176
Exact2.0-
1) Ref.[7]. Energy averaged in 7×106 step configurations2) Ref.[7]. Energy averaged in 7×107 step configurations


Figure 1. Energy variation during 1000 steps

エネルギーの時間ステップに対する変化をFigure 1に示す。ランダムに生成した初期配置から、100ステップ程度で平衡配置に達していることがわかる。それ以降はエネルギーが2.0周辺を最大0.3程度の誤差範囲で振動している。
フェルミ粒子の空間分布をFigure 2に示す。ここでは、2個のフェルミ粒子のx座標(x1およびx2)を、それぞれの軸として粒子の分布をプロットしている。基底関数の和の制限(- ¥< x2 < x1< ¥)のために、x2 < x1の領域にのみ粒子が分布する結果が得られる。図より、2個の粒子が同時に同じ座標を占めていないことが分かる。これは、フェルミ粒子が持つ反対称性の効果が正しくとりこまれていることを示している。また、粒子分布が最大となる点は、反対称性の効果により、ポテンシャルの極小点(0.0, 0.0)からややずれた地点(-0.7, 0.7)付近となっている。


Figure 2. Distribution of two Fermi particles along their coordinates, where maximum value is 1.0

1次元の系に対して要した計算時間は CPU:Intel Xeon 3.4GHz、メモリ2GBを搭載したマシンで400秒程度であった。

4. 2 3次元調和振動子

3次元調和振動子内の平衡スピンをもつ2個のフェルミ粒子の系に対して、BQMCシミュレーションで求めた基底状態のエネルギーをTable 2に示す。いずれの計算方法においても、基底状態のエネルギーを精度よく再現している。1次元の場合と異なり、エネルギーの精度はEgr1Egr2ともに同程度となっている。また、referenceの値よりも高い精度でエネルギーを得られた。

Table 2. Ground state energies of parallel spin two Fermi particles in a 3 dimensional harmonic potential (Hartree)
Energys
Egr13.500530.00122
Egr23.500510.000972
reference 1)3.4980.015
Exact3.5-
1) Ref.[6]. Energy averaged in 6000 step configurations

エネルギーの時間ステップに対する変化をFigure 3に示す。1次元系の計算よりも平衡に達する時間が長くなっている。平衡化に3000ステップ程度を要している。


Figure 3. Energy variation during 20000 steps

3次元の系に対する計算に要した時間は、1次元計算と同じマシンを用いて540秒程度となった。3次元の計算でも150秒程度(40%)の増加である。

4. 3 格子依存性

BQMCシミュレーションでは、計算領域を等間隔の格子に分割する。シミュレーション結果に対する格子依存性を確認するために、3次元調和振動子計算において、複数の格子間隔を用いた計算を行った。計算に用いた格子間隔は、b = 0.005, 0.01, 0.02, 0.03, 0.04, 0.06, 0.08, 0.10の8通りである。その際、tは、t = 3b2 の関係を保つように定めている。bt以外のパラメータについては、前節までと同等の条件を用いてBQMCシミュレーションを行った。エネルギーをEgr1で求めた結果をFigure 4に示す。


Figure 4. Mesh spacing dependence of the energy

格子間隔が0.02以下の計算では、いずれも3.50程度の正しい値を再現している。格子間隔が0.03以上では、エネルギーは厳密解よりも小さな値となり、格子間隔が広くなるほどその傾向が強くなる。

4. 4 配置数の変化

1次元、及び3次元調和振動子の計算を行った際の、各時間ステップにおける配置数の変化をFigure 5に示す。
1次元調和振動子の場合、配置の数は200から780程度の間で変化し、平均して約400個の配置が存在している。1次元のシミュレーションでは、100ステップごとに試行エネルギーETの値を更新したため、周期的に配置の数が変化する結果となっている。
3次元調和振動子の場合、配置の数は390から420程度の間で変化している。3次元の場合は、配置の数を保つために、1ステップごとに試行エネルギーETの値を更新した。このため1次元の場合に比べて配置数変化の少ない結果となっている。


Figure 5. Population variation during 20000 steps

BQMC法は、ポテンシャルの形をあらわに扱っているために、配置数の変化が大きいと考えられる。特に、ポテンシャルの値が発散する点を持つクーロン項で記述される原子・分子の系に適用する際には、配置数の変化や平均エネルギーの分散が、より大きくなることが予想される。分散を小さくする手法として、一般のMonte Carlo計算ではImportance SamplingやStratified Sampling等の手法が知られている。ガイド関数を必要としないBQMC法の利点を生かすためには、Stratified Samplingのような事前の分布を仮定しない手法を用いる必要があると考えられる。

5 まとめ

BQMC法を用いた計算により、並行スピンのフェルミ粒子からなる系の状態を効率よく求めることが可能であることが確認できた。この方法では、ガイド関数を仮定せずにシミュレーションを行うことができるため、さまざまな系への適用が期待できる。また、エネルギーの計算方法を変えることにより、より高い精度で基底状態のエネルギーを求めることができた。
今後は非相対論的な近似範囲で原子の計算を行い、方法の安定性を確認した上で分子に拡張することを試みる。

参考文献

[ 1] Schaefer H. F., Electronic Structure of Atoms and Molecules, Addison-Wesley (1972).
[ 2] Lester W. A. Jr., Recent Advances in Quantum Monte Carlo Method, World Scientific (1997).
[ 3] Anderson J. B., J. Chem. Phys., 63, 1499-1503 (1975).
[ 4] Reynolds P. J., Ceperley D. M., Alder B. J., and Lester W. A. Jr., J. Chem. Phys., 77, 5593-5603 (1982).
[ 5] Umrigar C. J., Nightingale M. P., and Runge K. J., J. Chem. Phys., 99, 2865-2890 (1993).
[ 6] Oksuz I., Arab. J. Sci. Eng., 9, 145-152 (1984).
[ 7] Oksuz I., Arab. J. Sci. Eng., 9, 239-249 (1984).
[ 8] Oksuz I., J. Chem. Phys., 81, 5005-5012 (1984).


Return