nano_exit

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

多重項を数える（２）

プログラミング 量子力学

以前に、軌道内の多重項を求めるスクリプトを書いた。
koideforest.hatenadiary.com

ここではそれらを用いた上で、さらに軌道間で多重項を合成する。

  • 軌道の記号と値を行き来する関数
def spec2l( spec ):
    if spec in { "s", "S" }:
        l = 0
    elif spec in { "p", "P" }:
        l = 1
    elif spec in { "d", "D" }:
        l = 2
    elif spec in { "f", "F" }:
        l = 3
    elif spec in { "g", "G" }:
        l = 4
    elif spec in { "h", "H" }:
        l = 5
    elif spec in { "i", "I" }:
        l = 6
    elif spec in { "k", "K" }:
        l = 7
    elif spec in { "l", "L" }:
        l = 8
    elif spec in { "m", "M" }:
        l = 9
    elif spec in { "n", "N" }:
        l = 10
    elif spec in { "o", "O" }:
        l = 11
    elif spec in { "p", "P" }:
        l = 12
    elif spec in { "q", "Q" }:
        l = 13
    elif spec in { "r", "R" }:
        l = 14
    elif spec in { "s", "S" }:
        l = 15
    elif spec in { "t", "T" }:
        l = 16
    elif spec in { "u", "U" }:
        l = 17
    else:
        l = -1
    return l

def l2spec( l ):
    if l == 0:
        spec = "S"
    elif l == 1:
        spec = "P"
    elif l == 2:
        spec = "D"
    elif l == 3:
        spec = "F"
    elif l == 4:
        spec = "G"
    elif l == 5:
        spec = "H"
    elif l == 6:
        spec = "I"
    elif l == 7:
        spec = "K"
    elif l == 8:
        spec = "L"
    elif l == 9:
        spec = "M"
    elif l == 10:
        spec = "N"
    elif l == 11:
        spec = "O"
    elif l == 12:
        spec = "P"
    elif l == 13:
        spec = "Q"
    elif l == 14:
        spec = "R"
    elif l == 15:
        spec = "S"
    elif l == 16:
        spec = "T"
    elif l == 17:
        spec = "U"
    else:
        spec = "NaN"
    return spec
  • 電子配置から軌道角運動量と電子数を得る関数
def n_l_num( orbital ):
    n = int( orbital[ 0 ] )
    l = spec2l( orbital[ 1 ] )
    num = int( orbital[ 2: ] )
    return [ n, l, num ]

def separate( configuration ):
    orb_eles = configuration.split( " " )
    while "" in orb_eles: orb_eles.remove("")    
    conf = []
    for oe in orb_eles:
        conf.append( n_l_num( oe ) )
    return cons
  • 多重項同士を合成する関数
def spin2_add( s1, s2 ):
    # s = 2S
    ss = []
    for s in range( abs( s1 - s2 ), s1 + s2 + 2, 2 ):
        ss.append( s )
    return ss

def om_add( om1, om2 ):  # orbital momentum addition
    oms = []
    for om in np.arange( abs( om1 - om2 ), om1 + om2 + 1 ):
        oms.append( int( om ) )
    return omg

def orb2sl( orb ):  # mp_sl_orb to mp_sl
    if len( orb ) == 2:
        sl = []
        s1 = orb[0][0] - 1
        s2 = orb[1][0] - 1
        l1 = orb[0][1]
        l2 = orb[1][1]
        for s in spin2_add( s1, s2 ):
            for l in om_add( l1, l2 ):
                sl.append( [ s + 1, l ] )
        return sl
    else:
        return "NaN"

def sl2slj( sl ):  # mp_sl to mp_slj
    # Caution: j = 2 * J + 1
    s = ( sl[0] - 1 ) / 2.
    l = sl[1]
    slj = []
    for ss in np.arange( - s, s + 1, 1. ):
        j = int( 2 * abs( l + ss ) + 1 )
        slj.append([ sl[0], sl[1], j ])
    return slj

def mp_sl_2_mp_slj( mp_sl ):
    mp_slj = []
    for mp in mp_sl:
        mp_slj += sl2slj( mp )
    return mp_slj
  • 多重項を記号で表す関数
