量子化学計算コード入門 (1): HeH+のHartree-Fock計算

目的

• 量子化学計算プログラムのコーディングを始めてみたい方のアシスト
• 計算化学業界のサステナビリティを指向し、技術伝承のための教育資料を作成

本稿の目標

• 「ザボの教科書」¹に記載されているシンプルな計算例 (HeH$^+$におけるRHF/STO-3G計算) に基づき、Hartee-Fock計算のエッセンスを理解
• 「ザボの教科書」に記載されているFORTRANコードをPythonコードに焼き直す

はじめに

(この項目は、量子化学研究におけるバックグラウンド的な内容を含みます。不要な方は読み飛ばしてください。)

国内の量子化学系の研究室に配属された学生さんは、大抵の場合、ザボ・オストランド著『新しい量子化学 上』 (業界では「ザボの教科書」と呼ばれる) という書籍を使って量子化学理論の基礎を勉強することになります。
「ザボの教科書」は現代の量子化学業界のバイブルのようなものです ²。

量子化学業界の研究内容は大きく二分されます。
1. 特定の分子や現象を理解するための応用的研究
2. 1の研究をするための理論・プログラム開発
業界にはGaussian³やGAMESS⁴など、汎用的な量子化学計算プログラムが既に存在するため、大抵は1の研究に取り組むケースが多いです。
しかしながら、1万分子を超えるような巨大系のダイナミクスを捉えたい、多光子吸収過程を緻密に理解したい、重元素化合物のNMRスペクトルを解析したい…などマニアックな研究に興味がある場合は、既存のプログラムで検討できない場合があります。そのときに、2の研究が必要になります。

2の理論・プログラム開発をする場合、先行研究で開発された理論や計算プログラムの拡張をするところから研究を始めることが多いと思います。
しかしながら、このようなプログラムは長きに渡って、秘伝のタレ的に継ぎ足しで作られており、コードの構成や可読性が決して良くない場合が多いです。
初学者にとっては独特のとっつきにくさがあるように感じます。このとっつきにくさを緩和させることが本稿の一つの目的であります。
決して、研究室配属された学生さんが楽をすることを目的としているわけではなく、導入を平易にして、かつ、基礎をしっかりと理解して、理論構築研究の本質的なところにより邁進していって欲しいという思いを持って執筆しています。

ところで、量子化学計算プログラムの多くは慣例的にFORTRAN言語で書かれていることが多いです。コーディングにおける独特な記法も多く存在します。
「ザボの教科書」P.265~に記載されているHatree-Fock計算のサンプルプログラムもFORTRANベースで書かれています。
FORTRAN言語はプログラミング初心者でも比較的理解しやすいと思いますが、時代の潮流に倣い、先のサンプルプログラムをPythonベースのコードに焼き直しつつ、その行間を解説していきたいと思います。

本文中に書ききれない「行間」(もとい、関連知識や慣習など) は脚注にまではみ出して付記しています。併せてご覧いただければと思います。

実装内容

• 対象: HeH$^+$分子 (最もシンプルな異核2原子分子)
• そもそもの仮定 (ハミルトニアン): Born-Oppenheimer近似下での非相対論的な取り扱い
• 計算方法 (波動関数理論): 制限Hartree-Fock (Restricted Hartree-Fock; RHF)
• 基底関数: STO-3G
• 入力パラメーター (ユーザーが指定するべきもの)
  • 分子構造: 核間距離 $R = 1.463\ \text{a.u.}$
  • 原子核の種類 (本質的には原子番号ないし核電荷Zを指定): 原子A = He ($Z_A = 2.0$), 原子B = H ($Z_B = 1.0$)
  • 電荷: $+1$
• 原子軌道 (AO) 積分の計算は簡単のためスキップ

電荷について、ザボの教科書のコード中では明示的に言及されていませんが⁵、量子化学計算を実行する際は必ず指定が必要になります。
原子核AがHe ($Z_A = 2.0$)、原子核BがH ($Z_B = 1.0$)であることから、電荷的に中性であるHeH分子の場合、電子数は $3$ です。これに対して、今回の対象系はHeH$^+$分子であることから、系全体の電荷は $+1$ です。したがって、HeH$^+$分子においては、電子数 $N_\text{elec} = 2$ です。
RHFでは制限スピン軌道を採用しており、共通した1種類の軌道に対して電子を2個ずつ詰めていくことを仮定します。すなわち、ユニークな占有軌道 (occupied orbitals) の数は $N_\text{occ} = 1$ となります。

多くの量子化学計算プログラムでは、計算の初期段階において、インプットファイルから計算する分子に含まれる元素種、構造、電荷、スピン多重度などを読み込み、上述した電子数 $N_\text{elec}$ や占有軌道数 $N_\text{occ}$ などの基本的な情報を変数に格納するところからコードが始まります。

