nano_exit

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

水素分子イオンの固有値方程式を行列代数で解くことについて(重なり有り)

前回、重なり積分がゼロの時について考察した。
koideforest.hatenadiary.com

今回は、より一般的な s \neq 0の場合について、複数の解法を考える。

解きたい行列方程式は、

\displaystyle
\begin{pmatrix}
\alpha & \beta \\
\beta & \alpha
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
=
\epsilon
\begin{pmatrix}
1 & s \\
s & 1
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
\\
\displaystyle
H {\bf c}
=
\epsilon S {\bf c}

 A {\bf x} = \lambda B {\bf x}固有値方程式を一般化固有値方程式と呼び、 B = Iの場合を標準固有値方程式と呼ぶ。
つまり、 s = 0 \, ( S = I )の時は標準固有値方程式だった訳である。

  • 単純に移項する。

これは以前の記事でも扱った方法。
koideforest.hatenadiary.com
今の問題の場合には、これが一番簡単。

\displaystyle
( H - \epsilon S ) {\bf c}
= O
\\
\displaystyle
\begin{pmatrix}
\alpha - \epsilon & \beta - \epsilon s \\
\beta - \epsilon s & \alpha - \epsilon
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
= O
\\
\displaystyle
\begin{vmatrix}
\alpha - \epsilon & \beta - \epsilon s \\
\beta - \epsilon s & \alpha - \epsilon
\end{vmatrix}
= 0
\\
\displaystyle
(\alpha - \epsilon)^2 - (\beta - \epsilon s)^2 = 0
\\
\displaystyle
\alpha - \epsilon = \pm ( \beta - \epsilon s )
\\
\displaystyle
( - 1 \pm s ) \epsilon = - \alpha \pm \beta
\\
\displaystyle
\therefore
\epsilon = \frac{ \alpha \mp \beta } { 1 \mp s } = \frac{ \alpha \pm \beta } { 1 \pm s }
ただし、未知数 \epsilonが非対角成分にも現れるため、数値計算的なアルゴリズム(つまり機械的な方法)に帰着させ難い。
別な言葉で言えば、標準化させずに解くため、解法の指針が立たず、一般化させ難い。

  •  S^{-1}の利用。

最も安直な標準化の方法。

\displaystyle
S^{-1} H {\bf c} = \epsilon {\bf c}
\\
\displaystyle
\therefore
H_s {\bf c} = \epsilon {\bf c} \quad ( H_s \equiv S^{-1} H )

二次元行列の逆行列は、次のように書ける。

\displaystyle
  M =
  \begin{pmatrix}
  a & b \\
  c & d
  \end{pmatrix}
\\
\displaystyle
M^{-1} =
\frac{1}{ad-bc}
  \begin{pmatrix}
  d & -b \\
  -c & a
  \end{pmatrix}
掃き出し法を使っても良いし、余因子を使っても示せる。

これにより、 S^{-1}は、

\displaystyle
S^{-1} = 
\frac{1}{1-s^2}
\begin{pmatrix}
  1 & -s \\
  -s & 1
\end{pmatrix}

したがって、

\displaystyle
S^{-1} H = 
\frac{ 1 }{ 1 - s^2 }
\begin{pmatrix}
\alpha - s \beta & \beta - s \alpha \\
\beta - s \alpha & \alpha - s \beta
\end{pmatrix}
\\
\displaystyle
\begin{vmatrix}
\frac{\alpha - s \beta}{ 1 - s^2 } - \epsilon & \frac{\beta - s \alpha}{ 1 - s^2 } \\
\frac{\beta - s \alpha}{ 1 - s^2 } & \frac{\alpha - s \beta}{ 1 - s^2 } - \epsilon
\end{vmatrix}
= 0
\\
\displaystyle
\frac{\alpha - s \beta}{ 1 - s^2 } - \epsilon = \pm \left( \frac{\beta - s \alpha}{ 1 - s^2 } \right)
\\
\displaystyle
\epsilon = \frac{\alpha - s \beta}{ 1 - s^2 } \mp \left( \frac{\beta - s \alpha}{ 1 - s^2 } \right)
  = \frac{ ( 1 \pm s ) \alpha \mp ( 1 \pm s ) \beta }{ 1 - s^2 }
\\
\displaystyle
\quad
  = \frac{ \alpha \mp \beta }{ 1 \mp s }
  = \frac{ \alpha \pm \beta }{ 1 \pm s }
しかし、逆行列をまともに扱うため、アルゴリズムは単純だが、行列が大きくなるとシンドイ。

  • コレスキー分解の利用。

一般化固有値方程式 A {\bf x} = \lambda B {\bf x}の話に戻る。
 Bがエルミート行列の時、上三角行列 Vを用いて、 Bを次のように表すことが出来る(コレスキー分解)。

\displaystyle
B = V^{\dagger} V
特に、 Bが実対称行列であれば、

\displaystyle
B = V^{T} V
この時、

