2017-01-18

条件反射シリーズ：ガンマ関数

数学固体物理量子力学

$x$ と $e^x$ が一緒に入っている積分：
ガンマ関数に持って行けないか、疑ってみましょう。

ちなみにガンマ関数の定義は、
$\Gamma (n) = \int_0^{\infty} x^{n-1} e^{-x} dx$

$x$ の指数が $n-1$ だったか $n$ だったかややこしいし、 $e$ の肩がプラスだったかマイナスだったか忘れるけど、とにかく

$x$ と $e^x$ の組み合わせ => $\Gamma (n)$

と思っておけば見通しは良いはず。
ガンマ関数を使いそうと思ってからググっても遅くはないはず。

例：Sommerfeld展開で見かけるあいつ

$\displaystyle \int_0^{\infty} \frac{ x^{n-1} }{ e^x + 1 } dx$

ガンマ関数にまとめるために、 $e^x$ を $e^{-x}$ の形に持って行けるように変形します。

$\displaystyle \int_0^{\infty} \frac{ x^{n-1} e^{-x} }{ 1 + e^{-x} } dx$

$e^{-x}$ の形にするために、分数を級数展開します。

$\displaystyle \frac{ 1 }{ 1 - ( - e^{-x} ) } = \sum^{\infty}_{m=1} ( - e^{-x} )^{m-1}$

これにより $e^{-mx}$ とまとまりますが、前に掛かっている $x$ と指数の肩を合わせるために、 $y = mx$ の積分変数の変換を行うと、 $dx = dy / m$ に注意して、

$\displaystyle \sum^{\infty}_{m=1} \frac{ (-1)^{m-1} }{ m^n } \int_0^{\infty} y^{n-1} e^{-y} dy \displaystyle = \sum^{\infty}_{m=1} \frac{ (-1)^{m-1} }{ m^n } \Gamma( n )$

という感じで、ガンマ関数にまとめることが出来ました。めでたし。
それで、前に掛かってる奴は、見るからにゼータ関数 $\zeta(n) = \sum^{\infty}_{m=1} 1 / m^n$ を含んでいるので、この形に持って行くと、

$\displaystyle \sum^{\infty}_{m=1} ( \frac{ 1 }{ m^n } - 2\frac{ 1 }{ (2m)^n } ) \Gamma( n ) \displaystyle = ( 1 - \frac{ 1 }{ 2^{n-1} } ) \zeta( n ) \Gamma( n )$

となり、具体的な $n$ の値に対して関数は求まっているので、計算が出来るということになります。

追記：

$\displaystyle \int_0^{\infty} \frac{ x^{n-1} }{ e^{x} + 1 } dx \displaystyle = ( 1 - \frac{ 1 }{ 2^{n-1} } ) \zeta( n ) \Gamma( n )$
から分かる通り、ゼータ関数とガンマ関数は繋がる。
積分の分母が $+1$ じゃなくて、 $-1$ のとき、分数の級数展開にはマイナスは出て来ないから、もの凄くストレートに

$\displaystyle \int_0^{\infty} \frac{ x^{n-1} }{ e^{x} - 1 } dx \displaystyle = \zeta( n ) \Gamma( n )$

となる。

2016-12-31

周期的境界条件の位相の和の関係

量子力学

前回、Tight-Bindingで使った位相の和の関係についてまとめておく。
koideforest.hatenadiary.com

$\displaystyle \frac{ 1 }{ N } \sum_k e^{ i k R_i } = \delta_{ R_i 0 } \\ \displaystyle \frac{ 1 }{ N } \sum_i e^{ i k R_i } = \delta_{ k 0 }$
これらを証明する。

簡単のため、格子定数 $a$ の一次元格子を考える。
格子点 $N$ 個を含む範囲 $L = N a$ （今の場合は長さと表現できる。2,3次元では箱とよく呼ばれる）を設定し、この $L$ が周期的に繰り返しているとする。ここに割と自分は嵌まった。もともと格子定数 $a$ で既に周期的であるが、わざわざそれよりも大きい入れ物で周期性を考えるのである。この $L$ の中で、格子点は $i = 0, 1, \cdots, N-1$ が含まれるとする。 $0$ と $N$ が同じ値を取るとするのが周期的境界条件である。