また、今回実装するプログラムでは、簡単のためAO積分の計算をスキップし、計算対象であるHeH$^+$に対応する積分値を読み込むだけとします。これは以下の二つの理由に依ります。
1. 本稿ではあくまでもHF計算の流れを理解することに主眼を置くため。
2. 実用的な積分計算アルゴリズムは非常に複雑であり、理論開発系の研究においても大抵の場合、既存のルーチンで計算した積分値を読み込むだけのケースが多いため。

計算のフロー

1. 分子情報・計算設定の読み込み
2. AO積分の計算
3. 初期電子密度の読み込み
4. SCFループ開始
5. Fock行列の計算
6. HFエネルギーの計算
7. Fock行列の対角化
8. 密度行列の計算
9. 収束判定 (収束していれば10へ; 未だであれば4に戻る)
10. ポスト処理 (計算結果の保存など)

プログラムの解説

ここでは上述した計算のフローに従い、Pythonコードに実装した内容を順に解説していきます。
本コードでは今後の拡張も見据えて、ザボの教科書のサンプルコードよりも多少汎用性の高い (= 他の系にも展開しやすいような) 実装をしています。また、FORTRANの実装と異なり、クラスを使っています。
コードの全文はAppendixとして本稿の末尾に付記します。

メインルーチン

本コードでは、szaboというクラスを作り、この中に種々の操作を入れ込みます。
汎用量子化学計算ソフトのようにインプットファイルを読み込むのではなく、szaboクラスの引数として分子情報や計算設定を指定することにします。
また、積分計算とRHF計算は、それぞれ対応するメソッド (integralsおよびRHF) を明示的に呼び出します。

calc1 = szabo(r=1.4632, zeta1=2.0925, zeta2=1.2400, za=2.0, zb=1.0, doDebugPrint=False)
calc1.integrals()
calc1.RHF()

分子情報・計算設定の読み込み

このコードで初めに実行されるのは、szaboクラスのコンストラクタ__init__内の処理です。

class szabo():

    def __init__(self, r, zeta1, zeta2, za, zb, doDebugPrint):
        # Settings from input
        self.r = r  # Bond length
        self.zeta1 = zeta1  # Exponent of basis function 1
        self.zeta2 = zeta2  # Exponent of basis function 2
        self.za = za  # Atomic number of atom A
        self.zb = zb  # Atomic number of atom B
        self.doDebugPrint = doDebugPrint
        self.numElec = 2  # Number of electrons
        self.numOcc = self.numElec // 2  # Number of occupied orbitals
        self.numAO = 2  # Number of AOs
        self.occupancy = np.zeros(self.numAO)
        self.occupancy[:self.numOcc] = 2.0  # Here, only the first orbital is occupied by two electrons

引数になっている、計算の入力パラメーターは以下の6種類です。
- r: 核間距離 (原子単位で指定)
- z1, z2: 基底関数1 or 2の軌道指数 (本実装では不使用)
- za, zb: 原子A or Bの核電荷 (原子番号)
- doDebugPrint: 積分値・行列値を出力するかどうかのフラグ

量子化学計算の重要なパラメーターである電子数・軌道に関する情報もここで指定しておきます。
- numElec: 電子数 ($N_\text{elec}$)
- numOcc: 占有軌道数 ($N_\text{occ}$)
- numAO: AOの次元数 ($N_\text{AO}$)
- occupancy: 占有数

今回の計算は、2つのSTO-3G関数を使った、2電子系のRHF計算であるため、上記の設定となります
ザボのコードに記述されていない情報として、ここでは占有数occupancyの指定も一緒にしておきます。
occupancyはnumAO次元⁶のベクトルであり、占有軌道に対応する要素にはその占有数、仮想軌道に対応する要素にはゼロを代入しています。
HF計算では、MO係数から密度行列を構築する際に、この情報が必要となります (詳細は後述)。

AO積分の計算

積分計算 (ここでは計算せず、計算値を読み込むのみ) には、integralsメソッドを呼びます。
積分には1電子積分と2電子積分があるので、それぞれを独立したメソッド (one_electron_integralsおよびtwo_electron_integrals) にしています。

def integrals(self):
    self.one_electron_integrals()
    self.two_electron_integrals()

1電子積分

one_electron_integralsメソッドでは、まず、各行列に対応する変数を宣言 (メモリ領域を確保) します。
- 1電子ハミルトニアン⁷ $\textbf{H}^\text{core}$: H
- 重なり積分 $\textbf{S}$: S
- 正準直交化のための行列 $\textbf{X}$⁸: X

def one_electron_integrals(self):
    self.H = np.zeros((self.numAO, self.numAO))  # Core Hamiltonian
    self.S = np.zeros((self.numAO, self.numAO))  # Overlap matrix
    self.X = np.zeros((self.numAO, self.numAO))  # X-matrix for canonical orthogonalization

宣言した行列に対して、種々の行列要素を計算・代入していきます。
ここで、Sxxは重なり積分、Txxは運動エネルギー積分、Vxxxは核-電子引力積分の要素を示しています。

