nano_exit

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

特殊相対論の計量テンソルと4元ベクトルの内積

計量テンソル g_{ \mu \nu }を以下のように定義。

 \displaystyle
\begin{align}
g_{ \mu \nu } = \left\{
\begin{array}{rl}
 1, & \mu = \nu = 0 \\
 -1, & \mu = \nu = 1,2,3 \\
 0, & \mu \neq \nu
\end{array}
\right.
\end{align}

それで以下の4元ベクトルの内積が、ローレンツ変換で不変であることが特殊相対論の特徴であった。

 \displaystyle
\begin{align}
g_{ \mu \nu } x'{}^\mu x'{}^\nu = g_{ \mu \nu } x^\mu x^\nu
\end{align}

で、これはアインシュタインの規約を用いているから実際には和である。
しかし、頭ではわかっているのだが、心がどうも胡散臭く感じているので、全部展開してみた、というのがこの記事の目的。


\begin{align}
g_{ \mu \nu } x^\mu x^\nu
= g_{ 0 0 } x^0 x^0
 + g_{ 1 0 } x^1 x^0
 + g_{ 2 0 } x^2 x^0
 + g_{ 3 0 } x^3 x^0 \\
 + g_{ 0 1 } x^0 x^1
 + g_{ 1 1 } x^1 x^1
 + g_{ 2 1 } x^2 x^1
 + g_{ 3 1 } x^3 x^1 \\
 + g_{ 0 2 } x^0 x^2
 + g_{ 1 2 } x^1 x^2
 + g_{ 2 2 } x^2 x^2
 + g_{ 3 2 } x^3 x^2 \\
 + g_{ 0 3 } x^0 x^3
 + g_{ 1 3 } x^1 x^3
 + g_{ 2 3 } x^2 x^3
 + g_{ 3 3 } x^3 x^3 \\
= g_{ 0 0 } x^0 x^0 + g_{ 1 1 } x^1 x^1 + g_{ 2 2 } x^2 x^2 + g_{ 3 3 } x^3 x^3 \\
= x^0 x^0 - x^1 x^1 - x^2 x^2 - x^3 x^3 \\
= (ct)^2 - (x)^2 - (y)^2 -(z)^2
\end{align}

まぁ、やる前から当たり前だが、やった後だともっと当たり前に感じることが出来るようになった。
要は、

 \displaystyle
\begin{align}
g_{ \mu \nu } x^\mu x^\nu = g_{ \mu \mu } x^\mu x^\mu
\end{align}

であることが既に明らかであるから、愚直に全部展開することまでする必要はなかった。

一般性を保つならば、記号を分けておいた方が良いことは良いが、一度は(必要とあれば何度でも)アホなことをするのは大事だと思った次第である。

ASE: 構造クラスターを得るスクリプト

半径Rの構造クラスターをASEで作る。
例として単体アルミニウムを使う。

import math
from ase.spacegroup import crystal
from ase.visualize import view

# parameters
LC = 4.05
R = 7.
NSL = math.ceil( 2. * R / LC )

# make bulk structure
aluminium = crystal( 'Al', [ ( 0, 0, 0 ) ], spacegroup=225,
                     cellpar=[ LC, LC, LC, 90, 90, 90 ] )
al_sc = aluminium.repeat( ( NSL, NSL, NSL ) ) # super cell

# cut to make sphere
v_com = al_sc.get_center_of_mass()
al_sc.translate( v_com )

al_sc.wrap() # take outside atoms into unit cell by periodic boundary condition

distances = al_sc.get_all_distances()[0] # distances from the atom at the origin
del al_sc[ [ i for i, dis in enumerate( distances ) if dis > R ] ]

al_sc.translate( - v_com )

view(al_sc)

ここでは原点の原子を中心として構造クラスターを得るようにしている。
unit cell の重心は atoms.get_center_of_mass() ですぐ求まるし、距離もatoms.get_all_distances()で丸ごと求まる。
atoms.wrap()は unit cell の外側にある原子を周期的境界条件で内側に入れ込むというもの。
atoms.translate( xyz ) だけでは、ただ並進移動するだけで、unit cell の外側に原子がはみ出たままになる。
この状態でwrap をしないと、del で不要な原子を除いた時に形がおかしくなる。
最後は一応、元の原点に戻るようにしている。

ASE: 窒素分子のCu(111)の二原子層における吸着構造の最適化

Python module のASE を使って、吸着構造を最適化。
ここでは簡単な有効媒質理論を使うが、calculatorを選ぶことで、AbinitとかVASPとかも使えるらしい。
Introduction: Nitrogen on copper — ASE documentation
元サイトのチュートリアルに少し追加して、吸着した時に結合距離がどう変わるかを見る。
吸着エネルギーを出しているが、ここではあまり関係ない。

from ase import Atoms
from ase.calculators.emt import EMT
from ase.constraints import FixAtoms
from ase.optimize import QuasiNewton
from ase.build import fcc111, add_adsorbate
from ase.visualize import view

h = 1.85 # initial adsorption distance
d = 1.10 # initial bond distance

# make Cu structure
slab = fcc111('Cu', size=(4, 4, 2), vacuum=10.0)
view( slab )

# calculate Cu energy
slab.set_calculator(EMT())
e_slab = slab.get_potential_energy()

# make N2 structure
molecule = Atoms('2N', positions=[(0., 0., 0.), (0., 0., d)])
molecule.set_calculator(EMT())
##e_N2 = molecule.get_potential_energy()

# optimise N2 structure
dyn = QuasiNewton(molecule, trajectory='N2.traj')
dyn.run(fmax=0.05)
print( molecule.get_all_distances() )

# calculate N2 energy
e_N2 = molecule.get_potential_energy()

# make a whole structure
add_adsorbate(slab, molecule, h, 'ontop')
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

# optimize adsorption
constraint = FixAtoms(mask=[a.symbol != 'N' for a in slab])
slab.set_constraint(constraint)
dyn = QuasiNewton(slab, trajectory='N2Cu.traj')
dyn.run(fmax=0.05)

print('Adsorption energy:', e_slab + e_N2 - slab.get_potential_energy())
view( slab )
print( [ dis for dis in slab.get_all_distances()[-2] if dis < 3. ] ) #  distance Cu-N1, N1-N1, N1-N2

add_adsorbateで足された原子群はリストの後ろから入る。
距離を出力する部分は、ちょっとサボって、3オングストローム以下だけ出力するようにした。
本当は、短い距離順にソートして、そこから三つ持ってくるのが良いと思う。

結果:
[1.8499999999999996, 0.0, 0.99809885325797332]
[2.0610952626456882, 0.0, 1.0480140019810147]

1.85は適当に決めた吸着距離だからどうでもよい。
N-Nの距離に注目すると、吸着前は 0.998 だったが、吸着後は 1.048 に伸びることが予言された。
「Cuに近い方の窒素原子の方がより引っ張られるから、結合が伸びる」というイメージを支持するものと思われる。

今はN2分子がCu表面に対して垂直だが、横にした時どうかとか、そういう比較をし始めるときには、Adsorption energyが必要になってくるだろう。

Path operator と Coherent Potential Approximation (CPA) (誤り)

#このページは誤っています。

Path operator を次のように用意しておく。

 \downarrow誤り(一般的なPath operatorの定義ではない。)

\displaystyle
\tau = G + GtG + GtGtG + ...
\\
\displaystyle
= G \left( 1 - tG \right)^{-1} = \left( G^{-1} - t \right)^{-1}
\\
\displaystyle
= t^{-1} \left[ \left( 1 - tG \right)^{-1} - 1 \right]= \left( t^{-1} - G \right)^{-1} - t^{-1}
 \uparrow誤り

 Gはサイト間の移動、 tはサイト上での繰り込まれた相互作用(site T matrix)を表すとする。
後のために tを一応、site potential  vを使って定義しておく。

\displaystyle
t = v + v g v + v + v g v g v + ... = \left( v^{-1} - g \right)^{-1}
 gはサイト内の伝搬を表すと定義するが、後の議論には出てこないのでここでは割愛。

サイトと tが一対一対応すれば良いが、例えば合金のように一つのサイトに対して複数の原子がある確率を持って占有する場合、 tをsite potentialの平均から作る近似(仮想結晶近似:VCA)は粗過ぎる。
何故なら、本来は固定された tのセットを全ての可能な配置に対して求め、それを配置数で割るという平均の取り方をするべきであり、複数種類の tが同時にそのサイトに重みを持って存在している訳ではない。

したがって、配置に対して、すなわちpathの取り方に対して、ある種の平均操作を行うのが妥当だと思われる。
そのような平均操作によって得られたpath operatorを \bar{ \tau }と定義する。
仮に \bar{ \tau }が求めまっていると仮定する。
この \bar{ \tau }を与えるような \bar{ t }を上でpath operatorを用意したように、次の様に定義する。

 \downarrow誤り(一般的なPath operatorの定義に沿っていない)

\displaystyle
\bar{ \tau } = G + G \bar{ t } G + ... = \left( G^{-1} - \bar{ t } \right)^{-1}
\\
\displaystyle
\therefore \; \bar{ t } = G^{-1} - \bar{ \tau }^{-1}
\\
\displaystyle
or
\\
\displaystyle
\bar{ t }^{-1} = \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1}
= G G^{-1} \left( G^{-1} - \bar{ \tau }^{-1} \right)^{-1} \bar{ \tau }^{-1} \bar{ \tau }
= G \left( \bar{ \tau } - G \right)^{-1} \bar{ \tau }
 \uparrow誤り