def value2symbol( v ):
    return "{:d}{:s}{:d}".format( v[0], l2spec( v[1] ), v[2] )


例として、Ce {}^{3+}イオンの L_{2,3} edgeの終状態について多重項を求めてみる。

# input
configuration = "2p5 4f1 5d1"  # Ce3+ L2,3 edge
orbs = separate( configuration )

# multiplet within each orbital
mp_sl_orb = []
for orb in orbs:
    n   = orb[0]
    l   = orb[1]
    num = orb[2]
    conf = []
    for i in range( num ):
        conf = add_electron( conf, l )
    mp_sl_orb.append( slater_configurations( conf, l, num ) )

# multiplet including all orbitals
mp_sl = []
for mp1 in mp_sl_orb[0]:
    for mp2 in mp_sl_orb[1]:
            mp_sl += orb2sl([ mp1, mp2 ])
if len( mp_sl_orb ) > 2:
    for i in range( len( mp_sl_orb ) )[ 2: ]:
        old = copy.deepcopy( mp_sl )
        mp_sl = []
        for mp1 in old:
            for mp2 in mp_sl_orb[ i ]:
                mp_sl += orb2sl([ mp1, mp2 ])    

# add total angular momentum j (j = 2J + 1)
mp_slj = mp_sl_2_mp_slj( mp_sl )
print( [ value2symbol( mp ) for mp in mp_slj ] )


結果、多重項がこれでもかと出てくる。

['2S2', '2S2', '2P2', '2P4', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '2P2', '2P4', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '2H10', '2H12', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '2H10', '2H12', '2I12', '2I14', '2S2', '2S2', '2P2', '2P4', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '4S4', '4S2', '4S2', '4S4', '4P2', '4P2', '4P4', '4P6', '4D2', '4D4', '4D6', '4D8', '4F4', '4F6', '4F8', '4F10', '4G6', '4G8', '4G10', '4G12', '2P2', '2P4', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '2H10', '2H12', '4P2', '4P2', '4P4', '4P6', '4D2', '4D4', '4D6', '4D8', '4F4', '4F6', '4F8', '4F10', '4G6', '4G8', '4G10', '4G12', '4H8', '4H10', '4H12', '4H14', '2D4', '2D6', '2F6', '2F8', '2G8', '2G10', '2H10', '2H12', '2I12', '2I14', '4D2', '4D4', '4D6', '4D8', '4F4', '4F6', '4F8', '4F10', '4G6', '4G8', '4G10', '4G12', '4H8', '4H10', '4H12', '4H14', '4I10', '4I12', '4I14', '4I16']

これを手でやる気にはなれない。

ここでは多重項のみを求めたが、実際にエネルギーを求めるには、ちゃんと波動関数の組み合わせまで再現する必要がある。
その場合には、おそらく他の方法が必要だが、とりあえずは目的が達成出来たと思う。

量子揺らぎと不確定性原理

量子力学

量子揺らぎと不確定性原理について言及してある記事を下記サイトで見つけた。
第一原理計算入門 密度汎関数法 理解への道

「量子揺らぎ」と聞くと、何だかよくわからないが、要は「状態が混ざる」ということである。

古典的な意味の平均（期待値） < x >および分散 (\Delta x)^2は、確率 pを用いて、
\displaystyle < x > = \sum_i x_i p_i \\ \displaystyle (\Delta x)^2 = < ( x - < x > )^2 > = \sum_i ( x_i - < x > )^2 p_i \\ \displaystyle = \sum_i ( x^2_i - 2 x_i < x > + < x >^2 ) p_i \\ \displaystyle = \left( \sum_i x^2_i p_i \right) - 2 < x > \left( \sum_i x_i p_i \right) + < x >^2 \left( \sum_i p_i \right) \\ \displaystyle = < x^2 > - < x >^2 \\ \displaystyle \left(\because \sum_i p_i = 1 \right)
と書ける。

これを量子力学的に書けば、演算子 \hat{A}の固有値 aとすると、
\displaystyle < A >_{ \varphi } = \sum_i a_i | < a_i | \varphi > |^2 \\ \displaystyle ( \Delta A_{ \varphi } )^2 \equiv < A^2 >_{ \varphi } - < A >^2_{ \varphi }

