読者です 読者をやめる 読者になる 読者になる

nano_exit

量子力学、固体物理、fortran、python、etc

水素原子の波動関数の数値計算

束縛状態の波動関数の計算に抵抗がありまくりんぐなので、作ってみました。
普段は連続状態を計算していて、境界条件って言っても解析解につなぐだけなので、「正しいエネルギー(固有値)を見つける」というプロセスに心のハードルが無限の井戸型ポテンシャル状態でした。

コードは以下のサイトを参考に、pythonで書きました。
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/radschro.html#S_Orbital

自分で分かっていなかったポイントを整理すると、

  • 原点周辺も無限遠周辺も、漸近解を初期値に使う。厳密な初期値かどうか気にするよりも、大雑把でも物理的要請を満たしている方がよっぽど大事。
  • 漸近解を用いるのであれば、ルンゲクッタだろうが予測子修正子だろうが、初期値はただ漸近解の値を使えば良いから、運用上の違いはそんなに気にしなくて良い。もちろん誤差は違うだろうが。
  • 最初に与える試行固有値は、やっぱり解析解をある程度は参照する必要があると思われる。例えば節の数が全然合わないときには真のエネルギー近傍に無理やり持って行ったりとかは知らないと出来ないなぁと。
  • 試行固有値が真の固有値かどうかの判定として、原点近傍と無限遠近傍の両方から計算して、両者の解があるところで一次微分までで繋がるかどうかを見ることが使える。で、そのようなある点で微分不連続を与えるようなポテンシャルはデルタ関数に違いないということを考え、デルタ関数の分を一次摂動としてエネルギーから差っ引く。これで真のエネルギーに近づくでしょうという発想。
  • 本来は真の波動関数があって、そこにデルタ関数があるから波動関数が滑らかにならず、真の波動関数で一次摂動を計算したものを差っ引くべきなんだが、ここで現在得られている不連続な波動関数を真の波動関数と読みかえ、それで近似的に一次摂動を計算している。なかなかアクロバットやなぁと感心する。
  • 節を数えるのも結構苦手意識があったけど、単純に前後の波動関数の値を掛けて負になったら一個節があると思えば十分だった。
  • 波動関数の規格化が長方形のやつを使っていて、普段スプライン積分で中身が半分ブラックボックスで使っている身としてはカルチャーショックだった。多分自分で一からやろうとしたら、「これで良いの?」と自問自答して前に進めなかったと思う。数値計算には漢気が大事だと学んだ。

自分で作ると、やっぱり感慨深いものがある。
試しに3s軌道が初期値からどのように修正されるかを見てみる。
f:id:koideforest:20170128201714p:plain

以下、コードを示す。

import numpy as np
import matplotlib.pyplot as plt
import prettyplotlib as ppl

z     = 1.
dx    = 0.0625
r_num = 255 + 1
r_far = r_num - 1
r_o   = 160
output_r_max = 64.0

r, u = [], []
for i in range( r_num ):
  r.append( np.exp( dx * ( i - r_o ) ) / z )
  u.append( - z )

def P( i, l, E ):
  return 2. * r[i] * u[i] - 2. * ( r[i]**2 ) * E + l * ( l + 1 )

def Pzero( l, E ):
  for i in range( r_far, -1, -1 ):
    if P( i, l, E ) < 0.: break
  return i

def Calc_dG( i, G, F, l, E ):
  return F

def Calc_dF( i, G, F, l, E ):
  return P( i, l, E ) * G + F

def HammingA( i_fit, l, E ):
  node     = 0
  Gpc, Fpc = 0., 0.
  G, F     = np.zeros( r_num ), np.zeros( r_num )
  dG, dF   = np.zeros( r_num ), np.zeros( r_num )
  for i in range( 0, 4 ):
    G[i]  = r[i] ** ( l + 1 )
    F[i]  = ( l + 1 ) * G[i]
    dG[i] = Calc_dG( i, G[i], F[i], l, E )
    dF[i] = Calc_dF( i, G[i], F[i], l, E )
  for i in range( 4, i_fit + 1 ):
    Gp = G[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dG[i-1] - dG[i-2] + 2. * dG[i-3] )
    Fp = F[i-4] + ( 4.0 * dx / 3.0 ) * ( 2. * dF[i-1] - dF[i-2] + 2. * dF[i-3] )
    Gm = Gp - ( 112. / 121. ) * Gpc
    Fm = Fp - ( 112. / 121. ) * Fpc
    dGm = Calc_dG( i, Gm, Fm, l, E )
    dFm = Calc_dF( i, Gm, Fm, l, E )
    Gc = ( 9. / 8. ) * G[i-1] - ( 1. / 8. ) * G[i-3] + ( 3. * dx / 8. ) * ( dGm + 2. * dG[i-1] - dG[i-2] )
    Fc = ( 9. / 8. ) * F[i-1] - ( 1. / 8. ) * F[i-3] + ( 3. * dx / 8. ) * ( dFm + 2. * dF[i-1] - dF[i-2] )
    Gpc = Gp - Gc
    Fpc = Fp - Fc
    G[i] = Gc + ( 9. / 121. ) * Gpc
    F[i] = Fc + ( 9. / 121. ) * Fpc
    dG[i] = Calc_dG( i, G[i], F[i], l, E )
    dF[i] = Calc_dF( i, G[i], F[i], l, E )
    if G[i] * G[i-1] < 0. : node += 1
  return node, G, F