\displaystyle
A {\bf x} = \lambda V^{T} V {\bf x}
\\
\displaystyle
\left( V^{T} \right)^{-1} A {\bf x} = \lambda V {\bf x}
\\
\displaystyle
\left( V^{T} \right)^{-1} A V^{-1} V {\bf x} = \lambda V {\bf x}
\\
\displaystyle
A' {\bf x}' = \lambda {\bf x}'

ただし、

\displaystyle
A' \equiv \left( V^{T} \right)^{-1} A V^{-1}
\\
\displaystyle
{\bf x}' \equiv V {\bf x}
と定義した。

一般に、 Sは実対称行列であるため、コレスキー分解が使える。

\displaystyle
V =
\begin{pmatrix}
v_{11} & v_{12} \\
0 & v_{22}
\end{pmatrix}
\\
\displaystyle
S = V^{T} V =
\begin{pmatrix}
v_{11} & 0 \\
v_{12} & v_{22}
\end{pmatrix}
\begin{pmatrix}
v_{11} & v_{12} \\
0 & v_{22}
\end{pmatrix}
\\
\displaystyle
\quad
=
\begin{pmatrix}
v_{11}^2 & v_{11} v_{12} \\
v_{11} v_{12} & v_{12}^2 + v_{22}^2
\end{pmatrix}
=
\begin{pmatrix}
1 & s \\
s & 1
\end{pmatrix}
\\
\displaystyle
\therefore
V = 
\begin{pmatrix}
1 & s \\
0 & \sqrt{ 1 - s^2 }
\end{pmatrix}
\\
\displaystyle
V^{-1} =
\frac{ 1 }{ \sqrt{ 1 - s^2 } }
\begin{pmatrix}
\sqrt{ 1 - s^2 } & -s \\
0 & 1
\end{pmatrix}
\\
\displaystyle
\left( V^T \right)^{-1} =
\frac{ 1 }{ \sqrt{ 1 - s^2 } }
\begin{pmatrix}
  \sqrt{ 1 - s^2 } & 0 \\
  -s & 1
\end{pmatrix}

したがって、

\displaystyle
\left( V^T \right)^{-1} H V^{-1}
=
\frac{1}{1-s^2}
\begin{pmatrix}
  \sqrt{ 1 - s^2 } & 0 \\
  -s & 1
\end{pmatrix}
\begin{pmatrix}
  \alpha & \beta \\
  \beta & \alpha
\end{pmatrix}
\begin{pmatrix}
  \sqrt{ 1 - s^2 } & -s \\
  0 & 1
\end{pmatrix}
\\
\displaystyle
\quad
=
\frac{1}{1-s^2}
\begin{pmatrix}
  \sqrt{ 1 - s^2 } & 0 \\
  -s & 1
\end{pmatrix}
\begin{pmatrix}
  \sqrt{ 1 - s^2 } \alpha & \beta - s \alpha \\
  \sqrt{ 1 - s^2 } \beta & \alpha - s \beta
\end{pmatrix}
\\
\displaystyle
\quad
=
\frac{1}{1-s^2}
\begin{pmatrix}
  ( 1 - s^2 ) \alpha & \sqrt{ 1 - s^2 } ( \beta - s \alpha ) \\
  \sqrt{ 1 - s^2 } ( \beta - s \alpha ) & ( 1 + s^2 )\alpha - 2 s \beta
\end{pmatrix}
\\
\displaystyle
\begin{vmatrix}
  \alpha - \epsilon & \frac{ \beta - s \alpha }{\sqrt{ 1 - s^2 }} \\
  \frac{ \beta - s \alpha }{\sqrt{ 1 - s^2 }} & \frac{ ( 1 + s^2 )\alpha - 2 s \beta}{ 1 - s^2 } - \epsilon
\end{vmatrix}
= 0
\\
\displaystyle
\left( \epsilon - \alpha \right) \left( \epsilon - \frac{ ( 1 + s^2 )\alpha - 2 s \beta}{ 1 - s^2 } \right) - \frac{ ( \beta - s \alpha )^2 }{ 1 - s^2 } = 0
\\
\displaystyle
\epsilon^2 - \frac{ ( 1 - s^2 ) \alpha + ( 1 + s^2 )\alpha - 2 s \beta}{ 1 - s^2 } \epsilon + \frac{ ( 1 + s^2 )\alpha^2 - 2 s \alpha \beta}{ 1 - s^2 } - \frac{ s^2 \alpha -2 s \alpha \beta + \beta^2 }{ 1 - s^2 } = 0
\\
\displaystyle
\epsilon^2 - 2 \frac{ \alpha - s \beta}{ 1 - s^2 } \epsilon + \frac{ \alpha^2 - \beta^2 }{ 1 - s^2 } = 0
\\
\displaystyle
\therefore
\epsilon =  \frac{ \alpha - s \beta}{ 1 - s^2 } \pm \sqrt{ \left( \frac{ \alpha - s \beta}{ 1 - s^2 } \right)^2 - \frac{ \alpha^2 - \beta^2 }{ 1 - s^2 } }
\\
\displaystyle
\quad
  =  \frac{ \alpha - s \beta }{ 1 - s^2 } \pm \frac{ \sqrt{ ( \alpha - s \beta )^2 - ( 1 -s^2 ) ( \alpha^2 - \beta^2 ) } }{ 1 - s^2 }