\varphiが離散的で、規格化が問題なく定義出来る場合には、 <>_{ \varphi } = 1であるため、古典的な場合と同じ方法で分散を導出できる。
しかし、 \varphiが連続量（例えば位置や運動量）の場合、デルタ関数規格化を採用していると <>_{ \varphi } = \inftyとなるため、分散を上記の様に定義するしかない。それでも、分散としての機能を果たしているので、問題は無い。

不確定性原理は、位置演算子 \hat{x}と運動量演算子 \hat{p}_xが交換しないというものであり、ここから位置と運動量は同時に決めることが出来ないということが導かれる。
\displaystyle [\hat{x} , \hat{p}_x ] = \hat{x} \hat{p}_x - \hat{p}_x \hat{x} = i \hbar \neq 0

一般に、異なる演算子が交換する場合には、固有状態を同時に両方の固有値を持つことが出来る。
\displaystyle [\hat{A} , \hat{B} ] = 0 \rightarrow | a, b >
つまり、 aと bを同時に決定出来る。（ただし、aとbが独立という意味では必ずしもないことに注意）
むしろ、同時に指定しなければ、完全に状態を決定したとは言えない。
決定出来るという意味は、分散を見るとより明らかだろう。
\displaystyle < A >_{ a_0, b_0 } = \sum_i a_i | < a_i, b_i | a_0, b_0 > |^2 = \sum_i a_i \delta_{ i 0} = a_0 \\ \displaystyle ( \Delta A_{ a_0, b_0 } )^2 \equiv < A^2 >_{ a_0, b_0 } - < A >^2_{ a_0, b_0 } = a_0^2 - a_0^2 = 0
したがって、分散はゼロになり、 b_0状態を指定しても、 a_0が分散無く求まることがわかる。

一方で、演算子が交換しない場合には、両方の固有値を同時に持つ状態が存在しない。
\displaystyle [\hat{C} , \hat{D} ] = f \neq 0 \rightarrow \hat{C} \hat{D} = \hat{D} \hat{C} + f \\ \displaystyle if \quad | \varphi > = | c, d > \\ \displaystyle \hat{C} \hat{D} | \varphi > = \hat{C} \hat{D} | c, d > = d \hat{C} | c, d > = c d | c, d >, \\ \displaystyle ( \hat{D} \hat{C} + f ) | \varphi > = ( c d + f ) | c, d > \neq c d | c, d > \\ \displaystyle \therefore | \varphi > \neq | c, d >

この時、片方の演算子の固有状態を、もう片方で展開して表そうとすると、複数の状態の混ぜ合わせる必要が出てくる。
\displaystyle | c > = \sum_i | d_i > < d_i | c > = \sum_i c_{d_i} | d_i > \displaystyle \\ c_{d_i} \equiv < d_i | c >
c_{d_i} は展開係数であり、仮に cと dが同時に決まるのであれば、展開係数は一つの状態でしか値を持たないが、演算子が交換しない場合には複数の状態が混ざるように展開係数が定まる。
これが量子揺らぎである。

つまり、固有状態でないもので期待値を取ると、分散はゼロでないということである。（同時固有状態でないと言っても、その状態の固有値で分散を取ればゼロである）
展開係数を用いれば、
\displaystyle < A >_{ \varphi } = \sum_i a_i | c^{\varphi}_{ a_i } |^2 \\ \displaystyle ( \Delta A_{ \varphi } )^2 = \sum_i a^2_i | c^{\varphi}_{ a_i } |^2 - \left( \sum_i a_i | c^{\varphi}_{ a_i } |^2 \right)^2

仮に二状態のみが寄与するとすると、
\displaystyle | c^{\varphi}_{ a_1 } |, | c^{\varphi}_{ a_2 } | < 1 \\ \displaystyle | c^{\varphi}_{ a_1 } |^2 + | c^{\varphi}_{ a_2 } |^2 = 1 \\ \displaystyle ( \Delta A_{ \varphi } )^2 = \left( a^2_1 | c^{\varphi}_{ a_1 } |^2 + a^2_2 | c^{\varphi}_{ a_2 } |^2 \right) - \left( a_1 | c^{\varphi}_{ a_1 } |^2 + a_2 | c^{\varphi}_{ a_2 } |^2 \right)^2 \\ \displaystyle = \sum^2_{i=1} a^2_i ( | c^{\varphi}_{ a_i } |^2 - | c^{\varphi}_{ a_i } |^4 ) - 2 a_1 a_2 | c^{\varphi}_{ a_1 } |^2 | c^{\varphi}_{ a_2 } |^2
このように、複数状態が混ざると分散が有限になり、量子力学的な状態に揺らぎが生じている。