# OEI values of HeH+ (R = 1.4632 a.u.) with STO-3G
S12 = 0.450770e+0
T11 = 2.164313e+0
T12 = 0.167013e+0
T22 = 0.760033e+0
V11A = -4.139827e+0
V12A = -1.102912e+0
V22A = -1.265246e+0
V11B = -0.677230e+0
V12B = -0.411305e+0
V22B = -1.226615e+0

# Form core Hamiltonian
self.H[0][0] = T11 + V11A + V11B
self.H[0][1] = T12 + V12A + V12B
self.H[1][0] = self.H[0][1]
self.H[1][1] = T22 + V22A + V22B

# Form overlap matrix
self.S[0][0] = 1.0
self.S[0][1] = S12
self.S[1][0] = S12
self.S[1][1] = 1.0

# Form X-matrix for canonical orthogonalization
self.X[0][0] = 1.0 / np.sqrt(2.0 * (1.0 + S12))
self.X[1][0] = self.X[0][0]
self.X[0][1] = 1.0 / np.sqrt(2.0 * (1.0 - S12))
self.X[1][1] = -self.X[0][1]

量子化学計算で扱う多くの行列 (正確には何らかの物理量に対応する行列) は、基本的にエルミート行列 $\textbf{A} = \textbf{A}^\dagger$ になっていることを心に留めておきましょう⁹。
ここでは実数の範囲しか扱わないので、対称行列 $\textbf{A} = \textbf{A}^\text{T}$ となり、A[1][0] = A[0][1] です。上記コード中でも、HとSはこの仕様を満たしています。
なお、行列 $\textbf{X}$ のような複数の「足」を繋ぐ変換行列については、この限りではないので注意してください¹⁰。

ついでに、フラグdoDebugPrintがTrueであれば、デバッグプリントとして各行列要素の値を露に出力することにします。

# Print out integral values
if self.doDebugPrint:
    print('The S array \n{}\n'.format(self.S))
    print('The X array \n{}\n'.format(self.X))
    print('The H array \n{}\n'.format(self.H))

2電子積分

2電子積分についても1電子積分と同様の流れで実装します。
two_electron_integralsメソッドでは、まず、配列の宣言から始めます。

def two_electron_integrals(self):
    self.TT = np.zeros((self.numAO, self.numAO, self.numAO, self.numAO))

ここではAO次元数numAOがたかだか2であり、2電子積分も全部で $2^4=16$ 個しか存在しないので、コードの簡素化・理解の平易化のためにメモリをケチらず、愚直に4インデックスの配列を宣言しています。
しかしながら、2電子積分の総数は $O(N_\text{AO}^4)$ なので、AO次元数が増加するとすぐにメモリがパンクしてしまいます。したがって、実際の汎用量子化学計算プログラムでは、
- 原理的に値が同じになる要素は計算・保存しない (詳細は本節にて後ほど補足)
- 複数のインデックスを一つにまとめて管理
- そもそも2電子積分をまとめて計算・保存せずに、個々のFock行列要素を計算するときに必要となる2電子積分を都度計算する¹¹
などのテクニックが用いられます。

続いて、各積分値を読み込みます。

# TEI values of HeH+ (R = 1.4632 a.u.) with STO-3G
V1111 = 1.307152e+0
V2111 = 0.437279e+0
V2121 = 0.177267e+0
V2211 = 0.605703e+0
V2221 = 0.311795e+0
V2222 = 0.774608e+0

# Form matrix of two-electron integrals
self.TT[0][0][0][0] = V1111
self.TT[1][0][0][0] = V2111
self.TT[0][1][0][0] = V2111
self.TT[0][0][1][0] = V2111
self.TT[0][0][0][1] = V2111
self.TT[1][0][1][0] = V2121
self.TT[0][1][1][0] = V2121
self.TT[1][0][0][1] = V2121
self.TT[0][1][0][1] = V2121
self.TT[1][1][0][0] = V2211
self.TT[0][0][1][1] = V2211
self.TT[1][1][1][0] = V2221
self.TT[1][1][0][1] = V2221
self.TT[1][0][1][1] = V2221
self.TT[0][1][1][1] = V2221
self.TT[1][1][1][1] = V2222

上述した通り、$N_\text{AO}^4 = 2^4 = 16$ 個の要素を代入しています。
ただし、ここでユニークな積分値が6個しかないことには留意しておきましょう。
上述した通り、2電子積分 (正確には電子反発積分¹²) のうち、いくつかは原理的に同じ値になります。2電子積分にはインデックスの交換に対する対称性が存在し、いわゆる化学者の記法において、
$$(\mu \nu | \lambda \sigma) = (\mu \nu | \sigma \lambda) = (\nu \mu | \lambda \sigma) = (\nu \mu | \lambda \sigma) = (\sigma \lambda| \mu \nu) = (\sigma \lambda| \nu \mu) = (\lambda \sigma| \mu \nu) = (\lambda \sigma| \nu \mu)$$
となります。

ついでに、フラグdoDebugPrintに応じて、積分値を4重ループで画面に出力します。

