量子化学計算コードを自作してみた記録 その１

緒言

量子化学計算に関わっていると、（計算手法などを新たに開発される方は別として）とりあえず普段行う計算は、GaussianやGAMESSといったパッケージのコードを使って行われることが多い。これらのコードは十分にテストされていて、しかも速いので化学や材料の研究にはこれらを用いるのが懸命と思われるが、これらのソフトがよくできすぎているために、量子化学計算はブラックボックスになりがちである。
　その理由の一つは、二電子積分をはじめとする電子積分のプログラム化が難しく、この部分がボトルネックとなってしまってプログラム作成に挫折してしまうためだと思われる。
　しかし、非常に遅い計算コードでも、一度は書いた方が良いかと思い、最も基本的な手法であるHartree Fock法を自作して見た。一度自分で書いてみると、積分以外のところも案外わかっていないものだと気づくし、また、Gaussianなどのコードで多数用意されているオプションが急にスッとわかるようになってきたりする。

　ちなみに、私は「量子化学計算の計算手法を考える人」ではないので、ここに出てくる各論（効率のよい積分方法、対角化アルゴリズム、SCFの際の密度行列の混ぜ合わせ、など）に関して深いことは知らない。もし間違いがあったら指摘してもらえると大変うれしい。
　また、計算手法の導出そのものを解説するわけではないので、Hartree Fock法そのものや、「そもそもSCF(Self-Consistent Field; 自己無撞着場)とはなんぞや」については、参考文献などで勉強してほしい。あくまで、これらの式とコードを結びつけることが本記事の目標である。

方針

ここでは、制限Hartree Fock法(RHF)を実装する。
ところで、プログラム作成の一般的な順序としては、

1. 基底関数（短縮ガウス関数）を扱うクラスを準備
2. 1. で準備した基底関数を使用して、積分を評価する関数を作る
3. 2.で作成した積分評価関数を呼び出す形で、HF方程式を解く関数を実装する。具体的には、フォック行列F、重なり行列Sを生成し、SCFを解かせる。

というものになるかと思う。しかしながら、手順２で出てくる積分評価の関数を書くのが恐らくもっとも難しく時間がかかる。そこで、本記事では、下記の要領で進める。

1. 始めにHF方程式を解く関数を実装する。各種積分の値は、Szaboの本に載っている計算値（答え）を利用し、ハードコードしておく。
2. その後で気が向けば積分を評価する関数を書いていく。

コード

完成したコードは、ここにあげてある。
また、積分も全て実装したC++コードは、ここにある。ただし、遅いので実用にはならない。

前提＆準備

プログラミング言語は、Pythonを使う。
余談だが、Pythonで書かれた高機能な量子化学計算ツールとして、PyQuanteやPySCFといったものがある。私も以前はPyQuanteのソースコードリーディングから始めた。

使用ライブラリとしては、下記のものを使う。

