計算科学技術特論B 第14回

計算科学技術特論B (2016)
第14回

1 /74

大規模MD並列化の技術1
名古屋大学 工学研究科
附属計算科学連携教育研究センター
特任講師
安藤 嘉倫

2016/7/21

講義の主旨計算科学技術特論B 第14回
2 /74

・分子動力学計算を例に, 三次元トーラスネットワーク上での高並列対応ソフトウェアのコーディング技術を紹介→ 京コンピュータ, FX10 などで有効な技術→ 他ソフトにも有効な技術. 講義内容にインスパイアされた並列化
ードを自作, 公表した場合は文献[1]を引用下さい.

コ・MODYLAS[1] での実装例をもとに, 並列 (MPI, OpenMP,
SIMD) の効率化技術をメインに, 演算の効率化技術についても紹介
・コードの詳細には極力立ち入らずコーディング概念を説明→ 詳細はライセンスに同意の上 www.modylas.org での公開ソースコー
ドをダウンロードして確認下さい.
[1] Y.Andoh, N.Yoshii, K. Fujimoto, K.Mizutani, H.Kojima, A.Yamada, S.Okazaki, K.Kawaguchi, H.Nagao, K.Iwahashi, F. Mizutani,
K.Minami, S. Ichikawa, H.Komatsu, S. Ishizuki, Y.Takeda, and M.Fukushima, J. Chem. Theory Comput., 9, 3201 (2013).

計算科学技術特論B 第14回

目次

3 /74

・分子動力学 (MD) 法・ 分子動力学計算の並列化特性
・並列化技術 1

データ構造・並列化技術 2

MPI・並列化技術 3

OpenMP, SIMD

第一回

第二回

計算科学技術特論B 第14回

分子動力学法(1) 基礎方程式
ニュートンの運動方程式
2

4 /74

rN

d ri
"
mi 2 = Fi = ! ! (r1, r2 ,..., ri ,..., rN )dt
"rii = 1,…, N

時間に対する座標の二階常微分方程式

N = 6	
t : 時刻	ri : 原子 i の座標	mi : 原子 i の質量	r1

Fi : 原子 i に作用する力	! (r N ) : 系のポテンシャルエネルギー	
N : 系に含まれる原子数

r3

F1
F2
r2

F5
r5

F4

F3

r4

r6

F6

計算科学技術特論B 第14回

分子動力学法(2) 運動方程式の数値積分
dr (t)vi (t) = idt

ベルレ法による差分化＋

Fi (t) d 2 ri (t)
=
mi
dt 2

dri (t)1 d 2 ri (t) 2
1 Fi (t) 2
ri (t + !t) = ri (t) +
!t +
!
t
+!
=
r
(t)
+
v
(t)
!
t
+
!t +!
i
i
dt
2! dt 22! m
dri (t)1 d 2 ri (t) 2
1 Fi (t) 2
ri (t ! !t) = ri (t) !
!t +
!
t
!!
=
r
(t)
!
v
(t)
!
t
+
!t !!
i
i
2
dt
2! dt
2! m

ri (t + !t) = 2ri (t) + ri (t ! !t) +ベルレの差分式

Fi (t) 2!t + O(!t 4 )m+Δtri (t + !t)t=t

Taylor 展開 -Δt

現在は「速度ベルレ法」
による差分化が標準的+Δt

ri (t)

Taylor 展開 +Δt

Fi (t + !t)

Fi (t)+Δt

5 /74

t = t + Δt

計算科学技術特論B 第14回

分子動力学法(3) 座標に対する境界条件

6 /74・周期境界条件

・自由境界条件

イメージセル

イメージセル

イメージセル

イメージセル

イメージセル基本セル

真空中での孤立系 (NVE)

イメージセル

イメージセル

イメージセルバルクでの孤立系 (NVE)計算科学技術特論B 第14回

分子動力学法(4)原子間相互作用7 /74

古典分子動力学法
古典近似された原子間相互作用関数 (力場) ! (r N )各項は以降のスライド参照! (r N ) =

2
2
K
(b(r
,
r
)
!
b
)
+
K
(
"
(r
,
r
,
r
)
!
"
)
" b i j 0 " " i j k 0
bonds

angles+ " K# #$1+ cos(n# (ri , rj , rk , rl ) ! $ )%& +dihedrals" K (% (r , r , r , r ) ! % )%impropers

12
61
# -'
%
*
'
*
/
/
q
q
'
'
,, ! )),, 2 + i j 5+ " 44& .))4 /( r(ri , rj ) + ( r(ri , rj ) + / r(ri , rj ) 5
nonbonds
3$ 0&

Lennard-Jonesd 2 ri"mi 2 = Fi = ! # (r N )dt"riCoulomb

i

j

k

l

0

2

計算科学技術特論B 第14回

分子動力学法(4)原子間相互作用8 /74

古典分子動力学法
古典近似された原子間相互作用関数 (力場) ! (r N )各項は以降のスライド参照! (r N ) =