\\
\displaystyle
\quad
  =  \frac{ \alpha - s \beta }{ 1 - s^2 } \pm \frac{ \sqrt{ s^2 \alpha^2 - 2 s \alpha \beta + \beta^2 } }{ 1 - s^2 }
\\
\displaystyle
\quad
  =  \frac{ \alpha - s \beta }{ 1 - s^2 } \pm \frac{ s \alpha - \beta }{ 1 - s^2 }
\\
\displaystyle
\quad
  =  \frac{ ( 1 \pm s )( \alpha \mp \beta ) }{ 1 - s^2 }
\\
\displaystyle
\quad
  =  \frac{ \alpha  \mp \beta }{ 1 \mp s }
  =  \frac{ \alpha  \pm \beta }{ 1 \pm s }
上三角行列が使えるため、数値計算上は得が多いが、手計算上は複雑になる。
コレスキー分解によるハミルトニアン行列の変形が、何か物理的に分かりやすい解釈を与えるかとホンノ少し期待したが、そんなことはなかった。

今は水素分子イオンなので2x2行列で簡単だったが、そのうち金属水素をtight-bindingで計算した時の S行列の影響について調べてみたい。

水素分子イオンの固有値方程式を行列代数で解くことについて(重なり無し)

行列代数を用いて固有値方程式を解くことを、具体的に考えてみる。

簡単のため、重なり積分 S S=0として考える。

\displaystyle
\begin{pmatrix}
\alpha & \beta \\
\beta & \alpha
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
=
\epsilon
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
\\
\displaystyle
H {\bf c} = \epsilon {\bf c}

ハミルトニアン行列 Hはエルミート行列なので、ユニタリー行列 Uを用いて、対角行列 \Lambdaに変換することが出来る。

\displaystyle
U H U^{\dagger} = \Lambda 
\\
\displaystyle
U^{\dagger} = U^{-1} \quad \left( U^{\dagger} U = U U^{\dagger} = I \right)
\\
\displaystyle
\Lambda = 
\begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix}
 I単位行列である。

したがって、

\displaystyle
H U^{\dagger} U {\bf c} = \epsilon U^{\dagger} U {\bf c}
\\
\displaystyle
U H U^{\dagger} U {\bf c} = \epsilon U {\bf c}
\\
\displaystyle
\Lambda U {\bf c} = \epsilon U {\bf c}
\\
\displaystyle
\Lambda \tilde{\bf c} = \epsilon \tilde{\bf c} \quad \left( \tilde{\bf c} \equiv U {\bf c} \right)

式で単に追っていくと、「ふーん」という感じなのだが、よく見てみると、

\displaystyle
\Lambda \tilde{\bf c} = \epsilon \tilde{\bf c}
\\
\displaystyle
\begin{pmatrix}
\lambda_1 & 0 \\
0 & \lambda_2
\end{pmatrix}
\begin{pmatrix}
\tilde{c}_1 \\ \tilde{c}_2
\end{pmatrix}
=
\begin{pmatrix}
\epsilon & 0 \\
0 & \epsilon
\end{pmatrix}
\begin{pmatrix}
\tilde{c}_1 \\ \tilde{c}_2
\end{pmatrix}
となり、「 \lambda_1 = \lambda_2 = \epsilon?」というよくわからん状態になっている。
この辺り、何が起こっているのかを考えることにする。

固有値については、以前に求めた。
koideforest.hatenadiary.com
今回は、 S=0として考えているから、

\displaystyle
\lambda_{\pm} = \alpha \pm \beta

これに対して、「オリジナル」の固有値方程式を満たす固有ベクトルは、

\displaystyle
{\bf e}_{\pm}
=
\frac{1}{\sqrt{2}}
\begin{pmatrix}
1 \\
\pm 1
\end{pmatrix}
\quad ( c_2 = \pm c_1 )
\\
\displaystyle
\because
H {\bf e}_{\pm}
=
\frac{1}{\sqrt{2}}
\begin{pmatrix}
\alpha & \beta \\
\beta & \alpha
\end{pmatrix}
\begin{pmatrix}
1 \\
\pm 1
\end{pmatrix}
=
\frac{1}{\sqrt{2}}
\begin{pmatrix}
\alpha \pm \beta \\ \beta \pm \alpha
\end{pmatrix}
=
(\alpha \pm \beta) \frac{1}{\sqrt{2}}
\begin{pmatrix}
1 \\ \pm 1
\end{pmatrix}
=
\lambda_{\pm} {\bf e}_{\pm}
 1/\sqrt{2}は規格化因子である。

