第一原理計算と密度汎関数理論

f:id:bobchan1915:20181112145120p:plain

第一原理計算とは

第一原理計算は「非経験的電子状態計算」とも呼ばれている経験に基づかない（パラメーターがない）、第一原理である量子力学を基本法則に立脚した計算手法で物性値を知るための計算手法である。

量子力学とは、私たちが中学の頃に学んだ力学では扱えないような微視系を扱う現代物理学の分野である。

力学では動きにくさを規定する慣性質量である $m[kg$ ]、力 $F[N$ ]、加速度 $\alpha$ [m/ $s^2$ ]をもってして以下のようなニュートンの運動方程式を基礎原理であった。

\begin{align}
m\alpha = F
\end{align}

しかし、ド・ブロイ波というものが提唱した波動性と粒子性を結びつける考えが量子力学の視点だった。

質量 $m$ の粒子が速さ $v$ で運動する場合、以下の式で示される波長 $\lambda$ に相当する波であると見なせることを示した。つまりは粒子であり、波でもあるということだ。

\begin{align}
\lambda = \frac{h}{mv}
\end{align}

しかしながら、これは日常的な生活ではもちろん考えもしないことだろう。例えば野球ボールを想定していよう。野球ボールは148.8 $g$ であり、日本人最速記録をもっている大谷翔平選手の球速は165 $\frac{km}{h}$ である場合、野球ボールの波動性はどのようなものだろう。

物質波の波長を計算すると以下のようになる。

\begin{align}
\lambda = 9.72 × 10^{-31}[m]
\end{align}

である。波長があまりにも短いので私たちの日常の生活では見ることができないのである。

このような微視的な粒子を扱う量子力学ではシュレディンガー方程式を基礎方程式とする。これは粒子に対応する基礎法則があるならば、波動に関する波動方程式が物質波にもあるのではないかという疑問から始まったものである。

以下はポテンシャル $V$ 下にある時間依存性があるシュレディンガー方程式である。
ここで $\psi$ は波動関数と呼ばれ、粒子の量子的な状態を記述するために用いられる。

\begin{align}
\displaystyle i \hbar = - \frac{\hbar^2}{2m}\frac{\partial^2{\psi}}{\partial{x^2}} + V\psi
\end{align}

かなり簡単に導出をする。まず任意の関数は以下のように書くことができる。

\begin{align}
\psi(x,t) = A\cos{\left(\frac{x}{\lambda}-vt\right)}
\end{align}

この式を偏微分して、オイラーの公式と合わせて時間に依存したシュレディンガー方程式を導出することができるので、余力がある人は試してみてください。

時間に依存しないシュレディンガー方程式は以下のようになる。

\begin{align}
\frac{d^2 \psi}{d x^2} + \frac{2m}{\hbar^2}[E-V]\psi = 0
\end{align}

ここでハミルトニアンと呼ばれる演算子 $H$ を導入すると、以下のようになる。

\begin{align}
H\psi = E \psi
\end{align}

話を戻して、これに対して「経験的電子状態計算」というものは計算パラメーターの選定や経験的な部分の移植性など間違って計算をしてしまうケースがあるが、第一原理計算はモデルの任意性を避けた手法になっていて、経験の浅い人でも計算を正しく行うことができる。

しかし、第一原理計算も多体問題というものがある。電子が多すぎると複雑すぎて、電子状態を決定することができない。

例えば、N個電子を持つ多電子系のポテンシャル項 $V$ は以下の通りである。

\begin{align}
V = \underbrace{\sum_{i=1}^N{V(\rm r_i)}}_{ポテンシャル項} + \underbrace{\sum_{i=0}^N \sum_{j < i}U(r_i, r_j)}_{電子同士の相互作用項}
\end{align}

このため、何らしかの近似が必要となってくるが、その手法については後ほど解説をする。

第一原理計算には大きく分けて二種類の手法があり、「ハートリー・フォーク法」と「密度汎関数理論」である。

近似的方法

これら二つを解説するために変分原理というものを理解しなければならないので、軽く解説する。

大学までの初等量子力学では、水素原子までしか習わないところが多いかもしれないが、私たちの体を構成する物質をはじめ、多くの電子を含み、一般的には原子核より１０倍ぐらいの電子を含んでいる。

これらはシュレディンガー方程式では解けないが、近似的方法を使って求めることができ、その一つが変分原理である。

