各軌道のエネルギーを見ると、2s軌道と2p軌道、3s軌道と3p軌道の縮退が解けており、それぞれ、2s軌道の方が2p軌道よりも、3s軌道の方が3p軌道よりも低いエネルギーを持つことがわかります。これは、動径方向の波動関数を見ると、2p軌道よりも2s軌道の方が、また、3p軌道よりも3s軌道の方が内側に局在しており、内側の軌道に入っている電子からの反発力を受けやすくなるためだと解釈できます。

この議論により、電子が主量子数nが同じ軌道に入る時に、s、p、dの順番で入るという事実は、水素原子のモデルでは説明ができませんでしたが、電子間相互作用を考えることで説明ができるということがわかりました。

f:id:AiRy:20190615001748p:plain — C原子の電子軌道の動径方向波動関数

軌道	エネルギー
1s	-9.944
2s	-0.501
2p	-0.199
3s	-0.006
3p	0.004

C原子の電子軌道のエネルギー

2019-06-02

密度汎関数法入門(3)-計算例(実装)

計算物理(量子力学) 計算物理(プログラム)

はじめに
計算対象
プログラムの流れ

はじめに

この記事は、密度汎関数法(DFT)についての紹介記事で、

密度汎関数法入門(1)-理論 - Airy's blog

密度汎関数法入門(2)-計算方法 - Airy's blog

の続きになります。

この記事では、DFTを自分で実装して遊んでみたいという方や、DFTの現実系への適用例を知りたいと言う方向けに、比較的手軽な実装例を紹介します。

プログラムの中身はここで公開しています。

GitHub - Airy07/DFT_of_atoms

記事では、このプログラムに沿って話をしていますが、中身を見なくても内容は理解できると思います。

なお、ここで公開したプログラムを、このページ/Githubのページに言及することなく再配布、編集して公開することは禁止いたします。(こんなことをわざわざ言うほど役に立つ計算ができるプログラムではありませんが念のため)

続きは

密度汎関数法入門(4)-計算例(結果) - Airy's blog

です。

計算対象

ここで紹介するプログラムでは、原子の電子軌道の計算を行います。計算においては以下の事柄を仮定しています。

1s, 2s, 2p, 3s, 3p軌道の順番に電子が詰まっていく。3d以上の軌道は考えない。
全ての軌道は軌道角運動量演算子 $\hat{L}^2$ の固有状態であり、s軌道は $l=0$ 、p軌道は $l=1$ に対応する。つまり、各軌道の波動関数は $\hat{L}^2$ の固有値 $\hbar^2 l(l+1)$ に対応する角度方向の関数 $Y_l(\theta, \phi)$ と動径方向の関数 $R(r)$ の積として $\Psi(r, \theta, \phi) = R(r)Y_l(\theta, \phi)$ とかける。
基底状態における電子密度は球対称と仮定する。
計算にはLDA汎関数の一つであるPZ81を用いる。

プログラムの流れ

計算の大まかな流れは前の記事で紹介した密度汎関数法の計算方法の通りですが、球対称性のために一部の計算が簡単になっています。このプログラムの大まかな流れは以下のようになります。

入力

入力から原子番号Zを受け取り、各軌道に入る電子の数(配列N)を計算します。3p軌道までしか考えていないので、原子番号Zは1~18までしか計算できません。

初期状態の設定

各軌道の初期状態として、原子核の電荷が+Zの水素様原子の解を用いています。この初期状態の関数の定義はparameters.hppの中に書かれています。

ここから、Kohn-Sham equation