固有ベクトル {\bf e}が直交規格化されていることから、 Hを対角化するユニタリー行列を以下のように作れる。

\displaystyle
U =
\begin{pmatrix}
{\bf e}_+^{\dagger} \\
{\bf e}_-^{\dagger}
\end{pmatrix}
=
\frac{1}{\sqrt{2}}
\begin{pmatrix}
1 & 1 \\
1 & -1
\end{pmatrix}
\\
\displaystyle
U U^{\dagger}
=
\begin{pmatrix}
{\bf e}_+^{\dagger} \\
{\bf e}_-^{\dagger}
\end{pmatrix}
\begin{pmatrix}
{\bf e}_+ &
{\bf e}_-
\end{pmatrix}
=
\begin{pmatrix}
{\bf e}_+^{\dagger} {\bf e}_+  & {\bf e}_+^{\dagger} {\bf e}_- \\
{\bf e}_-^{\dagger} {\bf e}_+  & {\bf e}_-^{\dagger} {\bf e}_-
\end{pmatrix}
=
\begin{pmatrix}
1  & 0 \\
0  & 1
\end{pmatrix}
\\
\displaystyle
U H U^{\dagger}
=
\frac{1}{2}
\begin{pmatrix}
1 & 1 \\
1  & -1
\end{pmatrix}
\begin{pmatrix}
\alpha & \beta \\
\beta  & \alpha
\end{pmatrix}
\begin{pmatrix}
1 & 1 \\
1  & -1
\end{pmatrix}
=
\begin{pmatrix}
\lambda_+  & 0 \\
0  & \lambda_-
\end{pmatrix}

したがって、

\displaystyle
\Lambda \tilde{\bf c} = \Lambda U {\bf c}
=
\begin{pmatrix}
\lambda_+ & 0 \\
0 & \lambda_-
\end{pmatrix}
\begin{pmatrix}
\frac{c_1 + c_2}{\sqrt{2}}
\\
\frac{c_1 - c_2}{\sqrt{2}}
\end{pmatrix}
=
\epsilon 
\begin{pmatrix}
\frac{c_1 + c_2}{\sqrt{2}}
\\
\frac{c_1 - c_2}{\sqrt{2}}
\end{pmatrix}
この式が成り立つためには、 \tilde{\bf c}の成分のうち唯一つが残って他はゼロになれば良い。
そして、その条件こそが c_2 = \pm c_1であり、これは固有ベクトルが満たす条件と等しい。
この事実は、以下のように Uの作り方から明らかなのだが、普段の問題を解く手続きとしては、固有値固有ベクトルを求めて終わってしまう( Uを作るところは飛ばしてしまう)ので、目にする機会は少ない気がする。


\displaystyle
\tilde{\bf e}_+ = U {\bf e}_+
=
\begin{pmatrix}
{\bf e}_+^{\dagger}  {\bf e}_+ \\
{\bf e}_-^{\dagger} {\bf e}_+
\end{pmatrix}
=
\begin{pmatrix}
1 \\ 0
\end{pmatrix}
\\
\displaystyle
\tilde{\bf e}_- = U {\bf e}_-
=
\begin{pmatrix}
{\bf e}_+^{\dagger}  {\bf e}_- \\
{\bf e}_-^{\dagger} {\bf e}_-
\end{pmatrix}
=
\begin{pmatrix}
0 \\ 1
\end{pmatrix}
また、このことからも明らかなように、固有値方程式を解いて得られる固有ベクトル {\bf e}のことであり、 \tilde{\bf e}ではない。

変換が多いと、どの表示の時の値を求めているのか、よくわからなくなるので、注意が必要に思う。

調和振動子のハミルトニアンを生成消滅演算子で表す。

生成消滅演算子の係数がいつも天下りだったので、自分で求めてみることにした。

調和振動子ハミルトニアン

\displaystyle
H = \frac{p^2}{2m} + \frac{m\omega^2}{2}x^2

ちなみに、古典軌道の調和振動子について復習しておくと、

\displaystyle
m\ddot{x}_c = F(x_c) = - k x_c
\\
\displaystyle
\ddot{x}_c = - \frac{k}{m} x_c = - \omega^2 x_c \quad \left( \omega = \sqrt{\frac{k}{m}}\right)
\\
\displaystyle
F(x_c) = - \frac{ d V(x_c) }{ d x_c } \rightarrow V(x_c) = \frac{1}{2} k x_c^2 = \frac{m\omega^2}{2} x_c^2
復習終了。

ハミルトニアンを平方完成するような演算子 a, a^{\dagger}を作りたい。
見たところ、 a = \alpha p + \beta x とおいて、二乗した時に交差項( xp or  px)が出て来なければ良さそうである。
ただし、ハミルトニアンはエルミートであることが要請されるので、 a^{\dagger} aの形で二乗を表現することにする。
とりあえず試してみると、