仮に求まっているとした平均的な場(と呼んで良いかは微妙だが、ここでは有効場と呼ぶことにする)に、注目している相互作用 t_{\alpha}が埋まっているとする。
その場合のpath operator  \bar{ \tau }_{\alpha}は次の様に求められる*(厳密になんでそうなるかは不勉強)。
 \downarrow誤り

\displaystyle
\bar{ \tau }_{\alpha}
= \bar{\tau} \left[ 1 + \left( \bar{ t }^{-1} - t^{-1}_{\alpha} \right) \bar{\tau} \right]
 \uparrow誤り
 \downarrow訂正

\displaystyle
\bar{ \tau }^{nn}_{\alpha}
= \bar{\tau}^{nn} \left[ 1 + \left( t^{-1}_{\alpha} - \bar{ t }^{-1} \right) \bar{\tau}^{nn} \right] ^{-1}
= \left[ \left( \bar{\tau}^{nn} \right)^{-1} + \left( t^{-1}_{\alpha} - \bar{ t }^{-1} \right) \right] ^{-1}
 nはサイト番号を表す。 \tau \neq \tau^{nn} に注意(前者は行列、後者はその要素)
 \uparrow訂正

これの意味は、次の変形を踏まえると分かり易くなる。
 \downarrow誤り(一般的なPath operatorの定義に沿っていない)