変分原理

たとえば以下のようなシュレディンガー方程式があるとする。

\begin{align}
\hat{H}\psi_0 = E_0 \psi_0
\end{align}

があり、ここで添字の $0$ は基底状態の物質量を表す。

ここで左側から $\psi_0^{*}$ をかけて全空間について積分すると以下のように基底状態のエネルギーが求まる。

\begin{align}
E_0= \frac{\int{\psi_0^{*} \hat{H} \psi_0 d\tau}}{\int{\psi_0^{*} \psi_0 d \tau}}
\end{align}

ここで $d\tau$ は適当な体積要素である。また、この基底状態の波動関数がわからなくても他の波動関数 $E_{\phi}$ とし、その時のエネルギー $E_{\phi}$ とすると以下のように求まる。

\begin{align}
E_{\phi} = \frac{\int{\phi_0^{*} \hat{H} \phi_0 d\tau}}{\int{\phi_0^{*} \phi_0 d \tau}}
\end{align}

となり、これは以下の式が成立する。

\begin{align}
E_{\phi} \geq E_0
\end{align}

これを変分原理という。また、この時に使用した $E_{\phi}$ を試行関数といって、その関数が何かしら基底状態の波動関数 $\phi_0$ を近ければ当然 $E_{\phi} = E_0$ となる。

しかしながら、どのように基底状態の波動関数 $\phi_0$ に試行関数を近づければいいのかはわからないので変分パラメーターというものをとる。

つまり、以下のように多変数パラメータをとって

\begin{align}
E(\alpha, \beta, \gamma,...) \geq E_0
\end{align}

のようにするわけであり、あるパラメータの変化 $d\alpha, d\beta...$ に対してのエネルギー変化を小さくなるように $\alpha, \beta...$ を定める。

ハートリー・フォーク法

ようやく変分原理を理解したので、ハートリー・フォーク法に取り組むことができる。この節ではスピン軌道 $\chi$ という波動関数 $\psi$ とスピン状態を含んだものを使う。

まず、原子核の取り扱いが重要となってくる。電子に対して非常に原子核は重いので、原子核は電子の運動につられず静止していると考えて、原子核の問題と電子の問題を分離してわかりやすくする。

これをボルン・オッペンハイマー近似と呼ぶ。つまり、電子の問題を考えている段階では、原子核のことについては考える必要がなくなるのである。

例えば、ヘリウムに対するシュレディンガー方程式を見てみよう。

以下のようになる。ここでHartree単位である。

\begin{align}
\left[ - \frac{1}{2} \nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{| \rm r_1 - \rm r_2|} \right] \chi(\rm x_1, \rm x_2) = E \chi(\rm x_1, \rm x_2)
\end{align}

ここで添字1, 2は電子1, 2の便宜的に区別するためのものである。左辺の第二項がクーロン相互作用になるのだが、この項のせいで微分方程式を変数分離で解析的に解くことができない。

そこで、その項を無視した、以下の式

\begin{align}
\left[ - \frac{1}{2} \nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{2}{r_1} - \frac{2}{r_2} \right] \chi(\rm x_1, \rm x_2) = E \chi(\rm x_1, \rm x_2)
\end{align}

で、また波動関数 $\chi(\rm x_1, \rm x_2)$ が以下のように変数分離すると仮定する。これをHatree積と呼ぶ。

\begin{align}
\chi(\rm x_1, \rm x_2) = \phi_1(\rm x_1)\phi(\rm x_2)
\end{align}

と変数分離をして考える。また、電子はフェルミ粒子であるので、２個の座標を入れ替えると符号が反転するはずである。その要請を満たすためには、いかのように式を改良しなければならない。

\begin{align}
\chi(\rm x_1, \rm x_2) = \frac{1}{\sqrt{2}}[\phi_1(\rm x_1) \phi_2(\rm x_2) - \phi_2(\rm x_1) \phi_1(\rm x_2)] \tag{1}
\end{align}

試しに $\rm x_1$ と $\rm x_2$ を入れ替えてみると要請を満たしているのがわかる。

\begin{align}
\chi(\rm x_2, \rm x_1) &= \frac{1}{\sqrt{2}}[\phi_1(\rm x_2) \phi_2(\rm x_1) - \phi_2(\rm x_2) \phi_1(\rm x_1)] \\
&= \frac{1}{\sqrt{2}}[-\phi_1(\rm x_1) \phi_2(\rm x_2) + \phi_2(\rm x_1) \phi_1(\rm x_2)] \\
&= - \chi(\rm x_1, \rm x_2)
\end{align}