位置と運動量で言えば、 < {\bf r} | {\bf k } > \propto \exp( i {\bf k } \cdot {\bf r} ) であるから、運動量固有状態には無限の数の位置固有状態が含まれていることになり、位置揺らぎが無限となる。

「位置と運動量が同時に定まらない」だけで済めばまだ良いが、実際はもっとこれが広範囲に影響している。
考えてみれば当たり前だが、「位置と運動量が同時に定まらない」＝「位置と運動エネルギーが同時に定まらない」と等価である。
\displaystyle [ \hat{H}_0, \hat{\bf r} ] \neq 0
そして、運動エネルギーと位置が交換しないということは、位置に依存したポテンシャルも同時に定まることがない。
\displaystyle [ \hat{H}_0, V( {\bf r} ) ] \neq 0

これにより、電子間相互作用は互いの電子の位置に依存しているため、多電子ハミルトニアン H_mは一電子ハミルトニアン H_1と交換しないことがわかる。
\displaystyle [ \hat{H}_1, \hat{H}_m ] \neq 0
電子間相互作用によって、多電子固有状態は、複数の電子配置（一電子固有状態の積）が混ぜ合わさったものになる。
混ぜ合わせる前と比べて計算結果が改善するため、混ぜ合わせたことによる効果を配置間相互作用と呼び、配置間相互作用を電子相関と呼ぶことが多い（ただし分野に依る）。
したがって、（最後、議論が粗いが）電子相関は（配置間相互作用無しの状態からの）量子揺らぎと言える。

密度汎関数理論（ DFT）では強配位子場しか計算出来ない？

固体物理 量子力学

学生の時に教授から言われた「DFTでは強配位子場しか計算出来ない」というフレーズをふと思い出した。
その時は何を言われているのかよくわからなかったが、今は「何も工夫しなければまぁそうだろう」と思う。

いくつか鍵となる概念がある。

  • DFTは、Kohn-Sham方程式（一電子方程式）を経由する限り、前提として一電子描像であり、電子相関の弱い時に有効。
  • 電子相関の最も強い系は、原子（束縛状態）である。
  • 電子相関の最も弱い系は、自由電子（連続状態）である。
  • 電子相関とは、（電子間クーロン力ではなく）パウリの排他原理による電子"状態"の避け合いであり、多電子系で発生する。（必ずしも電子位置そのものの避け合いではない）
  • 多重項分裂は、電子相関によって発生する。
  • 配位子場分裂は、結晶場分裂に軌道混成が加わって運動エネルギーが改善されたもので、基本的には結晶場分裂と同様に考えられる。
  • 結晶場分裂は、周囲に配置した点電荷による中心原子における電子軌道のエネルギー分裂であり、一電子状態の時にも発生する。
  • 強配位子場：配位子場（結晶場）分裂 > 多重項分裂
  • 弱配位子場：配位子場（結晶場）分裂 < 多重項分裂

DFTはあくまで「指定したN電子系」においての「基底状態電子密度を再現するような有効一電子軌道（Kohn-Sham軌道）」を教えてくれるだけである。
例として、(MnA {}_6) {}^{2-6n}（A {}^{n-}:アニオン）を考える。Mn {}^{2+}の開殻電子は(3d)^5である。
立方対称の結晶場で一電子軌道が t_{2g}と e_gに割れていて、強配位子場ではスピンが一個だけ残る (t_{2g})^5を、弱配位子場ではスピンが全部生きている (t_{2g})^3 (e_g)^2を、それぞれ基底状態とする。
Mn-Aの距離が近ければ強配位子場、遠ければ弱配位子場的になると考えられる。
相関ポテンシャルを入れずに計算すれば、無限にMn-Aが遠くない限り、 t_{2g}と e_gの分裂は発生し、強配位子場的な電子配置が実現する。
もちろん、無限にMn-Aが遠い、つまりはMnだけの状態で計算すれば、3d軌道は縮退して全部スピンが残る弱配位子場的になる。

