【プログラミングに向けて】量子化学計算のアルゴリズムを解説

量子化学計算で一般的なものに閉殻のハートリー・フォック法（RHF法）があります。

このRHFは最も簡単な量子化学計算であるものの、基底関数の準備や電子積分、繰り返し計算による分子軌道の更新など、量子化学計算の基本的なアルゴリズムが詰まっています。

また、量子化学計算のスタンダードであるDFT法のアルゴリズムも電子間相互作用の取り扱い以外はRHF法とほぼ同じなので、DFT法に興味がある方も量子化学の基本として内容を押さえておいてください。

量子化学の方程式とそれを解くアルゴリズム

Schrodinger方程式を分子に適応させたときに、解くべき方程式はハートリー・フォック方程式$\boldsymbol{FC}=\boldsymbol{SC\epsilon}$です。

各変数$\boldsymbol{F}$や$\boldsymbol{C}$などは行列です。$\boldsymbol{F}$はフォック行列と呼ばれ、分子が持つ電子系のエネルギーを表現します。$\boldsymbol{C}$は分子軌道、$\boldsymbol{S}$は重なり積分、$\boldsymbol{\epsilon}$は軌道エネルギーです。

この方程式を解くアルゴリズムとして次のようなものが一般的です。

このアルゴリズムに則れば簡単な分子の場合、200行ほどプログラミングすることで量子化学計算を行うことができます。

【ステップ１】電子積分の準備

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

基底関数の選択が終わったら、1電子積分（重なり積分、運動エネルギー積分、核引力積分）と2電子積分を計算する必要があります。この電子積分はガウス関数を利用している場合（STO-3Gや6-31G*など）、簡単に解析解を求めることができます。

【ステップ２】初期値となる分子軌道の生成

量子化学に限らず、シミュレーションを行う上では、初期値を設定する必要があります。量子化学のシミュレーションでは初期値として、一電子積分から計算できるハミルトニアン$\boldsymbol{H}$を利用して分子軌道の初期値をつくることが多いです。

ハミルトニアンは既に計算している運動エネルギー行列$\boldsymbol{T}$と核引力行列$\boldsymbol{V}$の和$\boldsymbol{H}=\boldsymbol{T}+\boldsymbol{V}$で与えられ、

行列方程式$\boldsymbol{H}\boldsymbol{C}=\boldsymbol{S}\boldsymbol{C}\epsilon$を解けば初期値となる分子軌道$\boldsymbol{C}$を得ることができます。

この方程式の解となる分子軌道$\boldsymbol{C}$は$\boldsymbol{H}^\prime\boldsymbol{W}=\boldsymbol{W}\epsilon$を満たす固有ベクトル$\boldsymbol{W}$を用いると$\boldsymbol{C}=\boldsymbol{XW}$で得ることができます。（ただし、$\boldsymbol{U}^\dagger\boldsymbol{SU}=\boldsymbol{M}$(対角行列)を満たす$\boldsymbol{U}$を使って、$\boldsymbol{X}=\boldsymbol{UM}^{-\frac{1}{2}}$と$\boldsymbol{H}^\prime=\boldsymbol{X}^\dagger\boldsymbol{HX}$とします。）

【ステップ３】繰り返し計算開始

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

分子軌道を（初期値とはいえ）得られたことで電子間の相互作用を加味したフォック行列$\boldsymbol{F}$と電子状態のエネルギー$E$を計算できます。

フォック行列$\boldsymbol{F}$は占有軌道$\phi_i=\sum^\text{NBasis}_{\mu}C_\mu^i\chi_\mu$から得られる密度行列$P_{\lambda\sigma}=\sum_i^{\text{occ}}C_\lambda^iC_\sigma^i$を利用することで$F_{\mu\nu}=H_{\mu\nu}$$+\sum_{\lambda\sigma}P_{\lambda\sigma}\{(\mu\nu|\sigma\lambda)-\frac{1}{2}(\mu\lambda|\sigma\nu)\}$で計算でき、$E=\frac{1}{2}\sum_{\mu\nu}^\text{Nbasis}P_{\mu\nu}(H_{\nu\mu}+F_{\nu\mu})$によって電子状態のエネルギーを計算することができます。

RHF計算ではこの電子状態のエネルギー$E$が変化しなくなるまで繰り返します。繰り返し計算開始直後（ループ1周目）では、このエネルギー$E$の変化量$\Delta{E}$はまだわからないので、計算されているフォック行列$\boldsymbol{F}$を利用して、新しい分子軌道の組$\boldsymbol{C}$を計算しましょう。

【ステップ４】新しい分子軌道を計算する

電子状態のエネルギー$E$がまだ下がりそうな場合は、フォック行列$\boldsymbol{F}$を利用して新しい分子軌道$\boldsymbol{C}$を$\boldsymbol{F}\boldsymbol{C}=\boldsymbol{S}\boldsymbol{C}\epsilon$を解くことで得ます。

これは先ほど$\boldsymbol{H}\boldsymbol{C}=\boldsymbol{S}\boldsymbol{C}\epsilon$を解いた時と同じように、$\boldsymbol{F}^\prime\boldsymbol{W}=\boldsymbol{W}\epsilon$を満たす固有ベクトル$\boldsymbol{W}$を用いると$\boldsymbol{C}=\boldsymbol{XW}$で得ることができます。（ただし、$\boldsymbol{F}^\prime=\boldsymbol{X}^\dagger\boldsymbol{FX}$で、$\boldsymbol{X}$は先ほどと同じ定義です。）

SCF計算ではこの新しい分子軌道$\boldsymbol{C}$をもとに密度行列$\boldsymbol{P}$を再度作り、電子状態のエネルギー$E$を計算します。この手続きを複数繰り返すことでエネルギーの低い電子状態を与える分子軌道を探すのです。

まとめ

このページでは、量子力学の基本であるハートリー・フォック方程式（RHF法）のアルゴリズムについて解説しました。

このアルゴリズムに則れば簡単な分子の場合、200行ほどプログラミングすることで量子化学計算を行うことができますので、皆さんも一度チャレンジしてみてください。