\displaystyle
a^{\dagger} a
  = (\alpha^* p + \beta^* x) (\alpha p + \beta x)
  = |\alpha|^2 p^2 + |\beta|^2 x^2 + \alpha \beta^* x p + \alpha^* \beta p x

ハミルトニアンと見比べて、

\displaystyle
  |\alpha|^2 = \frac{1}{2m}
\\
\displaystyle
  |\beta|^2 = \frac{m\omega^2}{2}
だと具合が良さそうである。

交差項を消すには、正準交換関係 [x,p] = i \hbarの利用を考える。
そのためには、

\displaystyle
  \alpha^* = - \alpha
\\
\displaystyle
  \beta^* = \beta
つまり、 \alphaは純虚数 \betaは実数であれば良さそうである。

この条件を課すと、

\displaystyle
  \alpha \beta^* x p + \alpha^* \beta p x = \alpha \beta [x,p] = i \hbar \alpha \beta
と表せる。

ここで、 a^{\dagger}, aが生成消滅演算子として働くように、 [a, a^{\dagger}] = 1の条件を課すと、

\displaystyle
[a, a^{\dagger}] = - 2 i \hbar \alpha \beta = 1
\\
\displaystyle
\therefore
2 \hbar (\Im\alpha) (\Re\beta) = 1
つまり、 \alpha \betaは同符号である。

 2 \hbar (\Im\alpha) (\Re\beta) = 1を一旦忘れて、これまでの情報を整理すると、

\displaystyle
\alpha = \frac{i}{\sqrt{2m}}
\\
\displaystyle
\beta = \sqrt{ \frac{m\omega^2}{2} } = \frac{m\omega}{\sqrt{2m}}
と出来そうだが、 2 \hbar (\Im\alpha) (\Re\beta) = 1を思い出すと、このままでは 2 \hbar (\Im\alpha) (\Re\beta) = \hbar \omega \neq 1である。

そこで、 \alpha, \beta \hbar \omegaで割ったもので再定義すれば上手く行きそうである。
これは、ハミルトニアン全体を \hbar \omegaで括り、その中身を平方完成することに対応する。

\displaystyle
H = \hbar \omega \left( \frac{p^2}{2 \hbar m \omega} + \frac{m^2\omega^2}{2 \hbar m \omega} x^2 \right)
  = \hbar \omega \left( a^{\dagger} a - i \hbar \alpha \beta \right)
\\
\displaystyle
\alpha = \frac{i}{\sqrt{2 \hbar m \omega } }
\\
\displaystyle
\beta = \frac{m\omega}{\sqrt{2 \hbar m \omega }}
\\
\displaystyle
\therefore
\\
\displaystyle
a = \frac{1}{\sqrt{2 \hbar m \omega }}( i p + m \omega x )
\\
\displaystyle
a^{\dagger} = \frac{1}{\sqrt{2 \hbar m \omega }}( - i p + m \omega x )

この時、 - i \hbar \alpha \beta = 1/2となり、めでたく、

\displaystyle
H = \hbar \omega \left( a^{\dagger} a + \frac{1}{2} \right)
と求まる。

水素分子イオンの各種一電子積分(プロット)

前回までに、水素分子イオンの結合・反結合軌道についてまとめた。
koideforest.hatenadiary.com
koideforest.hatenadiary.com

今回は、各値をプロットしてその挙動を確かめる。

f:id:koideforest:20190215181143p:plain
 Jの方が R \uparrowに対する減衰が速いため、より近距離で働くことが分かる。

f:id:koideforest:20190215181237p:plain

  •  \beta = T - 2J + S/Rの各成分

f:id:koideforest:20190215181221p:plain
 Rが遠いところでは、 Tはちょびっと負になっている。
しかし、安定化に大きく貢献しているのは、やはり Jである。
 Jは、重なった成分が各原子核からクーロン力で引っ張られている寄与なので、イメージし易い。
 S/Rは相変わらずどう解釈するべきかわからない。

  • 結合・反結合エネルギー

f:id:koideforest:20190215181254p:plain
ちゃんとプロットしてみると、結合エネルギーは常に反結合エネルギーよりも安定である。
approxは、 E^{\pm}_{approx} = \alpha \pm \betaと近似したもので、 Sを無視した形になっている。
基底関数的に言えば、 Sは基底が過完備だった時に対する補正のようなものである。
そのため、何か基底関数系を適当に作った時に補正が如何に重要であるかを端的に見れるのは、教育的に思う。
ちなみに、 E^{\pm}_{approx}の大小関係が逆転する点は、 \beta \ge 0の点であり、これは Sを無視した時の固有値方程式が

\displaystyle
\begin{pmatrix}
\alpha -E & \beta \\
\beta & \alpha -E
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
=O
であることに因る。

水素分子イオンの結合・反結合軌道

前回、水素分子イオンにおけるハミルトニアンの行列要素を求めた。
koideforest.hatenadiary.com

今回は、それらが求まっているとして、結合・反結合軌道がどのようなロジックで得られるかを紹介する。