問題は、Mn-Aの距離が有限の時に、相関ポテンシャルを入れた有効一電子軌道は「 t_{2g}と e_gが混ざった完全スピン偏極状態」が実現するか？ということである。
もしこれが計算出来るならば、「DFTでも弱配位子場を計算出来る」ということになるだろうが、普通のDFTでは結晶場（配位子場）の分裂の方が大きくて、相関ポテンシャルで全部のスピンを片側に寄せることは出来ないだろう。

その意味で、Uパラメータや、電子状態に拘束条件を課す等、何かしら工夫をしない限り、弱配位子場（の基底状態）は計算出来ないということになるだろう。

三つの散乱波動関数

量子力学

前回、位相シフトと T行列の関係についてまとめた。
koideforest.hatenadiary.com

散乱波動関数は、ざっくり３つに分けられると思うので、それぞれまとめてみた。
結局は、位相シフト \delta_lで全部繋がる。

  • 散乱振幅

\displaystyle \psi_{\bf k}( {\bf r} ) = \frac{ 1 }{ ( 2 \pi )^{3/2} } e^{ i {\bf k} \cdot {\bf r} } + f_k( \theta )\frac{ e^{ i k r } }{ r } \\ \displaystyle \cos \theta = \frac{ {\bf k} \cdot {\bf r} }{ k r } \\ \displaystyle f_k( \theta ) = 4 \pi \sum_{L} f_l( k ) Y_{L}( \hat{\bf r} ) Y^*_{L}( \hat{\bf k} ) \\ \displaystyle f_l( k ) = \frac{ e^{ 2 i \delta_l } - 1 }{ 2 i k }

  • 位相シフト

\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_L i^l e^{ i \delta_l } ( \cos \delta_l \, j_l( kr ) - \sin \delta_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} )

  • Green関数（Lippman-Scwinger方程式解）