2
2
K
(b(r
,
r
)
!
b
)
+
K
(
"
(r
,
r
,
r
)
!
"
)
" b i j 0 " " i j k 0
bonds

angles+ " K# #$1+ cos(n# (ri , rj , rk , rl ) ! $ )%& +dihedrals" K (% (r , r , r , r ) ! % )%iimpropers

12
61
# -'
%
*
'
*
/
/
q
q
'
'
,, ! )),, 2 + i j 5+ " 44& .))4 /( r(ri , rj ) + ( r(ri , rj ) + / r(ri , rj ) 5
nonbonds
3$ 0&

Lennard-Jones

Coulomb

分子内 (intra-molecule, bonded) 相互作用 +分子間 (inter-molecule, nonbonded) 相互作用	

j

k

l

0

2

計算科学技術特論B 第14回

分子動力学法(5)分子内 bonded 相互作用

K b (b(ri , rj ) ! b0 )2

improper dihedral	

bond	

i	
i	

j	
angle	

dihedral	

l	

i	
j	

k	

j	

9 /74

k	

i	

k	
l	

K! "#1+ cos(n! (ri , rj , rk , rl ) ! " )$%j	
K! (! (ri , rj , rk , rl ) ! ! 0 )2

K! (! (ri , rj , rk ) ! ! 0 )2
コーディングイメージ

do n=1,nbonds
do n=1,nangles
b=b(ri, rj)
theta=thteta(ri, rj, rk)
φbond,ij の計算
φangle,ijk の計算
φbond=φbond+φbond,ij
φangle=φangle+φangle,ijk
Fi, Fj の計算
Fi, Fj, Fk の計算
f(i)=f(i)+Fi
f(i)=f(i)+Fi
f(j)=f(j)+Fj
f(j)=f(j)+Fj
enddo 	f(k)=f(k)+Fk
enddo 	

do n=1,ndihedrals
do n=1,nimpropers
phi=phi(ri, rj, rk, rl)
psi=psi(ri, rj, rk, rl)
φdihedral,ijkl の計算
φimproper,ijkl の計算
φdihedral=φdihedral+φdihedral,ijkl
φimproper=φimproper+φimproper,ijkl
Fi, Fj, Fk, Fl の計算
Fi, Fj, Fk, Fl の計算
f(i)=f(i)+Fi
f(i)=f(i)+Fi
f(j)=f(j)+Fj
f(j)=f(j)+Fj
f(k)=f(k)+Fk
f(k)=f(k)+Fk
f(l)=f(l)+Fl
f(l)=f(l)+Fl
enddo 	enddo 	

計算科学技術特論B 第14回

分子動力学法(6)分子内・分子間 nonbonded 相互作用10 /74

