## SCF 法
SCF (Self-consistent Field) は、シュレーディンガー方程式を計算機で解く場合に用いられている数値解法の一種です。この方法では、方程式に解の初期値を与え、反復的に解くことによって真の解に限りなく近い近似解を探索します。
セルフコンシステントセルフコンシステント は、求めるべき解が自分自身を含むような問題、あるいはそのような問題に帰着させる解析手法のことである。得られる解が与えられる解の候補と一致しなければならないためセルフコンシステントとか、自己無撞着とか自己整合などと呼ばれる。撞着とは整合性がなく矛盾することを指し、無撞着であることは矛盾がなく整合することを意味する。多体問題を1体問題に近似させる場合、たとえば量子力学におけるハートリー＝フォック方程式や、統計力学における平均場近似（分子場近似）などはセルフコンシステントな方程式を解く問題としてしばしば取り上げられ、変分問題において重要な概念である。セルフコンシステント方程式の解が厳密に求まる場合はそれほど多くないが、セルフコンシステント方程式に解候補を与え、新たな解候補を作ることによって得られる解候補の列が収束するような場合、適切な解候補を初期値として選ぶことでより良い近似解を得られる。このような方法をセルフコンシステント法と呼ぶ。
ここでは、以下の行列方程式を SCF 法で解くコードを NumPy で実装したいと思います。
は 正方行列、 は 対角行列、 は以下の式で与えられる 正方行列です：
方程式は一見ふつうの固有値問題のように見えて、 が解である に依存して決まるためそのまま解けないという構造になっています。
ヒント
の共役転置 を求めてみる。 に対し、 であることから、
となることを利用すれば、
であることがわかる。すなわち、 がエルミートであれば もまたエルミートであり、ユニタリー行列 が存在して、次のように対角化が可能となる：
よって、初期解 に任意のエルミート行列を選べば、SCF 解 は を対角化するようなユニタリー行列に収束することがわかる。
## 方針
まず初期解 を与え、初期行列 を定めて固有値問題を解きます。得られた解 を使用して再び を計算し、同様に固有値問題を解きます。この手続きを、真の解になるべく近い値が得られるまで繰り返します。
実際には、 回目の試行で得られた解を とするとき、
となったときに計算をカットオフします（変分原理）。ただし、 は収束の閾値です。得られた解が妥当かどうかは、初期解 に依存します。
## 実装
import numpy as np
### 所与
まず、初期解 と収束の閾値 を与えます。ここでは、
および とします。
C_0 = np.ones([2, 2]) # 初期解
th = 1e-6 # 収束の閾値
### 計算関数
の計算と、その固有値・固有ベクトルを求める関数を作成します。固有値・固有ベクトルは np.linalg.eig() を使用して求めることができます。
def M(C):
return 2 * np.identity(2) + (C.T @ C) ** 3 @ C
def get_eig(M):
(val, vec) = np.linalg.eig(M)
return (val, vec)
### SCF ループ
SCF ループを行う関数を作成しましょう。初回の計算を行う部分は以下のように実装できます。
def scf_loop(M_0, C_0, th):
M_0 = M(C_0)
(val, vec) = get_eig(M_0)
tmp = (val, vec)
count = 0
tmp は、ループごとの結果を格納する一時変数です。以降のループは、while True 文を使って実装してみます。
def scf_loop(M_0, C_0, th):
# --- 初回の計算 ---
while True:
count += 1
M_new = M(tmp[1])
(val_new, vec_new) = get_eig(M_new)
diff = np.linalg.norm(tmp[0] - val_new)
print(f"ε = {diff}")
if diff < th:
print(f"SCF converged after {count} steps.")
print(f"C = {tmp[1]}")
break
else:
tmp = (val_new, vec_new)
diff = np.linalg.norm(tmp[0] - val_new) はステップ間の解の差分ノルム です。これが より小さければ結果を表示してループを終了し、そうでなければ tmp を新しい計算結果で更新します。
### 整理
最後に、処理を main 関数にまとめます。
def main():
M_0 = M(C_0)
scf_loop(M_0, C_0, th)
# --- 関数定義 ---
if __name__ == "__main__":
main()
## 完成したコード
import numpy as np
C_0 = np.ones([2, 2]) # 初期解
th = 1e-6 # 収束の閾値
def main():
M_0 = M(C_0)
scf_loop(M_0, C_0, th)
def M(C):
return 2 * np.identity(2) + (C.T @ C) ** 3 @ C
def get_eig(M):
(val, vec) = np.linalg.eig(M)
return (val, vec)
def scf_loop(M_0, C_0, th):
M_0 = M(C_0)
(val, vec) = get_eig(M_0)
tmp = (val, vec)
count = 0
while True:
count += 1
M_new = M(tmp[1])
(val_new, vec_new) = get_eig(M_new)
diff = np.linalg.norm(tmp[0] - val_new)
print(f"ε = {diff}")
if diff < th:
print(f"SCF converged after {count} steps.")
print(f"C = {tmp[1]}")
break
else:
tmp = (val_new, vec_new)
if __name__ == "__main__":
main()
## 実行
$ python scf.py
ε = 31.31685115084307
ε = 1.1260325006104934
ε = 1.0418242466972605
ε = 0.1790963391034448
ε = 0.029636692487829575
ε = 0.05159161960237516
ε = 0.041294139499986525
ε = 0.011897494606076985
ε = 0.00559038348319239
ε = 0.006815210309690685
ε = 0.0029511899187867033
ε = 0.0003422485813089852
ε = 0.0009304107303725144
ε = 0.0005892126941714231
ε = 0.00012700430815054
ε = 0.00010049462978659264
ε = 9.959467611813162e-05
ε = 3.7853619339729646e-05
ε = 6.70960240338517e-06
ε = 1.4359995706696292e-05
ε = 8.117894426856232e-06
ε = 1.3579272825675583e-06
ε = 1.68942799605633e-06
ε = 1.4404718789145576e-06
ε = 4.7650304526233737e-07
SCF converged after 25 steps.
C = [[ 0.90881315+0.j -0.24725334+0.33604234j]
[ 0.24725334+0.33604234j 0.90881315+0.j ]]
今回の場合、25回のループの後に以下の解が得られました：
np.linalg.inv() などを使用して逆行列 を求めてみましょう。（おおよそ） となっていることがわかると思います。したがって、 はユニタリー行列です。
つぎに、 を変えて実行してみましょう。今回得られた近似解に近い以下の行列を与えてみます：
Warn
コード上では虚数単位が j となることに注意が必要です。
C_0 = np.array([
[0.91, -0.25 + 0.34j],
[0.25 + 0.34j, 0.91]
])
$ python scf.py
ε = 0.008489633208246214
ε = 0.02819032609290374
ε = 0.015494468731080513
ε = 0.002390801048507048
ε = 0.0032483757748762044
ε = 0.0026966334167551264
ε = 0.000871229442128549
ε = 0.00026442103997421484
ε = 0.000410450971958739
ε = 0.00020414537123217554
ε = 2.459519218885558e-05
ε = 5.2239477591821945e-05
ε = 3.8366712776089546e-05
ε = 1.0582162654511265e-05
ε = 4.95267462076122e-06
ε = 6.129425438534219e-06
ε = 2.713318516361577e-06
ε = 3.011124498637939e-07
SCF converged after 18 steps.
C = [[ 0.9088131 +0.j -0.24725362+0.33604225j]
[ 0.24725362+0.33604225j 0.9088131 +0.j ]]
18 ステップで同じ解が得られました。このように、SCF の収束の速さは初期解に依存します。
## まとめ
SCF 法は著名な量子化学計算ソフトウェアでも用いられる重要な数値解法の一つですが、意味するところは比較的単純です。今回は、NumPy を使用して行列方程式を SCF 法で解くプログラムを実装しました。