計算科学技術特論B

大規模量子化学計算(2)小林 正人(北海道大学大学院理学研究院)
K-masato@sci.hokudai.ac.jp
2016/07/14
1

講義概要
 7/7 量子化学計算の概要と構成要素、高速化
 量子化学計算の目的と種類
 量子化学計算の手順、構成要素と高速化

 7/14 大規模系に適用するための量子化学計算法
 フラグメント分割に基づく方法
 フラグメント分子軌道(FMO)法 分割統治(DC)法 ラプラス変換MP2法
 2電子積分の密度フィッティング法
 MP2計算への応用(2)

量子化学計算にかかる時間と精度
分子の大きさの3乗に比例して計算時間増大
方法

Hartree-Fock
(HF)法

密度汎関数
理論(DFT)

MP2法
(摂動法)

CCSD法

CCSD(T)法

計算時間

O(N3)

O(N3)

O(N5)

O(N6)

O(N7)

近似レベル

電子相関理論

平均場理論

1000倍性能
の計算機で

×10

計算精度

定性的

×10

計算時間は
 精度の低い理論でもO(N3)

×4.0

×3.2

×2.7
正確

『京』をただ使うだけでは大きな分子を扱えない

 精度が上がるにつれて莫大に(3)

大規模量子化学計算手法
 並列化とプログラムの工夫で頑張る
 RSDFT (実空間密度汎関数理論)
 ProteinDF (タンパク質密度汎関数プログラム)

 密度行列を近似計算
 エネルギー最小化法, 密度行列purification法

 数学・アルゴリズムにより高速化
 積分のRI計算, Laplace変換MP2法

 系を分割して計算を簡略化

O(N)へ

 フラグメント分子軌道(FMO)法 エロンゲーション法
 分割統治(DC)法(4)

計算科学における分割統治(DC)法 マージソート (フォン・ノイマン, 1945)
 計算科学における最初の分割統治法
 n個のデータをソートするコスト: O(n log n) さまざまなDCアルゴリズム
 二分法 (求根, 探索) クイックソート カラツバ乗算法
 高速フーリエ変換 (FFT)

分割部分は高並列化が可能©Nuno Nogueira(な場合が多い)(5)

フラグメント分子軌道(FMO)法
 分子を(単結合で)切断
 結合に使われる電子は混成軌道を使って一方に寄せる (HOP)

 エネルギー等のプロパティ
は多体展開で求める
EFMO1 = ∑ EI (FMO1)
I

EFMO2=∑ E + ∑ ( E − E − E ) (FMO2)
I

I

I >J

IJ

I

J

EI: モノマーのエネルギー
EIJ: ダイマーのエネルギー
EIJK: トリマーのエネルギー

・・・

EFMO3= EFMO2 + ∑ ( EIJK − EIJ − EJK − EIK + EI + EJ + EK ) (FMO3)
I >J >K

 多体展開を打ち切り