\displaystyle
\bar{ t }^{-1} - t^{-1}_{\alpha}
=  \left( \bar{ v }^{-1} - g \right) -  \left( v^{-1}_{\alpha} - g \right)  
= \bar{ v }^{-1} - v^{-1}_{\alpha}
= \frac{ v_{ \alpha } - \bar{ v } }{ v_{ \alpha } \bar{ v } }
= \frac{ \Delta v }{ v_{ \alpha } \bar{ v } }
 \uparrow誤り

つまり、以下の様に、有効場を与える様なsite potential(決してただの平均ではない(と思う))を考慮したいsite potentialで置き換えることで、「埋まる」ことを表現している(なんで v_{\alpha} \bar{ v }で割るのかは不明。しかもこれだとただの一次摂動に思える)。
 \downarrow誤り(一般的なPath operatorの定義に沿っていない)

\displaystyle
\bar{ \tau }_{\alpha}
= \bar{\tau} + \bar{\tau} \frac{ \Delta v }{ v_{ \alpha } \bar{ v } } \bar{\tau}
 \uparrow誤り

有効場に埋まっている各相互作用、ここでは二元系として A, B、に対して、有効path operator \tau_A, \tau_Bが求まったが、これらが有効場を構成しているわけだから、この平均がもともとの有効場に一致すると考えれば、