# Print out integral values
if self.doDebugPrint:
    print('The TT array')
    for i in range(self.numAO):
        for j in range(self.numAO):
            for k in range(self.numAO):
                for l in range(self.numAO):
                    print('({}{}|{}{}) = {}'.format(i+1, j+1, k+1, l+1, self.TT[i][j][k][l]))
    print('')

SCF初期設定

いよいよRHFメソッドに入ります。
初めにSCFループの最大ステップ数maxIterと、密度行列に関する収束判定値tolScfConvを設定しておきます。

def RHF(self):
    # Settings
    self.maxIter = 25  # Maximum number of iterations
    self.tolScfConv = 1.0e-4  # Convergence criterion for density matrix

初期電子密度の読み込み

SCFループに入る前に、まずは密度行列 $\textbf{P}$¹³ の初期値を決めておきます。
ここで良い初期値を入力しておくと、SCF計算が速く収束することがあります。
今回は実装を簡単にするために、初期電子密度として、$\textbf{P}$ を零行列としておきます。
Fock行列要素 $F_{\mu \nu}$は、下式のように表せます。

\begin{align}
F_{\mu \nu} &= H_{\mu \nu}^\text{core} + G_{\mu \nu} \\
&= H_{\mu \nu}^\text{core} + \sum_{\lambda \sigma} P_{\lambda \sigma} \left[(\mu \nu|\sigma \lambda) - \frac{1}{2}(\mu \lambda|\sigma \nu) \right]
\end{align}

したがって、$\textbf{P} = \textbf{0}$ であれば、Fock行列 $\textbf{F}$ は、1電子ハミルトニアン $\textbf{H}^\text{core}$ に一致します。
このようにして初期軌道を構築する手法のことを、業界では「initial guessとしてcore Hamiltonianを使う」というように表現します。

さて、密度行列に対応する配列 P を確保して、ゼロで初期化します。
併せて、SCF計算に必要なその他の配列も確保しておきます。

# Clear arrays
P = np.zeros((self.numAO, self.numAO))  # Use core-Hamiltonian for initial guess (P = 0)
G = np.zeros((self.numAO, self.numAO))  # Two-electron part of Fock matrix
F = np.zeros((self.numAO, self.numAO))  # Fock matrix
Fprime = np.zeros((self.numAO, self.numAO))  # Transformed Fock matrix
E = np.zeros(self.numAO)  # Orbital energy matrix
C = np.zeros((self.numAO, self.numAO))  # MO coefficient marix
Cprime = np.zeros((self.numAO, self.numAO))  # MO coefficient marix in the orthonormal basis
oldP = np.zeros((self.numAO, self.numAO))  # Density matrix in the previous step

SCFループ開始

SCFループを開始します。
計算が収束するまで (あるいはカウンタiterがmaxIterに到達するまで)、後続の手順を繰り返します。

# SCF loop
for iter in range(self.maxIter):

Fock行列の計算

Fock行列の構築は、SCF前半部の肝となる手続きです。
先に示した式に従って、Fock行列要素を計算します。
まずは、Fock行列の2電子パートGを計算します。

# Form G
for i in range(self.numAO):
    for j in range(self.numAO):
        G[i][j] = 0.0
        for k in range(self.numAO):
            for l in range(self.numAO):
                G[i][j] += P[k][l] * (self.TT[i][j][k][l] - 0.5 * self.TT[i][l][k][j])

Gの計算後、先に構築しておいた1電子パートHと足し合わせ、Fに代入します。

# Form Fock matrix
F = self.H + G

HFエネルギーの計算

いよいよHF計算の出力値であるHFエネルギー $E_\text{HF}$ を計算します。

$$ E_\text{HF} = \frac{1}{2} \sum_{\mu \nu} P_{\mu \nu} \left( H_{\mu \nu} + F_{\mu \nu} \right) $$

この式を愚直に実装すると、以下の通りです。
行列の2つのインデックスに関する2重ループにより、各行列要素の足し算・掛け算をしていきます。

# Calculate energy
self.engHF = 0.0
for i in range(self.numAO):
    for j in range(self.numAO):
        self.engHF += 0.5 * P[i][j] * (self.H[i][j] + F[i][j])

よりPython的な実装をすると、下記のように表現できます。
行列の次元が大きくなると、Pythonの場合、行列要素をfor文で回す上記のコードよりも、組み込み関数を使う下記のコードの方が大分計算が速くなるようです。

self.engHF = (P * (self.H + F)).sum() * 0.5

Fock行列の対角化

SCF後半部の肝が、この対角化の手続きです。
まずは、変換行列Xを使って、AO基底で表されたFを正規直交基底に変換します。
式としては、
$$\textbf{F}^\prime = \textbf{X}^\dagger \textbf{F} \textbf{X} = \textbf{X}^\text{T} \textbf{F} \textbf{X}$$
となります。

# Transform Fock matrix to the orthonormal basis
Fprime = np.dot(self.X.T, np.dot(F, self.X))