計算の大幅な高速化[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method (CRC Press, 2009).

(6)

フラグメント分子軌道(FMO)法[1] フラグメントX [X = I (モノマー), IJ (ダイマー), …]の計算
𝑋𝑋 = 𝐻𝐻 𝑋𝑋 + 𝑉𝑉 𝑋𝑋 + 𝑃𝑃 𝑋𝑋 (μ, ν ∈ X)
�𝜇𝜇𝜇𝜇
 ハミルトニアン: 𝐻𝐻
𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇
𝜇𝜇𝜇𝜇
𝑋𝑋
 𝐻𝐻𝜇𝜇𝜇𝜇: フラグメントX自身のハミルトニアン𝑋𝑋
 𝑉𝑉𝜇𝜇𝜇𝜇: Xの外側の電子・原子核からの静電ポテンシャル
𝑋𝑋 = ∑
𝑉𝑉𝜇𝜇𝜇𝜇𝐾𝐾≠𝑋𝑋 ∑𝐴𝐴∈𝐾𝐾 𝜇𝜇 −𝑍𝑍𝐴𝐴
𝐫𝐫−𝐑𝐑 𝐴𝐴𝐾𝐾𝜈𝜈 + ∑𝜆𝜆𝜎𝜎∈𝐾𝐾 𝐷𝐷𝜆𝜆𝜆𝜆
Γ𝜇𝜇𝜈𝜈,𝜆𝜆𝜆𝜆𝑋𝑋 𝑃𝑃𝜇𝜇𝜇𝜇: Xに属していない混成軌道をプロジェクトアウト

� 𝑋𝑋 𝐂𝐂𝑋𝑋 = 𝐒𝐒 𝑋𝑋 𝐂𝐂𝑋𝑋 𝛆𝛆𝑋𝑋
 HF (KS)方程式: 𝐇𝐇
occ(𝑋𝑋)𝑋𝑋
 密度行列: 𝐷𝐷𝜇𝜇𝜇𝜇= 2 ∑𝑖𝑖

𝑋𝑋 𝑋𝑋
𝐶𝐶𝜇𝜇𝑖𝑖
𝐶𝐶𝜈𝜈𝑖𝑖電子数(とスピン)はあらかじめ指定が必要

密度行列Dの自己無撞着な決定(SCC)が必要[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method (CRC Press, 2009).

(7)

FMO計算に用いられる近似
 FMO2の計算時間: O(N2) [多体展開以外の近似なしの場合]
 カットオフ距離を使用して、いくつかの近似を導入
 RESPPC: フラグメント間静電ポテンシャルを2電子積分を用いずにMulliken電荷で近似
 RESDIM: ダイマー計算をあらわに実行せず、静電相互作用で近似
 RCORSD [post HF電子相関計算を実行する場合]: ダイマーの電子相関を計算せずに無視

 SCC計算はモノマーに対してのみ実行
計算時間を大幅に削減し、O(N)を実現[1] D.G. Fedorov and K. Kitaura eds., The Fragment Molecular Orbital Method (CRC Press, 2009).

(8)

FMO2計算手順
原子座標, フラグメントの情報
X = 1, N (モノマー)

モノマーSCCDK (K = 1-N)  HX

計算時間: O(N)RESPPCにより計算時間のプレファクターを削減
 並列化: DKの通信が必要→ RESPPCにより簡略化

HX  DX, EX
X
No

DK 収束?
Yes
DK (モノマー)  HX (X: ダイマー)
HX (X: ダイマー)  DX, E
End

ダイマー計算X

RESDIMによりO(N)化(静電相互作用計算のコストはnegligibly small)(9)スパコンを用いた生体分子計算
 インフルエンザの抗原抗体複合系の計算[1,2] 2350残基, 36000原子のMP2/MP3計算
 フラグメントMO法による大規模相互作用解析計算
赤: 黄色との引力相互作用
青: 黄色との斥力相互作用

阻害剤(タミフルなど)との相互作用計算にも応用[1] S. Tanaka et al., Annual Report of the Earth Simulator Center 2010, 187.
[2] A. Yoshioka et al., J. Mol. Graph. Model. 30, 110 (2011).

(10)

分割統治(DC)量子化学計算
 フラグメント(部分系)に分けて計算するのはFMOと同じ
 部分系の周辺をバッファ領域として計算に含め、環境の効果を考慮
 部分系の電子数はフェルミ準位を導入して自動決定(予め指定する必要なし)
SCFとpost-HF計算で取り扱いが大きく異なる
 分割統治SCF法 全系の密度行列を部分系の寄与の和で表現
 分割行列を導入

 分割統治post-HF電子相関計算
 部分系の相関エネルギーを部分系の分子軌道で表現(11)

分割統治量子化学計算のスキーム[1,2]
中央領域
(部分系)

局在化領域

バッファ領域

divide

部分系ごとの方程式を解く

部分系の分子軌道

共通のフェルミ準位を設定

全電子数を保存 バッファサイズにより

conquer

精度を制御
SCF

 電子が非局在化した全系の電子密度・エネルギー系の計算も可能[1] MK and H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics (2011), pp.97-127.
[2] W. Yang and T.-S. Lee, J. Chem. Phys. 103, 5674 (1995).

(12)

DC-HF/DFT法
 DC法: 密度行列を分割
Dµν =
DC

subsystem

P D
∑
α
α
µν

α
µν

Pα =

0

1/2

0

バッファ領域

1/2

1

1/2

中央領域
(部分系)

0

1/2

0

バッファ領域

局在化領域

 部分系のMOから部分系の密度行列を構築
MO(α )α≈ 2 ∑ f β (ε F − ε αp )Cµα p CναpDµνµ ,ν ∈ L (α ) : 部分系αの局在化
領域の基底

p

α

α

α

α α

占有数をFermi関数で 各部分系で決定 F C = S C ε→ Fermi準位は電子数の条件から決定

αne 2∑∑ f β (ε F − ε αp )∑ Pµνα Cµα p Cναp Sνµ
=
α

µν

p

1

fβ
0

 Fock行列とエネルギーの表式は、通常と同じ
core
DC 
+ ∑ DλσΓ µν ,σλ − 1 Γ µλ ,σν Fµν = H µν
2


λσ

EeDC

(

)

1 Tr  DDC H core + F 

2 

εF

(Hartree-Fock)(13)

DC-HF/DFT計算スキーム
core
DC 
+ ∑ DλσΓ− 1 Γ µλ ,σν 
全体のFock行列 Fµν = H µν
2 µν ,σλ


λσ

O(N2)

divideFMMなどでO(N)に各局在化領域のFock行列

O(N)

各局在化領域で対角化 Fα Cα = Sα Cα εαFermi準位 εF の決定

全電子数neを保存

α=
ne 2∑∑ f β (ε F − ε αp )∑ Pµνα Cµα p Cναp Sνµ
α

µν

i

各局在化領域の密度行列
conquerα≈ 2∑ f β (ε F − ε αp )Cµα p CναpDµνp

全体の密度行列 Dµν =DC

SCF

subsystem∑ααPµνα Dµν

エネルギー, プロパティ(14)

DC-HF計算: CPU時間
 最初のSCF cycleにかかるCPU時間[1]基底関数: 6-31G**
中央領域: C2H2(3) (1 unit)
バッファ: 6 units
Pentium4 / 3.20 GHz

 ポリエン鎖 CnHn+2
3000

Fock行列対角化

600

Fock行列生成(w/ FMM)

2500

Conventional HF
DC-HF
3.2

CPU time [sec]CPU time [sec]

800

O(n )

400
200

O(n1.3)0Conventional HF
DC-HF

2000

O(n1.6)

1500
1000

O(n1.1)

500
0

0

30

60

90

n

120

150

0

30

60

90

120

150

n

DC法によりほぼリニアスケーリングを達成[1] M. Kobayashi and H. Nakai, in Linear-Scaling Techniques in Computational Chemistry and Physics(2011), pp. 97-127.
(15)

Fermi準位の決定法
α
α
α
α α
=
n
f
−
P
C
C
2
(
ε
ε
)
 電子数保存の条件: e ∑∑ β F p ∑ µν µ p ν p Sνµ
α

µν

i

 εFを求める非線形方程式

α
α
α
α α
w
=
P
C
C
 部分系の軌道の重み p ∑ µν µ p ν p Sνµ を保存µν Fermi関数を計算する閾値を設定

1

 Brent法で求根 εFの上限と下限を設定(一方は0または前回の結果を利用)
0εF
 囲い込み、二分法、補間の組合せ−1 x  =f β ( x) exp  計算時間はO(N)だが、並列化は
 + 1kT b  容易ではない(繰り返しのたびに通信が必要)(16)

Density functional tight binding (DFTB)法[1] DFTB法: DFTをベースとした半経験的計算法
 DFTエネルギー
Tight binding近似： ρ = ρ0 + δρ
交換相関項に対して2次のTaylor展開: Exc[ρ0 + δρ]Charge dependent項に対してモノポール近似： δρA ≈ ΔqAF00Y00
DFT計算から決定するパラメータ DFTBエネルギー
Mulliken charge
Charge independent

Charge dependent

Repulsive energy

2電子積分不要!

 Fock行列
 Mulliken電荷を自己無撞着に決定 (SCC)[1] D. Porezag, Th. Frauenheim, Th. Köhler, G. Seifert, and R. Kaschner, Phys. Rev. B 51, 12947 (1995).

(17)

DC-DFTB計算のパフォーマンス
 計算時間のサイズ依存性[1] 水 Nwater 分子系

 計算機: Intel Xeon 1ノード (8コア)

計算時間がO(N)に 1ノードでも十分に高速 エネルギー勾配計算
もO(N)[1] H. Nishizawa, Y. Nishimura, MK, S. Irle, and H. Nakai, J. Comput. Chem. 37, 1983 (2016).

(18)

DC-DFTB: 並列化効率
 並列化による高速化率・・・ ポリエン鎖 C10000H10002

Buffer Subsystem(C2H2)nC2H2

Buffer(C2H2)n部分系: C2H2(3) (1 unit)
バッファ: 8 units
計算機: 京Acceleration

500ノード(4000コア)使用時の計算時間は2.2秒
 H0、S、γの計算やgradient計算の効率は90%以上

反発エネルギー計算は全計算
時間の1%未満
 ボトルネックは
SCCエネルギー計算
# nodes(19)

SCCエネルギー計算の並列化
 並列計算の手順

全部分系の軌道エネルギーを使ってFermi準位εFの求根(Brent法等を使用)

=
ne
・・・

MPI comm.
Set common Fermi level
MPI comm.
・・・

Δq

Δq

Δq

Δq

Δq

・・・

MPI comm.

f β ( ε − ε α ) wα
∑∑
α
F

p

p

p

逐次計算が繰り返し計算の中に存在
500ノード使用時の計算時間
 SCCエネルギー計算: 1.10 s
 うちFermi準位決定: 0.78 s

並列計算可能な新たなフェルミ準位
計算アルゴリズムの開発(20)

補間を用いた新たなεF計算アルゴリズム
 新アルゴリズムの手順:補間(内挿)法を利用
1. Fermi準位の推定範囲を指定
2. 各ノードで推定範囲内を等間隔にした各εFに対して、電子数を計算3. MPI_Reduceにより、εFの上限・下限を決定# electrons

Sum up# electrons# electrons

εF# electrons

εF

Correct No. electrons

εF

Approx. Fermi level

εF

4. 適切な桁まで求めたら、Spline補間でεFを決定(21)

Fermi準位の決定アルゴリズム: 比較 Fermi準位決定に要する時間のノード数依存性[1] 水4000分子系
 計算機: 京
部分系: 水1分子
バッファ: 6.0 Å

GSS (従来法)
 低ノード数の場合、第1サイクルはCRIより高速
CRI (内挿法)
 ノード数増加に伴って理想的に計算時間減少


全サイクル足すと、低ノード数でもCRIが高速

SCFが進むと内挿1回で十分な精度に[1] H. Nishizawa, Y. Nishimura, MK, S. Irle, and H. Nakai, J. Comput. Chem. 37, 1983 (2016).

(22)

DC-DFTB-K: 並列計算パフォーマンス
 計算時間の比較[1] 水256,000分子
部分系: 水1分子
バッファ: 6.0 Å 計算機: 京

Nnodes
640
1280
2560
5120
効率


Rep. 1-e int.
1.81
0.42
0.92
0.22
0.48
0.13
0.26
0.07
86%
74%

γ
SCC
Grad. Total
27.95 165.66 37.33 233.16
13.96 85.81 18.78 119.69
7.11 47.87
9.58 65.17
3.55 25.18
5.05 34.12
98%
82%
92%
85%

並列化効率: 1電子積分(0.5秒以下)を除き80%以上


特にγと勾配計算[O(N2)] は90%以上[1] H. Nishizawa, Y. Nishimura, MK, S. Irle, and H. Nakai, J. Comput. Chem. 37, 1983 (2016).

(23)

分割統治post-HF電子相関計算[1] 部分系の相関エネルギー DC-HF計算で求められる部分系(局在化領域)の軌道から相関エネルギーを計算
 中央領域だけの相関エネル
ギーの見積もりが必要

中央領域

エネルギー密度解析(EDA)[2]の利用
中央領域[1] M. Kobayashi, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007).
[2] H. Nakai, Chem. Phys. Lett. 363, 73 (2002).

バッファ領域

局在化領域(24)

分割統治post-HF電子相関計算
 相関エネルギー(Nesbetの定理)
occ vir∑∑ (ia | jb)[2t∆Eia , jb− tib , ja ]i , j a ,b

 部分系のΔEを部分系の軌道から計算∆Eαocc(α ) vir (α )

∑
i, j

α
α


−
µ
w
C
a
j
b
t
t
(
|
)[2
∑ occ ∑
ia , jb
ib , ja ]a ,bµ∈S (α )

α
µi

α

α α

EDA

S(α) : 部分系a の中央領域のAOα MP2の場合: tia , jb = −

局在化領域

α

α

α α

中央領域

(i a | j b )ε aα + ε bα − ε iα − ε αj

 全相関エネルギー = 部分系の相関エネルギーの和∆E ≈subsystem

∆E α
∑
α

DC-MP2法(25)

DC-MP2法の計算時間
 MP2計算にかかる時間

rb ＝ 6.0 Å
部分系: 1原子MP2/6-31G

 βストランドグリシンペプチド (Gly)n

Xeon (Paxville) 2.8 GHz (1CPU)

1800

O(n5.1)1600

CPU time [min]

1400
1200
1000

計算時間: 系の大きさに対してほぼ線形Canonical MP2
DC-MP2
 必要メモリ量: 系の大きさにほぼ非依存O(n1.5)


800
600
400
200

(最大部分系の大きさに依存)

0
0

5

10

15

20

25

30

n
[1] M. Kobayashi, Y. Imamura, and H. Nakai, J. Chem. Phys. 127, 074103 (2007).
(26)

DC-MP2計算の並列化
 部分系の相関エネルギー
 重複を防ぐため、中央領域

だけの相関エネルギーを算定
∆E

α

中央領域

occ vir

α
α
α α
α
α



−
µ
C
j
a
b
t
t
2
∑∑ ∑ µi ij ,ab ij ,ba 
i , j a ,b µ∈S (α ) 既存の高並列MP2アルゴリズムを利用可能 全系の相関エネルギー∆E ≈subsystem

∑
α

∆E α

 部分系ごとに独立の計算

2段階の並列化が可能!

中央領域

バッファ領域

局在化領域(27)

GDDIを用いた階層的並列処理
 Generalized distributed data interface (GDDI)
 GAMESSで多層並列処理を担うインターフェース 3種類のスコープ(MPIのコミュニケータ)を使用(1) DDI_WORLD: 全ノード = 1グループ (MPI_COMM_WORLD)
マスタ0

スレーブ
・・・

(2) DDI_GROUP: グループ分割して並列処理 (グループ内並列)
マスタ (N-1)

マスタ0

グループ内では通常の並列処理・・・
グループ0

グループ (N-1)(3) DDI_MASTERS: マスタ間で通信処理 (グループ間並列)(28)

DC-MP2の2段階並列アルゴリズム
部分系のエネルギー計算を各グループに割振り:
GDDI_SCOPE(DDI_GROUP)
グループ0

グループ (N-1)∆E α =1∆E α = N
部分系の計算はグループ内で並列化

部分系のエネルギーの足合わせ: GDDI_SCOPE(DDI_MASTERS)∆=
E

subsystem

∑
α

∆E α
グループ0

[1] M. Katouda, M. Kobayashi, H. Nakai, and S. Nagase, J. Comput. Chem. 32, 2756 (2011).

グループ (N-1)(29)

GDDI DC-MP2計算の擬似コード[1]1: 基底関数の数順に部分系をソート(ロードバランシングのため)
2: Call DDI_SCOPE(DDI_GROUP)
3: Call GDDICOUNT(-1,MYJOB)
4: EMP2TOT ← 0
5: Loop isub=1, nsub ; 部分系のループで並列化 (coarse-grain)
6:
Call GDDICOUNT(0,MYJOB)
7:
If (MYJOB=TRUE) Then
8:
EMP2 ← [MP2 correlation energy of isub subsystem](グループ内 [fine-grain] 並列化を利用)
9:
EMP2TOT ← EMP2TOT + EMP2
10:End If
11: End Loop
12: Call GDDICOUNT(1,MYJOB)
13: Call DDI_SCOPE(DDI_MASTERS)
14: Call DDI_GSUMF(EMP2TOT)
15: Call DDI_SCOPE(DDI_WORLD)[1] M. Katouda, M. Kobayashi, H. Nakai, and S. Nagase, J. Comput. Chem. 32, 2756 (2011).

(30)

2段階並列DC-MP2: 並列性能評価
 T2K-Tsukubaでの並列加速度

基底関数: 6-31G*
中央領域: AUTO (4 Å)
バッファ: 7 Å
NGROUP = Ncore / 16

(OpenMP化はしていない)– (Ala)40

512

512

Two-level
One-level

Two-level

Acceleration ratio

Acceleration ratio– βストランド (Ala)20

384

256

128

0

384

256

128

0
0

128

256

Ncore

384

512

0

128

256

384

512

Ncore

大規模系で特に高い並列計算効率を実現(31)

2段階並列DC-MP2: 「京」での性能評価
 京での並列加速度
 ポリエン鎖 C300H302
 DC-MP2/6-31G*

部分系: C2H2(3) (1ユニット)
バッファ領域: 左右8ユニット
NGROUP = Nnode / 18
SERIAL MP2アルゴリズム

 MPI+ARMCI/OpenMP hybrid

Nnode

Nthread*

FLOPS

計算時間 [s]α strong

72
144
288
576
1152

504
1008
2016
4032
8064

6.01%
5.80%
5.45%
3.99%
2.36%

4845
2478
1246
715
468

98%
97%
85%
65%

*各ノードでARMCIの通信スレッド立ち上がるため、1ノードにつき7スレッド利用
(32)

講義概要
 7/7 量子化学計算の概要と構成要素、高速化
 量子化学計算の目的と種類
 量子化学計算の手順、構成要素と高速化

 7/14 大規模系に適用するための量子化学計算法
 フラグメント分割に基づく方法
 フラグメント分子軌道(FMO)法 分割統治(DC)法 ラプラス変換MP2法
 2電子積分の密度フィッティング法
 MP2計算への応用(33)

Laplace変換MP2法[1,2]occ vir

 MP2エネルギー: ∆EMP2 =∑∑(ia | jb) [ 2(ia | jb) − (ib | ja ) ]i , j a ,bεi + ε j − εa − εb

 分母があるので、このままではO(N4)よりも小さくできない1=
 Laplace変換
x∞∫ exp(− xs)ds を利用
0 分子積分(ia|jb)もΓ になおす [O(N5)の積分変換を除去]∞∆EMP2 = − ∫ ∑ ∑ X µγ ( s )Yνδ ( s ) X κλ ( s )Yεσ ( s )Γγδ ,κε [2Γ µν ,λσ − Γ µσ ,λν ]ds
0

≡ −∫

∞

0

γδκε µνλσ∑ (µν | λσ )[2Γ µν ,λσ − Γ µσ ,λν ]dsµνλσ 積分(数値求積)が必要[1] M. Häser, Theor. Chim. Acta 87, 147 (1993).
[2] P. Y. Ayala and G. E. Scuseria, J. Chem. Phys. 110, 3660 (1999).

occ

X( s ) = ∑ eε i s Ci CiT
ivir

Y( s ) = ∑ e −ε a s Ca CTa
a(34)

Laplace変換MP2: 求積法
 最小二乗法で求積点を決定
 O(N4)の点数に対して実行するのは非効率

 Minimax法[1] 求積誤差の最大値を最小にする

 一般的な求積法
 Gauss-Laguerre Exponentialに減衰する[0, ∞]積分に有効

 Euler-Maclaurin法 台形公式の誤差を見積もる方法 Romberg積分

有限範囲への変数変換が必要 Euler-Maclaurin法の誤差への外挿法[1] A. Takatsuka, S. Ten-no, and W. Hackbusch, J. Chem. Phys. 129, 044112 (2008).

(35)

Laplace変換MP2: Euler-Maclaurin求積
 求積を台形公式で実行
1 τ
 k  1
≈
+
+
f
(
r
)d
r
f
f
(0)
f
(1)
[
]
∑
2
∫0 2
τ + 1  k =1 2  τ + 1  2 2
1

f 2 (r ) = e2 ( s )dsdr 積分範囲を有限にする変数変換が必要

 Euler-Maclaurin法による誤差の見積もり

f 2′ (0)f 2′′′ (0)f 2 (5) (0)f 2 (7) (0)
=
∆
−
+
−
+
2
4
6
8
12(τ + 1) 720(τ + 1) 30240(τ + 1) 1209600(τ + 1)

r = 0でヤコビアンが0になるような変数変換を利用
r 2 − (r 3 + r 4 + r 5 + 0.9r 6 ) 4  ds
s
0
=
=
(A)


(1 − r ) 2dr
 r =02
r 3 − 0.9r 4
d
s
d
s
2
0,=
0
=
+ r tan (π r / 2 ) =
(B) s

2
2

(1 − r )dr r =0 dr r =0(36)

Laplace変換MP2: 計算手順
 求積点ごとに以下を実行 (求積点: s, 重み: w)
X( s ) = ∑ e C C Y ( s ) = ∑ e
1. 行列X(s)とY(s)を求める
occ

vir

εi s

i

T
i

i

−ε a sCa CTaa

2. Schwarzのスクリーニングに用いる行列を求める
3. 各κεに対し、( µν | κε ) を求めてディスクに保存Γを計算し、( µδ=| κε ) ∑ X µγ Γγδ ,κε を足しこみγ( µν | κε ) = ∑ Yνδ ( µδ | κε ) を足しこみ3～6でSchwarzのδ
不等式を利用した
(
µν
|
λε
)
=
X
(
µν
|
κε
)
4.
を足しこみ
∑κ λκ
スクリーニング
5. ( µν | λσ ) = ∑ Yσε ( µν | λε ) を足しこみσ

6. [2Γ µν ,λσ − Γ µσ ,λν ] を求めてエネルギーに足しこみ(37)Laplace変換MP2: 計算時間
 MP2計算時間
 ポリエン鎖 CnHn+2
Computational time [hour]

250[Pentium4/3.0 GHz]

200

Canonical MP2
Laplace MP2150

計算時間の削減に成功
 nに対するスケー
リングも改善


100
50
0

0

30

60

n

90

120

150
(38)

Laplace変換MP2: 求積法の精度
 MP2相関エネルギーの求積法依存性
 ベンゼン/6-31G*

求積法
Gauss-Laguerre
Euler-Maclaurin (A)Euler-Maclaurin (B)Romberg (A)Romberg (B)Canonical MP2[Hartree]

求積点数
5
5
5
7
7

Ecorr
-0.733451-0.770241-0.784540-0.784643-0.783803-0.784761(diff.)(+0.051311)(+0.014521)(+0.000221)(+0.000118)(+0.000958)

Euler-Maclaurin (B)やRombergが良い結果
誤差解析の結果にも対応


[1] M. Kobayashi and H. Nakai, Chem. Phys. Lett. 420, 250 (2006).

(39)

2電子積分の密度フィッティング(RI近似)
 2電子積分 Γ µν λσ =∫∫ dr dr φµ (r )φν (r )r φλ (r )φσ (r )
,

1

2

1

1

−1
12

2

2

 メモリにストアすることは困難 (4階テンソル) 原子軌道の積 φµ (r )φν (r ) を補助基底関数で展開
φµ (r )φν (r ) ≈ ∑ d mµν ϕm (r ) ≡ θ µν (r )m

 誤差の自己反発積分を最小化するように決定∫∫ dr dr1

2Rµν (r1 ) Rµν (r2 )r12=Rµν (r ) ψ µ (r )ψ ν (r ) − θ µν (r )→ Min

d mµν = ∑ (m | n) −1 (n | µν )n∑

( µν | m)(m | n) (n | λσ ) まとめると Γ µν ,λσ =−1m,n2階・3階のテンソルの積和(40)

並列RI-MP2アルゴリズム[1]occ vir

 MP2エネルギー: ∆EMP2 =∑∑(ia | jb) [ 2(ia | jb) − (ib | ja ) ]εi + ε j − εa − εbi , j a ,b

 (m|n)は正定値 → Cholesky分解 ( m | n) =T
LL∑ ml ln (LAPACK)l

Cµ Cν Cλ Cσ Γ µν λσ を考慮∑µνλσ
Cµ Cν ( µν | l ) ≡ ∑ L (ia | l )(ia=| jb) = ∑ B B B
∑L ∑
µν

=
(ia | jb) 積分変換
ia
n

jb
n

i

j

−1
nl

ia
n

n

a

l

b

,

i

−1
nl

a

l

 並列RI-MP2アルゴリズム
1. (m|n)を計算し、Cholesky分解でL-1を計算して保存
2. (ia|l)を計算 (lに対して動的並列化)(ia|l)とBのデータ分散もia
3. Bn を計算 (iに対して静的並列化)
4. (ia|jb)を計算 (ijに対して静的並列化) 、MP2エネルギー積算[1] M. Katouda and S. Nagase, Int. J. Quantum Chem. 109, 2121 (2009).

(41)

RI-MP2法: 計算コスト

O

O
H
N
O

O

NH

 MP2計算時間とメモリ容量[1]HN

Elapsed time [min]

O
O

O

O
O

O
H
N

[Pentium4 640 (3.2 GHz)]

NH
O
O

メモリ
[MB/node]

スクラッチ
[GB]

RI-MP21206.1

32.4

Canonical

2401.4639.3



スーパーリニアスケーリングを達成
 RI近似による誤差は0.118 mEh



HN

O

O

# CPU

O

O

 バリノマイシン(右図) / 6-311G**
3500
3000
2500
2000
1500
1000
500
0

O

O

リソース量を大幅削減

FMO法やDC法と組み合わせて利用し、更なる高速化も可能[1] M. Katouda and S. Nagase, Int. J. Quantum Chem. 109, 2121 (2009).

(42)本日のまとめ(1) 「京」を使うだけでは大規模量子化学計算は不可能
 計算コストのオーダーを削減するさまざまな手法
 フラグメント分割法 (FMO, DC) フラグメントの計算結果を足し合わせて全体の結果を得る
 大規模並列化が可能

 Laplace変換MP2法
 摂動論で現れるエネルギー分母をLaplace変換で消去
 Schwarz不等式等のカットオフを利用してオーダー削減

 2電子積分の密度フィッティング法 (RI法) 4階のテンソルを3階以下のテンソルの積和で表現
 計算リソースを大幅削減 (時間オーダーは不変)(43)

本日のまとめ(2) 並列化を見据えたアルゴリズムの改善
 DC法におけるFermi準位決定 計算時間はごく短いが、DFTB法では問題に
 一見非効率な方法も、並列化した場合には良いことも 2段階並列アルゴリズム
 計算を大粒度で並列化し、その中身をさらに細粒度並列
 フラグメント分割計算では非常に有効 データ分散が可能なアルゴリズム
 うまく組めばスーパーリニアな並列効率が出る場合も(44)