• Numpy : 量子化学計算で固有値問題を解くのに使う。対角化なども、より適した技法(Davidson diagonaliationなど）があるが、現段階ではここはまず律速にはならないので考えない。

• Scipy : 当面はいらないのだが、最終的に二電子積分(ERI)や核引力積分(NAI)を評価する際、展開していくと、不完全ガンマ関数の評価の必要が出てくる。これを評価するのに用いる（その気があれば自分で書いても良いが、ここでは扱わない）。

それでは始めよう。

HF方程式の実装

Hartree Fock法は、つまるところローターン方程式（下記）に帰着する。

 FC = SC\epsilon 

　ここで、Fはフォック行列、CはFを対角化して得られる原子軌道の係数の行列、Sは基底関数の重なり積分の行列である。
　これらのうち、我々が計算機の上で行うことは、FとSの行列要素を全て求めてやり、このFを行列計算のライブラリなどを用いて対角化することである。しかしながら、「方針」のところでも述べたように、Fの中身を全て自前で計算するコードを書くのは大変な作業なので、ここではSzabo, Ostlundの本に書かれている「核間距離1.4オングストロームの水素分子をsto-3g基底関数系で解く時の値」を適宜利用することとする。
　まず、水素原子の基底関数は、以下の通り。

\phi_{H-1s} = 0.444635 g(0.109818)+ 0.535328 g(0.405771) + 0.154329 g(2.22766) \tag{3.221}

ここで、g(\alpha) は、$ \exp(-\alpha r^2) $に、規格化因子をかけたものである。

　積分計算を自作する代わりに用いる値（答え）は下記の通りで、これらは英語版だと160ページ付近に載っている。
　これらの積分は１回計算してしまえば、SCF全体に渡って使いまわせる。なお、式番号は、全て、Szabo, Ostlundの本の式番号である。

S = 
\begin{pmatrix}
 1.0 & 0.6593 \\
 0.6593 & 1.0 \\
\end{pmatrix} \tag{3.229} \\

H^{core} = 
\begin{pmatrix}
 -1.1204 & -0.9584 \\
 -0.9584 & -1.1204 \\
\end{pmatrix} \tag{3.233} \\ 

(\phi_1 \phi_2 | \phi_1 \phi_1 ) = (\phi_2 \phi_2 | \phi_2 \phi_2 ) = 0.7746 \\
(\phi_1 \phi_1 | \phi_2 \phi_2 ) = 0.5697 \\
(\phi_2 \phi_1 | \phi_1 \phi_1 ) = (\phi_2 \phi_2 | \phi_2 \phi_1 ) = 0.4441 \\
(\phi_2 \phi_1 | \phi_2 \phi_1 ) = 0.2970  \tag{3.235}

SCFループを書く

復習になるが、SCFループを解く流れは以下の通りだ。疑似コードで表す。

H, S = calc_H(), calc_S() # ループの最中、コアハミルトニアンH（運動エネルギー積分と、核引力積分の和）、重なり積分Sは使い回す。
X = symmetrize_matrix(S)  # X_adjcent() * S * X == E となるようなSを求める。
D = initial_guess()       # 初期密度行列。とりあえず0行列を返しておけばOK

loop:
   G = calc_G(D)    # 密度行列を用いて、行列G(クーロン反発、交換相互作用の和）を求める。
   F = H + G        # フォック行列の完成
   F_prime = X.adjcent * F * X
    
   eigenvalues, eigenvectors = F_prime.diagonalize() # eigenvalues が軌道エネルギー
   C = X*eigenvectors     # 上で求めたeigenvectorsは、変換されたF(つまり、F')に対して求めたものなので、逆変換してC得る。これが軌道の係数の行列となる。
   D_new = calc_D(C, n_occupied_orbital)

   E0 = calc_E0()
   if is_converged?(D_new, D)
       break
   else
       D = D_new
       continue  # ループの頭に戻る。

あとは、上の疑似コードを素直にコード化していくだけである。今回はH, Sは与えられる。また、calc_Gの中で本来はその都度評価する二電子積分もハードコードしてあるので、自作するのは、上記のうち、

• symmetrize_matrix()
• calc_G(D)
• calc_D()
• calc_E0()
• is_converged?()

だけだ。

重なり積分Sの対角化

FC = SCeを解くときに、Sが単位行列でないために、Fの固有値を求めても、Cが決まらない。そこで、このSの部分を消すために、

 X ^\dagger S X = E

となるXを求める必要がある。これには２種類のやり方があるが、今回はCanonical Orthogonalizationという方法を用いる。詳しい記述は、本を読んでもらいたいが、定義としては、

 X_{ij} = U_{ij}/s_j^\frac{1}{2} \tag{3.170}

であり、UはSの固有ベクトル、sは固有値である。

hf.py

def canonical_orthogonalization(S):
    # Szabo. pp.144 (3.169 - 3.172)
    l, U = np.linalg.eigh(S)
    row,col = S.shape
    X = np.zeros(S.shape)
    for i in range(row):
        for j in range(col):
            X[i,j] = U[i,j] / (math.sqrt(l[j]))
    return X

行列Gの計算

行列Gの定義は、

 G_{\mu v} = \sum_{\lambda,\sigma} P_{\lambda\sigma}[(\mu v | \sigma \lambda) - \frac{1}{2}(\mu \lambda|\sigma v)] \tag{3.154}
 

　式をみれば分かるように、ここで二電子積分が評価される。今回はハードコードしてあるが、二電子積分は単純に考えれば基底関数の数の４乗に比例していく。したがって、テーブルが膨大になってしまうため、最近はテーブルを作るのではなく、毎回計算し直す方法が主流のようである（ダイレクトSCF）。

今回ハードコードした値は、下記の関数で得られるが、このときの添字の交換について簡単に記述しておく。

hf.py

def eri_table(u,v,p,q):
    if u < v:
        u,v = v,u
    if p < q:
        p,q = q,p
    if (u+1)*(v+1) < (p+1)*(q+1):
        u,v,p,q = p,q,u,v
    
    # Szabo. pp162(p3.235)
    if (u,v,p,q) == (0,0,0,0) or (u,v,p,q) == (1,1,1,1):
        return 0.7746
    elif (u,v,p,q) == (1,1,0,0):
        return 0.5697
    elif (u,v,p,q) == (1,0,0,0) or (u,v,p,q) == (1,1,1,0):
        return 0.4441
    elif (u,v,p,q) == (1,0,1,0):
        return 0.2970
    else:
        # Never get here.
        print(u,v,p,q)
        raise

式の対称性から、(u,v,p,q)=(v,u,p,q)=(u,v,q,p)=(p,q,u,v)という関係がある。したがって、テーブルは必ずu < vかつp < qになるようにする。さらに、u*v > p*qになるようにしているのが、上記のコードで３個目のif文である。
（注: 添字をそれぞれ一つ増やしているのは、pythonでの配列のインデックスは0から始まるため、掛け算して0になると順序を決められないため）

密度行列Dの計算

密度行列Dの定義は、

 D_{pq} = 2 \sum_a^{N/2} C_{\mu a}C\ast_{va} \tag{3.145}

ここで、Nは電子の数である。したがって、RHFの場合は占有軌道の軌道係数だけが密度行列に入り込むことになる。

電子のエネルギーの計算

E_0 = 0.5 \sum_\mu \sum_v D_{\mu v}(H_{\mu v}^{core} + F_{\mu v}) \tag{3.184}

この時に使う密度行列は、上記で固有値計算をして得た密度行列を使う。エネルギーは、ボルン・オッペンハイマー近似下での、電子のみのエネルギーである。核同士の反発エネルギーは別途足す必要がある。
ちなみに、GaussianなどのソフトウェアでSCF収束時に出力される、エネルギー(SCF Doneでgrepしてみるとよい)は、核同士の反発エネルギーをこれに足し込んだものである。

収束の確認

　やり方はいろいろある。
　最も単純な実装は、エネルギーの変化が一定の値よりも小さくなったかどうか、という基準だろうか。
　多くのコードでは、密度行列の変化（差の二乗和の平均と、要素同士のもっとも大きな変化の差がともに閾値以下になった場合）を見ているかと思うので、ここではそれを使用した。
　収束していなければ、新しく出てきた密度行列を、保存しておき、次のループに入る。

最後に

今回は積分をハードコーディングしたが、気が向けば、これらも実装してようと思う。
また、今回作成したPythonのコードは実行すると、

**** iteration: 1 ****
E0:   0.0
RMSD: 0.363203624348
**** iteration: 2 ****
E0:   -1.8310386546
RMSD: 0.0

という感じの出力が得られるはずで、電子のエネルギーは、1.8310 (a.u.)程度に収束するはずである。

参考文献

Szabo, Attila, and Neil S. Ostlund. "Modern quantum chemistry: introduction to advanced electronic theory." NY: McGraw-Hill (1989).