\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_L i^l ( j_l( kr ) - 2 i k \, \tilde{t}_l( k ) \, h^{(1)}_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \tilde{t}_l( k ) = 2 \int dr' dr'' \, r'^2 r''^2 \, j_l( kr' ) t_l( r', r'' ) j_l( kr'' ) = - \frac{ e^{ 2 i \delta_l } - 1 }{ 2 i k } = - f_l( k )

散乱理論：位相シフトとT行列

量子力学

外側で値を持たないポテンシャルに対する外側の波動関数は、一般に次のように書ける。
\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_L i^l ( A_l \, j_l( kr ) + B_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad = \sqrt{ \frac{2}{\pi} } \sum_L i^l A_l( j_l( kr ) - \tan\delta_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \\ \displaystyle \left( \frac{ 1 }{ (2\pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r} } = \sqrt{ \frac{2}{\pi} } \sum_L i^l j_l( kr )Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} ) \right)
平面波の球面波展開と比較すれば、動径S-eq.の非正則解の線形結合になっていて、一般解の形になっているのがわかると思う。
-\infty \le \tan\delta_l \le \inftyであるため、この変換は厳密である。
（わざわざマイナスにしたのは、無限遠での漸近形における位相シフトの意味を分かりやすくするためのものだが、ここでは重要ではない。）

一方で、Lippman-Scwinger方程式を用いて解を表すこともできる。ここでは、Hartree原子単位を用いる。
\displaystyle \psi_{\bf k}( {\bf r} ) = \phi_{\bf k}( {\bf r} ) + \int d{\bf r}' \, G_0( {\bf r} - {\bf r}' ) V( {\bf r}' ) \psi_{\bf k}( {\bf r}' ) \\ \displaystyle \quad = \phi_{\bf k}( {\bf r} ) + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) \phi_{\bf k}( {\bf r}'' ) \\ \displaystyle \quad = \frac{ 1 }{(2 \pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r} } + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) \frac{ 1 }{(2 \pi)^{3/2} } e^{ i {\bf k} \cdot {\bf r}'' } \\ \displaystyle \quad = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( {\bf r} ) \\ \displaystyle \qquad \qquad \qquad + \int d{\bf r}' d{\bf r}'' \, G_0( {\bf r} - {\bf r}' ) T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( \hat{\bf r} ) \\ \displaystyle \qquad \qquad - 2 i k \sum_{L'} \int d{\bf r}' d{\bf r}'' \, j_{l'}( kr_< ) Y^*_{L'}( \hat{\bf r}_< ) h^{(1)}_{l'}( kr_> ) Y_{L'}( \hat{\bf r}_> ) \\ \displaystyle \qquad \qquad \qquad \qquad \qquad \qquad \times T( {\bf r}', {\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} )
自由電子Green関数の部分波展開は以前にまとめた。
koideforest.hatenadiary.com

{\bf r}がポテンシャルの範囲の外側の時、 {\bf r} > {\bf r}'が常に成り立つから、球対称ポテンシャルを仮定すると T行列は角運動量に対して対角的になり、
\displaystyle \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( {\bf r} ) \\ \displaystyle \qquad - 2 i k \sum_{L'} \sum_{L''} h^{(1)}_{l'}( kr ) Y_{L'}( \hat{\bf r} ) \int d{\bf r}' d{\bf r}'' \, j_{l'}( kr' ) Y^*_{L'}( \hat{\bf r}' ) \\ \displaystyle \qquad \qquad \qquad \qquad \times t_{l''}( r', r'') Y_{L''}( \hat{\bf r}' ) Y^*_{L''}( \hat{\bf r}'' ) j_l( kr'' ) Y_L( \hat{\bf r}'' ) \right) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad = \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) Y_L( \hat{\bf r} ) \\ \displaystyle \qquad - 2 i k h^{(1)}_{l}( kr ) Y_{L}( \hat{\bf r} ) \int dr' dr'' \, r'^2 r''^2 \, j_l( kr' ) t_{l}( r', r'') j_l( kr'' ) \right) Y^*_L( \hat{\bf k} ) \\ \displaystyle \quad \equiv \sqrt{ \frac{ 2 }{ \pi } } \sum_L i^l \left( j_l( kr ) - i k h^{(1)}_{l}( kr ) \, \tilde{t}_l( k ) \right) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} )

\tilde{t}_l(k)の中に「2倍」を含ませていることに注意。
したがって、 \delta_lと \tilde{t}_l( k )の関係は、 h^{(1)}_l = j_l + i n_lより、
\displaystyle A_l = 1 - i k \tilde{t}_l( k ) \\ \displaystyle - A_l \tan \delta_l = k \tilde{t}_l( k ) \\ \displaystyle \therefore \tan \delta_l = - \frac{ k \tilde{t}_l( k ) }{ A_l } = - \frac{ k \tilde{t}_l( k ) }{ 1 - i k \tilde{t}_l( k ) } = i \frac{ i k \tilde{t}_l( k ) }{ 1 - i k \tilde{t}_l( k ) }

答えをもう知っているのもあって、恣意的ではあるが、
\displaystyle \tan \delta_l = \frac{ \sin \delta_l }{ \cos \delta_l } = -i \frac{ e^{ 2 i \delta_l } - 1 }{ e^{ 2 i \delta_l } + 1 } = i \frac{ 1 - e^{ 2 i \delta_l } }{ 1 + e^{ 2 i \delta_l } } \\ \displaystyle \quad = i \frac{ 1 - e^{ 2 i \delta_l } }{ 2 - ( 1 - e^{ 2 i \delta_l } ) } = i \frac{ ( 1 - e^{ 2 i \delta_l } ) /2 }{ 1 - ( 1 - e^{ 2 i \delta_l } ) / 2 } = i \frac{ i k \tilde{t}_l( k ) }{ 1 - i k \tilde{t}_l( k ) } \\ \displaystyle \therefore i k \tilde{t}_l( k ) = \frac{ 1 - e^{ 2 i \delta_l } }{ 2 }

したがって、
\displaystyle \tilde{t}_l( k ) = - \frac{ e^{ 2 i \delta_l } - 1 }{ 2 i k }