各水素原子波動関数を基底関数として、その線形結合で水素分子イオンの波動関数を作るとする。

\displaystyle
\Psi = c_1 \psi_1 + c_2 \psi_2

この時に満たされるべきシュレーディンガー方程式は、

\displaystyle
H \Psi = E \Psi

決定したい係数が c_1, c_2と二つあるため、 c_1, c_2に関する方程式が二本必要。
そのために、左から \psi^*_1もしくは \psi^*_2を掛けて、位置変数で積分する。

\displaystyle
c_1 \alpha + c_2 \beta = E (c_1 + c_2 S)
\\
\displaystyle
c_1 \beta + c_2 \alpha = E (c_1S + c_2)

これを行列で表すと、

\displaystyle
\begin{pmatrix}
\alpha & \beta \\
\beta & \alpha
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
=
\begin{pmatrix}
E & ES \\
ES & E
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
\\
\displaystyle
\therefore
\begin{pmatrix}
\alpha - E & \beta - ES \\
\beta - ES & \alpha - E
\end{pmatrix}
\begin{pmatrix}
c_1 \\ c_2
\end{pmatrix}
=
O

 c_1, c_2が共にゼロでない解が得られる条件は、

\displaystyle
\begin{vmatrix}
\alpha - E & \beta - ES \\
\beta - ES & \alpha - E
\end{vmatrix}
=
0
\\
\displaystyle
(\alpha - E)^2 - (\beta - ES)^2 = 0
\\
\displaystyle
\alpha - E = \pm (\beta - ES)
\\
\displaystyle
\therefore
E^{\pm} = \frac{ \alpha \pm \beta }{ 1 \pm S}

 E = E^+の時、

\displaystyle
c_1 \left( \alpha - \frac{ \alpha + \beta }{ 1 + S} \right) + c_2 \left( \beta - \frac{ \alpha + \beta }{ 1 + S} S \right) = 0
\\
\displaystyle
c_1 \left( \frac{ \alpha S - \beta }{ 1 + S} \right) + c_2 \left( \frac{ - \alpha S + \beta }{ 1 + S} \right) = 0
\\
\displaystyle
c_1 \left( \alpha S - \beta \right) - c_2 \left( \alpha S - \beta \right) = 0
\\
\displaystyle
\therefore
c_2 = c_1

 E = E^-も同様にして、 c_2 = - c_1が示せる。

したがって、結合軌道 \Psi^+と反結合軌道 \Psi^-が得られた。

\displaystyle
\Psi^+ = c \left( \psi_1 + \psi_2 \right)
\\
\displaystyle
\Psi^- = c \left( \psi_1 - \psi_2 \right)

この時点では、エネルギーの評価をしていないので、どちらが安定化はわからないが、プロットしてみると結合エネルギーが常に反結合エネルギーよりも安定になる。

水素分子イオンの各種一電子積分

水素分子イオンの重なり積分、クーロン積分、交換積分を求める。

参考にしたPDF
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=7&cad=rja&uact=8&ved=2ahUKEwilj8WMrbvgAhVPBGMBHYacDWMQFjAGegQIAxAC&url=http%3A%2F%2Fwww.geocities.jp%2Fribake2006%2Fshiken-kankei%2Fshikepuri%2Fkozokagaku_komaba%2Fyasmin-2.pdf&usg=AOvVaw0Rmljxrb9VmhorUrBVIegl
https://hb3.seikyou.ne.jp/home/E-Yama//HuckelMOCalculation.pdf

 \displaystyle
\psi(r,\theta,\phi) = \frac{1}{\sqrt{\pi}} e^{-r}

  • 規格化の確認


\displaystyle
\int dr d\theta d\phi \, r^2 \sin\theta \, |\psi(r,\theta,\phi)|^2 
  = \frac{ 1 }{\pi} 4 \pi \int^{\infty}_0 dr \, r^2 e^{-2r}
  = \int^{\infty}_0 dr \, (2r)^2 e^{-2r}
\\
\displaystyle
\quad
  = \frac{1}{2} \int^{\infty}_0 dx \, x^2 e^{-x} \, (x = 2r)
\\
\displaystyle
\quad
  = \frac{1}{2} \left( \left. \frac{x^2 e^x}{-1} \right|^{\infty}_0 - \frac{2}{-1} \int^{\infty}_0 dx \, x e^{-x} \right)
\\
\displaystyle
\quad
  = \frac{1}{2} \left( 2 \left( \left. \frac{ x e^{-x} }{-1} \right|^{\infty}_0 - \frac{1}{-1} \int^{\infty}_0 dx \, e^{-x} \right) \right)
\\
\displaystyle
\quad
  = \left. \frac{ e^{-x} }{ -1 } \right|^{\infty}_0 = 1

  • 楕円座標系

原点を中心とし、z軸上に R離れた2点を固定する。
空間上のある一点を指定する際、その2点からの距離をそれぞれ r_1, r_2とすると、