また式1を行列式で表したものをスレータ行列式という。先ほどの例でもあったように、二つの電子が同じ波動関数を占めることができないというPauliの排他律を満たしていることがわかる。

\begin{align}
\chi(\rm x_1, \rm x_2) = \frac{1}{\sqrt{2}}\left|\begin{array}{cc} \phi_1(\rm x_1) & \phi_2(\rm x_1) \\ \phi_1(\rm x_2) & \phi_2(\rm x_2) \\\end{array}\right|
\end{align}

これを元のシュレディンガー方程式に代入すると以下のようになる。

\begin{align}
\left[ - \frac{1}{2} \nabla_1^2 - \frac{1}{2}\nabla_2^2 - \frac{2}{r_1} - \frac{2}{r_2} + \frac{1}{| \rm r_1 - \rm r_2|} \right] \phi_1(\rm x_1)\phi(\rm x_2) = E \phi_1(\rm x_1)\phi(\rm x_2)
\end{align}

これを左から $\phi^{*}(\rm r_2)$ を乗じて、そのあとに $\rm r_2$ で積分をすると以下のようになる。

\begin{align}
E^{\prime}=E-\int{d\vec{r}_{2}\phi^{\ast}\left( \vec{r}_{2}\right) \left( -\dfrac{1}{2}\nabla_{2}^{2}-\dfrac{2}{r_{2}}\right) \phi\left( \vec{r}_{2}\right)} \tag{2}
\end{align}

ここで規格化条件より

\begin{align}
\int \phi^{*}(\rm r_2)\phi(\rm r_2) d\rm r = 1
\end{align}

を利用している。

式2は自己無撞着問題の形をしている。つまり、解析的には解くことができず、数値的に解くしかない。

数値的な解法としては

初期波動関数を選ぶ
フォック演算子を使ってハートリー・フォック方程式の固有ベクトルと固有値を求める
エネルギーを計算する
一回前の波動関数とほぼ同等で、かつ、エネルギーもほぼ同等ならば終了
違うなら繰り返す

密度汎関数理論

密度汎関数理論の基本的な概念は相互作用のある多数の粒子の性質は基底状態の1電子密度 $n_o(r)$ の汎関数として考えることができるのである。

Hohenberg-Kohn(HK)の定理を基本とする密度汎関数理論が現れた。

HKと定理は以下の二つから構成される。