nonbonded = Lennard-Jones(LJ) + Coulomb(CL)	
12
6
qi q j
(!
$ !
$ ,*
*
"
"
&& ' ##&& 4! )##*+" r(ri , rj ) % " r(ri , rj ) % *.

r(ri , rj )r(ri , rj ) =(xi ! x j )2 + (yi ! y j )2 + (zi ! z j )2

分子内

1-4 scale	

分子間1-2 void	

コーディングイメージ

分子間do imol=1,nmol-1
do jmol=imol+1,nmol
do i=1,natom(imol)
do j=1,natom(jmol)
rij=rij(ri, rj)
φij の計算
φnonbond=φnonbond+φij
Fi, Fj の計算
f(i)=f(i)+Fi
f(j)=f(j)+Fj
enddo
enddo
enddo
enddo	

分子内 nonbond 項計算の注意点:

分子内do imol=1,nmol
do i=1,natom(imol)-1
do j=i+1,natom(imol)
rij=rij(ri, rj)
0 if 1-2,-3 voidx = s if 1-4 scale
1 elseφij の計算
φnonbond=φnonbond+x*φij
Fi, Fj の計算
f(i)=f(i)+x*Fi
f(j)=f(j)+x*Fj
enddo
enddo
enddo	

・1番目および2番目の隣接原子とは相互作用しない (1-2, 1-3 void)・3番目の隣接原子との相互作用は因子 s でスケールする (1-4 scale) [s は LJ, Coulomb べつ]・4番目以降の隣接原子とは通常の相互作用計算科学技術特論B 第14回

分子動力学法(7)! LJ (r) =126,(!
* # $ ! # $*&& ' ##&& 4" )##*+" r(ri , rj ) % " r(ri , rj ) % *.

LJ 相互作用のカットオフ

11 /74

if (r ≤ rcut) then	 else	
12

!LJ = !LJ (r)

!! $
# &
"r%

!LJ = 0カットオフ

r(ri , rj ) =12!! $ !! $# & '# &"r% "r%(xi ! x j )2 + (yi ! y j )2 + (zi ! z j )2

r > rcutを無視すること
で計算時間を削減i	

6

値・傾きともに
ほぼ 0

j	

rcut

"! %
!$ '
#r&

6

カットオフ半径rcut = 3~4σ	r	計算科学技術特論B 第14回

分子動力学法(7)・分子系での留意事項

コーディングイメージ(分子間)do imol=1,nmol-1
do jmol=imol+1,nmol
do i=1,natom(imol)
do j=1,natom(jmol)
rij=rij(ri, rj)if(rij > rcut) cycle
φij の計算
φnonbond=φnonbond+φij
Fi, Fj の計算
f(i)=f(i)+Fi
f(j)=f(j)+Fj
enddo
enddo
enddo
enddo	

LJ 相互作用のカットオフ

12 /74

LJ相互作用を計算しない
別の分子内の原子

LJ相互作用を計算しない同一分子内
の原子

j	
void	

rij	

scale	

i	

void	
void	

カットオフ半径 rcut	計算科学技術特論B 第14回

分子動力学法(8)!CL (r) =Coulomb 相互作用

13 /74qi q jr(ri , rj )

カットオフ
すると結果
に影響大

12

!! $
# &
"r%
r(ri , rj ) =値・傾きともに0 ではない(xi ! x j )2 + (yi ! y j )2 + (zi ! z j )2

i	

j	

1
r
"! %
!$ '
#r&

6

カットオフ半径rcut	

r	計算科学技術特論B 第14回

分子動力学法(8)Coulomb 相互作用

14 /74

カットオフしない Coulomb 相互作用の計算法
・自由境界条件下
すべての i, j 原子対を計算
多重極展開法, および高速多重極展開法 (FMM)・周期境界条件下
Ewald 法Particle Mesh Ewald (PME) 法
FMM + 多重極子のEwald法

今回の講義ではこれの並列化を説明FMM:
方法論の詳細は「CMSI計算科学技術特論A」 古典分子動力学法の高速化 (吉井) を参照

計算科学技術特論B 第14回

分子動力学法(8)Coulomb 相互作用
基本セル

15 /74l = L/2n

FMM での階層的セル分割
・基本セルは立方体 (一辺 L)・各辺を 2n 回に等分割
本講義では簡単のためこの2つを前提

L
用語の定義：
被相互作用点の含まれるサブセル

サブセル:
レベル0の分割セル (一辺 l = L/2n)
スーパーセル:
レベルmの分割セル (一辺 l×2m)
系には計8n 個のサブセル, レベルmには8n-m個のスーパーセル.

: レベル０計算科学技術特論B 第14回

分子動力学法(8)Coulomb 相互作用

16 /74l = L/2n

基本セル

FMM での階層的セル分割
・基本セルは立方体 (一辺 L)・各辺を 2n 回に等分割
本講義では簡単のためこの2つを前提

L
用語の定義：
被相互作用点の含まれるサブセル

サブセル:
レベル0の分割セル (一辺 l = L/2n)
スーパーセル:
レベルmの分割セル (一辺 l×2m)
系には計8n 個のサブセル, レベルmには8n-m個のスーパーセル.
サブセルはm=0のスーパーセル: レベル０: レベル１: レベル３: レベル４: レベル２

計算科学技術特論B 第14回

分子動力学法(8)FMM 基本演算

Coulomb 相互作用

17 /74多極子エワルド法
（無限遠からの寄与)

M2L (multipole to local)
・同じ階層
・well-separated条件を満たす相手セル

L2L

M2M(local to local)(multipole to multipole)

展開中心

展開中心
点電荷

多極子展開

p2M (particle to multipole)
p2M 	
  :	
  点電荷の多極子展開	
  	
  
M2M	
  :	
  多極子のマージ	
  
M2L 	
  :	
  多極子の局所展開 	
  	
  

点電荷

p2p
(particle to particle)

局所展開

L2p (local to particle)
L2L 	
  :	
  局所展開中心の移動	
  	
  
L2p	
  :	
  点電荷上の電場の評価	
  
p2p	
  :	
  粒子対計算 Hot spots

計算科学技術特論B 第14回

分子動力学法(8)Coulomb 相互作用
イメージセル

相互作用計算の要点
レベル０
・近傍サブセルに含まれる点電荷との p2p・サブセルAとの M2L・点電荷への局所展開 L2p
レベル1
・スーパーセルBとのM2L, L2L...以下同様 (入れ子構造)
最上レベル
・多重極子のEwald法, L2L
入れ子構造により,周期境界条件下では各レベル
103-53 = 875 個のスーパーセルとM2L

18 /74

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

A A A A A A A A A A
A A A A A A A A A A

B

B

B

B

B

A A A AA A

A A A A A A
A A A

B

B

B

B

B

A A
A A

A A A
A A A

B

B

B

B

B

A A
A A

A A A
A A A

B

B

B

B

B

A A A A A A A A A A
A A A A A A A A A A

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

B

基本セル
B
B

B

B

B

B

B

: レベル０: レベル１: レベル２: レベル３: レベル４: 粒子-粒子相互作用を計算計算科学技術特論B 第14回

分子動力学計算の並列化特性(1)19 /74

MD 計算の処理フロー
座標・速度, 力場パラメータの入力
座標よりポテンシャル・力・ビリアルの計算
cycle >106

分子内力
分子間力（LJ, CL粒子間）

Hot spot

分子間力（p2M, M2M, M2L, L2L, L2p） Hot spot
総演算量の9割以上

座標・速度の更新 (Verlet 法)
系の力学量・熱力学量 (H, Φ, K, T, P, ...)の計算, 出力
リスタート用 End
座標・速度の出力

計算科学技術特論B 第14回

分子動力学計算の並列化特性(1)20 /74

MD 計算の処理フロー
座標・速度, 力場パラメータの入力
座標よりポテンシャル・力・ビリアルの計算
cycle >106

分子内力
分子間力（LJ, CL粒子間）
分子間力（p2M, M2M, M2L, L2L,
MPIL2p）
並列化
座標・速度の更新 (Verlet 法)
系の力学量・熱力学量 (H, Φ, K, T, P, ...)の計算, 出力
リスタート用 End
座標・速度の出力

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 21 /74

d 2 rimi 2 = Fidt
i ごと独立な演算
(粒度N)
rcut
コーディングイメージ

do istep=1,Nstep
do i=1, N
ポテンシャルの計算
力の計算
enddo
座標 • 速度更新
力学 • 熱力学量計算
enddo

基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 22 /74

d 2 rimi 2 = Fidt
i ごと独立な演算(粒度N)
最も単純な MPI 並列化:
do istep=1,Nstep
do i=1+myrank, N, nprocs
ポテンシャルの計算
力の計算
enddo
call mpi_allreduce(ポテンシャル, 力)
座標 • 速度更新
力学 • 熱力学量計算
enddo
myrank : MPIプロセス番号 (0始まり)
nprocs : MPIプロセス数

rcut
myrank基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 23 /74

d 2 rimi 2 = Fidt
サブセル

i ごと独立な演算(粒度N)
空間ドメイン分割によるMPI 並列化:

rcut①基本セルを各辺にそって分割

基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 24 /74
2d ri
mi 2 = Fi
dt例) プロセスあたり 4 サブセル

