2018-08-25

第二量子化と場の演算子と波動関数

量子力学

以前、第二量子化についての記事を書いたが、誤って消してしまったので、再度投稿。

任意の波動関数 $\psi(x)$ を、基底関数 $u( x )$ で展開し、その展開係数を $c$ とする。
これに対応する場の演算子 $\hat{ \psi }(x)$ は、消滅演算子 $\hat{ a }$ を用いて、次のように書ける。

$\displaystyle \psi( x ) = \sum_{ \alpha } u_{ \alpha }( x ) \, c_{ \alpha } \\ \displaystyle \hat{ \psi }( x ) = \sum_{ \alpha } u_{ \alpha }( x ) \, \hat{ a }_{ \alpha }$

したがって、第二量子化は展開係数を量子化したものと言える。

基底関数 $u(x)$ は完備性により次の関係を満たす。

$\displaystyle \sum_{ \alpha } u^*_{ \alpha }( x' ) u_{ \alpha }( x ) = \sum_{ \alpha } < \alpha | x' > < x | \alpha > = \sum_{ \alpha } < x | \alpha >< \alpha | x' > \\ \displaystyle = < x | x' > = \delta( x - x' )$

これを用いて、場の演算子 $\hat{\psi}^{\dagger}(x')$ を真空状態 $|>$ に作用させた波動関数を考えると、
$\displaystyle < x | \hat{ \psi }^{\dagger}( x' ) |> = < x | \sum_{ \alpha } u^*_{ \alpha }( x' ) \, \hat{ a }^{\dagger}_{ \alpha } |> \\ \displaystyle = \sum_{ \alpha } u^*_{ \alpha }( x' ) < x | \alpha > = \sum_{ \alpha } u^*_{ \alpha }( x' ) u_{ \alpha }( x ) = \delta( x - x' )$

したがって、場の演算子 $\hat{\psi}^{\dagger}(x')$ は、本当に位置 $x'$ にだけ粒子の存在確率を与えることがわかる。
波動性は一切無い。

「粒子を足す」という言葉にも色々あって、「エネルギー $\varepsilon$ の粒子を足す」とかだと波動性が発生していることになり、注意が必要だと思う。

2018-08-09

水素様原子波動関数 (python)

プログラミング量子力学

以前、pythonにおける水素様原子のためのLaguerre陪関数について確認した。
koideforest.hatenadiary.com

せっかくなので、pythonで規格化まで含めた水素様原子波動関数のスクリプト（原子単位系）を作った。

from scipy.special import assoc_laguerre
from math import factorial

# Laguerre function
def Lnl( n, l, x ):
    alpha = 2 * l + 1
    beta  = n - l - 1
    return assoc_laguerre( x, beta, alpha )

# radial hydrogen like atomic wavefunction in atomic unit
def Rnl( n, l, r, Z ):
    rho = 2. * r / float( n ) * Z
    norm = Z**1.5 * ( 2. / float( n ) )**1.5 \
        * np.sqrt( float( factorial( n - l - 1 ) ) \
        / float( 2 * n * factorial( n + l ) ) )
    return norm * rho**l * np.exp( -0.5 * rho ) * Lnl( n, l, rho )

以下、式の羅列。
水素原子におけるシュレーディンガー方程式の解 - Wikipedia

$\displaystyle \phi_{nlm}( {\bf r } ) = R_{nl}( \rho ) Y_{lm}( \hat{\bf r } ) \\ \displaystyle \rho = \frac{ 2r }{ n } \frac{ Z }{ a_0 } \\ \displaystyle R( \rho ) = N_R \rho^l u( \rho ) e^{ - \rho / 2 }$

$u( \rho )$ が満たすべき方程式

$\displaystyle \rho \frac{ d^2 }{ d \rho } u( \rho ) + ( ( 2 l + 1 ) + 1 - \rho ) \frac{ d }{ d \rho } u( \rho ) + ( n - l - 1 ) \frac{ d }{ d \rho } u( \rho )$

ちなみにエネルギーは $E_n = - \frac{ 1 }{ 2 n^2 } E_{Hartree} Z^2$

（一般）Laguerre微分方程式
 Laguerre polynomials - Wikipedia
$\displaystyle x \frac{ d^2 }{ d x^2 } L^{\alpha}_{\beta}( x ) + ( \alpha + 1 - x ) \frac{ d }{ d x } L^{\alpha}_{\beta}( x ) + \beta L^{\alpha}_{\beta}( x ) = 0$

したがって、
$\displaystyle u(\rho) = L^{2l+1}_{n-l-1}(\rho) \\ R( \rho ) = N_R \rho^l e^{ - \rho / 2 } L^{2l+1}_{n-l-1}(\rho)$

規格化定数は以下の様になる（らしい）。
$\displaystyle N_R = \left( \frac{ Z }{ a_0 } \right)^{3/2} \left( \frac{ 2 }{ n } \right)^{3/2} \sqrt{ \frac{ ( n - l - 1 ) ! }{ 2n ( n + l ) ! } }$

2018-08-09

X線吸収の多重散乱理論でよく使いそうな関数 (python)

プログラミング量子力学

ほぼ自分用のメモ

from scipy.special import spherical_jn, spherical_yn, sph_harm
from sympy import *
from sympy.physics.wigner import clebsch_gordan, gaunt

# spherical Bessel function
def jl( l, rho ):
    return spherical_jn( l, rho )

# spherical Neumann function
def nl( l, rho ):
    return spherical_yn( l, rho )

# 1st kind of spherical Hankel function
def h1l( l, rho ):
    return spherical_jn( l, rho ) + 1j * spherical_yn( l, rho )

# 1st derivative of spherical Bessel function
def jlp( l, rho ):
    return spherical_jn( l, rho, True )

# 1st derivative of spherical Neumann function
def nlp( l, rho ):
    return spherical_yn( l, rho, True )

# spherical harmonics
def Ylm( l, m, theta, phi ):
    return sph_harm( m, l, phi, theta )

# Clebsh-Gordan coefficient
def CG( j1, m1, j2, m2, j3, m3 ):
    return float( clebsch_gordan( j1, j2, j3, m1, m2, m3 ) )

# Gaunt integral
def GLLL( l1, m1, l2, m2, l3, m3 ):
    temp = complex( gaunt( l1, l2, l3, m1, m2, -m3 ) )
    return (-1)**m3 * temp

# Coefficient for matrix element with multipole expansion
def ck( k, l1, m1, l2, m2 ):
    q = m2 - m1
    temp = GLLL( l1, m1, k, q, l2, m2 )
    return np.sqrt( 4 * np.pi / ( 2 * k + 1 ) ) * temp

2018-07-16

EXAFS振動のフーリエ変換について

固体物理プログラミング

「EXAFS振動をフーリエ変換する」というフレーズはX線吸収分光をやっていると腐る程見かけるが、実際にはArtemisやらLarchやらのソフトで流し込んで、実際に気にするところはパラメータのフィッティングだけということがほとんどな気がする。
ここでは「そもそもフーリエ変換するとはどういうことか？」について触れて行きたいと思う。

EXAFSについては、田渕先生の以下のスライドが非常に分かりやすい。
http://titan.nusr.nagoya-u.ac.jp/Tabuchi/BL5S1/lib/exe/fetch.php?media=nssr-long-201406.pdf

EXAFS振動 $\chi( k )$ を超簡略化して、一つのsin関数（かつ余計な項（位相シフト）を含まない）で与えられるとする。これは最も理想的なEXAFS振動である。
「理想的な」という意味は、フーリエ変換(波数 $k \rightarrow$ 距離 $R$ )すると、原子間距離（求めたい量） $R_0$ のところで（原理的に）厳密に鋭いピークを持つ。（位相シフトなどゴチャゴチャ入ると、 $R_0$ からピーク位置が少しズレる）

$\displaystyle \chi( k ) = \sin{( 2 k R_0 )}$

（"2"が付いているのは距離 $R_0$ を電子が「行って返って来る」ことに依る）

これをいきなりフーリエ変換する前に、通常のやり方（信号解析（時間 $t \rightarrow$ 周波数 $f$ ））を復習して、比較出来るようにする。
周波数 $f_0$ の信号 $s(t)$ がsin関数で次のように与えられているとする。

$\displaystyle s( t ) = \sin{( 2 \pi f_0 t ) } = \sin{( \omega_0 t ) }$

これをpythonを使ってフーリエ変換してみる。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

N = 600
delta_t = 1.0 / 600.0

min_t = 0.0
max_t = N * delta_t
t = np.linspace( min_t, max_t, N )

f0 = 10
s = np.sin( 2. * np.pi * f0 * t )

plt.plot( t, s, "-o" )
plt.savefig( "s.png" )
plt.close()

tilde_s = fft( s )
tilde_s = fftshift( tilde_s )
tilde_t = fftfreq( N, delta_t )
tilde_t = fftshift( tilde_t )

plt.plot( tilde_t, 1.0/N * np.abs( tilde_s ), "-o" )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

plt.plot( tilde_t, 2.0 * 1.0/N * np.abs( tilde_s ), "-o" )
plt.xlim( 0, 50 )
plt.grid()
plt.savefig( tilde_s.png )
plt.close()

f:id:koideforest:20180716015720p:plain
f:id:koideforest:20180716020903p:plain
f:id:koideforest:20180716020921p:plain

$f_0 = 10$ に対応したピークがキチンと得られている。
sin関数のフーリエ変換なので、負の振動数成分と対で現れる（ cos関数も同様）。
規格化は $1/N$ 。 $N$ はメッシュ数（今の場合、 $N = 600$ ）。
そのため、正の振動成分のみをプロットする場合、二倍しないと一見規格化が保たれていないように見える。

フーリエ変換後の $\tilde{t}$ の範囲は $[ - 1 / 2 \Delta_t, \, 1 / 2 \Delta_t$ ]（今の場合、 $\Delta_t = 1/600$ 、[ $-300, 300$ ]）。
これに対し、 $\Delta_{ \tilde{t} } = 1 / \Delta_t N$ （今の場合、 $\Delta_{ \tilde{t} } = 1$ ）。
したがって、幾ら細かく( $\Delta_t \ll 1$ )データを取っても、元データのデータ範囲 $N \Delta_t$ が変わらないとフーリエ変換後のメッシュは細かくならない。

単純にフーリエ変換すると、固有振動数 $f_0$ の位置にピークが得られたわけだが、これはつまり横軸が振動数（ $\tilde{t} = f$ ）に変換されており、角振動数 $\omega$ ではない。
$\omega$ に変換するには、以下の2パターンが考えられる。

1. $\omega_0 / 2 \pi$ の位置にピークが出て来ると解釈する。

$\displaystyle s( t ) = \sin( \omega_0 t ) = \sin \left( 2 \pi \frac{\omega_0}{2\pi} t \right )$
より $f_0 = \omega_0 / 2 \pi$ と思える。

2. 横軸 $\tilde{t}$ に $2\pi$ を掛けて、 $\omega$ に変えてしまう。（ $\omega = 2 \pi f$ ）

さて、ここでEXAFSの話に戻る。 $\chi( k ) = \sin( 2 k R_0 ) = \sin( R_0 ( 2 k ) )$ だったので、この形は $\sin( \omega_0 t )$ に近い。
やりたいことは、「 $\chi( k )$ をフーリエ変換して、 $R_0$ にピークが出る動径分布関数を得る」ことである。なので、

$k$ ではなく、 $2k$ でフーリエ変換する。
$R_0$ は $f_0$ ではなく $\omega_0$ に対応しているので、横軸 $\tilde{2k}$ に $2\pi$ を掛ける。

これを踏まえて、 $\chi(k)$ から動径分布関数を得る。

import numpy as np
from scipy.fftpack import fft, fftfreq, fftshift
from matplotlib import pyplot as plt

min_k = 2
max_k = 12
delta_k = 0.1
k = np.arange( min_k, max_k, delta_k )
N = len( k )

R0 = 3
k2 = 2. * k
delta_k2 = 2. * delta_k
chi = np.sin( R0 * k2 )
plt.plot( k, chi, "-o")
plt.savefig( "chi.png" )
plt.close()

tilde_chi = fft( chi )
tilde_chi = fftshift( tilde_chi )
tilde_k2 = fftfreq( N, delta_k2 )
tilde_k2 = fftshift( tilde_k2 )
R = 2. * np.pi * tilde_k2
plt.plot( R, 2.0 / N * np.abs( tilde_chi ), "-o" )
plt.xlim( 0, 10 )
plt.ylim( 0, 1 )
plt.grid()
plt.savefig( "tilde_chi.png" )
plt.close()

f:id:koideforest:20180716031948p:plain
f:id:koideforest:20180716032005p:plain

$R_0=3$ の付近にピークが得られているのがわかる。しかし、綺麗なピークとは言い難い。
これは $k$ の範囲を広げると改善される。
$k = [ 2, 24$ ]にすれば、
f:id:koideforest:20180716033511p:plain
f:id:koideforest:20180716033527p:plain

しかし、現実的には $k>12$ の実験データを得ることは難しい。
そこで、「窓関数が掛かっている」と見做して、ゼロ値のデータを足したものをフーリエ変換する。
f:id:koideforest:20180716034843p:plain
f:id:koideforest:20180716034922p:plain

ピーク強度が小さくなり、周りに小さなピーク構造が現れたが、その代わりにデータ点が増えて滑らかな関数を表すようになった。
これが、Artemis等のソフトで行われている有限フーリエ変換である。窓関数を掛けるのはデータの一部「のみ」でフーリエ変換するのではなく、追加でゼロデータを足しまくって無理矢理滑らかなピークにすることの宣言でもある。
これによって、一見、十分なデータ数があるように見えてしまい、「fittingパラメータを増やすためには広い $k$ の範囲と十分大きな $N$ が必要」と言われても全然しっくり来なくなってしまっているのである。

実態は、汚いピーク（少ないデータ点）を与えたフーリエ変換であり、それに対してfittingをすると思うと、「もっとたくさん（フーリエ変換後のピーク近傍の）データ点が欲しい！」と思えるのではないだろうか？

2018-05-30

メタン(CH4)の分子軌道とLiebフラットバンド

固体物理プログラミング量子力学

四つ足で御馴染みのメタンの分子軌道が、固体におけるLiebフラットバンドが現れる構造と同じということで、テンションが上がった。

以下の簡単な飛び移りだけを考えたハミルトニアン行列の固有値と固有ベクトルを求める。
対角成分（各原子上でのエネルギー）は簡単のためゼロにしてある。

import numpy as np
import sympy as sm

t = -1.0
h = np.zeros( ( 5, 5 ) )
for i in range( 5 ):
  for j in range( 5 ):
    if i == 0 and j != 0:
      h[ i ][ j ] = t
    elif i !=0 and j == 0:
      h[ i ][ j ] = t

# [ [ 0, t, t, t, t ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ],
#   [ t, 0, 0, 0, 0 ], ]

H = sm.Matrix( h )
H.eigenvects()

#1. -2.0, [  2.0,  1.0, 1.0, 1.0, 1.0 ]
#2.  0.0, [  0, -1.0, 1.0, 0, 0 ]
#3.  0.0, [  0, -1.0, 0, 1.0, 0 ]
#4.  0.0, [  0, -1.0, 0, 0, 1.0 ]
#5.  2.0, [ -2.0,  1.0, 1.0, 1.0, 1.0 ]

今、遷移行列 $t$ が負なので、1.の固有状態において節が無く広がっている（飛び移る）と安定になり、5.の固有状態において節が出来ると不安定になる様になっている。

Liebフラットバンドにおいて重要なのは、ゼロ固有値を持つのにヌルでない2., 3., 4.の状態である。
メタン的に言えば、単に「炭素上に節を持った状態」もしくは「非結合軌道」で終わってしまうところだが、もうちょっと言えば「水素上に波動関数が局在した状態」であり、「飛び移れるのに飛び移らない状態」が生じている。
今の計算は完全に有限系であるが、酸化物のReO3は正にLiebフラットバンドを持つ理想的な三次元構造を持っていて、バンド構造にもめちゃくちゃ平らなバンドが現れている。
でもパッとググったところ、誰も何も言っていないっぽい。まぁ普通にフェルミエネルギーがフラットバンドよりも上で金属だから、フラットバンド由来の物性が多分無いからだと思うが。。。

前に「飛び移りがあるのに飛び移らない状態があるのは普通なのか？」と悩んだことがあったので、自分の疑問が正当化された感じがした次第である。

2018-05-23

立体角積分 - Lebedev quadrature

数学プログラミング

数値計算で立体角積分をしたいとき、安直に二重積分するよりも、Lebedev求積法を使った方が効率が良い。
Lebedev quadrature - Wikipedia
http://people.maths.ox.ac.uk/beentjes/Essays/QuadratureSphere.pdf

効率が良いという意味は、メッシュ数が少なくても精度が出るという意味で、言い換えれば同じ精度を出すのに二重積分よりもずっと速く計算が出来るということである。

求積法は、つまるところ、積分を「上手く精度が出るようにメッシュ間隔等を調整しかつ各点に重み付けをした和」に変換したものである。
Lebedev法の場合には、

$\displaystyle \int sin \theta \, d\theta \, \int d\phi \, f(\theta,\phi) \approx 4 \pi \sum^N_{i=1} \omega_i \, f(\theta_i, \phi_i )$

というように $4 \pi$ は後で掛けることを想定している。

Lebedev求積法のFortranのフリーのソースコードとしては、以下のものが知られている。
http://ccl.net/cca/software/SOURCES/FORTRAN/Lebedev-Laikov-Grids/Lebedev-Laikov.F

subroutineの構成としては、

get_oh( ... )
LD0006( X, Y, Z, W, N )
LD0014( X, Y, Z, W, N )
...
LD5810( X, Y, Z, W, N )

という感じで、使うのはLDXXXX(...)の方である。get_ohは外から使うことは無い。

使い方のポイントを箇条書きすると、

引数の"X, Y, Z, W, N"は全てoutputであるため、こちらから何か値を用意する必要は無い。
ただし、"X, Y, Z, W"は配列サイズを事前に指定する必要がある
例えばLD0006(...)の"0006"はメッシュ数を表していて、これを使いたければ少なくとも6より大きい配列サイズの"X, Y, Z, W"を用意しておく。
"N"はメッシュ数であるためLD0006(...)ではN=6が返って来る。
"X Y Z"がLebedevメッシュであり、 $( \theta, \phi )$ に直したければ、"acos(z)"や"atan2(y,x)"で自分で変換する。
"W"が重み付けで、"sum( W )"はほぼ1になるようになっている。
どのLDXXXX(...)を使うべきかは、計算時間と計算精度との相談になる。

gfortranでコンパイルしてWの和が1になるかチェックする簡単なテストプログラム。

program check_Lebedev
implicit none
double precision, allocatable :: x(:), y(:), z(:), w(:)
! double precision, allocatable :: theta(:), phi(:)
integer :: n, n_lbdv

write( *, * ) 'number of mesh:'
read( *, * ) n
allocate( x(n), y(n), z(n), w(n) )
select case( n )
    case(    6 )
        call LD0006( x, y, z, w, n_lbdv )
    ...
    case( 5810 )
        call LD5810( x, y, z, w, n_lbdv )
    case default
        write( *, * ) ' n /= n_lbdv'
        stop
end select

write( *, * ) sum( w )

! allocate( theta( n_lbdv ), phi( n_lbdv ) )
! theta = acos( z )
! phi = atan2( y, x )

end program check_Lebedev

include "Lebedev-Laikov_CommentOutBy!.F"

リンク先のLebedev-Laikov.Fはコメントアウトが行頭の"c"で行われているため、sedか何かを使って"!"でコメントアウトするようにする必要がある（適当なコンパイルオプションを使えば、そんなことしなくても良いのかもしれないが）。

これで球面調和関数とかを掛けた計算とかが比較的簡単に出来るようになると思われる。

2018-05-08

Møller wave operator における加藤の連鎖律(?)

量子力学数学

N粒子系における二体ポテンシャル $V$ を、粒子ペアに番号を付けて、その間の相互作用の和（ $M=N(N-1)/2$ ）として定義する。

$\displaystyle \begin{align} V &\equiv \sum^{M}_{k=1} v_k = V_i + \bar{ V }_i = V_M = \bar{ V_0 }\\ V_i &\equiv \sum^{ i }_{k=1} v_k \\ \bar{ V }_i &\equiv \sum^{M}_{k= i + 1 } v_k \end{align}$

二体ポテンシャルとしたが、結局は相互作用に何かしら番号が付けられて、全体がその和として与えらえるのであれば、いろんな応用が考えられる（と思う）。
このような相互作用に対するN粒子系のGreen関数を以下のように定義。

$\displaystyle \begin{align} G &\equiv ( z - H_0 - V )^{-1} = G_M \\ G_i &\equiv ( z - H_0 - V_i )^{-1} \end{align}$

ここで次の量を考える。

$\displaystyle \begin{align} 1 + G \bar{ V }_i &= G G^{-1} \left( 1 + G \bar{ V }_i \right) = G \left( G^{-1} + \bar{ V }_i \right) = G \left( z - H_0 - V + \bar{ V }_i \right) \\ &= G \left( z - H_0 - V_i \right) = G \left( G_i \right)^{-1} \\ 1 + G_i v_i &= G_i \left( G_i \right)^{-1} \left( 1 + G_i v_i \right) = G_i \left( \left( G_i \right)^{-1} + v_i \right) = G_i \left( z - H_0 - V_i + v_i \right) \\ &= G_i \left( z - H_0 - V_{i-1} \right) = G_i \left( G_{i-1} \right)^{-1} \end{align}$

したがって、 $G = G_M$ に注意すると、

$\displaystyle \begin{align} G_M \left( G_i \right)^{-1} &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \left( G_i \right)^{-1} \\ &= G_M \left( G_{M-1} \right)^{-1} G_{M-1} \cdots \left( G_{i+1} \right)^{-1} G_{i+1} \left( G_i \right)^{-1} \end{align}$

すなわち、

$\displaystyle \begin{align} 1 + G_M \bar{ V }_i &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_{i+1} v_{i+1} \right) \end{align}$

これを $i = 0$ の場合に適用すれば、

$\displaystyle \begin{align} 1 + G V &= 1 + G \sum^{M}_{k=1} v_k \\ &= \left( 1 + G_M v_M \right) \left( 1 + G_{M-1} v_{M-1} \right) \cdots \left( 1 + G_1 v_1 \right) \end{align}$