定理1(存在定理): 外部ポテンシャル $V_{ex}(\rm r)$ 中で相互作用している粒子系に対しては、ポテンシャル $V_{ex}(\rm r$ は定数をのぞいて基底状態の１電子密度 $n_0(\rm r)$ によって一意に決まる。
定理2(変分原理): ハミルトニアンはエネルギーの定数倍のシフトを除けば完全に決まるので、あらゆる状態(基底と励起）に対して多体波動関数が決まることになる。つまり、汎関数であるエネルギーを最低にする電子密度は真の電子密度である。

ある点 $\rm r$ での電子密度 $n(\rm r)$ は、以下のように波動関数から求めることができる。

\begin{align}
n(\rm r) = 2 \sum_i{\phi^{*}_i\phi_i} \tag{3}
\end{align}

結論から言うと、大事な部分は波動関数 $\psi$ と密度 $\rho$ は外部ポテンシャル $v$ に通して一対一対応していて、エネルギーは電子密度の汎関数なので $E[n(\rm r)$ ]と書くことができることを証明した。

特に、第一定理では、変分原理における波動関数 $\psi$ を電子密度 $\rho$ に置き換えて計算できることをいっていて、これが新しいのである。

\begin{align}
\psi \iff v \iff \rho
\end{align}

定理2でもあるようにエネルギーは電子密度の汎関数である。つまり、以下のようになっている。

\begin{align}
E[(\psi_i)] = E_{known}[(\psi_i)] + E_{XC}[(\psi_i)]
\end{align}

ここで $E_{known}[(\psi_i)$ は単純な解析的な形で表せる項目で、そのほかは $E_{XC}$ に分けていて、これは交換相関汎関数といい、すべての量子力学的な効果が含まれている。しかし、交換相関汎関数の形はわからなく、シュレディンガー方程式を解いて、全電子波動関数を求めるのは絶望的です。

ここで登場するのがKohn-Sham方程式である。これは１電子系からなる連立方程式を解くことで電子密度を求められることを示した。

Kohn-Shamの方程式は以下のようになっている。

\begin{align}
\left[ -\frac{\hbar}{2m}\nabla^2 + V(\rm r) + V_H(\rm r) + V_{XC}(\rm r) \right] \psi_i(\rm r) = \epsilon_i \psi(\rm r)
\end{align}

表面的には、あまり変わらないように思える。しかし、ここでは全電子系にみられるような $\sum_i$ のような項がないのがわかる。最初の $V(\rm r)$ は $V_{known}$ に対応する。

第二項はHatree項と呼ばれるものである。

\begin{align}
V_{H}(\rm r) = e^2 \int{\frac{n(\rm r^{'})}{| \rm r - \rm r^{'}|}} d^3 r^{'}
\end{align}

このポテンシャルは、Kohn-Shamの方程式が扱う一つの電子と系の電子から作られた全電子密度との間の相互作用を表すものである。注意点としては一つの電子は、全電子の一部でもあるので自己相互作用を含んでいる。本来は、自己相関は存在するものではない。

これを補正するために最後の項 $V_{XC}$ に含められていて、これは交換相関エネルギーの汎関数微分で定義することができる。

\begin{align}
V_{XC}(\rm r) = \frac{\delta E_{XC}(\rm r)}{\delta n(\rm r)}
\end{align}

ここまできたが、結局はHatreeポテンシャルを求めるためには電子密度を知らなければならない。しかし式3にあるように、電子密度を知るには波動関数を知らなければならない。またもや自己無撞着問題の形をしている。アルゴリズムは以下の通りである。

試行電子密度を決める
電子密度がわかるとポテンシャルを知ることができるので計算する
Kohn-Shamの方程式を解いて、1電子波動関数を得る
1電子波動関数から電子密度を計算する
計算された電子密度と試行電子密度を比較する。異なれば試行電子密度に何かしらの更新をして再度繰り返す。一致すれば終了する。

しかし、大きな問題がある。それは未だ $E_{XC}$ の汎関数の形がわかっていないのである。大きな問題のため、一つサブセクションを設けることにする。

交換相関汎関数

Hohenberg-Kohnの定理で存在定理として証明されているが、どのような形かはわからないので、数々の研究が行われている。幸いに、厳密に解くことができる一様電子ガスと呼ばれる場合がある。一様電子ガスとは、あらゆる場所において電子密度 $n(\rm r)$ が一定であるようなものである。
ここでは二つのメジャーな近似方法について取り上げる。

1. 局所密度近似(LDA)

局所的に、一様電子ガスのように振る舞うように近似する。

2. 密度勾配近似(GGA)

電子密度の勾配 $\delta n(\rm r)$ を考慮に近似である。

密度汎関数理論計算の実際

何より自己無撞着問題は収束することが大事である。さもなければ計算は一生、終わらない。

まず空間に周期的に繰り返されるスーパーセルは格子ベクトル $a_i,a_j, a_k$ で決めることができる。これに対して、逆格子空間ベクトルで、実空間で定義した基本胞と同じ様に、逆格子空間でも定義することができ、これはブリュアンゾーンと呼ばれている。

固体物理学の中では基本的ではあるが、あまり馴染みがない人が多いかもしれない。

ここで、密度汎関数理論では各k点(逆空間の点)で積分することが多いのだが、計算コストが高いため既約ブリュアンゾーンと呼ばれるものを使って、対称性をうまく使いながら行う。しかし、k点を選ぶでも限度があるのでどのくらいk点をとるのかを決めなればならない。

その他にも切断エネルギーと呼ばれるものあるが、これらは、またの機会に説明をする。

密度汎関数理論でできないことは何

ここまで第一原理計算と密度汎関数理論を述べてきた。

しかし、実際の密度汎関数理論は一様電子ガスの場合をのぞいて、近似解しか求めることができないのである。それはHohenberg-Kohnの定理でも述べた通りである。

また、ファンデルワールス力が重要になる様な条件では、計算が合わないことが多い。

次回

次回は、強誘電体について解説しようかなと思います。

Kekeの日記

エンジニア、読書なんでも