【水素分子で学ぶ】DFT法ではどのように量子計算をするのか

DFT法は最も利用されている量子化学計算の手法です。GaussianやGAMESSを使ったことのある方はDFT計算を利用したことが一度はあるかと思います。

しかし、DFT法では具体的にどのような計算をしているのかよくわからない、または、HF法との違いをイメージできないまま使用されている方も多いのではないでしょうか？

このページでは、広く使用されているコーン・シャムによるDFT法について解説し、水素分子を用いた具体的な計算例をご紹介します。

DFTは最も使われている手法ですが、恣意性があります

量子化学計算で最も利用されている計算手法はDFT法（密度汎関数理論）ですが、DFTでは電子間反発の取り扱いや、計算方法がHF（ハートリーフォック法）と異なります。

HF法では与えられた原子軌道（基底）を使って全ての計算を解析的に行います。計算の過程で恣意性が入ることはありません。しかし、DFT法には計算を行う際に恣意性が入ります。主に2点をご紹介します。

	2電子相互作用	量子計算する際の恣意性
HF	ガウス積分を利用した解析解	なし
DFT	グリッドを利用した数値積分	あり

【恣意性その1】グリッドの導入

DFTでは電子間反発エネルギーを計算する際にグリッドを利用した数値積分を行うため、まず、グリッドの細かさを選ぶという恣意性が入ります。

【恣意性その2】汎関数にあるパラメーター

B3LYPやPBE0など、一般によく利用される汎関数には複数のパラメータが存在します。

例えば、B3LYPでは電子密度$\rho$だけで表した汎関数を使うと、化学結合を強く見積もりすぎるため、ハートリー・フォックの交換積分を混ぜ合わせています。

実際にどのような方程式を解いているか

では、実際にDFTではどのような方程式を解いているのでしょうか？

HF法での解法

まずは、比較としてHF法を考えましょう。HF法では次のような方程式を解き、分子軌道$\boldsymbol{C}$や軌道エネルギー$\boldsymbol{E}$を得ます。

$\boldsymbol{FC}=\boldsymbol{SCE}$、$\boldsymbol{F}=\boldsymbol{H}+\boldsymbol{J}+\boldsymbol{K}$

これら一つ一つは$N\times{N}$の行列になっており、例えば、フォック行列は$\boldsymbol{F}=\begin{pmatrix}F_{11}&\cdots\\\vdots&\ddots\end{pmatrix}$のように表されます（$N$は基底の数です）。

HF法のフォック行列は1電子の運動エネルギーと核ポテンシャルを合わせた$\boldsymbol{H}$、2電子間のクーロン積分$\boldsymbol{J}$、交換積分$\boldsymbol{K}$で与えます。

DFTでの解法

DFT法も解く方程式はHF法と同じです。つまり、$\boldsymbol{FC}=\boldsymbol{SCE}$を解きます。

しかし、DFTではHFとは異なり、分子軌道から得られる電子密度$\boldsymbol{\rho}$（やその微分$\boldsymbol{\nabla\rho}$）によって、交換ポテンシャル$\boldsymbol{v}_{\text{xc}}$の計算を行います。そして、フォック行列を$\boldsymbol{F}=\boldsymbol{H}+\boldsymbol{J}+\boldsymbol{v}_{\text{xc}}$として分子軌道の計算を行います。

基底関数が$N$個あれば、ハミルトニアン$\boldsymbol{H}$などを$N\times{N}$の行列で表現しました。密度$\boldsymbol{\rho}$も同じように$N\times{N}$の行列で表現します。

例えば最も簡単な交換相関ポテンシャルは$\boldsymbol{v}_{\text{xc}}=-(\frac{3}{\pi}\boldsymbol{\rho})^{1/3}$とするもので、グリッドを利用して計算した$\boldsymbol{\rho}=\begin{pmatrix}\rho_{11}&\cdots\\\vdots&\ddots\end{pmatrix}$から$\boldsymbol{{v}_{\text{xc}}}=\begin{pmatrix}{v}^{\text{xc}}_{11}&\cdots\\\vdots&\ddots\end{pmatrix}$を計算して、$N\times{N}$のフォック演算子行列$\boldsymbol{F}=\boldsymbol{H}+\boldsymbol{J}+\boldsymbol{v}_{\text{xc}}$を作成します。

H2分子でのDFT計算例

実際に水素分子（原子間距離$1.4\text{bohr}$）のDFT計算を行ってみましょう。

基底関数の準備

基底関数は原点の水素原子H1の1s軌道$\chi_1=0.37e^{-0.41(x^2+y^2+z^2)}$、もう一方の水素原子H2の1s軌道$\chi_2=0.37e^{-0.41((x-1.4)^2+y^2+z^2)}$の2つ用意します。

グリッドの用意

DFT計算を行うためのグリッドを用意しましょう。

今回は水素原子の動径$r$方向は$0\sim1.0\text{bohr}$を$100$分割、偏角$\theta$方向は$0\sim\pi$を12分割、$\phi$方向は$0\sim2\pi$を$24$分割することにしましょう。