この一次元格子上で平面波が満たすべき条件は、
$e^{ i k L } = e^{ i k 0 } = 1$ より、 $k = 2 \pi n / L = 2 \pi n / ( N a )$ （ただし、 $n$ は整数）となり、（結晶）運動量が定まる。 $L \rightarrow \infty \, ( N \rightarrow \infty )$ としない限り、運動量は離散化する。

ここで、この平面波を $L$ 内で和を取る。 $R_i = a n_i$ と出来るから、
$\displaystyle \sum_{ i } e^{ i k R_i } = \sum^{ N-1 }_{ n_i = 0 } e^{ i k a n_i } = ( 1 + e^{ i k a } + e^{ 2 i k a} + \cdots + e^{ ( N - 1 ) i k a } )$

この級数は、
$\displaystyle S = 1 + e^{ i k a } + e^{ 2 i k a} + \cdots + e^{ ( N - 1 ) i k a } \\ \displaystyle e^{ i k a } S = e^{ i k a } + e^{ 2 i k a} + \cdots + e^{ ( N - 1 ) i k a } + e^{ N i k a } \\$

ここで $e^{ N i k a } = e^{ N i ( 2 \pi / N a ) a } = 1$ より、

$\displaystyle S = 1 + e^{ i k a } + e^{ 2 i k a} + \cdots + e^{ ( N - 1 ) i k a } \\ \displaystyle e^{ i k a } S = e^{ i k a } + e^{ 2 i k a} + \cdots + e^{ ( N - 1 ) i k a } + 1 \\$

よって、 $( 1 - e^{ i k a } ) S = 0$ であり、 $k \neq 0$ のとき $S = 0$ 。
$k = 0$ のとき $S = 1 + 1 + \cdots + 1 = N$ となる。
これをまとめると、

$\displaystyle \frac{ 1 }{ N } \sum_{ i } e^{ i k R_i } = \delta_{k0}$

となる。kについても同様に、和の範囲を $N$ に制限すると、
$\displaystyle \sum_{ k } e^{ i k R_i } = \sum^{ N-1 }_{ n_k = 0 } e^{ i ( 2 \pi n_k / L ) R_i } = 1 + e^{ i (2 \pi R_i / L ) } + e^{ 2 i (2 \pi R_i / L ) } + \cdots + e^{ ( N - 1 ) i (2 \pi R_i / L ) } \\$

また、 $n_i$ は整数であるため、
$\displaystyle e^{ N i ( 2 \pi R_i / L ) } = e^{ N i ( 2 \pi n_i a / ( N a ) ) } = e^{ 2 \pi i n_i } = 1$

あとは、 $i$ の和と同様の級数の扱いで求まる。
この導出はよく忘れるので、個人的にやっとまとめられたなぁと感慨深くあります。

2016-12-30

第二量子化でTight-Binding

量子力学

※2019/12/09 修正
位置表示での対角項は、後で行うフーリエ変換（級数）しても対角的なので、省略する。
電子の移動のみをハミルトニアンで扱うと、

$\displaystyle H = \sum_{ i, j } ( t_{ i j } c_i^\dagger c_j + h.c. )$

これを、最近接のみに制限（近似）する。
$< i, j >$ を最近接のペアを表すとすれば、

$\displaystyle H = \sum_{ < i, j > } ( t_{ i j } c^\dagger_i c_j + h.c. )$