myranki ごと独立な演算(粒度N)
空間ドメイン分割によるMPI 並列化:

rcut①基本セルを各辺にそって分割②複数のサブセルを各プロセスに均等割り当て基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 25 /74
2d ri
mi 2 = Fi
dt例) プロセスあたり 4 サブセル

myranki ごと独立な演算(粒度N)
空間ドメイン分割によるMPI 並列化:

rcut①基本セルを各辺にそって分割②複数のサブセルを各プロセスに均等割り当て③ myrank の所持するサブセル内
原子と近傍 other ranks の所持するサブセル内原子との間で粒子対型の相互作用計算

other ranks
基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 26 /74
2d ri
mi 2 = Fi
dt例) プロセスあたり 4 サブセル

myranki ごと独立な演算
(粒度N)

jcell

空間ドメイン分割によるMPI 並列化:

rcut

icell

コーディングイメージ

do icell(myrank)
do jcell_list(myrank or otherranks)
ポテンシャルの計算
力の計算
enddo
座標 • 速度更新
力学 • 熱力学量計算
enddo

other ranks
基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 27 /74
2d ri
mi 2 = Fi
dt例) プロセスあたり 4 サブセル

myranki ごと独立な演算
(粒度N)

jcell

空間ドメイン分割によるMPI 並列化:

rcut

icell

コーディングイメージ

do icell(myrank)
do jcell_list(myrank or otherranks)
ポテンシャルの計算
力の計算
enddo
座標 • 速度更新
力学 • 熱力学量計算
enddo

作用反作用 Fij = -Fji
により除くことのできる
jcell 範囲

other ranks
基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(2) プロセス分割 28 /74
計算機 (ノード)

通信ケーブルmyrank

トーラスネットワーク

空間ドメイン分割によるMPI 並列化:
ノード形状とプロセス形状とを一致させることを前提に,
トーラスネットワークとの相性が非常に良いother ranks
基本セル
イメージセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 29 /74
座標・速度, 力場パラメータの入力
座標よりポテンシャル・力・ビリアルの計算
cycle >106

通信①(隣接)

分子内力
分子間力（LJ, CL粒子間）

通信②(隣接・近接)
分子間力（p2M, M2M, M2L, L2L, L2p）
通信③(隣接)座標・速度の更新 (Verlet 法)

系の力学量・熱力学量 (H, Φ, K, T, P, ...)の計算, 出力
通信④(全ノード間・総和型)

リスタート用 End
座標・速度の出力

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 30 /74
例) プロセスあたり 4 サブセル

通信①(隣接)サブセル

l

分子内力, 分子間力の計算に必要な相手原子座標 rj, rk, rl の通信.特徴:・近傍プロセスとのみ通信
・データサイズは小さい例) 40 原子/サブセル では, 
8 byte(real*8)*3(xyz)*40 
= 1 KB/サブセル

rcut