ザボの教科書では3重ループで行列積を計算するサブルーチンをコーディングしていますが、ここではNumPyのdot関数を使って行列積を計算しました。
2次元の行列積ではあまり気になりませんが、100次元くらいから、計算時間の差が明白に見えてくると思います。
汎用的な量子化学計算コードでも、行列積などの基本的な計算は出来合いのサブルーチンを使うことが普通です。
FORTRANやCで書かれたプログラムの場合、BLAS (Basic Linear Algebra Subprograms) という線形代数ライブラリに実装されているサブルーチン (行列積の場合はdgemmなど) を使うことが一般的です。

続いて、正規直交基底に変換されたFock行列Fprimeを対角化します。
対角化も出来合いのプログラムを使うことが一般的です。
FORTRANやCの場合は、LAPACK (Linear Algebra PACKage) というライブラリに実装されたdsygvなどのサブルーチンを使って対角化をします。
ここでは、NumPyのlinalg.eig関数を使って対角化をします。
linalg.eig関数の出力値として、軌道エネルギー ${\bf \varepsilon}$ に対応する固有値を配列Eに、正規直交基底で表されたMO係数行列 $\textbf{C}^\prime$ に対応する固有ベクトルを配列Cprimeに格納します。

# Diagonalize transormed Fock matrix
E, Cprime = np.linalg.eig(Fprime)

対角化で得られる固有値・固有ベクトルは順番が滅茶苦茶な可能性があるので¹⁴、軌道エネルギーの昇順にソートをしておきます。
ここでも、ザボの教科書の実装とは異なり、NumPyのargsort関数を使っています。

# Sort eigenvalues & eigenvectors
eigen_id = np.argsort(E)
E = E[eigen_id]
Cprime = Cprime[:,eigen_id]

対角化のステップにおける最後の処理は、固有ベクトルの逆変換です。

$$ \textbf{C} = \textbf{X}^\dagger \textbf{C}^\prime = \textbf{X}^\text{T} \textbf{C}^\prime$$

Fock行列の対角化は行列 $\textbf{X}$ を用いて正規直交基底において行いました。
つまり、得られる固有ベクトル $\textbf{C}^\prime$ も正規直交基底での表現 (言い換えれば、$\textbf{C}^\prime$ は正規直交基底 → MO基底の変換行列) となります。
そこで、再び $\textbf{X}$ (AO基底 → 正規直交基底の変換行列) を用いて $\textbf{C}^\prime$ を逆変換し、AO基底 → MO基底を繋ぐように (本来のMO係数 $\textbf{C}$ の形に) 直してあげます。

# Back-transform eigenvectors to the AO basis
C = np.dot(self.X, Cprime)

ここでもNumPyのdot関数を使って行列積の計算をしました。

密度行列の計算

ザボの教科書では、RHFにおける密度行列の計算式は
$$ P_{\mu \nu} = 2 \sum_i^\text{{occ}} C_{\mu i} C_{\nu i}^\ast $$
と表されていました。(ここで、${ \text{occ} }$ は占有軌道の集合を表すものとします。)
上式において、"$2$" という係数は、「RHFにおける各占有軌道は2個ずつ電子が占められている」という仮定に基づくものです。
したがって、占有数 $\textbf{n}$ を用いることで、密度行列の計算式をより一般化することができます。$\textbf{n}$ はMO次元のベクトルであり、占有軌道に対してはその占有数が、仮想軌道に対しては0が代入されているものとします。

$$ P_{\mu \nu} = \sum_i^\text{{occ}} C_{\mu i} n_i C_{\nu i}^\ast $$

この一般化された計算式を使うことで、RHF以外のHF計算 (UHFなど) や非整数占有数を用いた計算にも対応することができます。
本稿ではこの一般化された計算式にしたがって、密度行列の計算を実装してみます。
また、ザボの教科書における3重ループによる計算ではなく、Pythonらしくブール配列によるマスクとNumPyのdot関数を利用した実装を試みます。

# Form new density matrix
oldP = np.copy(P)  # Save the present density matrix before creating new one
mask = self.occupancy>0
P = np.dot(C[:,mask] * self.occupancy[mask], C[:,mask].T)

収束判定

SCFループにおける最終ステップとして、収束判定をします。
収束判定条件は色々あります。一般に、
1. 前ステップと現ステップの密度行列の差が閾値以下かどうか
2. 前ステップと現ステップのHFエネルギーの差が閾値以下かどうか
3. Brillouin条件 ($\textbf{F}\textbf{D}\textbf{S} - \textbf{S}\textbf{D}\textbf{F} = \textbf{0}$) を満たすかどうか
…などが収束の判断基準として用いられます。
ここではザボの教科書と同様に、1の密度行列を使った条件を実装します。
実装した内容を正確に書くと、「前ステップと現ステップにおける各密度行列要素の差の2乗和のルートの1/2」を計算し、deltaに代入しています。

# Calculate delta
delta = np.sqrt(((P - oldP) ** 2).sum()) * 0.5