\displaystyle
x_A \bar{ \tau }_A + x_B \bar{ \tau }_B = \bar{ \tau }
 x_A, x_B A, Bの割合である。

これで変数の閉じた自己無撞着方程式が得られた。
この方法でpath operatorを得る方法をCoherent Potential Approximation (CPA)と呼ぶ。

やっぱり、 \bar{ \tau }_{\alpha}の導出を理解しないと、わかった気になれんなぁ。。。 Path operatorの定義が間違っていたからわからんのは当たり前。

*Ref: V. Popescu et al., J. Synchrotron Rad. 6 (1999) 711.

グランドカノニカル平均からのKeldysh形式に移る前の導入の不満

参考文献:"Nonequilibrium Green's Functions Approach to Inhomogeneous Systems", Karsten Blazer and Michael Bonitz (Springer 2013).
該当箇所:2.1.1 Keldysh Contour

この手のもので自分が最初に詰まってしまうのは、

  • 時間に依存する場合、時間の異なるハミルトニアンの交換を仮定しているのかどうか

が曖昧な点である。

最終的には、一般化された順序演算子を導入して好き勝手に順番を弄れるようにするから何でも良いのだろうが、導入の部分だけを真面目に考えると、適用範囲が狭過ぎて「これ意味あるの??」と思ってしまった。でも導入だからあまり真面目に書かれていない気がする。

時間 tにおけるグランドカノニカル平均を次の様に定義。

\displaystyle
< \hat{A} > (t)
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta ( \hat{H}(t) - \mu \hat{N} ) } \hat{A}_S  \right\}
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ e^{ - \beta \hat{K}( t ) } \hat{A}_S  \right\},
\\
\displaystyle
Z_0 =  {\rm Tr } \left\{ e^{ - \beta \hat{K}( t_0 ) }  \right\}
\\
\displaystyle
\hat{K}( t ) = \hat{H}( t ) - \mu \hat{N}
ハットは演算子であることを表す(ハットの位置がズレて見えるのが不満で仕方がない)。

一般的な時間発展演算子は時間順序演算子 \hat{T}を使って表すと、
{ \displaystyle
\hat{U}( t_1, t_0 ) = \hat{T} {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right)
}

ハミルトニアンが時間に依存するが、 \left[ \hat{H}(t_1), \hat{H}(t_2) \right] = 0 の時(時間毎にエネルギーがwell defined)、
{ \displaystyle
\hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i }{ \hbar } \int^{t_1}_{t_0} dt \; \hat{H}( t ) \right)
}

最後に、ハミルトニアンが時間に依存しない場合、
{ \displaystyle
\hat{U}( t_1, t_0 ) = {\rm exp} \left( - \frac{ i ( t_1 - t_0 ) }{ \hbar } \hat{H}  \right)
}

ここで、 t_0における温度因子を次の様に定義。
{\displaystyle
\hat{U}( t_0 - i \hbar \beta, t_0 ) = {\rm exp} \left( - \beta \hat{K}( t_0 ) \right)
}
ここの取り扱いは \hat{K}( t ) = \hat{K}( t_0 ) として、ハミルトニアン(もしくは \hat{K})が時間に依存しないとした時の時間発展演算子に対応させたものになる。

 \hat{U}の中のハミルトニアン \hat{K}に読み替え、ハミルトニアン \hat{K})が時間に依存しなければ、
{\displaystyle
\hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0, t ) \hat{U}( t_0 - i \hbar \beta, t_0 ) = \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t )
}

これによって、
 
