量子化学計算の大規模化2
石村 和也(分子科学研究所 理論・計算分子科学研究領域)
ishimura@ims.ac.jp
2017年度計算科学技術特論A 第15回
2017年7月27日

本日の内容• 高速化・並列化例– 2次の摂動(MP2)計算 (MPI並列)– Hartree-Fock計算 (MPI/OpenMP並列)• 新たなオープンソースソフトウェアの開発– 大規模並列量子化学計算プログラムSMASH• メニーコア時代に向けた開発– OpenMP並列効率のさらなる向上– メモリ使用量削減– 通信時間削減• まとめ• 参考資料
1

高速化
•

•

•

•

演算量削減
– 演算量の少ない1,2電子積分計算アルゴリズムの開発– 効果的なカットオフの導入– 対称性の導入
近似の導入– FMO, DC, ONIOM, QM/MMなど分割法– ECP, Frozen core, 局在化軌道など化学的知見の利用
収束回数の削減– DIIS、SOSCF法などによるSCFサイクル数削減– 初期Hessian改良による構造最適化回数の削減
実行性能の向上– SIMD演算の利用– 時間のかかる演算回数を削減– データアクセスを効率的にしたアルゴリズム・プログラム開発– BLAS, LAPACKなど数学ライブラリの利用

2

並列化• ノード間(MPI)、ノード内(OpenMP)それぞれで並列化– 式の変形– 多重ループの順番の工夫• 均等な計算負荷分散– ノード間で分散、さらにノード内でも分散• 大容量データの分散– 京のメモリは1ノード16GB, 8万ノードでは1PB以上– ノード間では分散、ノード内では共有– 中間データ保存用としてディスクはあまり期待できない• 通信量、通信回数の削減– 並列計算プログラムのチューニングで最後に残る問題は通信関係(特にレイテンシ)が多い• 高速化と並列化の両立
3

2次の摂動(MP2)法•••

Hartree-Fock計算で分子のエネルギーの約99%を求めることができるが、定量的な議論を行うためには
残り1%の電子相関エネルギーが重要
MP2法は最も簡便な電子相関計算方法
積分変換(密行列-行列積)計算が中心

Hartree-Fock計算(mn|ls)計算 (O(N4))

𝑜𝑐𝑐 𝑣𝑖𝑟(mi|ls)計算 (O(N5))

𝑖𝑗(mi|lj)計算 (O(N5))

𝑎𝑖|𝑏𝑗 2 𝑎𝑖|𝑏𝑗 − 𝑎𝑗|𝑏𝑖
𝐸𝑀𝑃2 = ෍ ෍𝜀𝑖 + 𝜀𝑗 − 𝜀𝑎 − 𝜀𝑏𝑎𝑏 𝐴𝑂

𝑎𝑖|𝑏𝑗 = ෍ 𝐶𝜇𝑎 𝐶𝜈𝑖 𝐶𝜆𝑏 𝐶𝜎𝑗 (𝜇𝜈|𝜆𝜎)
分子軌道(MO) 𝜇𝜈𝜆𝜎
2電子積分

原子軌道(AO)
2電子積分

𝜀𝑖 : 軌道エネルギー, 𝐶𝜇𝑎 : 分子軌道係数𝑖, 𝑗: 占有軌道, 𝑎, 𝑏: 非占有軌道(ai|lj)計算 (O(N5))(ai|bj)計算 (O(N5))
MP2エネルギー計算 (O(N4))

4これまでのMP2並列計算アルゴリズム•

AOもしくはMOインデックス分散– AO:複数プロセスにある部分的なMO2電子積分の足し合わせ
総通信量がプロセス数に依存– MO:同じAO2電子積分を複数のプロセスで計算 or Broadcast
不完全な計算負荷分散 or 総通信量がプロセス数に依存
R. A. Whiteside, J. S. Binkley, M. E. Colvin, H. F. Schaefer (1987) (32CPU)I. M. B. Nielsen, C. L. Janssen (2000) (MPI + pthreads)••

Global arrays, ARMCI, DDI– 分散メモリを仮想的に共有メモリのように扱うライブラリ
前半AO、後半MOインデックス分散– 完全な計算負荷分散 and 総通信量がプロセス数に関わらずほぼ一定
J. Baker, P. Pulay (2002) (通信のためのデータソートに時間がかかる)

5新規MP2エネルギー並列計算アルゴリズム
K. Ishimura, P. Pulay, S. Nagase, J Comput Chem 2006, 27, 407.

•••

MPIのみで並列化、前半はAO、後半はMOインデックスを分散
シンプルなデータソーティング
ノード数にかかわらず、全体の通信量はほぼ一定
m (AO index) 各ノードに分散
do l,s
AO積分計算 (mn|ls) [n, mls] (all n)
第1変換(mi|ls) [i, mls] (all i)
第2変換(mi|lj)[ij, l, m] (i≧j)
end do l,s
第3変換(mi|bj)[b, ij](all b)(mi|bj)をディスクに書込み [b, ij, m]
end do mij (MO index) 各ノードに分散(mi|bj)をディクスから読込み + MPI_isend,irecv
第4変換(ai|bj)[b, a](all a,b)MP2エネルギー計算
end do ij
𝜇, 𝜈, 𝜆, 𝜎: 原子基底関数,
𝑖, 𝑗: 占有軌道, 𝑎, 𝑏: 非占有軌道

