nano_exit

基礎的なことこそ、簡単な例が必要だと思うのです。

モンテカルロ積分で水素原子のエネルギー期待値を確認

「スレーター型軌道の数値積分は大変」と呪文のように聞かされて続けて心にハードルが出来てしまったので、ここらで打破することを試みた。
前回、多次元のモンテカルロ積分をまとめた。
koideforest.hatenadiary.com
今回はスクリプトを組んで、どれくらいちゃんと動くのか見てみることにした。
目標は、非相対論水素原子の解析解を使って、ハミルトニアンの期待値を出す。
以下、スクリプトのポイントをまとめる。

  • 運動エネルギーを出すところの微分は、中心差分法を使った。簡単便利。
  • 積分範囲は、水素原子波動関数の値が0.01を下回る距離をrmaxと定義し、その範囲で水素原子の8分の1が入る箱を考えた。対称性から、その箱の中の積分値を8倍すれば全体の積分値が求まる。
  • そんなことをしたのは、乱数の範囲が [ 0, 1 ] の0始まりなので、全体を平行移動したりして水素原子全体を箱の中に無理矢理入れるのは余り得策ではなく感じたから。今回は素直に水素原子の原点を0にした。
  • 箱の中の全体で積分値を求めるから、乱数点に関する制限は必要ない。下手にif分を入れると遅くなるらしいので、そこまで一般化してスクリプトを作るのは余り良くないかもしれない。今回はそこの一般化は行っていない。

以下ソース。

import numpy as np
import time

def MonteCalro3d_ws( func, scales ): # ws: whole space
    M = 10**6
    vec_rand = np.random.rand( M * 3 ).reshape( M, 3 )
    # Caution: the number of samplings is still M, not M * 3
    vec_rand = np.array( scales ) * vec_rand
    tot = np.sum( [ func( v[0], v[1], v[2] ) for v in vec_rand ] )
    return np.prod( scales ) * tot / float( M )

def h1s( x, y, z ):
    vec = np.array( [ x, y, z] )
    r   = np.linalg.norm( vec )
    return 2. * np.exp( - r ) / np.sqrt( 4. * np.pi )

def d2func_dx2( func, x ):
    h = 1.0e-7
    return ( func( x + 2.*h ) - 2.*func(x) + func( x - 2.*h ) ) / (2.*h)**2

def integrand( x, y, z ):
    v = np.array( [ x, y, z] )
    r = np.linalg.norm( v )
    deriv2x = d2func_dx2( lambda x: h1s( x,     v[1], v[2] ), v[0] )
    deriv2y = d2func_dx2( lambda y: h1s( v[0],     y, v[2] ), v[1] )
    deriv2z = d2func_dx2( lambda z: h1s( v[0],  v[1],    z ), v[2] )
    ke = - ( deriv2x + deriv2y + deriv2z ) / 2.
    pot_n = - h1s( v2[0], v2[1], v2[2] ) / r
    return h1s( v[0], v[1], v[2] ) * ( ke + pot_n )

rmax = max( [ x_ for x_ in np.linspace( 0, 5, 1000 ) \
              if h1s( x_, 0., 0. ) / h1s( 0., 0., 0. ) > 0.01 ] )
scales = [ rmax, rmax, rmax ]

start = time.time()
result = 8. * MonteCalro3d_ws( integrand, scales )
print( "Hydrogen atom energy = {}".format( result )  )
end = time.time()
print( end - start )

結果:
-0.500215806919
148.051234961

  • 原子単位系は原子単位系でも、Hartree unitで計算しているので、理論値は-0.5でちゃんと計算できている。
  • ちなみにHartreeかRydbergかは運動エネルギーが 1/2されているかどうかで判断できる。
  •  1/2の無いRydbergの場合には、核ポテンシャルの項も二倍して理論値は-1になる。
  • 細かい値を忘れたが、ざっくり 1 Hartree ~ 27 eV, 1 Rydberg ~ 13 eV, 1 Hartree = 2 Rydberg である。最後の等式は厳密。
  • 値の安定性はサンプリング数 Mに依存するが、 M = 10^6であれば有効数字一桁は安心できると思う。
  • cythonとか使えばもっと速くなる。