\displaystyle
< \hat{A} > (t) = \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t, t_0 ) \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S  \right\}
\\
\displaystyle
= \frac{ 1 }{ Z_0 } {\rm Tr } \left\{ \hat{U}( t_0 - i \hbar \beta, t_0 ) \hat{U}( t_0, t ) \hat{A}_S  \hat{U}( t, t_0 ) \right\}

これは t_0 \rightarrow t \rightarrow t_0 \rightarrow t_0 - i \hbar \beta という時間発展を考えましょうということを表している。
でもこれ自身は、ハミルトニアンが時間に依存していない時の話だから、こんな仰々しいことしないで松原先生にお願いすることになると思われる。

多分大事なのは、もっと一般のハミルトニアンの時には、直接こうには書けないけれども、「後で順番整理するからとりあえず好きなように並べていいよ!」として同じような時間経路を組んだらどう計算出来るか、という流れなのだと思う。
まあ時間に依存しなきゃ順番なんて何でも良いから、その時に順序演算子は要らなくて、厳密にその形でかけるようになるのは当たり前と言えば当たり前ですが。。。

Formulationは触ったことありますが、実際に運用したことはないので、参考文献を読み進める次第です。

調和振動子と平面波の比較

零点振動している調和振動子と同じエネルギーを持つ平面波とで波動関数を比較する。

\displaystyle
l_0
= \sqrt{ \frac{ \hbar }{ m \omega } }
\\
\displaystyle
\xi
= \frac{ x }{ l_0 }
\\
\displaystyle
\psi_n( x )
= A_n H_n\left( \frac{ x }{ l_0 } \right) e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 }
= A_n H_n( \xi ) e^{ - \frac{ \xi^2 }{ 2 } }
\\
\displaystyle
A_n
= \sqrt{ \frac{ 1 }{ 2^n \; n! } } \; \sqrt{ \frac{ l_0 }{ \sqrt{ \pi } } }
\\
\displaystyle
H_0
= 1
\\
\displaystyle
E_0
= \frac{ \hbar \omega }{ 2 }
最後ら辺、サボりました。

ここで、 E_0 = E_kとして自由平面波のエネルギーを代入して、波数と有効距離 l_0を結びつけると、

\displaystyle
l_0
= \sqrt{ \frac{ \hbar }{ m } \frac{ \hbar }{ 2 E_k } }
= \sqrt{ \frac{ \hbar^2 }{ m } \frac{ 2 m }{ 2 \hbar^2 k^2_0 } }
= \frac{ 1 }{ k_0 }
綺麗に有効距離の逆数になる。

これで求めた k_0を持つ平面波( cos( k_0 x ) )を n=0調和振動子波動関数と比較する(規格化定数は無視)と、全体で一周期分をちょっと超えたかなくらいの波が調和振動子の中に入っていることがわかる。
f:id:koideforest:20180110005333p:plain

もう少し厳密に見てみると、 k_0を逆に固定して、 k_0 x_0 = \piとなる x_0における調和振動子(係数は無視)の値を求めると、

\displaystyle
e^{ - \frac{1}{2} \left( \frac{ x }{ l_0 } \right)^2 }
= e^{ - \frac{1}{2} ( k_0 x )^2 }
\rightarrow e^{ - \frac{1}{2} ( k_0 x_0 )^2 }
= e^{ - \frac{\pi^2}{2} } = 0.00719... \sim 0.01 \ll 1
したがって、十分小さいと見なせるところまでの範囲で、平面波一周期分が入っていることがわかる。

思っていたよりも調和振動子と平面波の対応関係がしっかりしていてビックリした。
原点近傍でポテンシャルがゼロに近くなるから、その近傍で自由解に近づくから当然といえば当然であるが。
自由解だけ波動関数があるような領域は、運動エネルギー(波動関数)ではなくポテンシャル項によってエネルギーが補填されているんだなぁとしみじみ感じた。

PythonのLaguerre陪関数について

Wikipediaには日頃からよくお世話になっているが、自分で確認するのは大事だなと感じた事件。

水素原子の波動関数の可視化とかカッコイイなぁと思い、ネットサーフィン中に以下のサイト様を訪問。
scipyの特殊関数 - 篠突く雨の日記