注意点:
・rcut ≤ n*l [今回の講義では n=2]・分子内結合原子 (特に dihedral)
を通信範囲内に収める

計算対象原子を含むセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 31 /74
通信①(隣接)l

データ通信が必要な
other ranks
63-23 = 208サ
ブセル
データサイズ~200 KB

分子内力, 分子間力の計算に必要な相手原子座標 rj, rk, rl の通信.特徴:・近傍プロセスとのみ通信
・データサイズは小さい例) 40 原子/サブセル では, 
8 byte(real*8)*3(xyz)*40 
= 1 KB/サブセル

rcut


注意点:
・rcut ≤ n*l [今回の講義では n=2]・分子内結合原子 (特に dihedral)
を通信範囲内に収める

計算対象原子を含むセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 32 /74
通信①(隣接)l

分子内力, 分子間力の計算に必要な相手原子座標 rj, rk, rl の通信.

データ通信が必要な
other ranks

周期境界条件による折り返し
特徴:・近傍プロセスとのみ通信
・データサイズは小さい例) 40 原子/サブセル では, 
8 byte(real*8)*3(xyz)*40 
= 1 KB/サブセル

rcut


注意点:
・rcut ≤ n*l [今回の講義では n=2]・分子内結合原子 (特に dihedral)
を通信範囲内に収める・周期境界条件による折り返し

イメージセル

基本セル

計算対象原子を含むセル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 33 /74
通信②(隣接・近接)レベル 0

多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrank計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 34 /74
通信②(隣接・近接)レベル 0

多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

M2M 対象セル

myrank

この場合
データ通信は必要なし

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 35 /74
通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	

レベル 0
プロセス内自セル位置(0,0), (1,0), (0,1), (1,1)


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrankM2L 対象セル

注意

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 36 /74
通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	

レベル 0
プロセス内自セル位置(0,0), (1,0), (0,1), (1,1)


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrank

M2L 対象セル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 37 /74
通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	

レベル 0
プロセス内自セル位置(0,0), (1,0), (0,1), (1,1)


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrank

M2L 対象セル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 38 /74
通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	

レベル 0
プロセス内自セル位置(0,0), (1,0), (0,1), (1,1)


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrank



	

	

M2L 対象セル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 39 /74
通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	

レベル 0
プロセス内自セル位置(0,0), (1,0), (0,1), (1,1)


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

myrank


	

	

	

データ通信が必要な
other ranks

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 40 /74
通信②(隣接・近接)レベル 1

多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

	

	

	

	

M2M 対象セル
データ通信が必要なother ranks
myrank計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 41 /74
レベル 1

通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

イメージセル基本セル

myrank


	

	

	

	

M2L 対象セル

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 42 /74
レベル 1

通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB
	

イメージセル基本セル

myrankデータ通信が必要な
other ranks

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 43 /74
レベル 1

通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KBy	

注意点 (上位階層):・軸ごとrank参照範囲が異なる
	

右の例では, x方向 -4 ~ +5	

y方向 -4 ~ +5	

イメージセル

+5

基本セル

+5

-4
myrank-4
データ通信が必要な
other ranksx計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 44 /74
レベル 1

通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KBy	

注意点 (上位階層):・軸ごとrank参照範囲が異なる
	

右の例では, x方向 -5 ~ +4	

y方向 -4 ~ +5	

イメージセル

+5

基本セル

+4

-5
myrank-4
データ通信が必要な
other ranksx計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 45 /74
レベル 1

通信②(隣接・近接)多極子情報の通信	

M2M	

: 周辺 8 スーパーセル 
M2L	

: 周辺 875 スーパーセル	


特徴:
・通信対象範囲が階層ごと異なる. 
上位階層ほど遠方.
・データサイズ例) nmax=4 では
16 byte(complex*8)*(4+1)2 = 0.4 KB/スーパーセル
階層ごと875*0.4KB=320KB

イメージセル基本セル

myrank
注意点:・軸ごとrank参照範囲が異なる
・周期境界条件による折り返し

データ通信が必要な
other ranks
周期境界条件による折り返し計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 46 /74
通信③(隣接)運動方程式の数値積分の結果, 原子が移動. これを各プロセスに再帰属させるための座標通信.

特徴:
・最近接 33-1=26 プロセスと通信
・データサイズは非常に小さい・myrank から出る粒子, myrank に
入ってくる粒子, 両方ある

myrankyx

データ通信が必要な
最近接 other ranks

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 47 /74
通信③(隣接)運動方程式の数値積分の結果, 原子が移動. これを各プロセスに再帰属させるための座標通信.

(Δx, Δy, Δz)= (-1, +1, 0)特徴:
・最近接 33-1=26 プロセスと通信
・データサイズは非常に小さい・myrank から出る粒子, myrank に
入ってくる粒子, 両方ある

myrank(Δx, Δy, Δz)= (+1, 0, 0)





注意点:
・x, y, z 軸方向への移動に加え,
斜め方向への移動が発生(Δx, Δy, Δz)= (0, -1, 0)

