nano_exit

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

多電子原子の波動関数(途中)

今日、朝起きてからずっと多電子原子波動関数の作成プログラムやっているけど、なかなかうまく行かない。
コードは引き続き下記のサイトを参照。
http://www-cms.phys.s.u-tokyo.ac.jp/~naoki/CIPINTRO/CIP/atom.html

とりあえずs軌道は大丈夫。
でもp軌道になった途端に、束縛状態がなくなってしまう。遠心力ポテンシャルに負けちゃう。

Liで結果を比較したけれども、イオン化エネルギーはバッチリ合っている。
でもトータルが一緒にならない。
全く同じコードを書いてるわけじゃないから、全く同じじゃないのは良いにしても、イオン化エネルギー合ってるのにトータルエネルギー合わないのはめっちゃ気持ち悪い。
てかp, d軌道を計算したい。。。
休みの間にケリがつかないのは気持ち悪いがしょうがない。

以下、コードです。

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

""" calculation condition """
dx    = 0.0625
r_num = 255 + 1
r_far = r_num - 1
r_o   = 160
output_r_max = 64.0

# 1s, 2s, 2p, 3s, 3p
z   = 3.
deg = [ 2., 1. ] # degeneracy = occupancy
n   = [ 1, 2 ] # principal quantum number 
l   = [ 0, 0 ] # partial wave
eps = 1e-6

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

""" functions """
# simga = r^2 | \phi |^2
# total sigma
def Calc_sigma( ): # Gset and deg are global values
  sig = np.zeros( r_num )
  for i_d, d in enumerate( deg ): # in this case, G[ r, ni ] is needed.
    NormCheck( Gset[ i_d ] )
    for i_r, gi in enumerate( Gset[ i_d ] ):
      sig[ i_r ] += d * ( gi**2 )
  sum = 0.
  for i_r, s in enumerate( sig ):
    sum += s * r[ i_r ] * dx
  print 'total electron', sum
  return sig

def TotalEnergy( TE0 ):
  corre = 0.
  for i, sig in enumerate( sigma ):
    corre += sig * r[i] * ( 2. * Calc_uh( i ) + Calc_uxc( i ) ) * dx
  return TE0 - 0.25 * corre

def Calc_uc( i ):
  return - z

def Calc_uh( i ):
  uh = 0.
  for j in range( r_num ):
    if j < i :
      uh += sigma[j] * r[j] * dx # dr = \delta e^{\delta * x } dx = \delta r dx
    else: # i < j
      uh += sigma[j] * r[i] * dx
  return uh

def Calc_rs( i ):
   return ( 6. * r[i] * r[i] / sigma[i] )**( 1. / 3. )

def Calc_beta( rs ):
  return 1. + 0.0368 * rs * np.log( 1. + 21. / rs )

def Calc_uxc( i ):
  rs   = Calc_rs( i )
  beta = Calc_beta( rs )
  return - 2. * beta * ( 3. / ( 32. * np.pi * np.pi ) * r[i] * sigma[i] )**( 1. / 3. )

def Calc_u( i ):
  uc  = Calc_uc( i )
  uh  = Calc_uh( i )
  uxc = Calc_uxc( i )
  return uc + uh + uxc

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

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

def Calc_dG( i, Gi, Fi, i_n ):
  return Fi

def Calc_dF( i, Gi, Fi, i_n ):
  return P( i, i_n ) * Gi + Fi

def HammingA( i_fit, i_n ):
  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[ i_n ] + 1 )
    F[i]  = ( l[ i_n ] + 1 ) * G[i]
    dG[i] = Calc_dG( i, G[i], F[i], i_n )
    dF[i] = Calc_dF( i, G[i], F[i], i_n )
  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, i_n )
    dFm = Calc_dF( i, Gm, Fm, i_n )
    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], i_n )
    dF[i] = Calc_dF( i, G[i], F[i], i_n )
    if G[i] * G[i-1] < 0. : node += 1
  return node, G, F

def HammingB( i_fit, i_n ):
  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], i_n )
    dF[i] = Calc_dF( i, G[i], F[i], i_n )
  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, i_n )
    dFm = Calc_dF( i, Gm, Fm, i_n )
    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], i_n )
    dF[i] = Calc_dF( i, G[i], F[i], i_n )
    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, norm

