2017-01-22

井戸型 vs クーロン

今更ながらIPythonデータサイエンスクックブックを購入した。もちろん私費である。

せっかくなので一次元のシュレーディンガー方程式でも解こうかと思った。
一階微分方程式に対しては、かくあきさんのサイトで紹介されているローレンツアトラクターのソースが簡潔で分かり易くて非常に参考になった。
Python NumPy SciPy : 1 階常微分方程式の解法 | org-技術

Hartree unitで、

$\phi''(x) = 2(E-V(x))\phi(x)$

井戸型ポテンシャルでは、

$V(x) = - 1 ( 0 < x < 1 )$
$V(x) = 0 ( otherwise )$

クーロンポテンシャルではゼロ割を防ぐために $10^{-7}$ を分母に足して、

$V(x) = - 1 / ( x + 10^{-7} )$

初期条件として、 $\phi(0) = 0, \phi'(x) = 1$ を用いた。~~理由は特にない。それっぽいだけ。~~実はいろいろ理由があったが、微妙に長くなるので割愛。軽く言及すると、井戸型は純粋に一次元の問題だが、クーロンポテンシャルは三次元の問題を球対称として動径部分のみの方程式に落とし込んでいて、 $\phi(x)= x \psi(x)$ を求めていることになる。

エネルギーを振って計算。

結果：
f:id:koideforest:20170122004507p:plain
f:id:koideforest:20170122004031p:plain

波動関数の大きさは規格化していないので置いといて、やはりクーロンポテンシャルの波の引きずり込み具合は素晴らしい。
Thomas-Fermiの遮蔽ポテンシャルも一瞬考えたが、Fermiエネルギーとかの設定がめんどいので今日はこの辺で。

井戸型ポテンシャルの場合のソースも一応曝しておく。

import numpy as np
import scipy.integrate as spi
import matplotlib.pyplot as plt
import prettyplotlib as ppl

def well( x ):
  a = 1.
  if -a < x < a: v = -1.
  else: v = 0.
  return v

""" Initialization """
q0 = np.zeros(2)
q0[0] = 0.
q0[1] = 1.

def f( q, x, e ):
  p, pdot = q[0], q[1]
  v = well( x )
  pdotdot = - 2 * ( e - v ) * p
  return np.r_[ pdot, pdotdot ]

x = np.linspace( 0., 10., 100 )
for e in np.linspace( 0.2, 1., 5 ):
  q = spi.odeint( f, q0, x, args=(e,) )
  ppl.plot( x[:], q[:,0], 'o-', mew=1, ms=8, mec='w', label='e={:4.2f}'.format(e) )
ppl.plot( x[:], map( well, x[:] ), 'o-', mew=1, ms=8, mec='w', label='v' )
ppl.legend()
plt.ylim( -1.0, 1.0 )
plt.show()

2017-01-21

ストークスの定理の最も簡単な例

固体物理数学量子力学

xy平面内でのみ回転している、大きさ1の回転を考える。
つまり、

$( \nabla \times \vec{A} )_x = 0$
$( \nabla \times \vec{A} )_y = 0$
$( \nabla \times \vec{A} )_z = 1$

ここから元のベクトルを復元すると、

$\vec{A} = \frac{1}{2}( - y, x, 0 )$

ただし、微分でしか定義されていないので、一意には決まらずに任意性が残る。
例えば以下のようなものも可能である。

$\vec{A} = ( 0, x, 0 )$

今回は上のものを考える。これをgnuplotで図示すると、

f:id:koideforest:20170121163708p:plain

グルグルとベクトルが渦を巻いているのが良くわかる。
原点から遠ざかるにつれて、そのベクトルの大きさは大きくなる。
実は、「最も簡単な回転しているベクトル」を想像したとき、これの全部のベクトルの大きさが1のver.を思いついたが、それは $( \nabla \times \vec{A} )_z$ が1にならず $1/\sqrt{x^2+y^2}$ になって原点で発散、外側で減衰していくので、個人的にはちょっとした発見があった。

これを用いてストークスの定理を考察する。
積分範囲を半径1の円に取る。そのため、線積分の方は、単位円の周りにそって一周する形になる。

面積積分の方は、何も考えずに

$\int_S ( \nabla \times \vec{A} )_z dS = \int_S dS = \pi$

と求まる。

線積分の方は体積素片を系に合わせて考察する必要があり、ここが読んでるだけではピンと来ないところではある。
単位円一周なので、 $r=1$ の極座標に対して $\theta$ を $0 \rightarrow 2\pi$ まで動かして積分すればよい。積分変数を $\theta$ に変換するのを忘れないように気を付けると、

$\displaystyle \oint_C \vec{A} d\vec{r} = \oint_C A_x dx + \oint_C A_y dy + \oint_C A_z dz \\ \displaystyle = \int^{2\pi}_{0} A_x( \theta ) \frac{ dx }{ d\theta } d\theta \displaystyle + \int^{2\pi}_{0} A_y( \theta ) \frac{ dy }{ d\theta } d\theta \displaystyle + \int^{2\pi}_{0} A_z( \theta ) \frac{ dz }{ d\theta } d\theta \\ \displaystyle = \int^{2\pi}_{0} ( - \frac{1}{2}{\rm sin}( \theta ) )( - {\rm sin}( \theta ) ) d\theta \displaystyle + \int^{2\pi}_{0} \frac{1}{2}{\rm cos}( \theta ) {\rm cos}( \theta ) d\theta \displaystyle = \frac{1}{2} \int^{2\pi}_{0} d\theta = \pi \\$

という感じで、面積積分と同じになることが確認できた。

ちなみに、rotationを計算したpythonのコードは以下のものを使用。

import numpy as np

def a( x, y ):
  x1 = 0.2 * x
  y1 = 0.2 * y
  ax = -0.5 * y1
  ay = 0.5 * x1
  aa = np.sqrt( ax**2 + ay**2 )
  s = '{:f} {:f} {:f} {:f} {:f}\n'.format( x1, y1, ax, ay, aa )
  return s

f = open( 'simple_rot.dat', 'w')
for x in range(6):
  for y in range(6):
    if np.sqrt( ( 0.2 * x )**2 + ( 0.2 * y )**2 ) > 1. : continue
    f.write( a(  x,  y ) )
    f.write( a( -x,  y ) )
    f.write( a(  x, -y ) )
    f.write( a( -x, -y ) )
f.close()

gnuplotのコマンドは以下の通り。

set terminal png
set output 'simple_rot.png'
set xr[-1.2:1.2]
set yr[-1.2:1.2]
set palette rgbformulae 31,13,10
pl 'simple_rot.dat' u 1:2:3:4:5 w vector lc palette
unset output

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を見てみて下さい（というか、本来はそうするのが当たり前なんだろうけれども。。。）。

nano_exit

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