2P096
大規模分子動力学計算を目指した DC-DFTB 理論の開発
1235（分子研 ，早大理工研 ，早大先進理工 ，JST-CREST4 ，名大院理 ）
1,22,3,4
西澤宏晃 ，中井浩巳，Stephan Irle5

Development of DC-DFTB method aimed at large scale molecular dynamics（IMS1 , Waseda Univ. RISE2 , Waseda Univ.3 , JST-CREST4 , Nagoya Univ.5 ）Hiroaki Nishizawa1,2 , Hiromi Nakai2,3,4 , Stephan Irle5【緒言】
密度汎関数強束縛 (DFTB) 法 [1, 2] は，密度汎10000.0
Total CPU time
関数理論 (DFT) に基づく半経験的手法である。2Construct H0 and S
8000.0原子間の距離とそれらの種類によって変動するパConstruct γDiagonalize
ラメータ [1] を用いており，AM1，PM3 などほか ].ec
Construct densitys[ 6000.0の半経験的手法に比べて高精度な結果を得ること emができる。また，積分を行うことなくすべての要 tiU 4000.0P
素を計算することが可能なため，その計算コスト C2000.0は比較的低い。本手法によりカーボンナノチュー
ブなどのナノスケールの系や溶液系などを取り扱0.0
うことで，それらの反応系を精密に解析し，その制
02000
4000
6000
8000
御を行うことが可能になると考えられる。しかしNumber of atoms
Fig. 1 に示すように，その計算コストは前因子はFig. 1 Computatinal time for energy
小さいものの，系の大きさ N に対して O(N 3 ) で
calculation of C2n H2n+2 obtained byあり，大規模系に対する分子動力学計算は困難でDFTB.
ある。エネルギー計算を行うための手続きごとに計算時間を調査すると，そのほとんどがハミルトニアンの対角化によるものであることが分かる。そこで本研究では，対角化に対する計算時間を削減するために，中井，小林らによって開発が進められてきた，大規模計算手法である分割統治 (DC) 法 [3] を DFTB 法に対して適用した DC-DFTB 法の開発を行った。分子動力学計算を行うためにエネルギー勾配計算も行ったので，これらの結果についても報告する。また，次世代スーパーコンピュータ「京」での実行を考慮し並列化も行っているた
め，併せて報告する。
【理論】
DFTB 法のエネルギーは次のように表される。

E0 =∑Dµν (hµν + Fµν ) + Erep =µν∑

1 ∑γAB ∆qA ∆qB + Erep2
atom0Dµν Hµν+µν(1)A,B0ここで，Hµν，γAB ，Erep は上述したパラメータにより求められる。Dµν は電子密度，∆qA は原子A の Mulliken 密度変化量である。DC-DFTB 法では全系を重なりのない部分系 α に分割し，周囲の効果を取り込むためにバッファ領域を加えて局在化領域とする。局在化領域に対して DFTB 方程式
α α αFα Cαi = ϵi S Ci(2)αααを解くことで，部分系 α の軌道係数 Cαi を得る。F ，ϵi ，S はそれぞれ部分系 α の Fock 行列，軌道エネルギー，重なり行列である。電子数保存の制約条件から共通の Fermi 準位を決定し，密度行列を次のように構築する。
DCDµν=∑ααpαµν Dµν =∑αpαµν∑

α αfβ (ϵF − ϵi )CµiCνi(3)i

ここで pαµν は分割行列，ϵF は Fermi 準位，fβ (x) は Fermi 関数である。∆qA が変化しなくなるまで式 (2) の方程式を繰り返し解き，エネルギーを求める。その後，得られた密度行列を用いてエネルギー勾配を計算する。エネルギー勾配計算には小林らの手法 [4] を用いた。

【結果と考察】
DC-DFTB 法をポリエチレン (C2n H2n+2 ) 分子に適用した。
初めに，Fig. 2 に DC-DFTB 法を用いて計算を行ったときの計算時間を示す。C2 H2 を 1 ユニットとし，部分系を 1 ユニット，バッファ領域を左右 8 ユニットとした。左図には比較のために通常のDFTB 法による結果も灰色で示した。図からも明らかなように，DC-DFTB 法により計算コストを大幅に削減することに成功している。また，右図に拡大図を示しているが，DFTB 法でボトルネックとなっていた対角化の計算コストが，O(N ) となっていることが確認できた。ほかの手続きによる計算時間も前因子が小さく，全体として O(N ) が達成できていることが分かる。
次に計算精度について述べる。Fig. 3 に DC-DFTB 法においてバッファサイズを変えたときの，通常法とのエネルギーの誤差を示した。部分系は 1 ユニットで変化しないため，その誤差は原子数に比例して増加していることが分かる。バッファ領域を左右 3 ユニットとしたとき，C2500 H2502 では30 mEh 以上の誤差があるが，バッファ領域を大きくするにつれてその誤差は減少する。左右 7 ユニットをバッファ領域とした場合，1 mEh 以下の誤差で全エネルギーを求めることができた。
最後にエネルギー勾配法について示す。Fig. 4 に部分系を 1 ユニット，バッファ領域を左右 8 ユニットとした DC-DFTB 法と通常の DFTB 法によるエネルギー勾配の計算時間を示した。式 (1)の右辺で表される第 1，2，3 項の計算時間をそれぞれ Normal, Gamma, Repulsive contributionとして示している。図からもわかるように，エネルギー勾配計算においてボトルネックとなるのはNormal contribution である。DC-DFTB 法ではこの項を部分系ごとに求めるため，計算時間が短縮され，さらに O(N ) が達成できていることが分かる。また，このときの平均絶対誤差は 10−5 であり，高い精度で計算を行うことが可能であることが確認できた。
これらの結果は，本手法を用いた大規模系に対する分子動力学計算への適用の可能性を示唆している。当日はカーボンナノチューブ，溶液系への適用に関しても報告を行う。
70.0

10000.0
Total CPU time
Construct H0 and S
Construct γ
Diagonalize
Determine the Fermi level
Construct density

8000.0.]ces 6000.0[eimtUPC 4000.0.]c 50.0es[ 40.0emit 30.0U
PCO(N1.00)

20.0

2000.00.0O(N1.09)

Total CPU time
Construct H0 and S
Construct γ
Diagonalize
Determine the Fermi level
Construct density

60.0O(N1.52)
O(N1.15)O(N1.54)

10.0
05000

10000

15000

0.020000

Number of atoms

05000

10000

15000

Number of atoms

20000

Fig. 2 Computatinal time for energy calculation of C2n H2n+2 obtained by DC-DFTB andconventional DFTB.
0.0

800.0]hEm[r -10.0or
reyg
re -20.0nE.]c 600.0
es[emti 400.0
UPC

500 atoms
1,000 atoms
2,000 atoms
3,000 atoms
4,000 atoms
5,000 atoms-30.0
2345678

Number of units for buffer region
Fig. 3 Buﬀer-size dependency of the energy
error in C2n H2n+2 obtained by DC-DFTB[in mEh ].

Normal contribution (DFTB)
Gamma contribution
Repulsive contribution
Normal contribution (DC-DFTB)
O(N1.00)

200.0

90.005000

10000

15000

Number of atoms

20000

Fig. 4 Computatinal time for energy gradient calculation of C2n H2n+2 obtained by
DC-DFTB and conventional DFTB.

[1] D. Porezag, Th. Frauenheim, Th. Köhler, Phys. Rev. B, 51 (1995) 12947.
[2] M. Elstner, D. Porezag, G. Jungnickel, J. Elsner, M. Haugk, Th. Frauenheim, S. Suhai, G. Seifert, Phys. Rev.
B, 58 (1998) 7260.
[3] M. Kobayashi, H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics: Methods andApplications (2011, Springer), pp. 97-127.
[4] M. Kobayashi, T. Kunisada, T. Akama, D. Sakura, H. Nakai, J. Chem. Phys., 134 (2011) 034105.