def NormCheck( vec ):
  sum = 0.
  for i_v, vi in enumerate( vec ):
    sum += r[ i_v ] * ( vi**2 ) * dx
  print 'Normalization check', sum  

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

""" program main """
dTE = 1e-5
Gset = [ np.zeros( r_num ) for i_n, ni in enumerate( n ) ]
sigma = [ 1e-12 for i in range( r_num ) ] # initialize the global sigma
u = np.array([ Calc_u( i ) for i in range( r_num ) ])
TE, TE_old = 0., 0.
while abs( dTE ) > eps:
  for i_n, ni in enumerate( n ):
    E = - 0.5 * ( z**2 ) / ( ni**2 ) - 0.05
    print '-' * 20
    print ' ' * 5, 'n=', ni, 'l=', l[i_n]
    print '-' * 20
    G, F = np.zeros( r_num ), np.zeros( r_num )
    dE  = 1e-5
    if l[ i_n ] == 1: 
      ppl.plot( r, [ P( i_r, i_n ) for i_r in range( r_num ) ] )
      plt.ylim( -5, 5 )
      plt.show()
    while abs( dE ) > eps:
      node, G, F, dE = Solver( i_n )
      E = E + dE
    Gset[i_n] = G
    sigma     = Calc_sigma()
    """ merge old and new versions of u """
    u = np.array([ 0.6 * ui + 0.4 * Calc_u(i) for i, ui in enumerate( u ) ])
    print deg[ i_n ], E
    TE += deg[ i_n ] * E
  dTE = TE - TE_old
  TE_old = TE
  TE = 0.
  print 'dTE', dTE
print 'Total Energy0', TE_old
TE = TotalEnergy( TE_old )
print 'Total Energy', TE
for i_gs, gs in enumerate( Gset ):
  ppl.plot( r, gs )
plt.show()
ppl.plot( r, map( Calc_u, range( r_num ) ) )
plt.show()

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

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

コードは以下のサイトを参考に、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()

ミクロカノニカルとカノニカル

ただの愚痴であるが、位相空間が苦手である。
大抵、ミクロカノニカル分布から入ると思うが、その説明に

が使われる。
その中で一番状態数の多い状態が平衡に対応するとし、変分をかまし、最後にマクスウェル-ボルツマン分布を得る。

正直、これだけよくわからんものをホイホイ導入されて、これをとにかく納得しろという方が無理。そもそもイメージ出来ないし。
あと、これで持って来られるとカノニカル分布との対応が非常にわかりにくい。
てかカノニカル分布と同じノリじゃダメなの?と思う。

カノニカル分布は

  • エネルギーで状態が指定される箱(系)を用意。
  • 箱を複数連結させる。これを集団(アンサンブル)と呼ぶ。
  • 箱の間でエネルギーのやり取り可能。
  • 集団全体でエネルギーは固定。
  • 集団のエネルギー = 箱1 + 箱2 + 箱3 +...
  • 全部の箱のエネルギーが一定とは全く限らないから、箱同士のエネルギーの奪い合いが、今、始まる!

的な感じで、

全体のエネルギーが保存する様に各箱にエネルギーを割り振ったときに、どんな風にエネルギーを割り振るのが一番可能性が高いか?

となる。
「箱の間のエネルギーのやり取り可能」と書きましたが、正直ここは僕はこの表現嫌いです。
「やり取りがあるからエネルギーが割り振れる」という立場からすると、必要と感じますが、やり取りを真面目に考えると箱のエネルギーを独立に足せないし、割り振っている間のエネルギー移動のプロセスは全く気にしないので、個人的には、「箱同士は独立に扱い、それらのエネルギー状態は確率的に発生する」とした方がしっくりくる。下手に物理的に扱うと変なことになる気がする。

ちなみに、全部の箱にエネルギーを均等に割り振る場合の数は「1通り」しかないから、ほとんど発生しない
箱1のエネルギーがちょっと高くて、その代わりに箱2のエネルギーがちょっと低いのは、その逆もあるし、箱2じゃなくて箱3が低くても良いから、こういう方が場合の数が多くなる。そして結果的にエネルギー分布が出来るのがわかる。
それでエネルギーの保存とか粒子数の保存とかの制限が入って、カノニカル分布と言われるものが得られる。