6

BLASルーチンの使い方•

BLAS level3のDGEMM(倍精度行列-行列積)を用いた演算と多次元配列のインデックス入れ替え(データの並び替え)– DGEMMを用いて転置した行列での演算– 通信におけるデータソーティングをシンプルにするため、ノード間で分散されているインデックスを外側に移動
𝐴𝑂– 例：第3積分変換 𝜇𝑖|𝑏𝑗 = ෍ 𝐶 𝜇𝑖|𝜆𝑗
𝜆𝑏

配列の並び
第2変換 (mi|lj)
第3変換 (mi|bj)

𝜆[ij, l, m][b, ij, m](i≧j)(all b)

ソースコードdo mcall dgemm(‘T’,’T’,...,C,...,(mi|lj),...,(mi|bj),..)
enddo

7

DGEMMでの行列転置• A(M,K), B(K,N)の場合、C=ABの計算でCの配列は(M,N)と(N,M)の2通り可能– 転置しない場合：DGEMM(‘N’,’N’,...,A,...,B,...,C,...)
N
M

C1

=

A

B

– 転置する場合：DGEMM(‘T’,’T’,...,B,...,A,...,C,...)
M

N

C2

=

B†

A†

8

MP2エネルギーテスト計算
計算機: Pentium4 3.0GHz
プログラム：GAMESS
分子: Taxol(C47H51NO14)
6-31G(d)6-311G(d,p)

基底次元数

1032

1484

占有軌道数

164

164

非占有軌道数8061258

1プロセス当たりメモリ使用量

0.67GB

0.96GB

1プロセス当たりディスク使用量

90GB/nproc

202GB/nproc

90GB

202GB

全体の通信量
CPU数

1

2

4

8

16

6-31G(d) (1032 AOs)
実時間 (hour)

10.2

5.08

2.54

1.31

0.64

Speed-up

1.0

2.0

4.0

7.8

15.8

6-311G(d,p) (1484 AOs)
実時間 (hour)

31.6

16.3

8.06

4.05

2.05

Speed-up

1.0

1.9

3.9

7.8

15.4

9

