1E18
非局所励起状態に対する分割統治(DC)型線形スケーリング法の開発（早大先進理工1，早大理工研2，JST-CREST3）
中井 浩巳1,2,3，吉川 武司1，野中 佑太郎1

Development of divide-and-conquer based linear-scaling method
for non-local excited states(Advanced Science and Engineering, Waseda Univ.1, RISE, Waseda Univ.2, CREST, JST Agency3)
Hiromi Nakai1,2,3, Takeshi Yoshikawa1, Yutaro Nonaka1【緒言】
通常の量子化学計算は、系が大きくなるにしたがって計算コストが飛躍的に増大し、計算が困難となる。この問題を克服するために、分割統治(DC)法(1991 年、Yang)、エロンゲーション(ELG)
法(1991 年、今村、青木)、分子仕立(MT)法(1994 年、Gadre)、フラグメント分子軌道(FMO)法(1999年、北浦)などの分割型計算理論が開発されてきた。我々のグループでも、DC 法とエネルギー密度解析(EDA)を組み合せることにより、MP2,CCSD,CCSD(T)などの電子相関計算の線形スケーリングを達成してきた。一方、励起状態計算に対する線形スケーリング法は発展途上であり、分割型計算理論の拡張が報告されているが、(A)部分系内の局所励起しか取り扱えない、(B)結合の切断に伴う非物理的な励起状態が現れるなどの欠点がある。そこで本研究では、非局所励起を取り扱うための線形スケーリング法の開発を検討する。

【理論】[1] DC-TDHF
ここでは、最も簡単な励起状態計算である時間依存ハートリー・フォック(TDHF)について示すが、時間依存密度汎関数理論(TDDFT)への拡張は容易である。TDHF 計算は、線形応答より導かれた固有値問題を解く場合が一般的である。DC 法では、部分系 s ごとに次式を解く。

 As
 s*
B0  X s Bs   X s s 1     sAs *   Y s 
0  1  Y (1)[2] DC-TDCPHF
周波数依存分極率は、外部振動場に対する応答であり、次の結合摂動(CP)方程式、TDCPHF を自己無撞着場(SCF)的に解くことにより求められる。そして、励起エネルギーは周波数分極率が発散する極として求められる。DC 法を適用すると、次の DC-TDCPHF 方程式が得られる。
(2)
F  , s C 0 , s  F 0 , s C  , s  S s C  , s  S s C  , s E 0 , s  S s C 0 , s E  , s
ここで、 F , C , S , E はそれぞれフォック行列、MO 係数行列、重なり行列、軌道エネルギーであり、添え字(0),(±)はそれぞれ無摂動および 1 次摂動を表す。全摂動密度行列 P 

は、部分系の密度行列に分割行列 p を作用させて求める。
sP     p s C  , s †n s C 0 , s  C 0, s † n sC  , s



(3)s

周波数依存分極率は、MO 表示の遷移双極子モーメント D を用いて次式から求められる    Tr P   D (4)

[3] DC-RPA
TDCPHF 方程式と等価な結果を与える表式として、乱雑位相近似(RPA)型の表式が知られている。
こ れ に DC を 適 用 す る と 、 部 分 系 に 対 す る 結 合 行 列 U   , s ( た だ し 、 C  , s  C 0 , sU  , s ,U  , s  U  , s †  0 )に対する 1 次方程式となる。

 I  A  B A  B  U    2A  B D2ssss

 ,sssss(5)ここで、 D は部分系の MO を用いた遷移双極子モーメントであり、 A , B は(1)式と等しい。
sss[4] DC-GF
外部振動場ごとに RPA 方程式を解く代わりに、グリーン関数を用いて書き換えられる。グリーン関数は無摂動状態((1)式と等価)の固有値・固有ベクトルを用いて表すことができる。DC 法を適用するとグリーン関数 G  , s および結合行列 U   , s は次式のように計算される。

G  , s   2 I s  As  B s As  B s  L  I  λ  L
1s2s 1ss†U  , s  2G  , s As  B s D s

ここで、 λ , L は A  Bssss

A  B  の固有値・固有ベクトルである。
s(6)
(7)

sCPU time [sec]【結果と考察】Table 1. Excitation energies (eV) of (H2N)C20H40(COOH) calculated by
表 1 は、置換ポリエチレ
TDHF, DC-TDHF, TDCPHF, and DC-TDCPHF methods.
ン(H2N)C20H40(COOH)(図 1)
TDHF
DC-TDHF
Main Configuration
TDCPHF DC-TDCPHF
Stateに対して、TDHF, DC-TDHF,1n(O)
*(C=O)
6.34(0.001)
6.31(0.004)
6.35
6.35
1 A"TDCPHF, DC-TDCPHF 計算
1n (N) *(N-H)
8.93 ( 0.009 )9.09 ( 0.009 )
8.84
8.84
2 A'
1により求めた励起エネル(C=O)
*(C=O)
9.28(0.209)
9.11(1.29)
9.25
9.25


3 A'
110.0210.02
ギーを示している。DC 法
2 A"  (C=O) * (C=O) 10.04 ( 0.005 ) 10.05 ( 0.074 )は図 1 のように 4 個の部分Fragment 2 Fragment 3Fragment 4Fragment 1系に分割し、バッファ領域として隣接 1 yユニットを考慮した。基底関数は
x6-31G**である。いずれの方法でも、
TDHF 計算の結果を比較的良く再現し
HOMO-3 n(O)ている。しかし、DC-TDHF 計算ではLUMO *(C=O)
1Figure1.
Fragmentsof(HN)C
H(COOH)and MOs related
220403 A’状態の振動子強度をかなり過大評with the first excited state.
価している。非物理的な状態は、
5000
DC-TDHF 計算では多数現れたが、DC-TDCPHF 計算
TDHF
では現れなかった。
4000
DC‐TDHF
次に、フッ化水素クラスター(HF)n に対する各手法
TDCPHF
3000
DC‐TDCPHFの計算時間を図 2 に示す。DC 法は HF モノマーを部GF分系とし、バッファ領域として隣接 1 ユニットを考慮
2000
DC‐GFした。基底関数は 6-31G である。従来法および DC 法
1000に対する TDCPHF, GF 計算では、振動電場 50 個に対する計算時間である。いずれの方法も従来法に比べ0010203040
DC 法により計算時間が短くなっている。特に、DC-GF
n法は非物理状態なども問題も解決され、同時にかなりFigure 2. CPU times of TDHF,
DC-TDHF, TDCPHF, DC-TDCPHF, GFの高速化も達成されていることがわかる。
and DC-GF of (HF)n clusters.