yx

データ通信が必要な
最近接 other ranks

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 48 /74
通信③(隣接)運動方程式の数値積分の結果, 原子が移動. これを各プロセスに再帰属させるための座標通信.

(Δx, Δy, Δz)= (-1, +1, 0)特徴:
・最近接 33-1=26 プロセスと通信
・データサイズは非常に小さい・myrank から出る粒子, myrank に
入ってくる粒子, 両方ある

イメージセル
周期境界条件による折り返し

myrank

基本セル(Δx,
Δy, Δz)= (+1, 0, 0)
注意点:
・x, y, z 軸方向への移動に加え,
斜め方向への移動が発生
・周期境界条件による折り返し(Δx, Δy, Δz)= (0, -1, 0)

yx

データ通信が必要な
最近接 other ranks

計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 49 /74
通信④(全ノード間・総和型)例) 運動エネルギー
系全体での量（ハミルトニアン, 温度, 圧力など）を計算するための通信
特徴:・データサイズは, real*8 の1変数	

・全プロセス間で総和をとる (allreduce)→ハードウェア機能の利用で高速化	

例えば, 京の TBI (Tofu barrier interface)~10,000ノード reduce 演算 10 µs
安島ら, FUJITSU. 63, 3, p. 260-264(05, 2012)

nprocs N

N

k
1
1
2
K = ! mivi = ! ! mivi2
i=1 2
k=1 i=1 2

kene

プロセスごとの和
wk_kenecall mpi_allreduce(wk_kene,kene,1,
mpi_double_precision,mpi_sum,
mpi_comm_world,ierr)計算科学技術特論B 第14回

分子動力学計算の並列化特性(3) 通信パターン 50 /74
通信④(全ノード間・総和型)例) 運動エネルギー
系全体での量（ハミルトニアン, 温度, 圧力など）を計算するための通信
特徴:・データサイズは, real*8 の1変数	

・全プロセス間で総和をとる (allreduce)→ハードウェア機能の利用で高速化			

例えば, 京の TBI (Tofu barrier interface)	nprocs N

N

k
1
1
2
K = ! mivi = ! ! mivi2
i=1 2
k=1 i=1 2

kenecall mpi_allreduce(wk_kene,kene,1,
mpi_double_precision,mpi_sum,
mpi_comm_world,ierr)~10,000ノード reduce 演算 10 µs

二重カウント安島ら, FUJITSU. 63, 3, p. 260-264(05, 2012)	

注意点:・rank 間での二重カウント
特にポテンシャルエネルギーとビリアル		

プロセスごとの和
wk_kenemyrankother rank

計算科学技術特論B
CMSI配信講義B 第14回
第10回

目次

/74

・分子動力学 (MD) 法・ 分子動力学計算の並列化特性
・並列化技術 1

データ構造・並列化技術 2

MPI・並列化技術 3

OpenMP, SIMD

第一回

第二回

CMSI配信講義B 第10回

MPI 並列化技術: データ構造

データ構造: MPI 並列性能 (および演算性能) に直結
MD 計算のように 1 サイクル当たりの経過時間が小さい [≤ O(10-3 sec)] 場合に特に重要.
システムが大規模・複雑化する将来, アプリによらず重要になる (?).

CPU

送信データ
の荷造り

これらを最小限に止めるような
データ構造

コピー

メモリー

ノードA

コピー

メモリー

データ転送
send通信用
バッファ

CPU

受信データ
の荷解き

recv通信用バッファ

ノードB[理想] プログラムを書く前に, ネットワーク構造に最適なデータ構造を考案する[現実] 途中で問題に直面して, 既存のプログラムをイチから書き直す

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
53 /74

プロセス分割した配列

通常のデータ構造

allocate( wk_x(n) )注) 実際には, 3次元のallocate( wk_xyz(3,n) )オリジナル配列
allocate( x(N) ) N: 原子数≈≈n = O(N/nprocs)

MPI プロセス
分割・・・
nprocs 個

頭詰め

頭詰め

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
54 /74転送された座標データの格納形態
wk_x(n)

荷解き&コピー

other ranksからの受信データ

荷解き&コピー

荷解き&コピー

荷解き&コピー

頭詰め

myrankの所持データ wk_x(1)MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
55 /74

新しく考案した座標のデータ構造: メタデータ構造
meta_x(n’)
n’ ~ n*margin

メタデータ：あるデータに関連する情報allocate( meta_x(n’) )+ここでは,あるデータ: 原子座標
関連する情報 : 原子の所属するサブセルの相対的な空間配置

myrankの所持データ

メタデータ管理配列

allocate( tag(nsbc, nsbc),
na_per_cell(nsbc, nsbc) )nsbc: 軸あたりのサブセル数
注) 3次元では,
allocate( tag(nsbc, nsbc, nsbc) )meta_x(1)MPI 並列化技術: データ構造 [1] 座標

56 /74

メタデータ
X-Y 平面上でのサブセル相対配置 = メタデータ(X,Y): セルの番地