$(-\frac{1}{2}\nabla^2 + \int \frac{n(r')}{|r-r'|}dr' + V_{ext}(r) + V_{xc}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

のSelf consistentな解を探す段階に入ります。

軌道の正規化

各軌道を正規化条件

$\int_0^\infty r^2 R^2(r)dr = \int_0^\infty |u(r)|^2 dr= 1$

が成立するように正規化します。

電子密度の計算

正規化した軌道を用いて電子密度

$n(r) = \sum_{n:軌道} N_n|R_n(r)|^2 = \sum_{n:軌道}N_n|\frac{u_n(r)}{r}|^2$

を計算します。ここで、 $N_n$ は軌道nに入っている電子の数を示しています。ここで、軌道とは1s, 2s, 2p, 3s, 3p軌道のことを指します。

ポテンシャルの計算

次は、Kohn-Sham方程式に現れるポテンシャルの計算をします。

外場ポテンシャル $V(r)$ は、中心の原子核によるクーロンポテンシャルですので、電子密度には依存しません。電子密度に依存する項は、電子間のクーロン相互作用 $\int \frac{n(r')}{|r-r'|}dr'$ とexchange correlation potential $V_{ext}(r)$ です。

ここで、電子密度は球対称と仮定しているので、電子間のCoulomb相互作用ポテンシャル

$V_{Coulomb}(r) = \int \frac{n(r')}{|r-r'|}dr'$

について、

$\frac{d^2}{dr^2} (rV_{Coulomb}(r)) = -4\pi r n(r)$

が成立します。したがって、電子間のCoulomb相互作用ポテンシャルは、面倒な積分を計算することなく、動径方向の微分方程式を適切な境界条件の元に解くだけで計算することができます。

また、exchange correlation potential $V_{ext}(r)$ についてはLDA汎関数を用いているので、先ほど計算した電子密度から直接計算できます。

電子軌道の計算

次に、動径方向の方程式を解きます。動径方向の方程式は、先ほど計算したポテンシャルを全てまとめて $V_{tot}(r)$ とおくと、

$\left(-\frac{1}{2}\frac{d^2}{dr^2}+V_{tot}(r) + \frac{l(l+1)}{2r^2} \right) u(r) = Eu(r)$

となります。ただし、 $u(r) = rR(r)$ です。これを、 $r$ が十分大きいところから数値的に積分し、 $r=0$ で $u(r)=0$ を満たすようなエネルギー固有値 $E$ を二分法で探索します。ここでは、二分法の初期条件を適切に工夫することで、各軌道角運動量 $l$ の状態に対して小さい順にエネルギー固有値を計算できるようにしています。

このプログラムでは、動径方向の方程式を解く関数はradial_eigensolver.hppで定義されています。

また、このプログラムではSCF equationを解くループで解が振動することを防ぐために、得られた動径方向の波動関数を正規化して、前のループで得られた波動関数との線型結合にするという操作を入れています。

収束条件の判定

得られた軌道エネルギーの変化が十分に小さくなったら、Self-Consistentな解が得られたと見なしてループを抜けます。収束条件が満たされない場合、得られた軌道を使って同じプロセスを繰り返します。

ループを抜けた場合、得られた解から得られる情報を計算します。

基底状態エネルギーの計算

得られた状態から、基底状態エネルギーを計算します。基底状態のエネルギーは軌道エネルギーの和ではないので、改めて計算する必要があります。

解の出力

得られた解を出力します。

実行など

$ make atom_DFT2

でコンパイルして

$ ./atom_DFT2 (原子番号)

で実行できます。

人が使いやすいように意識して書いたプログラムではないので、面倒で申し訳ありませんが、細かい設定を変えて遊んでみたいという方は、以下を参考にしてください。

現在は波動関数の計算結果などが大量に出力されるので、出力されるものを変えたい方はatom_DFT2.cppの中のmain関数の中身を変えてください(基本的に全ての出力はmainの中に書いてあるはずです)。
汎関数の形を変えたい方は、LDA-PZ81.hpp(現在使われているPZ81汎関数が書かれている)と同じようなheader fileを自分で作り、atom_DFT2.cppとradial_eigensolver.hppでincludeされている汎関数header fileを変更してください。
収束条件や最大繰り返し回数、動径方向の関数を計算する半径の刻み幅や最大半径などのパラメータを変更したい場合は、parameters.hppの中身を変更してください。

2019-05-24

密度汎関数法入門(2)-計算方法

計算物理(量子力学)

はじめに
Self-Consistent Field equation
Exchange Correlation Functional
- LDA(Local Density Approximation)
- GGA(General Gradient Approximation)

はじめに

この記事は、密度汎関数法についてのシリーズ記事で、

密度汎関数法入門(1)-理論 - Airy's blog

の続きになります。

この記事では、前回導出したKohn-Sham equationを数値的にどう解くかということを説明します。

続きは

密度汎関数法入門(3)-計算例(実装) - Airy's blog

密度汎関数法入門(4)-計算例(結果) - Airy's blog

です。

Self-Consistent Field equation

Kohn-Sham equation

$(-\frac{1}{2}\nabla^2 + \int \frac{n(r')}{|r-r'|}dr' + V_{ext}(r) + V_{xc}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

において、ポテンシャルの項のうち、Coulomb相互作用の項 $\int \frac{n(r')}{|r-r'|}dr'$

と、exchange correlation potential $V_{xc}(r)$ は基底状態の電子密度によって決まる項であり、外場ポテンシャル $V_{ext}(r)$ とは違って最初から形が決まっているものではありません。

したがって、Kohn-Sham equationを解いた時に、 $\psi_j(r)$ が、 $\psi_j(r)$ から求めた電子密度 $n(r)$ を使って求めた実効ポテンシャル下における一粒子方程式の解になっている必要があります。このような方程式を、self-consistent field equation(SCF equation)とよびます。

SCF equationであるKohn-Sham equationを数値的に解くには、まず、初期配位を仮定し、それを用いて電子密度を計算し、そこからCoulomb相互作用ポテンシャルと $V_{xc}(r)$ を計算し、得られたポテンシャル下での一粒子方程式を解きます。この段階では、Coulomb相互作用ポテンシャルと $V_{xc(r)}$ は最初に仮定した配位を用いて計算したものであり、一粒子方程式を解いて得られた解とはconsistentではないと考えられます。そこで、現在の配位から電子密度を計算、電子密度から実効ポテンシャルを計算、得られた実効ポテンシャルを用いて一粒子方程式を解く、というプロセスを状態が収束するまで繰り返すことで、self-consistentな解を得ることができます。*1

Exchange Correlation Functional

次に、exchange correlation potential $E_{xc}[\tilde{n}]$ について考えてみましょう。Kohn-Sham equationを解くには、電子密度[\tilde{n}]から $V_{xc}(r)$ を計算できる必要があり、そのためには $E_{xc}[\tilde{n}]$ の形がわかっている必要があります。しかし、 $E_{xc}[\tilde{n}]$ は[\tilde{n}]の関数形によって決まる汎関数であり、[\tilde{n}]の簡単な微分や積分によって書ける保証はないので、それを決定するのは非常に困難です。また、 $E_{xc}[\tilde{n}]$ の定義は

$E_{xc}[\tilde{n}] = F[\tilde{n}] -T_s[\tilde{n}] - \frac{1}{2}\int \frac{\tilde{n}(r)\tilde{n}(r')}{|r-r'|} dr dr'$

であり、全体から主要な項を引いた残りといった風になっているので、解析的な計算によって $E_{xc}[\tilde{n}]$ を決めることは難しいでしょう。

実際には、DFTは数値計算のための理論ですので、数値計算がしやすいような $E_{xc}[\tilde{n}]$ の近似式が用いられます。代表的な $E_{xc}[\tilde{n}]$ の近似としては、以下のようなものがあります。それぞれの詳細や、具体的な関数形については

常田貴夫著「密度汎関数法の基礎」
Naito, T., Akashi, R., & Liang, H. (2018). Application of a Coulomb energy density functional for atomic nuclei: Case studies of local density approximation and generalized gradient approximation. Physical Review C, 97(4), 044319.

などを参照してください。一つ目の文献では、様々な種類の汎関数が物理的正当性や再現性などの観点から詳しい説明がついて紹介されています。二つ目の文献では、汎関数の具体形がパラメータ込みで簡潔にまとめられているので、自分でDFTを実装してみたいときなどにおすすめです。

LDA(Local Density Approximation)

$V_{xc}(r)$ が考えている点の電子密度 $\tilde{n}(r)$ だけで決まるとする近似をLocal Density Approximation(LDA)とよびます。この時、exchange correlation functionalは

$E_{xc}[\tilde{n}] =\int \tilde{n}(r)\epsilon_{xc}(\tilde{n}(r)) dr$

という形に書くことができます。

GGA(General Gradient Approximation)

LDAに、 $\nabla \tilde{n}(r)$ による補正を取り入れた近似をGGAとよびます。

*1:このプロセスはいつも収束するとは限りません。これを収束させるには、良い初期状態を用いて計算することや、場合によってはこの他に状態を収束させるような工夫をする必要があります。

2019-05-20

密度汎関数法入門(1)-理論

計算物理(量子力学)

はじめに
量子多体系の計算の難しさ
概要
Hohenberg-Konhの定理
Kohn-Sham equation

はじめに

先日の五月祭で密度汎関数法(Density Functional Theory, DFT)に関する展示をしたのですが、我ながらなかなかの力作に仕上がったので、ここで解説記事を書くことにしました。

密度汎関数法は、相互作用する量子多体系の問題を、その基底状態の粒子密度を再現するような相互作用のない問題に書き換えることで解く方法で、化学や固体物理における電子状態計算などによく用いられる手法です。

このシリーズ記事では、

最初にDFTの基礎理論と、数値計算の手法の簡略な紹介
現実の系への応用例として、原子の電子軌道の計算の説明

について書きます。特に、この記事では最初の基礎理論の紹介を行います。

なお、以下では電子の状態を考えることを仮定して話をするので、他の量子多体系を考えたい場合はCoulomb相互作用を仮定している部分を一般の相互作用に書き換えて議論を追ってみてください。

また、この記事を書く上で

Kohn, W. (1999). Nobel Lecture: Electronic structure of matter—wave functions and density functionals. Reviews of Modern Physics, 71(5), 1253.
常田貴夫著「密度汎関数法の基礎」

などを参考にしました。

続きは

密度汎関数法入門(2)-計算方法 - Airy's blog

密度汎関数法入門(3)-計算例(実装) - Airy's blog

密度汎関数法入門(4)-計算例(結果) - Airy's blog

です。

量子多体系の計算の難しさ

DFTについて紹介する前に、なぜDFTが必要なのかということを説明します。

例えば、古典系においてNewtonの運動方程式に従って運動するN粒子の問題を考えると、この場合は各時間におけるN個の粒子の位置と運動量を保持しておいて、それをNewtonの運動方程式に従って時間発展させれば良いことがわかります。

それに対し、量子力学でN粒子系の状態を決めるにはN体波動関数 $\Psi(r_1, r_2, \cdots, \r_N)$ の関数形を決めてやる必要があります。数値的に関数形を保持するには、空間内のすべての位置に対応する関数の値を保持しておくことはできないので、空間内の有限の領域を離散化して、各格子点での関数の値を覚えておくことになります。ここで、 $x,y,z$ 軸をそれぞれ[n]分割したとすると、一つの位置がとりうる格子点の個数は $n^3$ であり、N粒子波動関数の引数の組み合わせは $n^{3N}$ 通りになることがわかります。これは、N体波動関数の自由度がNの指数関数で発散することを示しており、したがって波動関数を直接扱う方程式であるSchroedinger方程式は数値計算には向かない方程式であるということになります。このように、基礎方程式を直接数値的に計算できず、何らかの工夫をしなければならないということが古典と量子の場合の大きな違いになります。

いかに説明するように、DFTはこの問題を回避して問題を解ける形にする工夫の一つです。

概要

DFTの基本的な理論は、以下の二つの定理からなります。

定理1(Hohenberg-Kohnの定理)

量子多体系の全ての情報は、基底状態の粒子密度に含まれる。

定理2(Kohn-Sham equation)

Coulomb相互作用で相互作用する電子の基底状態の電子密度は、相互作用のない1粒子Schroedinger方程式

$\left(-\frac{1}{2}\nabla^2 + \int dr' \frac{n(r')}{|r-r'|} + V_{ext}(r) + V_{xc}(r) \right) \psi_j(r) = \epsilon_j \psi_j(r)$

を解いて得られる電子密度に等しい。

以下では、それぞれの定理について証明を含めてより詳しく説明します。

Hohenberg-Konhの定理

Hohenberg-Kohnの定理の正確なステートメントは以下の通りです。

Hohenberg-Kohnの定理

あるポテンシャルの下での電子ガスの基底状態の電子密度がわかっているとき、基底状態の電子密度がそれと等しくなるようなポテンシャルは定数のシフトを除いて一位に定まる。

証明

今、考えている基底状態の電子密度を $n(r)$ とし、基底状態の電子密度が $n(r)$ となるような二つの異なるポテンシャル $V_1(r), V_2(r)$ があったとする。

ここで、Hamiltonianのうち、外場ポテンシャル以外の項(運動エネルギーと電子館相互作用)を $H_0$ とおき、

$H_1 = H_0 + V_1$

$H_2 = H_0 + V_2$

とする。 $H_1, H_2$ の基底状態をそれぞれ $\Psi_1, \Psi_2$ とおくと、

$( \Psi_1, H_1\Psi_1) \lt (\Psi_2, H_1\Psi_2)$

より

$(\Psi_1, H_0\Psi_1) + \int dr n(r)V_1(r) \lt (\Psi_2, H_0\Psi_2) + \int dr n(r)V_1(r)$

同様に

$( \Psi_2, H_2\Psi_2) \lt (\Psi_1, H_2\Psi_1)$

より

$(\Psi_2, H_0\Psi_2) + \int dr n(r)V_2(r) \lt (\Psi_1, H_0\Psi_1) + \int dr n(r)V_2(r)$

辺々足して

$(\Psi_2, H_0\Psi_2) + (\Psi_1, H_0\Psi_1) \lt (\Psi_2, H_0\Psi_2) + (\Psi_1, H_0\Psi_1)$

となるが、これは矛盾である。したがって、題意が示された。

ただし、証明の途中で基底状態が縮退していないことを仮定した。

Kohn-Sham equation

Coulomb相互作用で相互作用する電子の基底状態の電子密度は、以下のKohn-Sham equation

$\left(-\frac{1}{2}\nabla^2 + \int dr' \frac{n(r')}{|r-r'|} + V_{ext}(r) + V_{xc}(r) \right) \psi_j(r) = \epsilon_j \psi_j(r)$

を解いて得られる電子密度に等しい。ただし、未定義の量 $V_{xc}(r)$ については以下の証明の過程で定義する。

証明

最初に

$T_s[\tilde{n}] = \min_{\Psi \to \tilde{n}} [ (\Psi| T |\Psi)]$

とする。ここで、 $T$ は運動量演算子、 $\min_{\Psi \to \tilde{n}}$ は、電子密度が[\tilde{n}]となるような状態の上でminを取るという意味である。これは、先ほどのHohenberg-Kohnの定理の議論をnon-interatcing electron gasについて行うと $T_s$ が電子密度 $\tilde{n}$ の汎関数として定義可能であることがわかる。

次に、

$F[\tilde{n}] = \min_{\Psi \to \tilde{n}} [ (\Psi| H |\Psi)] = T_s[\tilde{n}] + \frac{1}{2}\int \frac{\tilde{n}(r)\tilde{n}(r')}{|r-r'|} dr dr' + E_{xc}[\tilde{n}]$

として、 $F[\tilde{n}]$ とexchange correlation functional $E_{xc}[\tilde{n}]$ を定義する。ここで、 $T_s、F$ が外場ポテンシャルの形に依存しないuniversalな汎関数であるということに注意してください。

このとき、Hamiltonian Hの基底状態エネルギーは

$E_{gnd} = min_\tilde{n} (F[\tilde{n}] + \int \tilde{n}(r)V_{ext}(r) dr ) = min_\tilde{n} (T_s[\tilde{n}] + \frac{1}{2}\int \frac{\tilde{n}(r)\tilde{n}(r')}{|r-r'|} dr dr' + E_{xc}[\tilde{n}]+ \int \tilde{n}(r)V_{ext}(r) dr)$

である。基底状態の電子密度を $n(r)$ とすると、変分法より

$\left (\frac{\delta T_s[\tilde{n}] }{\delta \tilde{n}(r)} + \int \frac{\tilde{n}(r')}{|r-r'|} dr' +\frac{\delta E_{xc}[\tilde{n}] }{\delta \tilde{n}(r)} + V_{ext}(r) \right)_{\tilde{n} = n}=0$

を得る。

ここで一度、相互作用のない電子ガスを考えると、

$\frac{\delta T_s[\tilde{n}] }{\delta \tilde{n}(r)} + V_{ext}(r) = 0$

を満たすような電子密度を再現するには、外場の中で運動する一粒子のSchroedinger方程式

$(-\frac{1}{2}\nabla^2 + V_{ext}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

を解けば良いことがわかっている。

したがって、相互作用のある場合では、基底状態の電子密度を正しく再現するには、実効的な一粒子の方程式

$(-\frac{1}{2}\nabla^2 + \int \frac{n(r')}{|r-r'|}dr' + V_{ext}(r) + V_{xc}(r))\psi_j(r) = \epsilon_j \psi_j(r)$

を解けば良いことがわかる。ただし、exchange correlation potential $V_{xc}(r)$ は

$V_{xc}(r) = \frac{\delta E_{xc}[\tilde{n}]}{\delta \tilde{n}(r)}$

として定義される。

ここで得られた一粒子の方程式をKohn-Sham equationとよびます。

ここで、相互作用するN粒子の問題が相互作用しないN粒子の問題に書き換えられました。相互作用しないN粒子の問題は、実質的にはN個の一粒子問題と考えることができるので、問題を解ける形に落とし込むことができました。

また、DFTが優れている点は、上で示した二つの定理が厳密であるという点です。DFTではN体波動関数を扱わないので系の情報の一部を捨てていることになりますが、Hartree-Fock法などとは異なり、系の状態空間を制限するというタイプの近似が使われていません。したがって、exchange correlation functional $E_{xc}$ の形を良くすれば、原理的にはいくらでも厳密解に近い結果が得られることになります。

2019-02-25

Hartree-Fock method (1)-理論

計算物理(量子力学)

はじめに-シリーズについて
はじめに-(1)理論について
Slater行列式とは
Hartree-Fock equationの導出
いくつかの手法-UHFとRHFなど
Restricted Hartree-Fock equation
Fockエネルギーの解釈-Koopmans' theorem

はじめに-シリーズについて

このシリーズ記事では、量子化学でよく使われるHartree-Fock法を扱います。

Hartree-Fock法は、多電子系の多体波動関数をSlater行列式(一電子状態の反対称積)に制限し、変分法によってエネルギー期待値を最小化するSlater行列式を求めることで基底状態を近似しようという方法です。

多体系で相関を無視する近似はよく用いられる手法です(Ising modelの平均場近似など)が、現実系の計算をする際に良い近似になるのか、というのが気になると思うので

最初に、Hartree-Fock equationの導出と解釈について述べる。
次に、Hartree-Fock法の計算方法/手法についてまとめる。
最後に、実装方法と計算例を紹介する。

という構成にしようと思っています。理論自体はそれほど難しくないのですが、実装の細かい部分にかなり手間がかかるので、特に実装方法について詳しく述べたいです。

注意2) 続きはまだ書いてないです。書いたらここにリンクを貼っていきます。

なお、この記事は、

Thijssen著「Computaional Physics」
Szabo & Ostlund著「Modern Quantum Chemistry」

を参考にしています。

はじめに-(1)理論について

HF法とは、多体波動関数をSlater行列式で近似し、変分法によってエネルギー期待値を最小化することで、多電子系の基底状態を近似的に計算する方法です。

(1)-理論では、Hatree-Fock equationの導出と、それについての物理的解釈について述べます。

なお、これ以降、 $x = (r, s)$ で空間座標 $r$ とスピン $s$ の組を表すことにします。

Slater行列式とは

電子はフェルミオンであるので、波動関数は粒子の交換に対して反対称でなくてはならない。

最も単純な多粒子系の状態は、各粒子が完全に独立な一粒子状態を取る状態

$\Psi({x}_1,{x}_2,\cdots,{x}_N) = \psi_1({x}_1) \psi_2({x}_2)\cdots \psi_N({x}_N)$

であるが、これは波動関数の反対称性を満たさないので、多電子系の波動関数としては不適切である。反対称性を満たす最も簡単な多粒子系の波動関数は、これを反対称化した

$\Psi({x}_1,{x}_2,\cdots,{x}_N) = \frac{1}{\sqrt{N!}}\sum_{P\in S_n}\mathrm{sgn}(P) \psi_{P(1)}({x}_1) \psi_{P(2)}({x}_2)\cdots \psi_{P(N)}({x}_N)$

となる。これをSlater行列式という。ここで、 $\psi_i$ は互いに直行する規格化された一粒子状態である。

Hartree-Fock equationの導出

Hartree-Fock equationの導出については多くの文献があるので、ここでは概略だけにとどめます。必要であれば、

Thijssen 「Computational Phyisics」
Szabo, Ostlund 「Modern Quantum Chemistry」(混乱が起きない記法を用いているので、非常にわかりやすいです)
また、以下の記事もとてもわかりやすいと思います。
ハートリーフォック方程式の導出 - 物理とか

などを参照してください。

$N_A$ 個の原子核と $N$ 個の電子からなる系を考える。

原子核の運動を無視するBorn-Oppenheimer近似を用いると、原子核の影響は電子が感じる静的なポテンシャルとして扱えば良い。

電子の運動エネルギー
電子と原子核のCoulomb相互作用によるエネルギー
電子間相互作用のエネルギー

を考慮すると、考えている系のHamiltonianは

$H = \sum_{i=1}^{N}\left( -\frac{1}{2}\nabla_i^2 - \sum_{A=1}^{N_A} \frac{Z_A}{|r_i- R_A|} \right) + \frac{1}{2} \sum_{i\neq j}^{N} \frac{1}{|r_i - r_j|}$

と書ける。

前述のSlater行列式について、このHamiltonianのエネルギー期待値を計算すると、

$E = \sum_{i=1}^N [\psi_i|h|\psi_i ] + \frac{1}{2}\sum_{i, j}^{N}\{[ \psi_i\psi_i|\psi_j\psi_j] - [\psi_i\psi_j|\psi_j\psi_i] \}$

となります。ここで、

$[\psi_i|h|\psi_j] = \int dx \psi_i^*(x)(-\frac{1}{2}\nabla^2 + \sum_A \frac{Z_A}{|r - R_A|})\psi_j(x)$

$[ \psi_i\psi_j|\psi_k\psi_l] = \int dx_1 dx_2 \psi_i^*(x_1)\psi_j(x_1)\frac{1}{|r_1 - r_2|}\psi_k^*(x_2)\psi_l(x_2)$

です。さらに、変分法を用いて、

制約 $[\psi_i|\psi_j] = \delta_{ij}$ のもとで

$\frac{\delta E}{\delta \psi_i} = 0$

となる条件を求めると、一粒子状態の満たすべき方程式

$f\psi_i(x) = \epsilon_i\psi_i$

ただし

$f\psi_i(x) = h\psi_i(x) + \sum_{j=1}^{N} (\int dx_2\psi_i(x)\frac{1}{|r - r_2|}\psi_j^*(x_2) \psi_j(x_2) - \int dx_2 \psi_j(x)\frac{1}{|r - r_2|}\psi_j^*(x_2) \psi_i(x_2))$

が得られる。これをHartree-Fock equationとよぶ。

新たに演算子 $J,K$ を

$J\psi(x)= \sum_{j=1}^{N} \sum_{j=1}^{N} \int dx_2\psi(x)\frac{1}{|r - r_2|}\psi_j^*(x_2) \psi_j(x_2)$

$K\psi(x)= \sum_{j=1}^{N} \sum_{j=1}^{N} \int dx_2\psi_j(x)\frac{1}{|r - r_2|}\psi_j^*(x_2) \psi(x_2)$

と定義すると、

$f = h + J + K$

と書き直せる。 $K$ はnon-localな演算子である。また、Hartree-Fock equationを解いて得られる一粒子状態をHartree-Fock stateという。

いくつかの手法-UHFとRHFなど

ここまでの議論では、考える近似的な多電子系の基底状態の満たす条件はSlater行列式で書ける事だけであった。しかし、場合によっては計算を簡単にするために、さらに条件を強めることがある。

特に、各HF stateが空間波動関数とスピン波動関数の積

$\psi_{2i-1}(x) = \phi(r) \alpha(s)$

$\psi_{2i}(x) = \phi(r)\beta(s)$

と書けると仮定して計算する手法をRestricted Hartree Fock method(RHF)とよぶ。

なお、 $\alpha(s),\beta(s)$ は、それぞれスピンup、スピンdownのスピン波動関数である。ここで考えているHamiltonianは時間反転(スピン反転)の対称性を持つので、これは不自然な仮定ではない。また、考えている近似的な基底状態として一般のSlater行列式を認める場合をUnrestricted Hartree-Fock method(UHF)とよぶ。

全ての非占有準位にspin upの電子とspin downの電子が一つずつ入っているclosed-shellの状態(H $_2$ など、多くの安定な分子はclosed-shellである)ではRHFは比較的良い近似になることが知られているが、open shell(一部のイオンなど。電子の数が奇数の場合は必ずopen-shellになる)の場合は電子対を形成していない電子のスピンが他の準位の分裂を起こすので、UHF的な取り扱いをする必要がある。

今後、この記事ではRestricted closed-shell Hatree-Fockの場合を中心に説明します。

Restricted Hartree-Fock equation

偶数個の電子を含む系について、Restricted closed-shell HFの仮定をHF equationに代入すると、Restricted Hartree-Fock equation

$f\phi_i(r) = \epsilon_i \phi(r)$

が得られる。ただし

$f\phi_i(r) = h\phi_i(r) + \sum_{j=1}^{N/2} (2\int dr_2 \phi_i(r)\frac{1}{|r - r_2|}\phi_j^*(r_2) \phi_j(r_2) - \int dr_2 \phi_j(r)\frac{1}{|r - r_2|}\phi_j^*(r_2) \phi_i(r_2))$

であり、また、基底状態のエネルギーは

$E = 2\sum_{i=1}^{N/2} (\phi_i|h|\phi_i ) + \frac{1}{2}\sum_{i, j}^{N/2}\{ 2(\phi_i\phi_i|\phi_j\phi_j) - (\phi_i\phi_j|\phi_j\phi) \}$

で与えられる。ここで、

$(\phi_i|h|\phi_j) = \int dr \phi_i^*(r)(-\frac{1}{2}\nabla^2 + \sum_A \frac{Z_A}{|r - R_A|})\phi_j(r)$

$( \phi_i\phi_j|\phi_k\phi_l ) = \int dr_1 dr_2 \phi_i^*(r_1)\phi_j(r_1)\frac{1}{|r_1 - r_2|}\phi_k^*(r_2)\phi_l(r_2)$

である。

前に登場した $[]$ の記号では空間積分に加えてスピンについての和を取っている( $\int dx$ )に対して、ここで $()$ で書いているものは空間内の位置に関する積分だけを含むので注意してください。これらは主に化学で用いられる記法らしいのですが、考えている空間を明確にしているのでとてもわかりやすいと思います。

Fockエネルギーの解釈-Koopmans' theorem

Hartree-Fock equationを解くと、通常は無限個のFock演算子の固有状態が得られ、基底状態は最もエネルギー状態の低いN個の状態によるSlater行列式で与えられる。

この時、occupied stateの $\epsilon_j$ は系のイオン化エネルギーの、unoccupied stateの $\epsilon_j$ は系の電子親和力の第一近似になる。これをKoopmans' theoremという。

Koopmans' thoeremの証明は以下の通りである。

まず、 $\psi_a$ がoccupied stateの場合を考える。イオン化が起こり、波動関数が[\phi_{2a-1}(x) = \psi_a(r) \alpha(s)]で表される電子が系から出て行ったとき、残った電子の状態が他のoccupied stateのSlater行列式で書けると仮定する(この仮定はHF近似の下でも厳密ではない)と、イオン化エネルギーは

$\Delta E = \sum_{i=1}^N [\psi_i|h|\psi_i ] + \frac{1}{2}\sum_{i, j}^{N}\{[ \psi_i\psi_i|\psi_j\psi_j] - [\psi_i\psi_j|\psi_j\psi_i] \}$

$- \sum_{i\neq 2a-1}^N [\psi_i|h|\psi_i ] + \frac{1}{2}\sum_{i, j\neq a}^{N}\{[ \psi_i\psi_i|\psi_j\psi_j] - [\psi_i\psi_j|\psi_j\psi_i] \}$

$= (\phi_a|h|\phi_a) + \sum_{i=1}^{N/2}\{2(\phi_a\phi_a|\phi_i\phi_i) - (\phi_i\phi_j|\phi_j\phi_i) \}$

$= \epsilon_a$

と計算できる。同様にして、unoccupied stateの[\epsilon_j]が電子親和力になることも証明できる。

実際は、occupied stateの[\epsilon_j]はイオン化エネルギーの妥当な見積もりを与えるのに対して、unoccupied stateの[\epsilon_j]は電子親和力の良い近似にはならない場合が多いらしいです。

2019-02-10

phase shiftと量子力学の散乱の計算について

計算物理(量子力学)

はじめに
仮定と用いる公式
計算方法
具体的な問題設定
計算結果
まとめ

はじめに

量子力学で散乱理論をやった時に、

式は追えたけど、何が起こっているかがわからなかった。
Born近似を用いた計算、低エネルギー極限、高エネルギー極限の場合以外の例を知らず、実際の現象に対するイメージがわからなかった。
特にphase shiftの意味がわからなかった。

などと納得できない部分が多かったが、計算機を用いた散乱理論の計算と可視化をすることで、散乱理論の使い方がある程度はっきりしてきたので

一般のポテンシャルに対して散乱断面積を計算する方法
phase shiftを用いた散乱断面積の式の有用性

についてまとめる。

なお、この記事における式や計算方法はThijssen著 Computational Physics第2章の内容を参考にしています。また、途中式や主張が間違っている/不正確という点があれば教えていただけるとありがたいです。

仮定と用いる公式

細かい公式の導出は

Thijssen著 Computational Physics 2.3節
J.J. Sakurai著 Modern Quantum Mechanics

などの教科書を参照して下さい。ここでは、大まかなストーリーと計算物理的に重要な式の列挙にとどめます。

3次元の球対称ポテンシャル $V(r)$ によって質量 $m$ の粒子が散乱される状況を考える。

このとき、入射方向周りの回転対称性より、波動関数は

$\psi(\boldsymbol{r}) = \sum_{l=0}^\infty A_l \dfrac{u_l(r)}{r} P_l(\cos\theta)$

と表すことができる。Schroedinger方程式

$\left(-\dfrac{\hbar^2}{2m}\nabla^2 + V(r)\right)\psi(\boldsymbol{r}) = E\psi(\boldsymbol{r})$

に代入すると、 $u_l(r)$ のみたす方程式は

$\dfrac{\hbar^2}{2m}\dfrac{d^2}{dr^2}u_l(r) = \left( V(r) + \dfrac{\hbar^2 l(l+1)}{2mr^2} - E \right)u_l(r)$

と計算できる。 $r\to\infty$ で、 $\frac{1}{r^2},V(r)$ はともに0に収束するので、 $u_l(r)$ は $r\to \infty$ の極限で正弦曲線へ漸近することがわかる。この漸近形をの位相を用いて $l$ 次のphase shift $\delta_l$ が

$u_l(r) \simeq \sin(kr - \frac{l\pi}{2} + \delta_l)$

と定義される。

この時、微分散乱断面積と全散乱断面積はphase shiftを用いて

$\dfrac{d\sigma}{d\Omega} = |\sum_{l=0}^\infty (2l+1)e^{i\delta_l} \sin(\delta_l)P_l(\cos\theta)|^2$

$\sigma_{tot} = \dfrac{4\pi}{k^2}\sum_{l=0}^\infty (2l+1)\sin^2\delta_l$

と表されることが知られている。

計算方法

$u_l(r)$ が原点で特異的な振る舞いをしないことを要請すると、方程式

$\dfrac{\hbar^2}{2m}\dfrac{d^2}{dr^2}u_l(r) = \left( V(r) + \dfrac{\hbar^2 l(l+1)}{2mr^2} - E \right)u_l(r)$

は定数倍を除いて数値的に解くことができる。 *1

ここで、 $V(r)$ が十分小さくなる領域まで微分方程式の解を積分し、漸近形

$\dfrac{u_l(r)}{r} \propto \cos(\delta_l)j_l(kr)-\sin(\delta_l)n_l(kr)$

を用いるとphase shiftが計算できる。*2

具体的な問題設定

ここでは、Thijssen著 Computational Physics2章に載っていた例である、Lennard-Jonesポテンシャルを仮定したKr原子によるH原子の散乱を考える。

$V(r)= \epsilon \left(\left(\dfrac{\rho}{r}\right)^{12} - \left(\dfrac{\rho}{r}\right)^6 \right)$

where

$\epsilon = 5.9$ meV

$\rho =0.357 \nu$ m

計算結果

まず、「遠心力の効果」を取り入れた有効ポテンシャル

$F(l,r) = V(r) + \dfrac{\hbar^2 l(l+1)}{2mr^2} - E$

を異なる $l$ に対してプロットすると、次の図のようになる。これを見るとわかる通り、 $l$ が大きくなると、遠心力の効果がポテンシャルの効果よりも強くなり、相対的にポテンシャルの効果が弱まると考えられる。

f:id:AiRy:20190210020447p:plain — 異なる $l$ に対する有効ポテンシャルのプロット

また、 $l=0,8$ に対して動径方向の波動関数と $j_l(kr)$ (の定数倍)をプロットしたものを以下の図に示した。

これより、 $l$ が大きくなると、動径方向の波動関数の立ち上がりが遅くなることがわかる。これは、 $j_l(x)\sim x^l (x\sim 0)$ であることや、 $l$ が大きくなると原点付近の有効ポテンシャル $F(l,r)$ が大きくなり、波動関数が大きい値を持つことができないことと関係している。したがって、 $l$ が大きい時、波動関数は $r\sim 0$ でほぼ0なので、相互作用ポテンシャルの形状に依存しなくなる。

f:id:AiRy:20190210020719p:plain — $l=0$ の場合の動径方向の波動関数と $j_0(kr)$ との比較

f:id:AiRy:20190210020739p:plain — $l=8$ の場合の動径方向波動関数と $j_8(kr)$ の定数倍との比較

以上の議論より、 $l\to \infty$ で[\delta_l\to 0]が成立することが定性的に確かめられた。

最も手計算でよく用いられるBorn近似は、実質的には摂動論なので、 $n$ 番目の項は $\left(\frac{V}{E}\right)^2$ 程度にと考えられるなる。よって、Born近似は入射粒子の運動エネルギーに比べて相互作用が強い場合に収束が遅くなる、もしくは発散する。それに対して、phase shiftを用いた展開の場合は、 $\delta_l$ が遠心力ポテンシャルに対する相互作用に対する相互作用ポテンシャルの大きさに対応すると考えると、 $l$ に関する展開項は $l$ が大きくなると0に収束する。よって、phase shiftを用いた展開は通常の摂動論とは全く別のアイデアに基づいており、*3収束性という点から数値計算に適した非常に有用な式であると言える。