def HammingB( i_fit, l, E ):
  node     = 0
  Gpc, Fpc = 0., 0.
  G, F     = np.zeros( r_num ), np.zeros( r_num )
  dG, dF   = np.zeros( r_num ), np.zeros( r_num )
  kapper   = np.sqrt( np.abs( 2. * E ) )
  for i in range( r_far, r_far - 4, -1 ):
    G[i]  = np.exp( - kapper * r[i] )
    F[i]  = - kapper * r[i] * G[i]
    dG[i] = Calc_dG( i, G[i], F[i], l, E )
    dF[i] = Calc_dF( i, G[i], F[i], l, E )
  for i in range( r_far - 4, i_fit - 1, -1 ):
    Gp = G[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dG[ i + 1 ] - dG[ i + 2 ] + 2. * dG[ i + 3 ] )
    Fp = F[ i + 4 ] - ( 4.0 * dx / 3.0 ) * ( 2. * dF[ i + 1 ] - dF[ i + 2 ] + 2. * dF[ i + 3 ] )
    Gm = Gp - ( 112. / 121. ) * Gpc
    Fm = Fp - ( 112. / 121. ) * Fpc
    dGm = Calc_dG( i, Gm, Fm, l, E )
    dFm = Calc_dF( i, Gm, Fm, l, E )
    Gc = ( 9. / 8. ) * G[ i + 1 ] - ( 1. / 8. ) * G[ i + 3 ] - ( 3. * dx / 8. ) * ( dGm + 2. * dG[ i + 1 ] - dG[ i + 2 ] )
    Fc = ( 9. / 8. ) * F[ i + 1 ] - ( 1. / 8. ) * F[ i + 3 ] - ( 3. * dx / 8. ) * ( dFm + 2. * dF[ i + 1 ] - dF[ i + 2 ] )
    Gpc = Gp - Gc
    Fpc = Fp - Fc
    G[i] = Gc + ( 9. / 121. ) * Gpc
    F[i] = Fc + ( 9. / 121. ) * Fpc
    dG[i] = Calc_dG( i, G[i], F[i], l, E )
    dF[i] = Calc_dF( i, G[i], F[i], l, E )
    if G[i] * G[ i + 1 ] < 0. : node += 1
  return node, G, F

def Normalization( GA, FA, GB, FB, i_fit ):
  for i in range( i_fit ):
    G[i] = GA[i]
    F[i] = FA[i]
  coeff = GA[ i_fit ] / GB[ i_fit ]
  for i in range( i_fit, r_num ):
    G[i] = GB[i] * coeff 
    F[i] = FB[i] * coeff
  sum = 0.
  for i in range( r_num ):
    sum += ( G[i]**2 ) * r[i] * dx # dr = e^x dx = r dx
  norm = 1. / np.sqrt( sum )
  for i in range( r_num ):
    G[i] *= norm
    F[i] *= norm
  return G, F

def Solver( n, l, E ):
  true_node = n - l - 1
  i_fit = Pzero( l, E )
  nodeA, GA, FA = HammingA( i_fit, l, E )
  LogA = FA[ i_fit ] / ( r[ i_fit ] * GA[ i_fit ] )
  nodeB, GB, FB = HammingB( i_fit, l, E )
  LogB = FB[ i_fit ] / ( r[ i_fit ] * GB[ i_fit ] )
  node = nodeA + nodeB
  G, F = Normalization( GA, FA, GB, FB, i_fit )
  dE0 = - 0.5 * ( G[ i_fit ] ** 2 ) * ( LogB - LogA )
  if not true_node == node:
    dE = float( true_node - node ) / ( n * n * n )
  elif dE0 > 0.01:
    dE = dE0 * 0.5
  else:
    dE = dE0
  E = E + dE
  print 'true_node', true_node, 'node', node
  print 'dE0', dE0, 'dE', dE
  print 'E', E
  ppl.plot( r, G )  
  if not E < 0:
    print 'E', E
    print ' not a bound state'
    exit()
  return node, G, F, dE

eps = 1e-6
for n in range( 1, 5 ):
  E = - 0.5 * ( z**2 ) / ( n**2 ) - 0.05
  for l in range( 4 ):
    if l >= n: continue
    print 'n=', n, 'l=', l
    G, F = np.zeros( r_num ), np.zeros( r_num )
    dE  = 1e-5
    while abs( dE ) > eps:
      node, G, F, dE = Solver( n, l, E )
      E = E + dE
    plt.xlim( 0, 50 )
    plt.ylim( -0.8, 0.8 )
    plt.show()