これにより、展開係数 A_lも求まり、
\displaystyle A_l = 1 - i k \tilde{t}_l( k ) = 1 + \frac{ e^{ 2 i \delta_l } - 1 }{ 2 } = \frac{ e^{ 2 i \delta_l } + 1 }{ 2 } = e^{ i \delta } \cos \delta_l \\ \displaystyle \therefore \, \psi_{\bf k}( {\bf r} ) = \sqrt{ \frac{2}{\pi} } \sum_L i^l e^{ i \delta_l } ( \cos \delta_l \, j_l( kr ) - \sin \delta_l \, n_l( kr ) ) Y_L( \hat{\bf r} ) Y^*_L( \hat{\bf k} )

注意として、部分散乱振幅 f_l(k)（これを t_l(k)と書く流儀もある）とは逆符号になっている。
\displaystyle f_l( k ) = \frac{ e^{ 2 i \delta_l } - 1 }{ 2 i k } = - \tilde{t}_l( k )
散乱振幅 - Wikipedia

そのうち、散乱波も加えてまとめたい。

多重項を数える（１）

プログラミング 量子力学

前回、binary表現を使って、電子配置を作ることを試みた。
koideforest.hatenadiary.com

次のステップとして多重項を数え上げるため、電子配置を作るのと同時に角運動量とスピンも一緒にリストで保存するように改良した。

def add_electron( configurations0, l ):
    configurations = []
    if configurations0 == []:
        for m in range( -l, l + 1 ):
            for ms in ( 0, 1 ):
                configurations.append([ binary_configuration( l, m, ms ), ms - 0.5, m ])
    else:
        for conf0 in configurations0:
            m0 = check_highest_m( conf0[ 0 ], l )
            for m in range( -l, l + 1 ):
                if m < m0: continue
                for ms in ( 0, 1 ):
                    if   check_unoccupied( conf0[ 0 ], l, m, ms ) \
                       & check_unpair( conf0[ 0 ], l, m, ms ):
                        configurations.append([ 
                            conf0[ 0 ] | binary_configuration( l, m, ms ),
                            conf0[ 1 ] + ( ms - 0.5 ),
                            conf0[ 2 ] + m
                        ])
    return configurations

角運動量とスピンの組から、多重項を抜き出す方法として、スレーターの方法（正式名称なのかは不明）を用いる。

  1. 角運動量およびスピンが最大の多重項から探索を始める。
  2. 必要な要素が全て揃った場合に多重項が存在すると判定。
  3. 判定後、それらの要素をリストから消す。
  4. リストが空になるまで繰り返す。

スピンが半整数の時と整数の時とあるので、ループの設定に多少苦労した。

def slater_s2l( s2l, l, n_ele ):
    lmax = l * n_ele
    smax = 1/2 * n_ele
    s_range = 1/2 if n_ele % 2 == 1 else 1
    multiplet = []
    for s_inv in range( int( smax + s_range ) ):
        for l_inv in range( lmax + 1 ):
            if not s2l: break
            multiplet_exist = True
            while multiplet_exist == True:
                for mms in np.arange( - ( smax - s_inv ), ( smax - s_inv ) + 1, 1 ):
                    for mm in range( - ( lmax - l_inv ), ( lmax - l_inv ) + 1 ):
                        if not [ int( 2 * mms ), mm ] in s2l:
                            multiplet_exist = False
                if multiplet_exist == True:
                    multiplet.append([ int( 2 * ( smax - s_inv ) ) + 1, ( lmax - l_inv ) ])
                    for mms in np.arange( - ( smax - s_inv ), ( smax - s_inv ) + 1 ):
                        for mm in range( - ( lmax - l_inv ), ( lmax - l_inv ) + 1 ):
                            s2l.remove( [ int( 2 * mms ), mm ] )
    return multiplet, s2l

def slater_configurations( configurations, l, n_ele ):
    s2l = [ [ int( 2 * conf[1] ), conf[2] ] for conf in configurations ]
    multiplet, s2l_rest =  slater_s2l( s2l, l, n_ele )
    if not s2l_rest == []:
        print( 'multiplet still exists in')
    return multiplet

(nd)^3の電子状態に対して、適用してみる。

l = 2
n_ele = 3