MP2エネルギー微分計算式
 x oallxx 2  H pq   2( pq | kk  ( pk | qk  Ppq( 2)
pq 
k

all

E

x
MP2

oall

 2 S W
x
ij

i

( 2)
ij

vall

I   2 S W I   4 Saix Wai( 2) I 

j

a

oall

S W
x
ij

( 2)
ij

a

ab

oall

II    S W II   4 S W II    SijxWijx III   2 tijab (ia | jbx
( 2)
ab

x
ai

i

( ja | kb, P( 2)  
Dab
jk

i

oall vall

x
ab

ab

Pij( 2)   tikab

( 2)
ab

b

vall

ij

k

oall

x
ab

iJ

tikab
k

( 2)
ai

a

ij

(Ja | kb, P( 2) 
i   Jabab

tijacij

ij(ib | jc, P( 2) 
Dijbc

c

aB

ab

tijacijc(iB | jca  B

Wij( 2) I    tikab ( ja | kb, Wab( 2) I    tijac (ib | jc, Wai( 2) I    t abjk (ij | kb
k

ab

ij

c

jk

b

II   P ( i   j , W II   P ( a   b , W II   P  , W III    Ppq(2) Apqij
MO

( 2)
ij

W

( 2)
ij

( 2)
ab

( 2)
ab

( 2)
ai

( 2)
ai
i

( 2)
ij

pq
oall

vall

jk

bc

bc
Lai   Pjk( 2) 4(ia | jk   2(ik | aj    Pbc( 2) 4(ia | bc  2(ib | ac  2 Na  t abjk (ij | kl  2 Ni  tij (ab | jc
jk

b

j

bc

10

新規MP2エネルギー微分並列計算アルゴリズム
K. Ishimura, P. Pulay, S. Nagase, J. Comput. Chem., 2007, 28, 2034.

do l (積分変換計算)
AO積分計算 (mn|ls
第1変換 (mn|lj)
第2変換 (mi|lj)
第3変換 (ai|lj)
ディスク書込み (ai|lj)
Enddo ldo s (MP2ラグランジアン計算)
ディスクから読込み t s
ns
計算とディスク書込み t
AO積分計算 (mn|ls
第1変換 (mn|js)
各項計算 L1mn, 2 , L4mienddo saijij

do a(MP2エネルギー、アンプリチュード計算) CPHF計算
ディスク読込みと送受信 (ai|lj)
do m (AO積分の微分計算)
第4変換 (ai|bj), (ai|kj)
H mx n , S mx n
微分項計算
MP2エネルギーとアンプリチュード計算
tns
tns
ディスク読込みと計算
smlij
各項計算 t , P , P ,W I ,W I ,W I , ts
微分項計算 (mn|ls)xt
送受信とディスク書込み
enddo menddo a
ab
ij

( 2)
ij

( 2)
ab

( 2)
ij

( 2)
ab

( 2)
ai

a
ij

a
ij

𝜇, 𝜈, 𝜆, 𝜎: 原子基底関数,

𝑖, 𝑗, 𝑘: 占有軌道, 𝑎, 𝑏: 非占有軌道

11MP2エネルギー微分計算アルゴリズムのポイント•••

全ての計算の負荷分散+データ分散+高速化+通信量削減
ノード数にかかわらず、全体の通信量はほぼ一定
演算量、メモリ使用量が少なくなるように式を展開– 例えば 𝑖𝑎|𝑏𝑐 (𝑖:占有軌道、 𝑎, 𝑏, 𝑐:非占有軌道)の積分変換はコストがかかるので、全てを変換せずに計算– 積分変換をしないために出てくる分子軌道係数𝐶の処理を、あらかじ
め行う(𝑃𝜆𝜎 ← 𝑃𝑏𝑐 )𝐴𝑂

𝐿′𝑎𝑖 = ෍ 𝑃𝑏𝑐 𝑖𝑎|𝑏𝑐 = ෍ ෍ 𝐶𝝀𝑏 𝐶𝝈𝑐 𝑃𝑏𝑐 𝑖𝑎|𝝀𝝈
𝑏𝑐

𝑏𝑐 𝝀𝝈
𝐴𝑂

= ෍ 𝑃𝝀𝝈 𝑖𝑎|𝝀𝝈
𝝀𝝈
𝐴𝑂

𝐴𝑂

= ෍ 𝐶𝝂𝑎 ෍ 𝑃𝜆𝜎 𝑖𝝂|𝜆𝜎𝝂𝑃𝜆𝜎 = ෍ 𝐶𝜆𝑏 𝐶𝜎𝑐 𝑃𝑏𝑐
𝑏𝑐

𝜆𝜎

𝑋𝑖𝝂
12

MP2エネルギー微分テスト計算
計算機: Pentium4 3.0GHz
プログラム：GAMESS
分子: Taxol(C47H51NO14)
6-31G

6-31G(d)

基底次元数

660

1032

占有軌道数

164

164

非占有軌道数434806

1プロセス当たりメモリ使用量

0.78GB

1.84GB

1プロセス当たりディスク使用量

147GB/nproc

426GB/nproc

147GB

426GB

全体の通信量

CPU数

1

2

4

8

16

32

実時間 (hour)

17.4

8.09

3.82

1.92

0.97

0.53

Speed-up

1.0

2.1

4.5

9.0

17.9

33.0

実時間 (hour)

31.1

15.1

7.57

3.86

2.05

Speed-up

2.0

4.1

8.2

16.1

30.4

6-31G (660基底)

6-31G(d) (1032基底)

13

FMO-MP2法
D. Fedorov, K. Ishimura, T. Ishida, K. Kitaura, P. Pulay, S. Nagase, J. Comput. Chem. 2007, 28, 1476.

•••

開発したMP2ルーチンとFMOルーチンの組み合わせ
分子：Trp-cage設計蛋白質TC5b (C98H150N27O29)
GAMESSプログラムに実装
基底関数：6-31(+)G* 2480基底
6-311(+)G* 3246基底
FMO計算：1フラグメントに2残基
メモリ使用量：FMO-MP2 数百MB/プロセス
分割無しMP2 4.7GB/プロセス

14

FMO-MP2計算結果
FMOn (n体相互作用まで計算)
MP2とのエネルギー誤差(kcal/mol)
6-31(+)G*
6-311(+)G*

FMO2
2.1
-2.9

FMO3
-0.2
-0.5

計算時間(hour)
6-31(+)G*
6-311(+)G*
•
•

FMO2

FMO3

3.2
7.6

22.7
56.8

分割無し
MP2
20.4
46.5

困難であった巨大分子FMO-MP2計算の精度確認が可能になった
FMO-MP2の高速化・高並列化
15

Hartree-Fock法
HartreeFock計算

原子軌道の線形結合係数(分子軌道係数)を求める

原子軌道Gauss関数

FC  εSC
Fock行列F: Fock行列, C: 分子軌道係数
S: 基底重なり行列, : 分子軌道エネルギー

Fmn  H mn   Cli Csi 2(mn | ls   (ml | ns 
i ,l ,s

原子軌道(AO)2電子積分

1𝜇𝜈|𝜆𝜎 = න 𝑑𝒓1 න 𝑑𝒓2 𝜙𝜇 𝒓1 𝜙𝜈 𝒓1
𝜙 𝒓 𝜙 𝒓𝑟12 𝜆 2 𝜎 2

𝜙𝜇 𝒓1 : 原子軌道Gauss関数

初期軌道係数C計算

AO２電子反発積分計算+
Fock行列への足し込み (O(N4))
Fock行列対角化 (O(N3))
分子軌道C収束
計算終了

分子軌道
収束せず

16

MPI/OpenMP並列アルゴリズム1
K. Ishimura, K. Kuramoto, Y. Ikuta, S. Hyodo, J. Chem. Theory Comp. 2010, 6, 1075.

Fock行列Fmn  H mn   Cli Csi 2(mn | ls   (ml | ns 
i ,l ,s

原子軌道(AO)2電子積分

MPI/OpenMPハイブリッド並列化後
並列化前

do m=1, nbasis
do n=1, m
do l=1, m
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo

!$OMP parallel do schedule(dynamic,1) reduction(+:Fock)do m=nbasis, 1, -1<-- OpenMPによる振り分け
do n=1, m
mn=m(m+1)/2+n
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc<-- MPIランクによる振り分け
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)
17

MPI/OpenMP並列アルゴリズム2
•
•

•
•

開発当初は、GAMESSプログラムに実装
1番目のループをOpenMP、3番目のループをMPIで並列化– MPIランクはプロセスに固有の値なので、 MPIランクによるプロセス間の分散箇所が決まれば、OpenMP並列が後でも先でも、各プロセスに割り当てられる仕事は変わらない
MPIランクとmod計算だけで、IF文を使わずにプロセス間分散– MPIランクによる分散までの演算はすべてのプロセスが実行
MPI通信はOpenMP領域外で実行!$OMP parallel do schedule(dynamic,1) reduction(+:Fock)do m=nbasis, 1, -1<-- OpenMPによる振り分け
do n=1, m
mn=m(m+1)/2+n
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc<-- MPIランクによる振り分け
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)

18

MPI/OpenMP並列アルゴリズム3•••

OpenMPの分散を最外ループで行うことにより、スレッド生成や分散などのオーバーヘッド(余分なコスト)を削減
演算量の多いインデックスから動的に割り振ることでプロセス内の負荷分散を均等化
すべての変数について、スレッド間で共有するか(shared)、別々の値を持つか(private)を分類– privateにすべきcommon変数は、threadprivateを使わずにサブルーチンの引数に変更– 引数の数を減らすため、x、y、zなどのスカラ変数をxyz配列に書き換え!$OMP parallel do schedule(dynamic,1) reduction(+:Fock)do m=nbasis, 1, -1<-- OpenMPによる振り分け
do n=1, m
mn=m(m+1)/2+n
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc<-- MPIランクによる振り分け
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)19

初期軌道計算などの高速化・並列化••

初期軌道計算のハイブリッド並列化
軌道の射影の高速化・並列化
1
11



1
11

C1  S S12C2 C S S S12C2
t
2

t
12



1

2

C1: 初期軌道係数
行列-行列積の多用で C2: 拡張Huckel法による軌道係数
S11: 実際の計算で用いる基底の重なり積分
高速化・並列化
S12: 拡張Huckel法で用いた基底と実際の計算で用いる基底との重なり積分
D. Cremer and J. Gauss, J. Comput. Chem. 7 274 (1986).

•

SCF計算の途中は対角化をせずに、Newton-Raphsonを基にした
Second-Order SCF法を採用– 対角化はSCFの最初と最後の2回のみ

20

Hartree-Fockエネルギーテスト計算条件
計算機: Cray XT5 2048 CPUコア(Opteron 2.4GHz, Shanghai, 8コア/ノード)
コンパイラ: PGI fortran compiler-8.0.2
BLAS,LAPACKライブラリ: XT-Libsci-10.3.3.5
MPIライブラリ: XT-mpt-3.1.2.1 (MPICH2-1.0.6p1)
プログラム: GAMESS
分子: TiO2クラスター(Ti35O70)(6-31G, 1645 functions, 30 SCF cycles)21

TiO2クラスター計算結果
1400.0MPI/OpenMP
ハイブリッド化

1200.0

全計算時間の
並列加速率
並列加速率

1000.0

初期軌道計算改良

800.0
600.0
400.0
オリジナルGAMESS(Flat
MPI)
オリジナルGAMESS (MPIのみ)

200.0

改良版 (MPI MPI)
のみ)
本研究(Flat
改良版 (MPI/OpenMP)
本研究(MPI/OpenMP)

0.0
0

512

1024
1536
CPUコア数

2048

Table 全計算時間(秒)と並列加速率(カッコ内)
CPUコア数オリジナル
MPIのみ
GAMESS
改良版

MPIのみ

改良版

MPI/OpenMP

16
18176.4(16.0)
18045.6(16.0)
18121.6(16.0)

256
1368.6(212.5)1241.2(232.6)1214.6(238.7)1024
527.6(551.2)428.7(673.5)381.1(760.8)2048
383.5(758.3)273.7(1054.9)
234.2(1238.0)

22

TiO2クラスター計算の解析• Fock行列計算では、MPI/OpenMPが最も速く並列加速率も高い• オリジナルの初期軌道計算時間は、2048コアで全体の40%弱を占める
Table Fock行列計算時間(秒)と並列加速率(カッコ内)
CPUコア数オリジナル
MPIのみ
GAMESS
改良版

MPIのみ

改良版

MPI/OpenMP

16
17881.8(16.0)
17953.5(16.0)
17777.6(16.0)

2561175.2(243.5)1175.2(244.4)1150.4(247.3)1024
334.0(856.6)360.0(797.9)316.4(899.0)2048
188.6(1517.0)
203.1(1414.4)
174.8(1627.2)

Table 初期軌道計算時間(秒)と全計算に占める割合(カッコ内)
CPUコア数オリジナル
MPIのみ
GAMESS
改良版

MPIのみ

改良版

MPI/OpenMP

16
166.2
(0.9%)
20.2(0.1%)
18.6(0.1%)

256
143.6(10.5%)
18.6(1.5%)
13.2(1.1%)

1024
143.6(27.2%)
18.9(4.4%)
13.6(3.6%)

2048
143.8(37.5%)
19.2(7.0%)
13.8(5.9%)
23

Hartree-Fock計算ハイブリッド並列化のまとめ••

使用CPUコア数が多くなると、ハイブリッド並列化の効果が顕著に出た
元々1%以下の初期軌道計算が、2048コアでは約4割を占めた– すべての計算を高速化・並列化する必要がある•

OpenMP導入のため、GAMESSの大幅な書き換えを行った– 多くのCommon変数をサブルーチンの引数に変更したため、Hartree-Fock
及びDFT計算しか実行できないコードになった– common変数が多く使われているプログラムにOpenMPを導入して、計算時間削減と並列化効率向上を両立させるのは大変✓ GAMESSに実装しても、特定の計算しか実行できないためそのソース
コードを公開することができない✓ MPIとOpenMPのハイブリッド並列をさらに進めるには、設計段階からしっかりと考慮した新たなプログラムが必要

24

新たなオープンソースソフトウェアの開発
•
•

•
•

キーワードはシンプル(実行方法、入力形式、ライセンス、開発)
MPI/OpenMPハイブリッド並列化と高速計算アルゴリズムを組み込み、一つのプログラムでスカラ型CPU搭載計算機をカバー– 演算量削減、計算負荷とデータの均等な分散、通信の最適化– 今後ノード当たりのコア数はさらに増加– 京、おそらくポスト京も研究室レベルの計算機と基本的な構成は同じ
よく用いられるルーチンのライブラリ化・モジュール化
オープンソースライセンスで配布– ますます複雑になる計算機を理解した上で、最適な式、アルゴリズ
ム、近似を開発、選択してプログラムを作る必要があり、開発コスト削
減は不可欠– 複数の人が同じような開発やチューニングをしない仕組み作り– 技術・ノウハウの共有– 計算機科学との連携
25

SMASHプログラム
•
•
•
•
•
•
•
•
•
•
•
•

大規模並列量子化学計算プログラムSMASH (Scalable Molecular Analysis
Solver for High performance computing systems)
オープンソースライセンス(Apache 2.0)
http://smash-qc.sourceforge.net/
2014年9月1日公開
対象マシン：スカラ型CPUを搭載した計算機(PCクラスタから京コンピュータまで)
エネルギー微分、構造最適化計算を重点的に整備
現時点で、Hartree-Fock, DFT(B3LYP), MP2計算が可能
MPI/OpenMPハイブリッド並列を設計段階から考慮したアルゴリズム及びプ
ログラム開発 (Module変数、サブルーチン引数の仕分け)
言語はFortran90/95
1,2電子積分など頻繁に使う計算ルーチンのライブラリ化で開発コスト削減
電子相関計算の大容量データをディスクではなくメモリ上に分散保存
2014年9月からの約3年間で、ダウンロード数は1137
26

SMASH ウェブサイト

最新ソースコードと日本語マニュアル
27

AO2電子積分ルーチンのインターフェイス••

2電子積分ルーチン：データはすべて引数で受け渡し、ブラックボックス化
call int2elec(twoeri, exijkl, coijkl, xyzijkl, nprimijkl, nangijkl, nbfijkl, maxdim,
mxprsh, threshex)

twoeriexijklcoijklxyzijkl
nprimijkl
nangijkl
nbfijkl
maxdimmxprshthreshex

2電子積分計算値 (Output)
primitive関数の指数 (Input)
primitive関数の係数
xyz座標
primitive関数の数
軌道角運動量(s=0, p=1, d=2,...)
基底関数の数(s=1, p=3, d=5or6,...)
最大twoeriの次元数
最大primitive関数の数
exp(-x2)計算の閾値28

プログラム開発の効率化•••

2電子積分ルーチンはHartree-Fock,DFT計算以外の高精度電子相関計
算でも利用
2電子積分の微分は軌道角運動量の異なる2電子積分の線形結合になるため、微分計算で2電子積分ルーチンを再利用可能
適切な係数と軌道角運動量を引数で渡し、2電子積分ルーチンをcallした
後適切な要素へ結果を足しこむだけで実装完了
call int2elec(twoeri, exijkl, coijkl, xyzijkl, nprimijkl, nangijkl, nbfijkl, maxdim,
mxprsh, threshex)
例：(𝒑𝒙 𝑠|𝑠𝑠)の微分計算の場合
𝜕 𝑝𝑥 𝑠 𝑠𝑠 /𝜕𝑋𝑎 = 𝟐𝜶𝒂 𝒅𝒙𝒙 𝑠 𝑠𝑠 − (𝒔𝑠|𝑠𝑠)
𝜕 𝑝𝑥 𝑠 𝑠𝑠 /𝜕𝑌𝑎 = 𝟐𝜶𝒂 𝒅𝒙𝒚 𝑠 𝑠𝑠
𝜕(𝑝𝑥 𝑠|𝑠𝑠)/𝜕𝑍𝑎 = 𝟐𝜶𝒂 (𝒅𝒙𝒛 𝑠|𝑠𝑠)29

MPI/OpenMPハイブリッド並列アルゴリズム
Hartree-Fock計算

MP2計算!$OMP parallel do schedule(dynamic,1) &!$OMP reduction(+:Fock)do m=nbasis, 1, -1<- OpenMP
do n=1, m
mn=m(m+1)/2+n
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc <- MPIランク
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)do ml (AO index pair) MPIランクによる振り分け!$OMP parallel do schedule(dynamic,1)
do s
AO積分計算(mn|ls) (all n)
第1変換(mi|ls)(all i)
enddo!$OMP end parallel do
第2変換(dgemm) (mi|lj)(i≧j)
end do ml
do ij (MO index pair)
MPIランクによる振り分け
MPI_sendrecv (mi|lj)
第3変換(dgemm) (mi|bj)(all b)第4変換(dgemm) (ai|bj)(all a)MP2エネルギー計算
end do ij
call mpi_reduce(MP2エネルギー)••

膨大な中間データをすべてメモリ上に分散保存
BLAS level3(dgemm)で効率的な演算とスレッド並列化を実現30

B3LYPエネルギー並列計算性能• 10万コアで5万倍のスピードアップ、実行性能13%• 10万コアで360原子系のB3LYP計算が2分半• 行列対角化(LAPACK,分割統治法)3回分の時間は約35秒→ScaLAPACK、EigenExaなどプロセス並列化されているライブラリ導入が必要
60000
50000

Speed-up

40000
30000
20000
10000
0
0

24576

49152

73728

98304

計算機 : 京コンピュータ
分子 : (C150H30)2 (360原子)
基底関数 : cc-pVDZ(4500基底)
計算方法 : B3LYP
SCFサイクル数 : 16

CPUコア数
CPUコア数

6144

12288

24576

49152

98304

実行時間 (秒)1267.4

674.5

377.0

224.6

154.2
31

Hartree-Fockエネルギー1ノード計算性能• 1ノードでは、GAMESSよりHartree-Fock計算時間を最大40%削減• SMASHではS関数とP関数を別々に計算するため、SP型の基底では
GAMESSとの差は小さい
GAMESSとの比較• Xeon E5649 2.53GHz 12core、1ノード利用• Taxol(C47H51NO14)のHartree-Fock計算時間(sec)• 同じ計算条件(積分Cutoffなど)
基底関数
6-31G(d) (1032 functions)GAMESS
706.4

cc-pVDZ (1185 functions)2279.9SMASH
666.6
1434.332

B3LYPエネルギー微分並列計算性能• エネルギー計算と異なり対角化計算が無いため、並列化
効率はほぼ100%
16384.0

Speed-up

12288.08192.04096.00.0
04096 8192 12288 16384

計算機 : 京コンピュータ
分子 : (C150H30)
基底関数 : cc-pVDZ (2250 functions)
計算方法 : B3LYP

CPUコア数

Table B3LYPエネルギー1次微分計算時間(秒)と並列加速率(カッコ内)
CPUコア数微分計算全体1024
402.0(1024.0)

4096
101.2
(4067.7)8192
50.8(8103.3)

16384
25.5
(16143.1)
33

MP2エネルギー計算性能• 通信量と回数を削減し、1万CPUコアでも高い並列性能を実現

18432.0

計算機 : 京コンピュータ
分子 : C150H30 (180原子)
基底関数 : 6-31G(d) (2160基底)計算方法 : MP2

Ｓｐｅｅｄ－ｕｐ13824.09216.04608.00.0
04608

9216 13824 18432

CPUコア数

CPUコア数

4608

6912

9216

18432

MP2計算時間(秒)

152.5

105.7

83.4

46.9

Speed-up

4608.06648.28425.914983.434

構造最適化回数の削減• Redundant座標と力場パラメータを使い、初期ヘシアンを改良することで最適化回数が1/5から1/6になった• 2サイクル目以降のヘシアンはBFGS法で更新
Table B3LYP/cc-pVDZ構造最適化回数(初期構造HF/STO-3G)
Cartesian
Redundant
座標
座標
Luciferin(C11H8N2O3S2)63
11
Taxol (C47H51NO14)
203
40

Taxol (C47H51NO14)

Luciferin(C11H8N2O3S2)
35

メニーコア時代に向けた開発
•

•

•

•
•

OpenMP並列効率の向上– ノードあたりのコア数は大幅に増えると予測される– 場合によっては、演算量を増やしても同じデータへのアクセス競合やReductionを削減した方が計算時間削減につながる
使用メモリ量の削減– ノードあたりの演算性能向上に比べて、メモリ量はそれほど増えない可能性が高い
通信時間の削減– 計算と通信のオーバーラップが重要になる (富士通FX100では1ノードに32演
算コア+2アシスタントコア)– 通信レイテンシの向上はあまり望めないと予測される
SIMD長の増加
開発コストの削減– モジュール化、ライブラリ化の推進• 分野全体での共通ルーチンの共有によるコスト削減• 情報とノウハウ共有– Pythonを用いた効率的な開発
36

OpenMP並列効率のさらなる向上1•

振り分けるタスクの数を増やす (タスクの粒度を小さくする)例) Hartree-Fock計算
改良前

改良後

!$OMP parallel do schedule(dynamic,1) &!$OMP reduction(+:Fock)do m=nbasis, 1, -1
do n=1, m
mn=m(m+1)/2+n
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc <- MPIランク
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)!$OMP parallel do schedule(dynamic,1) &!$OMP reduction(+:Fock)do mn=nbasis*(nbasis+1)/2, 1, -1
mnからmとnを逆算
lstart=mod(mn+mpi_rank,nproc)+1
do l=lstart, m ,nproc <- MPIランク
do s=1, l
AO2電子積分(mn|ls)計算+Fock行列に足し込み
enddo
enddo
enddo
!$OMP end parallel docall mpi_allreduce(Fock)
37

OpenMP並列効率のさらなる向上2••

計算機：Fujitsu PRIMEGY CX2500 (Xeon E5-2697, 28コア/ノード)Intel VTune Amplifier2015を用いた解析– 計算条件：Luciferin(C11H8N2O3S2)、Hartree-Fock/cc-pVDZ(300基底)– 改良前は、8コアではスレッド待ち時間とオーバーヘッド(オレンジ色)はほとん
どないが、28コアでは全計算時間の約1/4まで増加– 改良後の待ち時間とオーバーヘッドは、数%程度まで減少

改良前
8コア

改良前
28コア

改良後
28コア
38

OpenMP並列効率のさらなる向上3•••

第1世代Xeon Phi(Knights Corner, 5110P, 60コア, 1011GFLOPS)
Taxol(C47H51NO14),cc-pVDZのHartree-Fock計算時間(sec)Taxol分子では、60スレッド以上で並列加速率と計算時間に差が出ている
スレッド数

15
30
60
120
240

改良前

改良後

計算時間 並列加速率

計算時間 並列加速率9034.64559.02391.92652.92888.6

15.0
29.7
56.7
51.1
46.9

9057.94565.82302.6
1592.91333.5

15.0
29.8
59.0
85.3
101.9

参考：Xeon (E5-2698 v3, 32コア, 972GFLOPS)
32スレッドで345.9sec

39

メモリ使用量削減1••

演算量を増やしてでも、使用メモリ量を削減
MP2エネルギー計算では、AO2電子積分を複数回計算することで削減可
𝑜𝑐𝑐 𝑣𝑖𝑟

𝐸𝑀𝑃2 = ෍ ෍𝑖𝑗𝑎𝑏𝑎𝑖|𝑏𝑗 2 𝑎𝑖|𝑏𝑗 − 𝑎𝑗|𝑏𝑖
𝜀𝑖 + 𝜀𝑗 − 𝜀𝑎 − 𝜀𝑏𝐴𝑂

𝑎𝑖|𝑏𝑗 = ෍ 𝐶𝜇𝑎 𝐶𝜈𝑖 𝐶𝜆𝑏 𝐶𝜎𝑗 (𝜇𝜈|𝜆𝜎)𝜀𝑖 : 軌道エネルギー𝐶𝜇𝑎 : 分子軌道係数𝑖, 𝑗: 占有軌道𝑎, 𝑏: 非占有軌道𝜇𝜈𝜆𝜎

do ml (AO index pair)
do s
AO積分計算(mn|ls) (all n)
第1変換(mi|ls)(all i)enddo
第2変換(dgemm) (mi|lj)(i≧j)
end do ml
do ij (MO index pair)
MPI_sendrecv (mi|lj)
第3変換(dgemm) (mi|bj)(all b)第4変換(dgemm) (ai|bj)(all a)MP2エネルギー計算
end do ij
call mpi_reduce(MP2エネルギー)

do ij-blockdo ml (AO index pair)
do s
AO積分計算(mn|ls) (all n)
第1変換(mi|ls)(partial i)
enddo
第2変換(dgemm) (mi|lj)(partial j)
end do ml
do partial-ij (MO index pair)
MPI_sendrecv (mi|lj)
第3変換(dgemm) (mi|bj)(all b)第4変換(dgemm) (ai|bj)(all a)MP2エネルギー計算
end do ij
end do ij-block
call mpi_reduce(MP2エネルギー)
40

メモリ使用量削減2•

全プロセス合計の使用メモリ量はO2N2/2 → O2N2/(2nij-block) (O:占有軌道数、N:基底関数次元数, nij-block:占有軌道分割数, (mi|lj)保存用配列)– 計算条件：Taxol (C47H51NO14)、MP2/6-31G(d) (O=164軌道、 N=970基底)– 計算機：Fujitsu PRIMEGY CX2500 (Xeon E5-2697, 28コア/ノード)– MP2エネルギー計算では、メモリ使用量削減のための追加演算コストは比較的小さい計算時間(秒)とメモリ使用量(GB)
占有軌道分割数

1

2

3

MP2計算時間

765.7

943.4

1121.5

618.8

803.6

981.2

101

51

34

そのうちAO2電子
積分と第1変換メモリ使用量41

メモリ使用量削減3•MP2エネルギー微分計算の並列加速率– 計算条件：C150H30、MP2/cc-pVDZ– 計算機：「京」– 使用ノード数が増えるほどトータルのメモリ量は増えて、占有軌道分割数が減るため、使用ノード数以上の並列加速率が得られた
計算時間(秒)と並列加速率
ノード数

192

384

768

1536

3

2

1

1

MP2微分実行時間

4253.6

2002.0

743.6

382.0

並列加速率

192.0

407.9

1098.22137.9占有軌道分割数

42通信時間削減1••

演算と通信のオーバーラップにより、通信時間を削減
ノンブロッキング通信MPI_Isend,MPI_Irecvを使用
𝑜𝑐𝑐 𝑣𝑖𝑟

𝐸𝑀𝑃2 = ෍ ෍𝑖𝑗𝑎𝑏𝑎𝑖|𝑏𝑗 2 𝑎𝑖|𝑏𝑗 − 𝑎𝑗|𝑏𝑖
𝜀𝑖 + 𝜀𝑗 − 𝜀𝑎 − 𝜀𝑏𝐴𝑂

𝑎𝑖|𝑏𝑗 = ෍ 𝐶𝜇𝑎 𝐶𝜈𝑖 𝐶𝜆𝑏 𝐶𝜎𝑗 (𝜇𝜈|𝜆𝜎)𝜀𝑖 : 軌道エネルギー𝐶𝜇𝑎 : 分子軌道係数𝑖, 𝑗: 占有軌道𝑎, 𝑏: 非占有軌道𝜇𝜈𝜆𝜎

do ml (AO index pair)
do s
AO積分計算(mn|ls) (all n)
第1変換(mi|ls)(all i)enddo
第2変換(dgemm) (mi|lj)(i≧j)
end do ml
do ij (MO index pair)
MPI_sendrecv (mi|lj)
第3変換(dgemm) (mi|bj)(all b)第4変換(dgemm) (ai|bj)(all a)MP2エネルギー計算
end do ij
call mpi_reduce(MP2エネルギー)

do ml (AO index pair)
do s
AO積分計算(mn|ls) (all n)
第1変換(mi|ls)(all i)enddo
第2変換(dgemm) (mi|lj)(i≧j)
end do ml
do ij (MO index pair)
MPI_isend, irecv (次の(mi|lj))
第3変換(dgemm) (mi|bj)(all b)第4変換(dgemm) (ai|bj)(all a)MP2エネルギー計算
MPI_Waitall (次の(mi|lj))
end do ij
call mpi_reduce(MP2エネルギー)

43

通信時間削減2•「京」96, 192ノードでC150H30 (cc-pVDZ, 2250基底)のMP2エネルギー計
算時間– 演算と通信をオーバーラップさせることで、通信時間を大幅に削減– 前述の占有軌道分割数は96ノードでは2、192ノードでは1MP2実行時間と通信時間(秒)ノード数
改良前 MP2実行時間(うち通信時間)
改良後 MP2実行時間(うち通信時間)

961089.99.7
1082.2
5.2

192
421.9
8.7
416.7
3.4

44

まとめ
•

•
•
•

•

並列アルゴリズム開発では、計算負荷分散、データ分散、通信コスト削
減、高速化を同時に考慮して行う必要があり、計算式から考え直す場
合もある
MPI/OpenMPハイブリッド並列は、現在の計算機では必要不可欠
並列化効率向上のためには、初期値計算も含めたすべての計算を並列化する必要がある
アルゴリズム・プログラム開発コストはますます上昇する可能性が高く、分野全体でのコスト削減と情報共有が重要になり、そのための手段としてオープンソースでの公開が挙げられる
これからのメニーコア時代の計算機を効率的に使うためには、特にOpenMP並列効率、通信時間の削減、メモリ使用量の削減が重要になると考えられる

45

参考資料
•

•
•

量子化学
– “分子軌道法” 藤永茂– “新しい量子化学 上・下” A.ザボ, N.S.オスランド著, 大野公男,阪井健
男,望月祐志訳– “Introduction to Computational Chemistry” Frank Jensen– “Molecular Electronic-structure Theory” Trygve Helgaker, Poul
Jørgensen, Jeppe Olsen
計算機– “プロセッサを支える技術” Hisa Ando
BLAS、LAPACK– http://www.intel.co.jp/jp/software/products/

46