最後に計算したdeltaが閾値tolScfConv未満であるかどうかチェックします。
見事、deltaが閾値を下回れば、収束メッセージを出してループを抜けます。
もし何らかのエラーによりSCFステップ数がmaxIter = 25を上回るようであれば、エラーメッセージを出してループを抜けます。

# Check for convergence
if delta < self.tolScfConv:
    print('\nSCF calculation converged!')
    break
elif iter + 1 >= self.maxIter:
    print('\nNo convergence in SCF...')

ポスト処理 (計算結果の保存など)

SCFループを抜け出した後は、計算結果をテキストファイルに出力したり、得られた各行列をストレージに出力したりします。
ここでは、画面にHFエネルギー $E_\text{HF}$ (engHF) と全エネルギー $E_\text{tot}$ (engTot) を出力するのみとします。
全エネルギーは先に求めたHFエネルギーに、核間反発エネルギーを足し合わせたものです。
ここではBorn-Oppenheimer近似のもと、原子核を固定された点電荷として扱っていることから、高校生的なCoulombの法則の公式に従って計算するのみです。
$$E_\text{tot} = E_\text{HF} + \frac{Z_A Z_B}{R_{AB}}$$

# Calculate nuclear repulsion energy
self.engNuc = self.za * self.zb / self.r
self.engTot = self.engHF + self.engNuc
print('-- Electronic energy = {} hartrees'.format(self.engHF))
print('-- Total energy      = {} hartrees'.format(self.engTot))

実行結果

6ステップで全エネルギー-2.860660 hartreeに収束しました。

Step #0: energy = 0.0, delta = 0.8828667658215907
Step #1: energy = -4.14186137727303, delta = 0.2791765149586448
Step #2: energy = -4.226490444584946, delta = 0.029661524244687706
Step #3: energy = -4.227521451919089, delta = 0.0023182354979531866
Step #4: energy = -4.227527794398263, delta = 0.00017439175521713346
Step #5: energy = -4.227527830309551, delta = 1.3078900094163124e-05

SCF calculation converged!
-- Electronic energy = -4.227527830309551 hartrees
-- Total energy      = -2.8606606897956093 hartrees

まとめ

ザボの教科書を参考に、量子化学計算の基本となるRHFのプログラムをPythonで作成しました。
このコードは結合長 1.4632 a.u. のHeH$^+$分子のSTO-3G計算を想定したものです。
他の系に適用したい場合は、電荷等のパラメーターや積分値を適切なものに変える必要があります。

今後の展望

どこまで記事を書けるかわかりませんが、興味のあることを列挙しておきます。
- 本稿の計算結果の解釈
- (今回スキップした) 積分値の計算
- UHFなど、他のHF計算への拡張
- post-HF法への展開
- DFTへの展開
- 構造最適化
- 各種property計算 (電荷解析など)

Appendix: コード全体

上記パーツを全てつなぎ合わせたコードを付記します。
ご参考になれば幸いです。

###############################################################
# Minimal basis STO-3G calculation on HeH+                    #
# From "Modern Quantum Chemistry" by A. Szabo & N. S. Ostlund #
# Coded by H. Hatori (@hatori_hoku) in 2019                   #
###############################################################

# Modules
import numpy as np

################################################################################