これをミクロカノニカルにそのまま転用すれば、

  • エネルギーで状態が指定される粒子を用意。
  • 粒子を箱の中に入れる。
  • 粒子間でエネルギーのやり取り可能
  • 箱でエネルギーは固定。
  • 箱のエネルギー = 粒子1 + 粒子2 + 粒子3 +...
  • 粒子の箱のエネルギーが一定とは全く限らないから、粒子同士のエネルギーの奪い合いが、今、始まる!

で良くない?
ここでもカノニカルのときと同じで、粒子は独立だけど確率的にエネルギー状態が発現して、箱でエネルギーが保存していなければならないから、全ての粒子のエネルギーが等しいのではなく、マクスウェル-ボルツマン分布に従う形になる。

逆に、粒子が独立じゃなくて箱のエネルギーを各粒子の和で書けないから、いっそ箱単位の扱いに拡張しようというのがカノニカル分布だと理解した。
例えば「クーロン相互作用のある電子の集団」を考えたとき、ハートリー・フォック近似で独立粒子近似をしたとしても、全エネルギーは各粒子の和にはならない。そのまま足すと相互作用を二重に数えてしまう。

何となく、バンド計算とかを齧っていると、
ミクロカノニカル → カノニカル
unit cell → supercell
な対応に近いなぁとも最近思う。
例えばunit cellに原子一個だけとかだと、原子一個解いてあと並進対称性でズルすればよいけど、不純物が入っているとかなんとかでsupercellにすると、その大きいセルの中で真面目に解かないといけなくなる。

一応、言及しておくが、等重率の原理は使っている。場合の数を数えるときにはすべて同じ重みで扱う。
更に言及しておくが、多分、田崎先生の統計力学に似たようなことが書いてあったと思う。
今、手元にないので、割と適当なことを言ってしまっていると思う。
でも、阿部大先生の本の1章読んでて、やっぱり位相空間から入るのは性に合わないと感じた。

条件反射シリーズ:近似:固有値の異なる波動関数の積

グロッソ(中)に載ってた。

 \phi^*_k( x ) \phi_{k'}( x ) = \int dx \phi^*_k( x ) \phi_{k'}( x ) / \int dx

その発想は無かった。
ただし、積分区間は有限でないとヤバいことになるので、問題は選ぶだろう。

井戸型 vs クーロン

今更ながらIPythonデータサイエンスクックブックを購入した。もちろん私費である。

https://images-fe.ssl-images-amazon.com/images/I/51UNyto6saL._SL75_.jpg

せっかくなので一次元のシュレーディンガー方程式でも解こうかと思った。
一階微分方程式に対しては、かくあきさんのサイトで紹介されているローレンツアトラクターのソースが簡潔で分かり易くて非常に参考になった。
Python NumPy SciPy : 1 階常微分方程式の解法 | org-技術

Hartree unitで、

 \phi''(x) = 2(E-V(x))\phi(x)

井戸型ポテンシャルでは、

 V(x) = - 1 ( 0 < x < 1 )
 V(x) = 0 ( otherwise )

クーロンポテンシャルではゼロ割を防ぐために 10^{-7}を分母に足して、

 V(x) = - 1 / ( x + 10^{-7} )

初期条件として、 \phi(0) = 0, \phi'(x) = 1を用いた。理由は特にない。それっぽいだけ。実はいろいろ理由があったが、微妙に長くなるので割愛。軽く言及すると、井戸型は純粋に一次元の問題だが、クーロンポテンシャルは三次元の問題を球対称として動径部分のみの方程式に落とし込んでいて、 \phi(x)= x \psi(x)を求めていることになる。

エネルギーを振って計算。

結果:
f:id:koideforest:20170122004507p:plain
f:id:koideforest:20170122004031p:plain
f:id:koideforest:20170122004036p:plain

波動関数の大きさは規格化していないので置いといて、やはりクーロンポテンシャルの波の引きずり込み具合は素晴らしい。
Thomas-Fermiの遮蔽ポテンシャルも一瞬考えたが、Fermiエネルギーとかの設定がめんどいので今日はこの辺で。

井戸型ポテンシャルの場合のソースも一応曝しておく。

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
import prettyplotlib as ppl

def well( x ):
  a = 1.
  if -a < x < a: v = -1.
  else: v = 0.
  return v

""" Initialization """
q0 = np.zeros(2)
q0[0] = 0.
q0[1] = 1.

def f( q, x, e ):
  p, pdot = q[0], q[1]
  v = well( x )
  pdotdot = - 2 * ( e - v ) * p
  return np.r_[ pdot, pdotdot ]