ここまでは系に依らない（この項だけを扱うのが正当化されるかどうかはもちろん系によるが、それはまた次元の違う近似の話である）。
ここで一次元格子に制限すると、 $< i, j > = ( i, i + 1), (i, i - 1)$ と具体化される。正直、ここまで具体化しないと全然良くわからん。
注意として、
$\displaystyle \left( \sum_i t_{ i, i+1 } c_i^\dagger c_{i+1} \right)^\dagger = \sum_i t^*_{ i, i+1 } c_{i+1}^\dagger c_i = \sum_i t^*_{ i - 1, i } c_{i}^\dagger c_{i-1} \left( = \sum_i t_{ i, i-1 } c_{i}^\dagger c_{i-1} \right)$
であるため、逆方向への飛び移りはエルミート共役に対応することがわかる。

したがって、 $h.c.$ は逆方向の飛び移りだから、
$\displaystyle H = \sum_i ( t_{ i i+1 } c_i^\dagger c_{i+1} + t_{ i i-1 } c_i^\dagger c_{i-1} )$

これを対角化したいというのが人間の性というものである。
位置基底が固有状態　=> 局在（動かない）
運動量基底が固有状態 => 非局在（動く）
という発想の元、フーリエ変換（級数）すると、

$\displaystyle c_i = \frac{ 1 }{ \sqrt{N} } \sum_{ \vec{k} } c_\vec{k} e^{ i \vec{k} \cdot \vec{R}_i }, \\ \displaystyle c_{ \vec{k} } = \frac{ 1 }{ \sqrt{N} } \sum_i c_i e^{ - i \vec{k} \cdot \vec{R}_i }$

変換の変換が元に戻るか確かめると、

$\displaystyle c_i = \frac{ 1 }{ \sqrt{N} } \sum_\vec{k} c_\vec{k} e^{ i \vec{k} \cdot \vec{R}_i }, \\ \displaystyle = \frac{ 1 }{ N } \sum_{ \vec{k} } \sum_j c_j e^{ i \vec{k} \cdot ( \vec{R}_i - \vec{R}_j ) }, \\ \displaystyle = \sum_j c_j \delta_{ ij }, \\ \displaystyle = c_i$

ここで、

$\displaystyle \frac{ 1 }{ N } \sum_{ \vec{k} } e^{ i \vec{k} \cdot ( \vec{R}_i - \vec{R}_j ) }, \\ \displaystyle = \delta_{ ij }$

