3D04
高精度な大規模系計算のための divide-and-conquer (DC)法の改良(早大理工) ○小林正人，赤間知子，中井浩巳【緒言】Hartree-Fock (HF)や DFT 計算は、Fock 行列の生成とその対角化の 2 つのステップが計算時間のボトルネックであり、ほとんどの時間が費やされる。これまでにこれらの計算の基底関数の数 N に対するオーダーを改善する様々な取り組みが行われてきた。前者に対しては高速多重極展開法などが実用的に用いられている。一方後者に対しては、全 Fock 行列の対角化に代わる方法のひとつとして、全系を部分系に分割して計算を行い、その結果を合わせて全系のエネルギーを求める divide-and-conquer (DC)法[1,2]が DFT 計算などに対して用いられている。
最近の DFT 計算の主流は、LDA や GGA 汎関数よりも精度の高い結果を得ることができる hybrid汎関数を用いた計算となっている。この計算を行うためには、新たに HF 交換項の計算を行う必要がある。最近当研究室で DC 法を HF 及び hybrid DFT 計算に適用した事例を報告した[3,4]。その結果、DC 法により得られるエネルギーと通常の計算によるエネルギーの差は、pure DFT 計算のときよりも HF 計算や hybrid DFT 計算のときの方が大きくなることが示された。本研究ではこの原因を突き止め、精度を向上させる手法を検討するために、DC 法を用いたときの密度行列、一電子項、Coulomb 項、HF 交換項、DFT 交換相関項などについて精査したので報告する。また DC法においては、通常の計算よりも SCF 収束性が悪化するという問題も生じていたが、DIIS を DC法における部分系の計算に拡張することによって収束性を改善する DC-DIIS 法を新たに考案したので、併せて報告する。
【DC 法における各エネルギー項の寄与と密度行列の精度】DC 法では全系の密度行列を部分系αの密度行列の和として表すことにより、計算コストのオーダーを O(N)まで落としている。
occ

ρ µν = 2∑ Cµi Cν i = ∑ pαµν ρ µν = ∑ ραµν , ραµν = 2pαµν ∑η (ε F − ε iα )Cµ i Cν iiαα(1)iここでε F は Fermi 準位、η (x)は階段関数である。部分系α の密度行列は、α の領域にその周囲の領域を buffer として含めた localization region の Fock 行列を作成して対角化を行うことにより決定する。計算に際しては、電子数 ne を保存させるため、ραを(2)式で近似して(3)式の制約条件を課す。

ραµν = 2pαµν ∑ f β (ε F − ε iα )Cµi Cν i(2)ine =1ραµν Sαµν∑∑2 α µν(3)ここで fβ (x)は Fermi 関数である。全系のエネルギーは、上記のスキームにより得られる全系の密度行列を用いて、通常の HF 及び DFT 計算と同様の手続きにより求める。DC 法の計算プログラ
ムを独自に作成し、GAMESS に実装した。部分系は全系を構成する各原子とし、buffer は原子を中心とした半径 rbuf の球状領域とした。
Fig. 1 に rbuf を変化させて DC-HF 計算を行ったときのポリアセチレン C180H182 分子の全エネルギ
ー及び一電子、Coulomb、HF 交換の各エネルギー項の、従来の HF 法からの誤差を示した。DC法では buffer のサイズが大きくなるにしたがって、各エネルギー項が一様に HF 計算のものに近づき、この系では rbuf =22Å では mhartree レベルで収束している。この系では HF 交換エネルギー
は Coulomb エネルギー、一電子エネルギーに比べて絶対値でそれぞれ 20 分の 1 以下、40 分の 1以下であるにもかかわらず、各エネルギー項の誤差は HF 交換項が他に比べて大きい。このことから DC 法による打ち切りが HF 交換項に大きく影響し、これを含む計算ではエネルギーの誤差が大きくなることが分かった。Fig. 2 にはこの計算を行ったときの密度行列要素の平均絶対誤差

12.0た。DC 法では rbuf =12Å でも密度行列を MAEで 10-5 の精度で再現できている。またグラフから、密度行列の再現性と全エネルギーの再現性には密接な対応があることがわかる。
【DC 法の SCF 収束性】一般に SCF 収束性を改善させる有効な方法として、DIIS 法[5]が知られている。DIIS 法では SCF が収束したときに以下で定義される誤差行列 ee = FρS − SρFenergy difference [mhartree](MAE)と全エネルギーの誤差を対数で表し

10.0

1 electron term

8.0

Coulomb
HF exchange

6.0

Total

4.0
2.00.010-2.01520
r buf [Å]2530-4.0-6.0(5)が零行列になることを利用し、過去の Fock

Fig. 1. Buffer-size dependence of energy component
differences in DC-HF calculation from conventional HF.

は、F に対応するρを近似により求めている
ため、SCF 解が正しく(5)式を満たさない。
そこで、DIIS は途中の精度まで用い、その後は他の手法で収束させなくてはならなかった[6]。しかし SCF 収束性を向上させるためには最終段階の加速が最も重要である。
我々は DIIS 法を部分系の表式に拡張する

10-1
total energy difference [hartree]

小自乗的に予測する。しかし DC 法において

10-1

10-2

density matrix

10-3-4案した。DC 法による SCF 解では、各部分系

10-3-4101010-510-510-610-610-710-710-810-810152025r buf [Å]

ことにより、最後まで DIIS を用いて収束性を向上させることの可能な DC-DIIS 法を考

10-2

total energy

30mean absolute error of density matrix

行列の線形結合によって次の Fock 行列を最35

Fig. 2. Total energy difference and MAE of density matrix
elements in DC-HF calculation from conventional HF.

ごとに SCF が収束しており、部分系間の電子移動がないと仮定すると、収束点において以下の式が成り立つ。

eα ≡ Fα ρα Sα − Sα ρα Fα = 0 ∀α = subsystems
して含まれるよりも外側の領域の電子状態は外場とみなすことができる。この外場との相互作用と部分系間の電子移動が、各部分系の SCF よりも早く収束するならば、(6)式のeαを用いて各部分系ごとに DIIS の手続きを行うことによって、SCF 収束を加速することが可能である。Fig. 3 には taxol 分子の DC-HF計算(rbuf = 4.5Å)を行ったときの各 SCF サイ
クルにおける密度行列要素の変化の最大を示した。通常の DIIS に DC-DIIS を用いることでさらに SCF を加速させることができた。

maximum change of density matrix

DC 法における部分系の計算では、buffer と(6)1.0(A) No DIIS

0.1(B) DIIS(C) DIIS+DC-DIIS

10-2
10-3
10-4-510

convergence criteria

10-610-7

turning on DIIS

0510

turning on DC-DIIS1520
SCF iteration2530

Fig. 3. Maximum change of density matrix element in SCF
iteration of DC-HF calculation (A) without DIIS, (B) withconventional DIIS, and (C) with DIIS and DC-DIIS.

[1] W. Yang, Phys. Rev. Lett. 66, 1438 (1991).
[2] W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).
[3] 赤間, 小林, 中井, 日本化学会第 85 春季年会, 1PC-013, 神奈川 (2005 年 3 月).
[4] 赤間, 小林, 中井, 日本コンピュータ化学会 2005 春季年会, 2P05, 東京 (2005 年 5 月).
[5] P. Pulay, J. Comput. Chem. 3, 556 (1982).
[6] D. S. Shaw and A. St-Amant, J. Theor. Comput. Chem. 3, 419 (2004).