x = np.linspace( 0., 10., 100 )
for e in np.linspace( 0.2, 1., 5 ):
  q = spi.odeint( f, q0, x, args=(e,) )
  ppl.plot( x[:], q[:,0], 'o-', mew=1, ms=8, mec='w', label='e={:4.2f}'.format(e) )
ppl.plot( x[:], map( well, x[:] ), 'o-', mew=1, ms=8, mec='w', label='v' )
ppl.legend()
plt.ylim( -1.0, 1.0 )
plt.show()

ストークスの定理の最も簡単な例

xy平面内でのみ回転している、大きさ1の回転を考える。
つまり、

 ( \nabla \times \vec{A} )_x = 0
 ( \nabla \times \vec{A} )_y = 0
 ( \nabla \times \vec{A} )_z = 1

ここから元のベクトルを復元すると、

 \vec{A} = \frac{1}{2}( - y, x, 0 )

ただし、微分でしか定義されていないので、一意には決まらずに任意性が残る。
例えば以下のようなものも可能である。

 \vec{A} = ( 0, x, 0 )

今回は上のものを考える。これをgnuplotで図示すると、

f:id:koideforest:20170121163708p:plain

グルグルとベクトルが渦を巻いているのが良くわかる。
原点から遠ざかるにつれて、そのベクトルの大きさは大きくなる。
実は、「最も簡単な回転しているベクトル」を想像したとき、これの全部のベクトルの大きさが1のver.を思いついたが、それは ( \nabla \times \vec{A} )_zが1にならず 1/\sqrt{x^2+y^2}になって原点で発散、外側で減衰していくので、個人的にはちょっとした発見があった。

これを用いてストークスの定理を考察する。
積分範囲を半径1の円に取る。そのため、線積分の方は、単位円の周りにそって一周する形になる。

面積積分の方は、何も考えずに

 \int_S ( \nabla \times \vec{A} )_z dS = \int_S dS = \pi

と求まる。

積分の方は体積素片を系に合わせて考察する必要があり、ここが読んでるだけではピンと来ないところではある。
単位円一周なので、 r=1極座標に対して \theta 0 \rightarrow 2\piまで動かして積分すればよい。積分変数を \thetaに変換するのを忘れないように気を付けると、

 
\displaystyle \oint_C \vec{A} d\vec{r} = \oint_C A_x dx + \oint_C A_y dy + \oint_C A_z dz \\
\displaystyle   = \int^{2\pi}_{0} A_x( \theta ) \frac{ dx }{ d\theta } d\theta
\displaystyle   + \int^{2\pi}_{0} A_y( \theta ) \frac{ dy }{ d\theta } d\theta
\displaystyle   + \int^{2\pi}_{0} A_z( \theta ) \frac{ dz }{ d\theta } d\theta \\
\displaystyle   = \int^{2\pi}_{0} ( - \frac{1}{2}{\rm sin}( \theta ) )( - {\rm sin}( \theta ) ) d\theta
\displaystyle   + \int^{2\pi}_{0} \frac{1}{2}{\rm cos}( \theta ) {\rm cos}( \theta ) d\theta 
\displaystyle   = \frac{1}{2} \int^{2\pi}_{0} d\theta = \pi
\\

という感じで、面積積分と同じになることが確認できた。

ちなみに、rotationを計算したpythonのコードは以下のものを使用。

import numpy as np

def a( x, y ):
  x1 = 0.2 * x
  y1 = 0.2 * y
  ax = -0.5 * y1
  ay = 0.5 * x1
  aa = np.sqrt( ax**2 + ay**2 )
  s = '{:f} {:f} {:f} {:f} {:f}\n'.format( x1, y1, ax, ay, aa )
  return s

f = open( 'simple_rot.dat', 'w')
for x in range(6):
  for y in range(6):
    if np.sqrt( ( 0.2 * x )**2 + ( 0.2 * y )**2 ) > 1. : continue
    f.write( a(  x,  y ) )
    f.write( a( -x,  y ) )
    f.write( a(  x, -y ) )
    f.write( a( -x, -y ) )
f.close()

gnuplotのコマンドは以下の通り。

set terminal png
set output 'simple_rot.png'
set xr[-1.2:1.2]
set yr[-1.2:1.2]
set palette rgbformulae 31,13,10
pl 'simple_rot.dat' u 1:2:3:4:5 w vector lc palette
unset output