を使用。
iの和の場合は、
$\displaystyle \frac{ 1 }{ N } \sum_{ i } e^{ i ( \vec{k} - \vec{k}' ) \cdot \vec{R}_i }, \\ \displaystyle = \delta_{ \vec{k} \vec{k}' }$

細かい話は省略。
この変換により、例えば上の一次元話でi+1からiへの移動の項を考えると、

$\displaystyle \sum_i t_{ i i+1 } c_i^\dagger c_{ i+1 }, \\ \displaystyle = \frac{ 1 }{ N } \sum_i t_{ i i+1 } \sum_{ k, k' } c_k^\dagger c_{k'} e^{ - i k R_{ i } + i k' R_{ i+1 } } = \frac{ 1 }{ N } \sum_i t_{ i i+1 } \sum_{ k, k' } c_k^\dagger c_{k'} e^{ i ( k' - k ) R_{ i } + i k' a }$

$a$ は隣のサイトの距離である。今、一次元の格子だから格子定数と言って良い。
$t_{ij}$ がサイト間の相対座標 $R_i - R_j$ にしか依らないとすると、 $i$ で和が取れるから、

$\displaystyle t \sum_{ k, k' } c_k^\dagger c_{k'} e^{ i k' a } \delta_{ kk' } = t \sum_{ k } c_k^\dagger c_{k} e^{ i k a }$

となり、kに対して対角的になる。つまり、動くものは運動量の固有状態になり、隣のサイトに行った分、位相がずれる。
ただし、無限小変位に対して不変ではないため、運動量 $\vec k$ は離散的であり、結晶運動量と呼ばれる。

結局、
$\displaystyle H = \sum_k ( t ( c_k^\dagger c_k e^{ i k a }+ c_k^\dagger c_k e^{ - i k a } ) ) = \sum_k (2 t {\rm cos}( k a ) c_k^\dagger c_k )$

よって、対角化が完了。
$t$ が絶対座標に依存したら並進対称性がないわけで、それはバンドにならんことがわかります。

2016-12-28

逐次代入法を用いた最もシンプルな例と収束因子：図のテスト

数学

前回、簡単な例で逐次代入法を考察した。koideforest.hatenadiary.com

図の使い方のテストもかねて、その内容を図にしてみた。

収束因子を使わない場合、 $x = 1-x$ は以下のようになる。
f:id:koideforest:20161228135911p:plain
ただの長方形をグルグルまわるだけで一向に収束しないのが分かる。

収束因子 $\alpha = 0.7$ を使うと、
f:id:koideforest:20161228135918p:plain
のように $y = x$ と $y = 1-x$ の交点にドンドン近づいていくのが分かる。

ちなみに、これらはgnuplotで出力しているが、直線の連続はデータ点を読み込ませていて、 $y = x$ と $y = 1-x$ はgnuplotをそのまま使っている。つまり、

pl 'iteration.dat' u 1:2 w l
rep x
rep 1-x

と入力している。
各ファイルの中身は、

iteration.dat:
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000
1.50000 1.50000
1.50000 -0.50000
-0.50000 -0.50000
-0.50000 1.50000

iteration_broyden.dat:
-0.50000 -0.50000
-0.50000 1.50000
0.90000 0.90000
0.90000 0.10000
0.34000 0.34000
0.34000 0.66000
0.56400 0.56400
0.56400 0.43600
0.47440 0.47440
0.47440 0.52560

である。
これらのデータはpythonで計算させて出力させていて、収束因子を使った方は、

#!/usr/bin/env python

def fn( x ):
    return 1.0 - x

f = open( 'iteration_broyden.dat', 'w' )

x     = -0.5
alpha = 0.7
for i in map( float, range(5) ):
    f.write( '{:10.5f} {:10.5f}\n'.format( x, x       ) )
    f.write( '{:10.5f} {:10.5f}\n'.format( x, fn( x ) ) )
    x = ( 1.0 - alpha ) * x + alpha * fn( x )

f.close()

というような感じになっている。

2016-12-28

MKL11.3：逆行列計算ルーチンの罠

プログラミング

自分の今の環境では、MKL11.3を使ってfortran90でプログラムを弄っている。違う環境では、ここで言うことは当てはまらないかも知れない。

逆行列を計算したい。私の願いはそれだけであった。

一般の行列に対して対応可能なsubroutineは、GETRF (LU分解)とGETRI(それを踏まえて逆行列計算)である。頭文字のGはGeneralのGであろう。両方セットで使う。引数は、変換したい行列a(*)、変換するときの情報が入っているipiv (integer)、変換が上手く行ったか教えてくれるinfo (integer)、の三つのみである。aは実数でも倍精度実数でも複素数でも倍精度複素数でも、GETRFやGETRIは内部でそれぞれの型に対応した関数に渡して処理してくれる。
実際、GETRFやGETRIはLAPACK95のモジュールを使用しているのだが、これはInterfaceの集まりであり、計算そのものを実行しているのはF95_PRECISIONの中のCGETRF（単精度複素数用）、ZGETRF（倍精度複素数用）といった、型がキッチリ指定されたやつである。使うときには、LAPACK95のみをuseすると宣言してGETRF, GETRIをただ使えば十分である。

罠からは外れるが、使うときにはlapack95.modがあれば良く、f95_precision.modは無くて大丈夫。
ややこしい設定ファイルは、Windowsだと/Qmklとか/QMkl:sequentialとか、Linuxだと-mklで通ると書いてあるが、僕の環境では通らなかったので、自力でファイルを持って来たりした。使ったファイルだけ列挙すると、
libiomp5md.lib
mkl_core.lib
mkl_intel_ilp64.lib
mkl_intel_thread.lib
mkl_lapack95_ilp64.lib
である。libiomp5md.libは mkl/lib/intel64_win/ の中ではなく、compiler/lib/intel64_win/ に入っているため、探すのに苦労した気がする。
後、ilp64のファイルを使う場合には、コンパイルオプションに /4I8 を追加する必要がある。
本来は、Makefileをキチっと作って使い回すものなんでしょうけれども、Makefileはまだ勉強中であります。

んで、何が罠だったかというと、自分がやりたいのはエルミート行列の逆行列計算なので、エルミート行列用のHETRFとHETRIが使いたかった。頭文字のHはもちろんHermitianのHであろう。これまでは適当にGETRF,GETRIを使っていて、せっかくaがエルミートだから速く計算できる方が良かろうということであるのだが、引数の並びに問題があった。

call GETRF( a, [ipiv], [info] )
call GETRI( a, [ipiv], [info] )

に対して、

call HETRF( a, [uplo], [ipiv], [info] )
call HETRI( a, ipiv, [uplo], [info] )

である。uploは、Bunch-Kaufman diagonal pivoting methodにおいて上三角行列と下三角行列どちらを使うかというもので、'U'か'L'が入る。[...]はOPTIONAL引数で、省略可能である。
この引数の順番の違いで、コンパイルが全然通らなくてだいぶ困った。多分エラーメッセージが良くない。
「この汎用サブルーチン呼び出しと一致する特殊サブルーチンがありません。」
って言われたら、真っ先に設定ファイルが足りてないことを疑うのではないのだろうか？設定ファイルの多さと、GETRF,GETRIの流れから引数を疑うまでになかなか至らなかった。。。
もしMKLで詰まったら、大元のlapack.f90を見てみて下さい（というか、本来はそうするのが当たり前なんだろうけれども。。。）。

2016-12-20

逐次代入法を用いた超シンプルな例と収束因子

数学

$x = 1-x$

を考える。もちろん答えは $x = 0.5$ である。
これを敢えて逐次代入法で解いてみよう。この方法は固体物理や第一原理計算で用いるSCF法と同じノリである。

最初に初期値を設定する。答えに近ければ近い程良いが、とりあえず $x_0 = 0.7$ とおく。これを右辺に代入したものが左辺を与えるとして、

$x_1 = 1 - x_0 = 0.3$
$x_2 = 1 - x_1 = 0.7$

と次々に項が求まる。が、よく見ると全然0.5に近づいておらず、このままだとずっと振動して $x$ が求まらない。

そこで収束因子 $\alpha$ を導入する。 $\alpha$ は新しい $x$ をどれくらい古い $x$ に混ぜるかというもので、例えば $x_1$ は次のようになる。

$x'_0 = 1 - x_0$
$x_1 = \alpha x'_0 + ( 1 - \alpha) x_0$

$\alpha$ が小さい程、古い $x$ の情報が残り、 $x$ の変化が緩やかになる。
具体的に値を入れてその振舞を見てみると、

$\alpha = 0.2$
$x_0 = 0.7$
$x'_0 = 0.3$
$x_1 = 0.2 \cdot 0.3 + 0.8 \cdot 0.7 = 0.06 + 0.56 = 0.62$
$x_2 = 0.2 \cdot 0.38 + 0.8 \cdot 0.62 = 0.072 + 0.496 = 0.568$

というように振動が抑えられ、0.5に近づいていくのが分かる。
では本当にこれを続けて0.5になるかどうかを解析的に解いてみよう。

$x_1 = \alpha ( 1 - x_0) + ( 1 - \alpha) x_0 = x_0 (1 - 2\alpha) + \alpha$
$x_2 = x_1 (1 - 2\alpha) + \alpha = x_0 ( 1 - 2\alpha)^2 + \alpha (1 - 2\alpha) + \alpha$
$\vdots$
$x_n = x_0 ( 1 - 2\alpha)^n + \alpha ( 1 + (1 - 2\alpha) + (1 - 2\alpha)^2 + \cdots + (1 - 2\alpha)^{n-1})$

よって、これを無限に繰り返すと、 $|1 - 2\alpha| < 1$ であれば、

$x_{\infty} = \alpha \frac{1}{1-(1-2\alpha)} = \frac{1}{2}$

となり、初期値にも収束因子にも依存せず答えがちゃんと求まることが示せた。
$|1 - 2\alpha| < 1$ 、つまり $0 < \alpha < 1$ という条件は、元々が新しい $x$ と古い $x$ の割合ということを思い出せば全く自然な範囲である。
ちなみに、最初に $x_0=0.5$ を代入すると一発で求まるし、 $\alpha=0.5$ でも一発で求まるので、問題によって適切な初期値と収束因子を選ぶことが重要であることがわかる。しかし、一般に収束因子に関しては物理的なものではないので類推は困難であり、経験的に適切だと思われる値を入れることが多いと個人的には思う。
最適な収束因子を求める問題とか数値計算の方でありそうな気がするが、不用意に飛び込むと全ての数値計算法に満足出来なくなりそうで怖いのでやめておこう。。。

2016-12-19

初等的な縮退の話

量子力学

サクライの章末問題に載っててフムフムとなった縮退の話。

演算子 $A_1, A_2$ が互いに交換せず( $[ A_1, A_2 ] \neq 0$ )、かつそれぞれがハミルトニアンと同時固有状態を作るとする $[ A_1, H ] = 0, [ A_2, H ] = 0$ 。
このときに、一般に縮退が存在することが証明出来る。

これを真っ向から挑む場合に、必要十分であることを言うのはなかなかしんどい気がする。ホワンとした証明しか思いつかん。
これは背理法を使うと結構スッキリ話を持って行けて気持ちいい。

背理法として、縮退が無いとまず仮定する。あるエネルギー固有状態 $|E_i>$ を $A_1$ の固有値 $a^{(1)}_i$ でラベル付けるか、 $A_2$ の固有値 $a^{(2)}_i$ でラベル付けるか好きな方を選択出来るが、 $[ A_1, A_2 ] \neq 0$ の前提条件により $a^{(1)}_i$ と $a^{(2)}_i$ を一つの状態に同時に指定することは出来ない。しかし、今の仮定では縮退が無いため、

$|a^{(1)}_i, E_i > = |a^{(2)}_i, E_i >$

である。これは次の式で見るように、 $[ A_1, A_2 ] \neq 0$ に抵触している。

$( AB - BA ) |a^{(1)}_i, E_i > = (AB-Ba^{(1)}_i) |a^{(1)}_i, E_i > \\ = (AB-Ba^{(1)}_i)|a^{(2)}_i, E_i > = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(2)}_i, E_i > \\ = (A a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > = (a^{(1)}_i a^{(2)}_i - a^{(2)}_i a^{(1)}_i) |a^{(1)}_i, E_i > \\ = 0$

よって、縮退が無いという仮定が誤りであり、実際には縮退していることが言える。

これの具体例で良さそうだなぁと思ったのが点群で、例えば $C_{3v}$ を考えて、 $A_1 = C_3$ 、 $A_2 = \sigma_v$ と対応付けられて、図を描くとわかるが、 $C_3 \sigma_v = \sigma''_v$ 、 $\sigma_v C_3 = \sigma'_v$ （鏡映面の異なる鏡映操作）となり、交換しない。加えて、系が $C_{3v}$ に属していればハミルトニアンは $C_{3v}$ の対称操作に対して変わらない。よって、エネルギー固有状態は縮退している。
実際、 $C_{3v}$ は二次元の既約表現 $E$ を持ち、二重に縮退する。もちろん、一次元の全対称既約表現 $A_1$ （どの対称操作に対しても不変）しかないような部分空間を用意するならば、話は別で縮退は無い（この部分空間においては $[ C_3, \sigma_v ] = 0$ が成立）。
でもやっぱ点群とかは図が無いと説明が絶望的に伝わらんなぁ。。。でもサクライだと角運動量 $L_z$ と $L_x$ で説明してて、交換しないのがただ数式的な感じがして個人的には縮退のイメージがピンとこなかったんだよなぁ。。。悩ましい。。。