conf = []
for i in range( n_ele ):
    conf = add_electron( conf, l )

multiplet = slater_configurations( conf, l, n_ele )
print( multiplet )

結果は、

[[4, 3], [4, 1], [2, 5], [2, 4], [2, 3], [2, 2], [2, 2], [2, 1]]

つまり、 {}^4F, {}^4P, {}^2H, {}^2G, {}^2F, {}^2D, {}^2D, {}^2P であり、正しい結果と一致している。

(nd)^3の時の多重項は、以下のサイトを参照。
LS多重項((nd)^3)の場合) [物理のかぎしっぽ]

後は、異なる軌道毎に多重項を出して、各多重項の合成が出来れば、任意の電子状態に対して多重項を求めることが出来る。

バイナリーと量子力学

プログラミング 量子力学

例えば2p軌道とかのp軌道の場合、一電子状態は以下のどれかに対応する。

(1,u), (1, d), (0, u), (0, d), ( -1, u), ( -1, d)

これを、binary形式、「電子がいる軌道=1」、「電子がいない軌道=0」とすれば、例えば、(1,u)と(1, d)が占有されているとき、

110000

と表すことが出来る。
軌道に電子を足す操作を、バイナリーの操作にそのまま置き換えることが出来るので、電子数と軌道を指定したときの全電子配置を機械的に求めることが出来る。
なので試しにpythonでやってみた。pythonでのbinary演算は以下のサイトを参照した。
Python ビット演算 超入門

以下、関数の説明。

  • binary_configuration：軌道角運動量 l、磁気量子数 m、スピン(up=1 or down=0)から値を作る。
  • check_unoccupied：与えた電子配置bにおいて、（軌道角運動量 lの）磁気量子数 m、スピン(up=1 or down=0)に電子がいなければTrueを返す。
  • check_unpair：与えた電子配置bにおいて、（軌道角運動量 lの）磁気量子数 mにdownスピンを置こうとしたときに、upスピンがいなければTrueを返す。
  • check_highest_m：与えた電子配置bにおいて、占有されている軌道の中で最大の磁気量子数 mを返す。
  • add_electron：与えた電子配置のリストの各要素に対して、電子を一つ加えたときに取り得る電子配置を求めて、新しいリストに格納する。

気をつけたところとしては、

  1. 同じ配置を数えないように、次に足す電子は必ず前の配置の電子の磁気量子数以上の軌道に入れるようにした。
  2. 電子のペアを作る時には、前の電子がdownスピンの時にのみ作るようにした。

以下ソースコード。

def binary_configuration( l, m, ms ):
    return 2**( int( 2 * ( m + l ) + ms ) )

def check_unoccupied( b, l, m, ms ):
    return b & bin_configuration( l, m, ms ) == 0

def check_unpair( b, l, m, ms ):
    if ms == 0:  # down spin
        return b & bin_configuration( l, m, ms + 1 ) == 0 # check up spin
    else:
        return True

def check_highest_m( b, l ):
    highest_m = -l
    for m in range( -l , l + 1 ):
        if b & ( 0b11 << 2 * ( m + l ) ) != 0:
            highest_m = m
    return highest_m

def add_electron( configurations0, l ):
    configurations = []
    if configurations0 == []:
        for m in range( -l, l + 1 ):
            for ms in ( 0, 1 ):
                configurations.append( binary_configuration( l, m, ms ) )
    else:
        for conf0 in configurations0:
            m0 = check_highest_m( conf0, l )
            for m in range( -l, l + 1 ):
                if m < m0: continue
                for ms in ( 0, 1 ):
                    if   check_unoccupied( conf0, l, m, ms ) \
                       & check_unpair( conf0, l, m, ms ):
                        configurations.append( conf0 | binary_configuration( l, m, ms ) )
    return configurations

l = 3  # 0 : s,  1 : p,  2 : d,  3 : f
number_electrons = 7
conf = []
for i in range( number_electrons ):
    conf = add_electron( conf, l )
#print( list( map( bin, conf ) ) )
print( len( conf ) )

試しにGd {}^{3+}の (4f)^7に対してやってみたところ、 {}_{14}C_7=3432がちゃんと返ってきた。
ここから多重項を求められると楽しそうに思う。