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

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の中身を変更してください。