実際のDFT計算では動径方向の領域半径は原子ごとに変化させ（Bragg-Slater半径）、空間分割も滑らかに行うファジーセル法などが主流です。

電子間反発がないとして初期分子軌道を準備する

初期分子軌道を得るため、まずは水素分子の2電子間に相互作用が働かないとして量子化学計算を行いましょう。

詳細は割愛しますが、用意した基底を使うと1電子ハミルトニアンは$\boldsymbol{H}=\begin{pmatrix}-1.06&-0.90\\-0.90&-1.06\end{pmatrix}$、重なり積分は$\boldsymbol{S}=\begin{pmatrix}1.0&0.66\\0.66&1.0\end{pmatrix}$が得られるので、$\boldsymbol{HC}=\boldsymbol{SCE}$を解くことで初期分子軌道$\boldsymbol{C}=\begin{pmatrix}0.55&1.22\\0.55&-1.22\end{pmatrix}$を得ることができます。

分子軌道をもとに密度行列を数値積分で得る

DFT計算を行うため、密度行列$\boldsymbol{\rho}$を計算していきましょう。電子は分子軌道$\phi_1(\boldsymbol{r})=0.55*0.37(e^{-0.41(x^2+y^2+z^2)}$$+e^{-0.41((x-1.4)^2+y^2+z^2)})$に2つ入っているので、電子密度は$\rho(\boldsymbol{r})=2\phi_1(\boldsymbol{r})^2$で得ることができます。

例えば、水素原子H1から$(r,\theta,\phi)=(0.4,\frac{\pi}{2},0)$の位置にあるグリッド点$\boldsymbol{r}_1$を考えてみましょう。この点はxyz座標で$(x,y,z)=(0.4,0,0)$にあたるので、$\phi(\boldsymbol{r}_1)=0.32$、$\rho(\boldsymbol{r}_1)=0.32^2=0.10$になります。

一方、基底関数の値は$\chi_{\text{H}1}(\boldsymbol{r}_1)=0.34$、$\chi_{\text{H}2}(\boldsymbol{r}_1)=0.24$なので、グリッド点$\boldsymbol{r}_1$による密度行列$\rho_{12}$への寄与は$\chi_{\text{H}1}*\rho*\chi_{\text{H}2}*r^2\sin\theta\Delta{r}\Delta{\theta}\Delta{\phi}$$=0.34*0.10*0.24*0.4^2$$*\sin\frac{\pi}{2}*\frac{2}{100}*\frac{\pi}{12}*\frac{2\pi}{24}=1.78*10^{-6}$になります。

この計算を全グリッド点で行い、足し合わせると$\rho_{12}=0.0566$になります。ほかの行列成分も同じようにして計算することで$\boldsymbol{\rho}=\begin{pmatrix}0.0726&0.0566\\0.0566&0.0726\end{pmatrix}$が得られます。

交換相関ポテンシャル行列を計算する

$\boldsymbol{\rho}$の表現行列が得られたので、$\boldsymbol{v}_{\text{xc}}=-(\frac{3}{\pi}\boldsymbol{\rho})^{1/3}$から$\boldsymbol{v}_{\text{xc}}$を計算しましょう。

よくある間違いですが、$\boldsymbol{\rho}^{1/3}=\begin{pmatrix}{\rho}_{11}^{1/3}&{\rho}_{12}^{1/3}\\{\rho}_{21}^{1/3}&{\rho}_{22}^{1/3}\end{pmatrix}$ではありません。

$\boldsymbol{\rho}^{1/3}$は$\boldsymbol{\rho}$の対角化を経由して計算する必要があります。

具体的には$\boldsymbol{\rho}^{1/3}=\begin{pmatrix}0.3787&0.1268\\0.1268&0.3787\end{pmatrix}$です。

この$\boldsymbol{\rho}^{1/3}$を利用すると$\boldsymbol{v}_{\text{xc}}=-(\frac{3}{\pi}\boldsymbol{\rho})^{1/3}=\begin{pmatrix}-0.3730&-0.1249\\-0.1249&-0.3730\end{pmatrix}$と求めることができます。

新しい分子軌道を得て、エネルギーを計算する

得られた$\boldsymbol{v}_{\text{xc}}$を利用して新しいフォック演算子$\boldsymbol{F}=\boldsymbol{H}+\boldsymbol{J}+\boldsymbol{v}_{\text{xc}}$を作り、$\boldsymbol{FC}=\boldsymbol{SCE}$を解くことで新しい分子軌道$\boldsymbol{C}$と電子エネルギーを計算します

この手続きをエネルギーが安定になるまで繰り返すことでDFT法による分子軌道として$\boldsymbol{C}=\begin{pmatrix}0.55&1.22\\0.55&-1.22\end{pmatrix}$を得ることができます。

まとめ

このページでは、広く使用されているコーン・シャムによるDFT法について解説し、水素分子を用いた具体的な計算例をご紹介しました

DFTは最も利用されている量子化学計算の手法ですので、具体的な計算方法や、また、HF法との違いをイメージできるようにしておきましょう。