\displaystyle
(x,y,z) \rightarrow (\xi, \eta, \phi)
\\
\displaystyle
\xi = \frac{r_1 + r_2}{R} \quad (1 \le \xi \le \infty)
\\
\displaystyle
\eta = \frac{r_1 - r_2}{R} \quad (-1 \le \eta \le 1)
\\
\displaystyle
dxdydz = \frac{R^3}{2^3} (\xi^2 - \eta^2) d\xi d\eta d\phi

水素分子イオンのハミルトニアンは以下のように書ける。

\displaystyle
H = H_1 - \frac{1}{r_2} + V_n = H_2 - \frac{1}{r_1} + V_n
\\
\displaystyle
H_1 = \frac{\nabla^2}{2} - \frac{1}{r_1}
\\
\displaystyle
H_2 = \frac{\nabla^2}{2} - \frac{1}{r_2}
\\
\displaystyle
V_n = \frac{1}{R}


\displaystyle
S \equiv \int dxdydz \, \psi^*_1(x,y,z) \psi_2(x,y,z)
\\
\displaystyle
\quad
  = \int^{\infty}_1 d\xi \int^{1}_{-1} d\eta \int^{2\pi}_0 d\phi \, \frac{R^3}{2^3} (\xi^2 - \eta^2)\frac{ e^{- r_1 - r_2} }{\pi}
\\
\displaystyle
\quad
  = \frac{R^3}{2^3}  \frac{1}{\pi} 2\pi \int^{\infty}_1 d\xi \int^{1}_{-1} d\eta\, (\xi^2 - \eta^2)e^{- R\xi}
\\
\displaystyle
\quad
  = 2 \frac{R^3}{2^3} \int^{\infty}_1 d\xi \, (2\xi^2 - \frac{2}{3})e^{- R\xi}
\\
\displaystyle
\quad
  = 2^2 \frac{R^3}{2^3} \left( \left. (\xi^2 - \frac{1}{3})\frac{e^{- R\xi}}{-R} \right|^{\infty}_1 - \frac{2}{-R} \int^{\infty}_1 d\xi \, \xi e^{- R\xi} \right)
\\
\displaystyle
\quad
  = 2^2 \frac{R^3}{2^3} \left( - ( 1 - \frac{1}{3})\frac{e^{- R}}{-R} - \frac{2}{-R} \left( \left. \xi \frac{e^{- R \xi}}{-R} \right|^{\infty}_1 - \frac{1}{-R} \int d\xi \, e^{- R\xi} \right) \right)
\\
\displaystyle
\quad
  = 2^2 \frac{R^3}{2^3} \left( \frac{2}{3}\frac{e^{- R}}{R} + \frac{2}{R} \left( - \frac{e^{- R} }{-R} + \frac{1}{R} \left( -\frac{e^{- R } }{-R}\right) \right) \right)
\\
\displaystyle
\quad
  = 2^3 \frac{R^3}{2^3} e^{- R} \left( \frac{1}{3} \frac{1}{R} + \frac{1}{R^2} + \frac{1}{R^3} \right)
\\
\displaystyle
\quad
  = e^{- R} \left( 1 + R + \frac{R^2}{3} \right)
指数関数的に減少。
 R\rightarrow0で規格化条件に一致。

  • (一電子)クーロン積分


\displaystyle
\alpha \equiv \int dxdydz \, \psi^*_1(x,y,z) H \psi^*_1(x,y,z)
\\
\displaystyle
\quad
  = \int dxdydz \, \psi^*_1(x,y,z) \left( H_1 - \frac{1}{r_2} + V_n \right) \psi_1(x,y,z)
\\
\displaystyle
\quad
  = E_{1s} - U + \frac{1}{R}
\\
\displaystyle
U \equiv \int dxdydz \, \psi^*_1(x,y,z) \frac{1}{r_2}\psi_1(x,y,z)
\\
\displaystyle
\quad
  = \frac{R^3}{2^3} \int d\xi d\eta d\phi \, (\xi^2 - \eta^2)\frac{2}{R(\xi - \eta)} \frac{e^{-R(\xi+\eta)}}{\pi}
\\
\displaystyle
\quad
  = \frac{R^2}{2} \int d\xi d\eta \, (\xi + \eta) e^{-R(\xi+\eta)}
\\
\displaystyle
\quad
  = \frac{R^2}{2} \int^{1}_{-1} d\eta \, \left( \frac{1}{R} + \frac{1}{R^2} + \frac{\eta}{R} \right) e^{-R(1+\eta)}
\\
\displaystyle
\quad
  = \frac{R^2}{2} \left(
    \left( \frac{1}{R} + \frac{1}{R^2} \right) \left( \frac{e^{-2R} - e^0}{-R} \right)
     + \frac{1}{R} \left( \frac{ e^{-2R} - (-e^0) }{-R} + \frac{1}{R}\frac{ e^{-2R} - e^0 }{-R} \right)
 \right) 