meta_x(n)Ymeta_x(n’)

計算科学技術特論B 第14回

myrankの
所持データ

(1,5)

(2,5)

(3,5)

(4,5)

(5,5)

(1,4)

(2,4)

(3,4)

(4,4)

(5,4)

(1,3)

(2,3)

(3,3)

(4,3)

(5,3)

(1,2)

(2,2)

(3,2)

(4,2)

(5,2)

(1,1)

(2,1)

(3,1)

(4,1)

(5,1)

meta_x(1)meta_x(1)X

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
57 /74

メタデータ構造

メタデータ管理配列

サブセルの格納原子数:
allocate( na_per_cell(1:5,1:5) )データへのアクセス方法

(1,5)

(2,5)

(3,5)

(4,5)

(5,5)

(1,4)

(2,4)tag(3,3)
(3,4)

(4,4)

(5,4)

(1,3)

(2,3)

(4,3)

(5,3)

(3,3)

na_per_cell(3,3)

myrank所持データ:
tag(3,3) ~tag(3,3)+na_per_cell(3,3)-1

meta_x(n)Y

サブセルの先頭原子番号:
allocate( tag(1:5,1:5) )(X,Y): セルの番地

(1,2)

(2,2)

(3,2)

(4,2)

(5,2)

(1,1)

(2,1)

(3,1)

(4,1)

(5,1)

meta_x(1)X

MPI 並列化技術: データ構造 [1] 座標
転送された座標データの格納形態

計算科学技術特論B 第14回
58 /74

メタデータ構造Y

meta_x(n)myrankの所持データ

meta_x(1)X

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
59 /74

メタデータ構造転送された座標データの格納形態

Y

meta_x(n)荷造り&コピー①x軸方向転送

myrankの所持データ

X軸方向, 他ランクからの転送データ

X軸方向, 他ランクからの転送データ

右詰め

左詰め

荷解き&コピー

荷解き&コピー

meta_x(1)X

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
60 /74

メタデータ構造転送された座標データの格納形態(X,Y): セルの番地
Y

meta_x(n)

荷造り
&コピー

②y軸方向転送

myrankの
所持データ

荷解き
&コピー

Y軸方向, 他ランクからの転送データ

Y軸方向, 他ランクからの転送データ

X軸方向, 他ランクからの転送データ

X軸方向, 他ランクからの転送データ

Y軸方向, 他ランクからの転送データ

Y軸方向, 他ランクからの転送データ

meta_x(1)X

MPI 並列化技術: データ構造 [1] 座標

計算科学技術特論B 第14回
61 /74

メタデータ構造データへのアクセス方法
myrank所持データ:

帯の部分のデータ:
tag(1,5) ~tag(5,5)+na_per_cell(5,5)-1
tag(1,4) ~tag(5,4)+na_per_cell(5,4)-1
tag(1,3) ~tag(5,3)+na_per_cell(5,3)-1

meta_x(n)Ytag(3,3) ~tag(3,3)+na_per_cell(3,3)-1

(1,5)

(2,5)

(3,5)

(4,5)

(5,5)

(1,4)

(2,4)tag(3,3)
(3,4)

(4,4)

(5,4)

(1,3)

(2,3)

(4,3)

(5,3)

(3,3)

na_per_cell(3,3)

tag(1,2) ~tag(5,2)+na_per_cell(5,2)-1(1,2)
tag(1,1)(2,2)(3,2)tag(5,1)
(4,2)

(5,2)

tag(1,1) ~tag(5,1)+na_per_cell(5,1)-1

(1,1)

(2,1)

(3,1)

(5,1)

(4,1)

meta_x(1)
連続アクセス例) 40原子x5サブセル=200原子 (5KB)→ L1キャッシュ(~32KB)に乗り, かつSIMD並列に十分な長さX
na_per_cell(5,1)

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子
下位階層の多極子

62 /74

下位の定義:「M2L範囲 < スーパーセル分割数」

10

allocate( wm0(:,16,16) )<オリジナル配列
myrank

プロセス
分割

レベル0

レベル0

16

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

63 /74

下位階層の多極子
allocate( wm0(:,16,16) )オリジナル配列
M2L 対象セル
myrankプロセス
分割

M2M 対象セル

レベル0でのM2M,
M2L演算に必要な
データ範囲

レベル0

レベル0

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子
下位階層の多極子

64 /74

配列を袖部付きで局所化

allocate( wm0(:,16,16) )allocate( wm_local0(:,10,10) )x方向: 10=4+2+4
y方向: 10=4+2+4

オリジナル配列
M2L 対象セル
myrank

+4

-4

プロセス
分割

+4

M2M 対象セル-4

レベル0でのM2M,
M2L演算に必要な
データ範囲

レベル0

レベル0

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子
下位階層の多極子

65 /74(X,Y): セルの番地
allocate( wm_local0(:,10,10) )Y