自分も早速コピーしてやってみて、ふむふむと散布図の便利さに舌鼓を打っていた。
ちょっとアレンジして、波動関数の正負で色を変えるようにしたら気付いてしまった。
「ん?1s波動関数に節があんぞ??」
波動関数の部分はscipy.special.assoc_laguerreで実装されていて、原因はこいつであることがわかった。
何か引数がおかしいのかと思い、Wikipediaで水素原子の解析解におけるLaguerre陪関数を確認。
水素原子におけるシュレーディンガー方程式の解 - Wikipedia

\displaystyle
\rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( 2l + 2 - \rho ) \frac{ d u( \rho ) }{ dx } + ( n - l - 1 ) u( \rho ) = 0 \\
\displaystyle
\rho \frac{ d^{ 2 } u( \rho ) }{ d\rho^{ 2 } } + ( ( 2l + 1 ) + 1 - \rho ) \frac{ d u( \rho ) }{ dx } + ( ( n + l ) - ( 2l + 1 ) ) u( \rho ) = 0 \\
\displaystyle
u( \rho ) = c L^{ 2 l + 1 }_{ n + l }( \rho ) \\
\displaystyle
L^m_k( \rho ) = \frac{ d^m }{ d\rho^m } e^{ \rho } \frac{ d^k }{ d\rho^k } ( e^{ -\rho } \rho^k ) \\
Legendreではlとmは直接打ち込めば良かったけど、Laguerreはそういえばややこしかったんだよなぁ、と遠い記憶が蘇ってきた。
しかし、可視化スクリプトの引数はこれと対応したものが入っている。
実は上付きと下付きは逆の表記になっているとかあるか??と思い、他も参照。

ラゲールの陪多項式 - Wikipedia

\displaystyle
\left( x \frac{ d^{ k + 2 } }{ dx^{ k + 2 } } + ( k + 1 - x ) \frac{ d^{ k + 1 } }{ dx^{ k + 1 } } + ( n - k ) \frac{ d^k }{ dx^k } \right) L^k_n( x ) = 0 \\
\displaystyle
L^k_n( x ) = \frac{ d^{ k } }{ dx^{ k } } L_n( x ) = \frac{ d^{ k } }{ dx^{ k } } \left( e^x \frac{ d^{ n } }{ dx^{ n } } ( x^n e^{ -x } ) \right)
上がkだったり下がkだったりしてめちゃくちゃややこしいが、とりあえずLaguerre多項式(下付き添字のみ)を上付き添字の分だけ微分しまくるのが陪関数であることがわかった。
その意味で、水素原子におけるシュレーディンガー方程式の解 - Wikipediaの上付き下付きの意味はラゲールの陪多項式 - Wikipediaと同じ。

scipy.special.assoc_laguerreのサイトに行ってみた。
scipy.special.assoc_laguerre — SciPy v1.1.0 Reference Guide
見る限り、スクリプトの引数の順番も間違っていなさそうだった。

しかし、

scipy.special.assoc_laguerre( x, n + l, 2 * l + 1 )

はn=1, l=0で -x + 1のプロットを示し、意味不明。

そして英語版のwikipediaに辿り着く。
Laguerre polynomials - Wikipedia

\displaystyle
xy'' + ( \alpha + 1 - x ) y' + \beta y = 0 \\
\displaystyle
L^{ ( \alpha ) }_0( x ) = 1 \\
\displaystyle
L^{ ( \alpha ) }_1( x ) = 1 + \alpha - x
これや、、、これやったんや。。。

つまり、scipy.special.assoc_laguerreを水素様原子の波動関数に使うには、
 
\displaystyle
L^{ \alpha }_{ \beta } = L^{ 2 l + 1 }_{ n - l - 1 }
ということになる。
確かにこれで1s (n=1, l=0)の時 L^1_0 = 1となり、めでたく波動関数の節が無くなる。

お手軽に特殊関数が使える様になったとはいえ、油断は禁物ですな。