今日、朝起きてからずっと多電子原子波動関数の作成プログラムやっているけど、なかなかうまく行かない。
コードは引き続き下記のサイトを参照。
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()