\\
\displaystyle
\quad
  = - \frac{1}{2} \left(
    \left( 1 + \frac{1}{R} \right) \left( e^{-2R} - 1 \right)
     + \left( e^{-2R} + 1 + \frac{ e^{-2R} }{R} - \frac{1}{R} \right)
 \right) 
\\
\displaystyle
\quad
  = - e^{-2R} \left( 1 + \frac{1}{R} \right) + \frac{1}{R}
\\
\displaystyle
\therefore
  -U = e^{-2R} \left( 1 + \frac{1}{R} \right) - \frac{1}{R}
 -Uは隣の原子核からの引力ポテンシャルによる寄与であり、 R \rightarrow 0 で水素原子内のクーロンポテンシャルの期待値に一致する。
ある意味、非常に量子力学っぽい(電子を古典的な点電荷と思うと、原子核に無限に近付ければクーロン力は負に発散するが、量子力学的に電子「雲」としての構造を考慮するため、発散が抑えられて適切に扱えている)。
よって、 R \ll 1におけるクーロン積分内の斥力は、原子核間反発に帰着する。
 R \gg 1では第一項は指数関数的に減衰し、第二項が原子核同士の斥力を打ち消すから、遠いところでは何も損得が起きないのは納得できる。


\displaystyle
\beta \equiv \int dxdydz \, \psi^*_1(x,y,z) H \psi^*_2(x,y,z)
\\
\displaystyle
\quad
  = \int dxdydz \, \psi^*_1(x,y,z) \left( H_2 + V_n - \frac{1}{r_1} \right) \psi_2(x,y,z)
\\
\displaystyle
\quad
  = \left( E_{1s} + \frac{1}{R} \right) S  - J
  = T + \frac{S}{R} - 2 J
\\
\displaystyle
T = E_{1s} S + J = \int dxdydz \, \psi^*_1(x,y,z) \frac{ \nabla^2_2 }{ 2 } \psi_2(x,y,z)
\\
\displaystyle
J \equiv \int dxdydz \, \psi^*_1(x,y,z) \frac{1}{r_1} \psi_2(x,y,z)
  = \int dxdydz \, \psi^*_1(x,y,z) \frac{1}{r_2} \psi_2(x,y,z)
\\
\displaystyle
\quad
  = \frac{R^3}{2^3} \int d\xi d\eta d\phi \, ( \xi^2 - \eta^2 ) \frac{ 2 }{ R ( \xi + \eta ) } \frac{ e^{ -R \xi } }{ \pi }
\\
\displaystyle
\quad
  = \frac{ R^2 }{ 2 } \int d\xi d\eta \, ( \xi - \eta ) e^{ -R \xi }
  = \frac{ R^2 }{ 2 } \int d\xi \, 2\xi e^{ -R \xi }
\\
\displaystyle
\quad
  = R^2 \left( \frac{ e^{-R} }{ R } + \frac{ e^{-R} }{ R^2 } \right)
  = e^{-R} \left( 1 + R \right)
\\
\displaystyle
\therefore
T = -\frac{S}{2} + J = -\frac{1}{2} e^{-R}\left( 1 + R + \frac{R^2}{3} \right) + e^{-R}\left( 1 + R \right)
\\
\displaystyle
\quad
  = \frac{1}{2} e^{-R}\left( 1 + R - \frac{R^2}{3} \right)
\\
\displaystyle
\beta = ( T - 2J ) + \frac{S}{R}
\\
\displaystyle
\quad
  = - \frac{3}{2} e^{-R}\left( 1 + R + \frac{1}{9} R^2 \right) + e^{-R} \left( \frac{1}{R} + 1 + \frac{1}{3} R \right)
電子が飛び移ることで、運動エネルギー Tが負(つまり得をする)になることがわかる。 Rが小さいと、運動エネルギー的には損するが、交換積分Jの寄与で得をする。
電子波動関数が重なって不安定になるのが原子核間反発( S/R)って、納得しにくい。波動関数の重なりと原子核間反発は直接的には繋がっていないから、この項はよくわからない。

わかっているつもりでも、やってみるとわかってないことがゴロゴロ出てきて、凹むが、再発見出来るのはそれはそれとして面白い。

matplotlibで二文字以上の変数を上付きにする。

matplotlibで普通に

from matplotlib import pyplot as plt
atomic_number = 29
plt.title( 'Cu$^{}$'.format( atomic_number ) )

とやると、「Cu ^29」となってしまう。

これは、内部で「$^29$」と書かれたと見做されているためである。
ベタ打ちするならば、「$^{29}$」で回避出来る。
しかし、今の場合は変数を入れたいので、「$^{{}}$」などが思い付くが、認識されずに上手く行かない。

そのため、これを避けるには、一文字ずつに分けるしかない、という結論に至った。

from matplotlib import pyplot as plt
atomic_number = 29
an = str( atomic_number )
plt.title( 'Cu$^{0}$$^{1}$'.format( an[0], an[1] )

もっと良い方法があれば知りたい。