量子化学計算のアルゴリズムを解説

2023年12月30日2024年11月18日

このページでは量子化学計算のアルゴリズムが分かるように説明します。1電子積分と2電子積分が求まっていれば、RHF法で分子軌道を得るアルゴリズムは意外にシンプルです。ここでは、Hartree-Fock方程式を具体的に解いていく手続きについて解説します。

解くべき方程式とアルゴリズム

早速ですが、Schrodinger方程式を化学分子に適応させたときに、解くべき方程式は次のHartree-Fock方程式です。

$\textbf{FC}=\textbf{SC}\epsilon$

各変数$\textbf{F,C,S,}\epsilon$は行列です。

この方程式を解くアルゴリズとして下図のようなものが考えられます。以下、順を追って説明します。

1電子積分と2電子積分の準備

量子化学計算を行う上では基底関数と呼ばれるものを準備する必要があります。基底関数は原子軌道$\varphi$のようなもので、分子軌道は原子軌道の線形和として表現できると考えます（$\phi=\sum_ic_i\varphi_i$）。基底関数の種類としては最も簡単なSTO-3Gから6-31G*、電子相関を取り込みやすいとされるcc-PVTZなど様々な種類のものが開発されています。

基底関数（原子軌道）の準備が整ったら、早速ですが、1電子積分と2電子積分の値を計算する必要があります。1電子積分や2電子積分の詳細は別ページをご覧ください。

初期値となる分子軌道の生成

一電子積分のみのハミルトニアン$\textbf{H}$計算

量子化学計算を行うためには初期値が必要です。この初期値として、一電子積分だけから計算できるハミルトニアン$\textbf{H}$を利用して分子軌道の初期値をつくることが多いです。ハミルトニアン$\textbf{H}$の各成分は次の形です。

$\textbf{H}_{\mu\nu}=(\mu|H|\nu)=(\mu|-\frac{1}{2}\nabla^2+\sum\frac{Z_A}{|\textbf{r}-\textbf{R}_A|}|\nu)$

このハミルトニアンは先ほどの1電子積分の結果を用いると各行列成分を計算することができます。

つまり、次の行列方程式を解けば初期値となる分子軌道$\textbf{C}$とエネルギー準位$\epsilon$を得ることが出来ます。

$\textbf{HC}=\textbf{SC}\epsilon$

行列の要素の形をあらわに書くと次のような$n\times{n}$の行列方程式になります。$n$は基底関数の個数に一致します。例えば水分子をSTO-3Gで計算すると、基底関数の個数は7つです（O：1s、2s、2px、2py、2pz　H：1sが2つ）。

$\left(
\begin{array}{cccc}
H_{11} & \ldots & H_{1n} \\
H_{21} & \ldots & H_{2n} \\
\vdots & \ddots & \vdots \\
H_{n1} & \ldots & H_{nn}
\end{array}
\right)$$\left(
\begin{array}{cccc}
C_{11} & \ldots & C_{1n} \\
C_{21} & \ldots & C_{2n} \\
\vdots & \ddots & \vdots \\
C_{n1} & \ldots & C_{nn}
\end{array}\right)$

$~~~~~~~=\left(
\begin{array}{cccc}
S_{11} & \ldots & S_{1n} \\
S_{21} & \ldots & S_{2n} \\
\vdots & \ddots & \vdots \\
S_{n1} & \ldots & S_{nn}
\end{array}\right)$$\left(
\begin{array}{cccc}
C_{11} & \ldots & C_{1n} \\
C_{21} & \ldots & C_{2n} \\
\vdots & \ddots & \vdots \\
C_{n1} & \ldots & C_{nn}
\end{array}\right)$$\left(
\begin{array}{cccc}
\epsilon_{11} & \ldots & \epsilon_{1n} \\
\epsilon_{21} & \ldots & \epsilon_{2n} \\
\vdots & \ddots & \vdots \\
\epsilon_{n1} & \ldots & \epsilon_{nn}
\end{array}\right)$

分子軌道$\textbf{C}$を得る手順

この方程式$\textbf{HC}=\textbf{SC}\epsilon$を解く方法を考えます。

まず重なり行列$\textbf{S}$を対角化する$\textbf{U}$を計算します。

$\textbf{U}^\dagger\textbf{SU}=\textbf{T}$

本ホームページのプログラムではこの計算を含め、対角化にはLAPACKを利用しています。
得られた$\textbf{U}$から変換行列${X}$を作ります。

$\textbf{X}=\textbf{UT}^\textrm{-1/2}$

$\textbf{X}^\dagger$と$\textbf{X}$を$\textbf{h}$にかけて、これを$\textbf{h}^\prime$とします。

$\textbf{h}^\prime=\textbf{X}^\dagger\textbf{hX}$

この$\textrm{h}^\prime$を\textbf{W}で対角化すると、実はこの固有値がエネルギー準位$\epsilon$です。$\textbf{h}^\prime$は対称行列（実数のエルミート行列）なので、$\textbf{W}は直交行列（実数のユニタリー行列）です（線形代数の知識です）。

$\textbf{W}^\dagger\textbf{h}^\prime\textbf{W}=\epsilon$

分子軌道$\textbf{C}$は対角化行列$\textbf{W}$と変換行列$\textbf{X}$の積で得ることが出来ます。

$\textbf{C}=\textbf{XW}$

$\textbf{C}=\textbf{XW}$が$\textbf{HC}=\textbf{SC}\epsilon$を満たす証明。

$\textbf{C}=\textbf{XW}$を$\textbf{HC}=\textbf{SC}\epsilon$に代入します。

$\textbf{HXW}=\textbf{SXW}\epsilon$

左右から$\textbf{X}^\dagger$と$\textbf{W}^\dagger$をかけます。

(左辺) =$\textbf{W}^\dagger\textbf{X}^\dagger\textbf{HXW}=\textbf{W}^\dagger\textbf{h}^\prime\textbf{W}=\epsilon$

(右辺)
$=\textbf{W}^\dagger\textbf{X}^\dagger\textbf{SXW}\epsilon$
$=\textbf{W}^\dagger\textbf{X}^\dagger\textbf{SUT}^{-1/2}\textbf{W}\epsilon$
$=\textbf{W}^\dagger\textbf{T}^{-1/2}\textbf{U}^\dagger\textbf{SUT}^{-1/2}\textbf{W}\epsilon$
$=\textbf{W}^\dagger\textbf{T}^{-1/2}\textbf{TT}^{-1/2}\textbf{W}\epsilon$
$=\textbf{W}^\dagger\textbf{W}\epsilon=\epsilon$

よって、$\textbf{C}=\textbf{XW}$とすることで、量子化学方程式$\textbf{HC}=\textbf{SC}\epsilon$の解とできる。

いよいよ繰り返し計算開始

初期条件としての$\textbf{C}$を得た後はいよいよ繰り返し計算（SCF計算）を行っていきます。

分子軌道$\textbf{C}$が与えられたら、密度行列$\textbf{P}$として$P_{\mu\nu}=2\sum{C}_i^{\mu}C_i^\nu$を導入することで、系のエネルギーと電子間相互作用も加味したハミルトニアン（Fock演算子）を計算することができます。密度行列のsumはすべての占有軌道の分子軌道に対してとります。