class szabo():

    def __init__(self, r, zeta1, zeta2, za, zb, doDebugPrint):
        # Settings from input
        self.r = r  # Bond length
        self.zeta1 = zeta1  # Exponent of basis function 1
        self.zeta2 = zeta2  # Exponent of basis function 2
        self.za = za  # Atomic number of atom A
        self.zb = zb  # Atomic number of atom B
        self.doDebugPrint = doDebugPrint
        self.numElec = 2  # Number of electrons
        self.numOcc = self.numElec // 2  # Number of occupied orbitals
        self.numAO = 2  # Number of AOs
        self.occupancy = np.zeros(self.numAO)
        self.occupancy[:self.numOcc] = 2.0  # Here, only the first orbital is occupied by two electrons

    def integrals(self):
        self.one_electron_integrals()
        self.two_electron_integrals()

    def one_electron_integrals(self):
        self.H = np.zeros((self.numAO, self.numAO))  # Core Hamiltonian
        self.S = np.zeros((self.numAO, self.numAO))  # Overlap matrix
        self.X = np.zeros((self.numAO, self.numAO))  # X-matrix for canonical orthogonalization

        # OEI values of HeH+ (R = 1.4632 a.u.) with STO-3G
        S12 = 0.450770e+0
        T11 = 2.164313e+0
        T12 = 0.167013e+0
        T22 = 0.760033e+0
        V11A = -4.139827e+0
        V12A = -1.102912e+0
        V22A = -1.265246e+0
        V11B = -0.677230e+0
        V12B = -0.411305e+0
        V22B = -1.226615e+0

        # Form core Hamiltonian
        self.H[0][0] = T11 + V11A + V11B
        self.H[0][1] = T12 + V12A + V12B
        self.H[1][0] = self.H[0][1]
        self.H[1][1] = T22 + V22A + V22B

        # Form overlap matrix
        self.S[0][0] = 1.0
        self.S[0][1] = S12
        self.S[1][0] = S12
        self.S[1][1] = 1.0

        # Form X-matrix for canonical orthogonalization
        self.X[0][0] = 1.0 / np.sqrt(2.0 * (1.0 + S12))
        self.X[1][0] = self.X[0][0]
        self.X[0][1] = 1.0 / np.sqrt(2.0 * (1.0 - S12))
        self.X[1][1] = -self.X[0][1]

        # Print out integral values
        if self.doDebugPrint:
            print('The S array \n{}\n'.format(self.S))
            print('The X array \n{}\n'.format(self.X))
            print('The H array \n{}\n'.format(self.H))

    def two_electron_integrals(self):
        self.TT = np.zeros((self.numAO, self.numAO, self.numAO, self.numAO))

        # TEI values of HeH+ (R = 1.4632 a.u.) with STO-3G
        V1111 = 1.307152e+0
        V2111 = 0.437279e+0
        V2121 = 0.177267e+0
        V2211 = 0.605703e+0
        V2221 = 0.311795e+0
        V2222 = 0.774608e+0

        # Form matrix of two-electron integrals
        self.TT[0][0][0][0] = V1111
        self.TT[1][0][0][0] = V2111
        self.TT[0][1][0][0] = V2111
        self.TT[0][0][1][0] = V2111
        self.TT[0][0][0][1] = V2111
        self.TT[1][0][1][0] = V2121
        self.TT[0][1][1][0] = V2121
        self.TT[1][0][0][1] = V2121
        self.TT[0][1][0][1] = V2121
        self.TT[1][1][0][0] = V2211
        self.TT[0][0][1][1] = V2211
        self.TT[1][1][1][0] = V2221
        self.TT[1][1][0][1] = V2221
        self.TT[1][0][1][1] = V2221
        self.TT[0][1][1][1] = V2221
        self.TT[1][1][1][1] = V2222

        # Print out integral values
        if self.doDebugPrint:
            print('The TT array')
            for i in range(self.numAO):
                for j in range(self.numAO):
                    for k in range(self.numAO):
                        for l in range(self.numAO):
                            print('({}{}|{}{}) = {}'.format(i+1, j+1, k+1, l+1, self.TT[i][j][k][l]))
            print('')

    def RHF(self):
        # Settings
        self.maxIter = 25  # Maximum number of iterations
        self.tolScfConv = 1.0e-4  # Convergence criterion for density matrix

        # Clear arrays
        P = np.zeros((self.numAO, self.numAO))  # Use core-Hamiltonian for initial guess (P = 0)
        G = np.zeros((self.numAO, self.numAO))  # Two-electron part of Fock matrix
        F = np.zeros((self.numAO, self.numAO))  # Fock matrix
        Fprime = np.zeros((self.numAO, self.numAO))  # Transformed Fock matrix
        E = np.zeros(self.numAO)  # Orbital energy matrix
        C = np.zeros((self.numAO, self.numAO))  # MO coefficient marix
        Cprime = np.zeros((self.numAO, self.numAO))  # MO coefficient marix in the orthonormal basis
        oldP = np.zeros((self.numAO, self.numAO))  # Density matrix in the previous step

        # SCF loop
        for iter in range(self.maxIter):

            # Form G
            for i in range(self.numAO):
                for j in range(self.numAO):
                    G[i][j] = 0.0
                    for k in range(self.numAO):
                        for l in range(self.numAO):
                            G[i][j] += P[k][l] * (self.TT[i][j][k][l] - 0.5 * self.TT[i][l][k][j])

            # Form Fock matrix
            F = self.H + G

            # Calculate energy
            self.engHF = (P * (self.H + F)).sum() * 0.5

            # Transform Fock matrix to the orthonormal basis
            Fprime = np.dot(self.X.T, np.dot(F, self.X))

            # Diagonalize transormed Fock matrix
            E, Cprime = np.linalg.eig(Fprime)

            # Sort eigenvalues & eigenvectors
            eigen_id = np.argsort(E)
            E = E[eigen_id]
            Cprime = Cprime[:,eigen_id]

            # Back-transform eigenvectors to the AO basis
            C = np.dot(self.X, Cprime)

            # Form new density matrix
            oldP = np.copy(P)  # Save the present density matrix before creating new one
            mask = self.occupancy>0
            P = np.dot(C[:,mask] * self.occupancy[mask], C[:,mask].T)

            # Calculate delta
            delta = np.sqrt(((P - oldP) ** 2).sum()) * 0.5

            # Print out information
            print('Step #{}: energy = {}, delta = {}'.format(iter, self.engHF, delta))
            if self.doDebugPrint:
                print('The G array \n{}\n'.format(G))
                print('The F array \n{}\n'.format(F))
                print('The F\' array \n{}\n'.format(Fprime))
                print('The C\' array \n{}\n'.format(Cprime))
                print('The E array \n{}\n'.format(E))
                print('The C array \n{}\n'.format(C))
                print('The P array \n{}\n'.format(P))

            # Check for convergence
            if delta < self.tolScfConv:
                print('\nSCF calculation converged!')
                break
            elif iter + 1 >= self.maxIter:
                print('\nNo convergence in SCF...')

        # Calculate nuclear repulsion energy
        self.engNuc = self.za * self.zb / self.r
        self.engTot = self.engHF + self.engNuc
        print('-- Electronic energy = {} hartrees'.format(self.engHF))
        print('-- Total energy      = {} hartrees'.format(self.engTot))