データへのアクセス方法(10,1)(10,10)myrank所持データ:
wm_local0(1:(nmax+1)2,5,5)
wm_local0(1:(nmax+1)2,6,5)
wm_local0(1:(nmax+1)2,5,6)
wm_local0(1:(nmax+1)2,6,6)(5,6) (6,6)

注）多極子にはもともとスーパーセルの相対配置情報が含まれているため, メ
タデータ管理配列は要らない(=厳密にはメタデータ構造ではない)(5,5) (6,5)
myrank所持データ(1,1)(10,1)

X

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子
下位階層の多極子

66 /74(X,Y): セルの番地
allocate( wm_local0(:,10,10) )Y

データへのアクセス方法
myrank所持データ:
wm_local0(1:(nmax+1)2,5,5)
wm_local0(1:(nmax+1)2,6,5)
wm_local0(1:(nmax+1)2,5,6)
wm_local0(1:(nmax+1)2,6,6)(10,1)

Y軸方向, 他ランクからの転送データ

Y軸方向, 他ランクからの転送データ

X軸方向, 他
ランクからの転送データ

x,y軸方向転送荷造り&コピー

other ranks所持データ:

荷解き
&コピー

(10,10)(5,6) (6,6)(5,5) (6,5)

X軸方向, 他
ランクからの
転送データ

Y軸方向, 他ランクからの転送データ

Y軸方向, 他ランクからの転送データ

wm_local0(1:(nmax+1)2,1:10,9:10)(1,1)
wm_local0(1:(nmax+1)2,1:10,7:8)
wm_local0(1:(nmax+1)2,1:10,5:6)
wm_local0(1:(nmax+1)2,1:10,3:4)
全範囲に連続アクセス
wm_local0(1:(nmax+1)2,1:10,1:2)(10,1)

X

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子
上位階層の多極子

67 /74

上位の定義:「M2L範囲 ≥ スーパーセル分割数」

10

allocate( wm1(:,8,8) )>オリジナル配列

プロセス
分割

レベル1

myrankレベル1

8

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

68 /74

上位階層の多極子
オリジナル配列
allocate( wm1(:,8,8) )M2L 対象セル

レベル1でのM2M,
M2L演算に必要なデータ範囲
myrankM2M 対象セル

基本セル
イメージセル

レベル1

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

69 /74

上位階層の多極子
オリジナル配列レベル1でのM2M,
M2L演算に必要なデータ範囲

allocate( wm1(:,8,8) )(10,10)+4

配列を袖部付きで局所化+4-5

allocate(wm_local1(:,10,10) )(6,6)

myrank所持データ:
wm_local1(:,6,6)

Y

x方向: 10=5+1+4
y方向: 10=5+1+4
M2M 対象セル

-5

レベル1
(1,1)

X

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

70 /74

上位階層の多極子
M2L 対象セル

レベル1でのM2M,
M2L演算に必要なデータ範囲

y 方向に -1 隣のプロセスを考えると,Y-1

myrankM2M 対象セル

レベル1
X

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

71 /74

上位階層の多極子
オリジナル配列レベル1でのM2M,
M2L演算に必要なデータ範囲

allocate( wm1(:,8,8) )+5

配列を袖部付きで局所化
allocate(wm_local1(:,10,10) )上位階層では, 各軸の +, - 袖部の値 (4 or 5) がプロセス位置に
依存する

+4

-5

x方向: 10=5+1+4
y方向: 10=4+1+5

(6,5)

Y

myrank所持データ:
wm_local1(:,6,5)(10,10)

M2M 対象セル

-4
レベル1
(1,1)

X

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

72 /74「M2L範囲 > スーパーセル分割数」

上位階層の多極子, 通信削減の工夫

10

allocate( wm1(:,4,4) )

>

4

≈

オリジナル配列+4

myrank

≈

≈

プロセス
分割

+4

-5
M2M 対象セル

レベル2

≈レベル2
-5

計算科学技術特論B 第14回

MPI 並列化技術: データ構造 [2] 多極子

73 /74「M2L範囲 > スーパーセル分割数」

上位階層の多極子, 通信削減の工夫

10

allocate( wm1(:,4,4) )>4≈
下階層のように袖部付き局所化
配列にすると、+4

・折り返し部分のデータを何度も冗長に通信→ 演算部であたかも袖部にあるかのようにアクセスすることで
通信を削減

≈

≈

もともとの配列 wm1(:,4,4) の全要素を通信すれば十分

myrank+4-5
M2M 対象セル

≈レベル2
-5

まとめ

計算科学技術特論B 第14回
74 /74・分子内相互作用と分子間 nonbonded 相互作用の計算, およびFMM を用いた Coulomb 相互作用の計算を含む分子動力学計
算について, その並列化特性を解説した.
・MPI 並列性能および演算性能を向上させるデータ構造について, 座標および多極子 (下位, 上位階層べつ) に解説した.
次回は, このデータ構造を基にした 3 次元トーラスネットワーク
上でのMPI 並列化技術, および OpenMP, SIMD 並列化技術を具体的に説明する.