################################################################################

# <Input parameters>
# r           : Bond length (in a.u.)
# zeta1       : Slater orbital exponent (function 1)
# zeta2       : Slater orbital exponent (function 2)
# za          : Atomic number of atom A
# zb          : Atomic number of atom B
# doDebugPrint: Flag to print out arrays

calc1 = szabo(r=1.4632, zeta1=2.0925, zeta2=1.2400, za=2.0, zb=1.0, doDebugPrint=False)
calc1.integrals()
calc1.RHF()

---

1. A. ザボ, N. S. オストランド著, "新しい量子化学 上", 東京大学出版会 (1987). ↩

2. 教授世代は「藤永先生の教科書」で勉強された方が多いようです。現在絶版になっており、たまに出版社が復刊と称して重刷してくださることがあります。藤永 茂, "分子軌道法", 岩波書店 (1980). ↩

3. 化学業界で最も使われている商用パッケージ。http://gaussian.com/ ↩

4. 業界ではGaussianに次いで著名なフリーの計算パッケージ。https://www.msg.chem.iastate.edu/GAMESS/ ↩

5. ザボの教科書では、HeH$^+$分子のみを計算することを想定しており、コード中では電子数 (占有軌道数) が固定値として入力されている。 ↩

6. 占有数はMOに関する情報 (電子が占有するのはAOではなくMO) なので、本当はMOの次元数 (≦ AOの次元数) だけの配列を用意すれば良い。しかしながら、MOの次元数 ≒ AOの次元数であることから、簡単のために、この手のベクトル・行列はAOの次元数だけメモリを確保することが多い。また、このように指定しておくことで、このメモリ領域を後でAO系の情報の保存にも使い回しできるようになる。 ↩

7. 1電子ハミルトニアン $\textbf{H}^\text{core}$ は、core Hamiltonianあるいはbare-nucleus Hailtonianとも呼ばれる。 ↩

8. ザボの教科書では $\textbf{X}$ と表記されているが、おそらく一般的には $\textbf{Q}$ と書かれることが多く、そのためこの行列は "Q-matrix" と呼ばれることもある。 ↩

9. 「観測可能な物理量 (observable) に対して、線形なエルミート演算子が各々存在する」という量子力学的な要請に由来します。 ↩

10. 正準直交化のための行列 $\textbf{X} \in \mathbb{R}^{N_\text{AO} \times N_\text{MO}}$ はAO基底を正規直交基底 (orthonormal basis) に変換するものです。MO係数行列 $\textbf{C}$ であれば、AO基底からMO基底への変換を担います。これらの変換行列は物理量に対応するエルミート行列とは異なり、一般に非対称行列であり、非正方行列になることもあります。たとえば、MO次元数 $N_\text{MO}$ ≦ AO次元数 $N_\text{AO}$ ですから、$\textbf{C}$ は非正方行列になり得ます。長方形状である行列 $\textbf{C}$ の行方向にAO情報、列方向にMO情報が対応します。また、業界用語で「$\textbf{C}$ はAO足からMO足への変換行列である」といった表現をすることがあります。独特ですが、基底のことを「足」と呼ぶ習慣があります。これは、行列要素の下付き添字で何基底なのか区別できることに由来するのだと思います。たとえば $\textbf{C}$ の場合、その要素 $C_{\mu i}$ の $\mu$ はAO足、$i$ はMO足を表す添字です。 ↩

11. このアルゴリズムはDirect SCFと呼ばれる。 ↩

12. Born-Oppenheimer近似下での非相対論的なハミルトニアンを取り扱う範囲では、2電子積分といえば電子反発積分 (electron-repulsion integral; ERI) に対応します。逆にいえば、Born-Oppenheimerを超えた取り扱い、相対論的な取り扱いをする場合には、電子反発積分以外の2電子積分を扱うことがあり、これらは電子反発積分と異なる数学的性質を持つ場合があることに注意してください。 ↩

13. 密度行列 (density matrix) は、記号 $\textbf{D}$ で表されることも多いです。しかしながら、ここではザボの教科書の表記に倣って $\textbf{P}$ としておきます。 ↩

14. 対角化の結果得られる固有値・固有ベクトルの格納順序は、どのサブルーチンを利用するのかに依存します。LAPACKのdsygvルーチンなどでは、固有値の昇順にソートされた結果が配列に格納されます。 ↩