計算科学技術特論Ａ

古典分子動力学計算法の高速化
2017年7月6日
名古屋大学大学院工学研究科
附属計算科学連携教育研究センター
吉井範行

内容
ＭＤ計算の要素
各要素の高速化、並列化
ＭＤ計算とは• MD計算の流れ、取り扱う原子数と周期境界条件、相互作用（力）、分子モデルと運動方程式• ＭＤ計算の例、歴史、現状
非並列のＭＤ計算• 短距離相互作用（bookkeeping, cell index法）• 長距離相互作用（クーロン相互作用） Ewald, PME法
• 効率化のための方法
拘束動力学（Δtを大きくする）
Multiple time step（Δtを大きくする）
On the flyによる最適化（収束計算をなくす）
高速化のために、キャッシュの有効利用、キャッシュミス
ＭＤ計算の並列化• 高速多重極展開法(FMM)
データ構造 メタデータ法、通信量の最小化、データ局在化による計算効率の向上
ベンチマークMODYLAS

分子動力学(MD)計算とは• 多数の原子や分子が集まってできた分子集団系• 原子、分子、もっと大きな塊の力を計算• 運動方程式を数値積分し、運動を追跡• 得られた原子、分子、もっと大きな塊の軌跡から熱力学、統計力学的性質を抽出

MD計算の流れ• 初期条件の設定• 分子間力の計算• 運動方程式の数値積分
時間刻みΔt• 初期状態の平衡化• 平衡化できたか？
• 軌跡、物理量の出力• 十分に統計はとれたか？
• 最終座標・速度の出力

取り扱う原子数と周期境界条件
・原子数100~100,000,000 個 （長距離相互作用のある場合）
対象系の規模や物理現象による。
・基本セルの一辺の長さ・周期境界条件
基本セル （立方体、直方体、平行六面体etc）
イメージセルを基本セルの周囲に配置。
→表面効果の除去
3次元のみでなく、2，1次元方向にイメージセルを配置
2次元：界面、平面膜、2枚の平板にはさまれた状態
1次元： 細孔、棒状分子
・ポテンシャルカット
Lennard-Jones相互作用
Ewald法の実空間部分
通常、原子直径の3倍(3σ)程度。

相互作用（力）•分子間相互作用（力）Vintermolecule (r N ) =2体力の和∑ ∑V (rij )i> jLennard-Jones（LJ）相互作用
クーロン相互作用•

VCoulomb (r ) =12
σ 6   σ VLJ (r ) = 4ε    −   
 r    r q1q2r

分子内相互作用（力）
分子内の原子を剛体的に取り扱わず、フレキシブルに運動できるとする。（フレキシブルモデル）
分子内ポテンシャルを伸縮、変角、torsionなどで表現する。

Vintramolecule = Vst + Vbe + Vtor
122
= ∑ k st ( r − r0 ) + ∑ k be (θ − θ 0 ) + ∑ kntor {1 + cos ( nφ − φ0 )}n 2分子間相互作用の精度がMD計算の質を決める

分子モデルと運動方程式
分子モデル

単原子
分子

質点

座標変数

運動方程式
運動方程式

質点の位置
座標

ニュートンの
第二法則

デカルト座
標

Ar, Na+質量中心の位質量中心の位
置
ねじれ運
動なし

剛体回転
子モデル

H2O, CO2,ベンゼン

分子ニュートンの
第二法則

オイラー角
四元数
四元数

オイラー方程式

原子の位置座標

ニュートンの第二法則

結合長の拘束結合長の拘束ラグランジュ
未定乗数法

一般化座標

拘束動力学
拘束動力学

多原子
分子
ねじれ運
動あり

原子の位置座標

C4H10、
高分子

フレキシブ
ルモデル

分子内ポテンシャル

ニュートンの
第二法則

デカルト座
標

MD計算の例

名古屋大学 岡崎研HP http://simulo.apchem.nagoya-u.ac.jp/animation1/

MD計算の例

MD計算の例

ＭＤ計算の例
対象
低分子液体、溶液

系・現象
水、水溶液、
溶媒和、化学反応

対象
タンパク質

固体
不均一系、界面

相転移、フォノン
ＤＮＡ、ＲＮＡ
結晶成長、吸着、トライ 生体膜
ボロジー、蒸発、凝縮、透過、分配、界面増強ス
ペクトル、界面張力

系・現象
折りたたみ、酵素、イオンチャ
ンネル、阻害剤、タンパク質複合体
脂質膜、ラフト、抗生物質、薬物透過、吸収両親媒性分子水溶液 球状ミセル、棒状ミセル、 ウイルス
ミセル、単分子膜、二 Ｌ膜、ＬＢ膜、ベシクル
分子膜、バイセル
溶液内化学反応、
化学反応
ガラス
電解反応、
自由エネルギー曲面
超臨界流体
分子系、混合系、クラス 金属、半導体
ターの生成

カプシドタンパク質

高分子

対流、ずり、核生成、流体力
学

分離膜、電解質膜、複雑 非平衡マクロ動力学系、スローダイナミクス

ガラス化、スローダイナミクス

電極、液体金属

ＭＤ計算の歴史
発表年
1953
1957
1957
1964
1969
1971
1971
1977
1977
1980
1980
1981
1984
1985
1986
1987
1988
1991

手法
MC
MC
MD
MD
MC
MD
MC
MD
MD
MD
MD
MD
MD
MD
MD
MD, MC
MD, MC
MD

対象
初めてのMC、剛体球
球状分子
初めてのＭＤ、剛体球
球状分子
水
水（回転の運動方程式）
イオン
拘束動力学
タンパク質分子圧力一定
圧力テンソル一定
経路積分、量子系
温度一定
第一原理ＭＤ
溶液内化学反応
FMM
汎用ポテンシャル
非断熱量子動力学

発表者
Metropolis et al.
Wood and Parker
Alder and WainwrightRahmanBarker and Watts
Rahman and Stillinger
Woodcock
Rychaert et al.
McCammon et al.
Andersen
Parrinello and Rahman
Chandler and Wolynes
能勢
Car and Parrinello
Hynes et al.
Greengard
Jorgensen et al.
Rossky et al.

MDの現状・汎用ポテンシャルセットの開発
対象系の分子のポテンシャルが簡単に構築できる・粗視化モデル
大規模系が取扱い可能に

プログラム名開発者

AMBER

Kollman

CHARMM

KarplusGROMACS

Hess

ポテンシャル名

開発者

AMBER

Kollman

TINKER

PonderCHARMM

Karplus

NAMD

Schulten

OPLS

Jorgensen

DESMONDShaw

GROMOS

Berendsen

LAMMPS

Plimpton

MODYLAS

Okazaki

GENESIS

SugitaMRABLE

Ikeguchi非並列のMD計算
シングルコアで高速化を進める
MD計算手法の選択
短距離力の計算（bookkeeping, cell index）
長距離力の計算（Ewald, PME, FMM）
拘束動力学の利用
Multiple time stepの導入
on the fly計算の利用• 高速化のための工夫

短距離相互作用 bookkeeping法• ある分子からの距離がrcよりも短い相手分子の番号すべてをリストアップ．
• リストの作成は何ステップかに1回毎.
固体のような流動性のないものなら更新不要• rcよりもΔrだけ長い距離のものまで余分にリス
トに載せる．
• リストを用いて相互作用を計算する際には，リストを複数のノードやスレッドに割り振って並列計算させるループ分割法が容易に実装可能．

L
9σ
12σ
15σ
18σ

N
600
1500
3000
5000

N’
150
150
150
150

L/rｌ
2.6
3.4
4.3
5.1

N’/N
0.25
0.10
0.05
0.03

短距離相互作用 bookkeeping法
k=0
do i=1,n-1
do j=i+1,n
i, j間の距離がrl以下
k=k+1
list(1,k)=i リスト作成
list(j,k)=j リスト作成
enddo
enddo
do i=1,k
リストをもとに相互作用計算
enddo
MPI並列なら相互作用の計算結果をノード間で通信
スレッド並列
通信不要• 数ステップに一度相互作用対のリス
ト(list)を更新する• 毎ステップリストを用いて相互作用計算する。

短距離相互作用 Cell index法• 帳簿を作成，更新する際のN(N-1)/2通りの距離計算が大きな負荷.大きなメモリを要する．
• L/ rlが4以上の場合には，一辺の長さがrlよりも大きくなる範囲内で基本
セルそのものをx，y， zそれぞれの方向に4分割以上の分割を行う．
• 注目している分子が属しているセルとそれに直接接しているセルあわせて27個の分割セル内にある分子との間だけでbookkeeping法を行う．
• 分割セルの一辺の長さを よりも短くする場合には，第二近接のセルまで考慮に入れる• Cell index法についても作成したリストを複数のノードやスレッドに割り振って相互作用の並列計算を行えばよい．

長距離相互作用 クーロン相互作用• 直接計算 O(N2)→計算方法はＬＪ相互作用と同じ
カットオフの悪影響大のため用いない。
• エワルド法 O(N2) あるいは
Particle Mesh Ewald法 O(NlogN)→小、中規模系向け• 高速多重極展開法（Fast Multipole Method: FMM) O(N)
→大規模系向け

長距離相互作用 Ewald法
Ewald法の要点
法の要点
・減衰の遅いクーロン相互作用を、減衰の速い関数と遅い関数に分ける。
・遅い関数はフーリエ変換することにより逆格子ベクトルを用いて評価。
・速い関数は実空間で評価。
VN (r N ) =Qi Q j1
1′
2 ∑∑∑j 4πε 0 | ri − r j + Ln |n i

（減衰の遅い1/r）
減衰の遅い ）＝ （減衰の速い関数）＋（減衰の遅い関数）

1 erf (α r ) erfc(α r )
=
+
r
r
r

Ln

減衰の
減衰の速い関数
速い関数 ・・・ この関数形で評価
erfc(α r )r

減衰の遅い関数
減衰の遅い関数 ・・・ フーリエ変換したものでerf (α r )評価rn 整数ベクトル(−1, 2, 3)

長距離相互作用 Ewald法• Ewald法のよるポテンシャル
VN = VN(1) + VN(2) + VN(3)VN(1) =

bookkeeping法
セルインデックス法

Qi Q j erfc(α | ri − r j + Ln |)
1
′
2 ∑∑∑
| ri − r j + Ln |実空間での相互作用
n ij 4πε 0

LJと一緒に計算する

VN(2) =22exp(− | G |2 /4α 2 )     ∑ Qi cos ( G ⋅ ri )  +  ∑ Qi sin ( G ⋅ ri )  ∑
2ε 0 L3 G ≠ 0| G |2 i   i1

Ｇについてのループ分割
VN(3) = −Qi2 α
∑ 4πε 0 πi

→ PMEによる高速化

Particle mesh Ewald（ＰＭＥ）法
逆格子ベクトルの項は
 | G |2 
exp  −
1
4α 2 
(2)
VN =Qi exp(iG ⋅ ri )∑ Q j exp(−iG ⋅ r j )∑3 ∑2
2ε 0 L G ≠ 0
|G |
i
j
=

1

H (G ) S (G ) S (−G )
∑
2ε L
0

3

G ≠0

S(G)を毎ステップ計算しなければならない。
S(G)を格子点uにおける補間を用いて表現すると
S (G ) = ∑ Sˆ (u) exp(iG ⋅ u)uとなる。ここでS(u)はS(G)のuにおける重み因子。これによりVN(2) =

1

Hˆ ∗ Sˆ (u) Sˆ (u)
∑
2ε L
0

3

u

S(G)の計算に補間法を利用（格子点で評価）
高速フーリエ変換（FFT）利用 O(NlogN)

 G2
exp  − 2  4α 
H (G ) =G2

S (G ) = ∑ Qi exp(iG ⋅ ri )i

効率化のための方法
Δtを大きくする• 拘束動力学（⇔剛体回転子）• Ｍｕｌｔｉｐｌｅ time step
収束計算をなくす• On the fly による最適化

拘束動力学• 拘束条件つきの運動方程式
拘束条件

運動方程式

ラグラジアン

拘束力のある場合
ラグランジュの未定乗数λオイラー・ラグランジュ方程式

運動方程式

運動方程式

Multiple time step
• 同じMD計算のなかで変化の速いものと遅いものがある場
合には，速い変化のほうに対応した小さなΔｔを選択しなければならない．
例
United atomのブチルアルコールのフレキシブルモデル

CH3-CH2-CH2-CH2-O-H∆t = 2fs∆t = 0.1fs• Hに対する運動方程式とほかの重い原子に対する運動方
程式を区別して，異なったΔtで解くことができれば非常に都合のよい（部分的にΔtを大きくできる）ことになる．この考え方に基づいた方法がMultiple time step法である．

Multiple time step
• 時間発展演算子法により、重い（遅い）運動と軽い（速い）運動とを分離して数値積分する。
• 時間発展演算子法
リウビル演算子 L

時間発展演算子

異なる自由度など2つの部分に分離できると

eiL∆t = eiL p ∆t /2 eiLr ∆t eiL p ∆t /2 + O(∆t 3 )Multiple time step
• 時間発展演算子法により、重い（遅い）運動と軽い（速い）運動とを分離して数値積分する。
リウビル演算子 L

時間発展演算子

eiL∆t ≈ eiLheavy ∆t /2 ei( Lr +Llight )∆t eiLheavy ∆t /2iLheavy ∆t /2≈eiLlight ∆t δ t 2 i ( Lr )δ t iLlight δ t 2 n iLheavy ∆t /2

{e

e

e

}e

On the fly による最適化•••

第一原理MD（Car-Parrinello法）
電荷揺らぎモデル（TIP4P-FQ、SPC-FQモデル）
分極モデル（点双極子モデル）

•
•
•
•
•

最適値を得るために繰り返し計算を必要する変数を含む
繰り返し計算を行わず、最適値近傍を動的に揺らがせる近似解で満足できる場合
その変数に仮想的な慣性量、運動エネルギーを与えて運動方程式を解く。
繰り返し計算がなくなるため高速化が可能。
運動エネルギーの相当する分だけ最適値からずれて、誤差が生じる。
分極モデル

U = U 0 + U pol
1
U pol = −∑ Ei ⋅ µi +µi ⋅ µi∑
2α i
iEi = Tij ⋅ µi

µi = α Ei

拡張系のラグラジアン
粒子と双極子モーメントの自由度がある。
1 N
1 N
2
L ( r , µ ) = ∑ mri + ∑ M µ i2 − U
2 i =12 i =1NN

双極子モーメントの運動方程式Mµ i =µ∂L ( r N , µ N ) = − i + Ei∂µα

高速化のために遅い演算（除算と組み込み関数）• 除算は他の四則演算より遅い。
• 組み込み関数も時間のかかる演算である。
• ホットスポット内ではできるだけこれらの演算を減らすように工夫する。
do i=1,nb(i)=a(i)/x
enddo

do i=1,nx(i)=a(i)/b(i)/c(i)enddo

y=1.d0/x
do i=1,nb(i)=a(i)*yenddo

do i=1,nx(i)=a(i)/(b(i)*c(i))
enddo

除算をループの外へ移動

除算を乗算に変更

do i=1,nx(i)=a(i)/b(i)y(i)=c(i)/d(i)enddo
do i=1,ntmp=1.d0/(b(i)*d(i))x(i)=a(i)*d(i)*tmp
y(i)=c(i)*b(i)*tmp
enddo
通分の利用

高速化のために遅い演算（組み込み関数）
組み込み関数のfortranコード化 環境によって速くなる場合とそうでない場合がある。両者の比較が必要。

• ex×eyex+y• logax + logay
logaxy• sinθ ×cosθ
0.5sin2θ• If(sqrt(x*x+y*y) < R)
If(x*x+y*y < R2)
y=sign(x,a)
If(a > 0.d0) theny=abs(x)else
y=-abs(x)endif

imax=0do i=1,n
imax = max(imax,m(i))
enddoimax=0do i=1,nif(m(i) > imax) imax=m(i)enddom=nint(a)If(a > 0.d0) thenm=a+0.5d0
else
m=a-0.5d0
endif高速化のために分子シミュレーションへの適用• Ｌｅｎｎａｒｄ－Ｊｏｎｅｓ相互作用121 1V =   − r r

6

r 6=r 2 *r 2 *r 2
r 12 =r 6*r 6
V=1.d0/r 12 -1.d0/ r6
r 6=r 2 *r 2 *r 2
r 12 =r 6*r 6
V=(1.d0 -1.d0*r6)/r 12
LJ相互作用への適用

キャッシュの有効利用
CPUの必要なデータに高速アクセス① データはメモリに格納② 必要に応じて隣接するデータと共にキャッシュ上に③ 演算で必要となるとレジスタにコピー④ 演算された結果は再びキャッシュに上書き⑤ キャッシュ上で用いられなくなるとメモリ上に上書きコピーされる．

キャッシュミス
• データ：レジスタ←キャッシュ←メモリ• レジスタ←キャッシュ：高速アクセス• キャッシュ←メモリ：←低速アクセス• 演算を行う時にデータがキャッシュ上にない
ときはメモリまで取りに行く。キャッシュミス。
• データがキャッシュに乗っているように工夫することが重要。

キャッシュミスを防ぐ• 1次元配列の変数：引数の順にキャッシュに格納• 1次元配列の変数へのアクセス→配列の引数について連続的にアクセスするように心がける。
• 2次元以上の配列の変数：Fortranでは配列左側の引数の順にキャッシュに格納。
（C言語は右側から）• 2次元以上の配列の変数へのアクセス→Fortran
では配列の左側の引数について連続的にアクセスするように心がける。 A(I,j)（C言語は右側）ＭＤ計算の並列化• ループ分割：計算負荷の大きなdo ループを並列化⇒並列化が容易⇒全ノード間で通信（MPI)⇒MPIの高並列になると効率が悪い• 領域分割：基本セルを複数の領域に分割⇒並列化が大変。
⇒通信は隣接領域間のみ。
⇒超並列でも効率がよい。
１
データ構造を大幅に変更
５
ほぼ最初から書き直し
９

２

３

４

６

７

８

１０ １１ １２

１３ １４ １５ １６

高速多重極展開法(FMM)• 多極子モーメントを用いた多重極展開

原点０の周りに電荷Qiが分布
QiとQとの相互作用を多重極展開で表す。
余弦定理
1
=
ri

1
r 2 − 2r ∆ri cos γ i + ∆ri2

テイラー展開

1

=

∆ri

∆r
r 1 − 2 cos γ i +  i 
r
 r 

1
1
3
5
= 1 − x + x 2 − x3 + ⋅ ⋅ ⋅ ⋅ ⋅
2
8
16
1+ x

2

x <1

1 1 ∆ri∆ri2∆ri3
2= + 2 cos γ i + 3 ( 3cos γ i − 1 ) + 4 ( 5cos3 γ i − 3cos γ i ) + ⋅ ⋅ ⋅ ⋅ ⋅ri r r2r
2r

高速多重極展開法(FMM)• QiとQとの相互作用の多極子モーメントを用いた表式
V =

1

QQ

i
∑
4πε 0
ri
i

 1 ∆r∆r 2
Q ∑ Qi  + 2i cos γ i + i3 ( 3cos2 γ i − 1 ) + ⋅ ⋅ ⋅ ⋅ ⋅ 
4πε 0 i2r
r r
 ∑ Qi

1
1
1
1
2

i
=
Q
+ 3 r ⋅ ∑ Qi ∆ri + 5 r ⋅ ∑ Qi 3∆ri ∆ri − ( ∆ri ) I ⋅ r + ⋅ ⋅ ⋅ 


4πε 0
r
2
r
r
i
i


=

1

{

Z r ⋅ µ r ⋅Q ⋅r=
Q  + 3 ++ ⋅ ⋅ ⋅ 
5
4πε 0  r

r
r
1

}

Z = ∑ Qii

µ = ∑ Qi ∆ri
iQi
3∆ ri ∆ri − ( ∆ri )2 I
i 2

Q=∑

Legendre多項式を用いた表式
1 1∆ri n=Pn ( cos γ i ) ri r ∑
 r 
n

1 Q∆ri nV =Qi ∑ Pn ( cos γ i ) r
4πε 0 r ∑in

展開のしかた・多極子モーメント
・Legendre多項式
・球面調和関数等々{}高速多重極展開法(FMM)• 空間分割

２次元：４分木構造
３次元：８分木構造

レベル

サブセル一辺の長さ

0

L

1

L/2

2

L/22

3

L/23

最も細かく分割したレベルのサブセルについて領域分割を適用

高速多重極展開法(FMM)• 粒子と分割されたサブセルとの相互作用

周期境界条件あり

自セルAの粒子iのポテンシャルܸ௜Vi = Vi d + Vi m

距離に応じてレベルの異なる多極子モーメントを利用する高速多重極展開法(FMM)• 多重極展開と局所展開

ܸ௜m

ܸ௜d
手順(d) より遠いサブセル群E,FについてもDと同様。M(2), M(1)を求め、M2L、L2LによってL(4)を求める。
手順(e) 第３隣接以遠無限遠に至るイメージセルにM(0)があるものとして、基本セル上に作るL(0)を求める。L2LによりL(4)を求める。
手順(f) L(4)を用いてポテンシャルを評価する。
これを効率的に行うために手順(c), (d), (e)で得られたLそれぞれのレベルで足し合わせてからより小さなサブセルにL2L変換する。

高速多重極展開法(FMM)
計算精度• 静電+Lennard-Jones相互作用
原子数10,125,000
個
体積
46.83
nm3
セル数 643 (=262,144) 個 (セルあたり39個)
V / 10-13 J

Fx / 10-10 N

Fy / 10-10 N

Fz / 10-10 N

correct-218735559-4.357177

6.380146-1.477348

４次-218735742-4.356659

6.380802-1.475152

８次-218735559-4.357201

6.380145-1.477340

エネルギー
力

展開
4次
８次
4次
８次

6～
～7桁
桁
8～9桁
4～
～5桁
桁
6～7桁

• 計算精度を展開次数で制御

データ構造メタデータ法原子のデータ構造• 原子座標、速度、力、セグメントデータ → メタデータ化• プロセスに属するサブセルを割り当てて固定 （図では４×４のサブセルを1つのプロセスに割り当て）• サブセル内の原子のデータをプロセスに局在化させる。グローバルに持たない→省メモリ。
• データ領域に空きを用意 → サブセル内の原子数の増減に対応。各サブセル内の原子数に合わせて要素数が変化。要素数情報もあわせて保持。
• メタデータ配列のまま演算も通信も行う。サブセルを表す3次元配列を持っているので、取り扱いが容易• バッファリングやindirect accessがない。
• サブセル単位のメタデータがキャッシュに載るようにサブセル内の原子数を調整。
• 原子の運動に伴って原子の帰属するサブセルを更新 。

Cell index Y・セグメントデータ
一塊の原子団のデータ
CH3、H2Oなど。
原子間距離拘束の情報

Cell index Z

通信量の最小化

Rank 2

Rank 3
Step-1
Step-1

Step-2

Step-2Cell index Y

自セルと相互作用する原子情報を集める。
• 通信回数と通信量の最小化• 隣接セルと通信して1次元（z軸）方向にデータを結合し棒状のデータを作る。
• 棒状データをy軸方向に通信し結合させて、平面状データにする。
• x軸方向に平面状データを通信し結合させて、立方体状のデータにする。
• K－スパコンの多重同時通信（各軸ごとに＋、－方向）

原子情報

Step-2

Step-2Cell index Z

Rank 1

Rank 0

Step-1
Step-1

データ局在化による計算効率の向上
キャッシュの利用効率の向上• L1キャッシュ－L2キャッシュ－メモリ• L2－メモリ間、L1－L2間のデータ通信の最小化
最重要ホットスポット• 2体間相互作用計算• FMMのM2L計算
Cell index Y

2体間相互作用計算• 3重ループ構造
do 5つのサブセルを指定するインデックス
Cell index Z
do 注目する1つのサブセル内の粒子
do 隣接する5つのサブセル内の原子相互作用計算
enddo
enddo
enddo
• これにより、5つのサブセルのデータが効率的に再利用される。
• ただし5つのサブセルの座標データをL1キャッシュに載せておくことが重要

MODYLAS
長距離力の厳密計算
実質的に FFT free !

並列化
・完全領域分割・多階層隣接通信化
多極子展開
多極子展開
i

i

×
×

大きな通信はすべて隣接通信化
Qi = ∑ qss∈iµi = ∑ qs rss∈iΘ i = ∑ qs rs rs

+

s∈i

×

×

局所展開
局所展開
局所展開

多極子展開
・
・
・

+

周期境界条件下における多極子に対するエワルド法

×

=

厳密計算
（精度は展開次数で制御）

ウイルスカプシドの分子動力学計算
全原子計算
マイクロ秒程度の軌跡
マイクロ秒程度の軌跡~１０００万原子系
１０００万原子系
・周期境界条件下・長距離相互作用

小中規模系の計算
・２～３万原子系
・Ｅｗａｌｄ法
領域分割
ＦＦＴ

京スパコン
専用計算機

１００ｍｓ／ｓｔｅｐ・大規模並列が不可欠
・通信もｍｓ・ＡＮＴＯＮ
５１２並列数µｓ／ｓｔｅｐｓ／ｓｔｅｐ・ＭＤＧＲＡＰＥ

従来のスパコン
従来のスパコン

隣接通信化が不可欠（ＦＦＴは不可）

1マイクロ秒：１０年
マイクロ秒：１０年（解析不可能）
秒：１０年（解析不可能）（実効性能 0.1ペタフロップス
0.1ペタフロップス）
ペタフロップス）

MODYLAS

ベンチマークテスト
PYP 水溶液、1000万原子系
Ncore
10000

100000T / ms1000

1000000
10000

10000

1000

1000

100

100

10

10

1

1

T / ms1000

acceleration

10000

Ncore
1000

10000

1000001000000

pairwise additive forces
FMM calculation
intramolecular forces

100

10

1

communication time
0.1

100

1000

10000

0.1
100000

Nnode

The measured overall calculation time
per one MD step(Dt) for the 10-million
atomistic system, the PYP solution, andthe acceleration ratio with respect to the
64-nodes calculation as a function ofdegree of parallelism, Nnode.

0.1

100

1000

Nnode

10000

100000The measured direct calculation time of the pairwise
additive forces(circle), the FMM calculation(triangle),
the intramolecular forces(diamond), and thecommunication time(brown) per one MD step(Dt)for the 10-million atomistic system, the PYP solution,
as a function of degree of parallelism, Nnode.

今日のまとめ
ＭＤ計算の要素
各要素の高速化、並列化
ＭＤ計算とは• MD計算の流れ、取り扱う原子数と周期境界条件、相互作用（力）、分子モデルと運動方程式• ＭＤ計算の例、歴史、現状
非並列のＭＤ計算• 短距離相互作用（bookkeeping, cell index法）• 長距離相互作用（クーロン相互作用） Ewald, PME法
• 効率化のための方法
拘束動力学（Δtを大きくする）
Multiple time step（Δtを大きくする）
On the flyによる最適化（収束計算をなくす）
高速化のために、キャッシュの有効利用、キャッシュミス
ＭＤ計算の並列化• 高速多重極展開法(FMM)
データ構造 メタデータ法、通信量の最小化、データ局在化による計算効率の向上
ベンチマークMODYLAS

参考文献
ＭＤ全般• D. Frenkel, B. Smit，“Understanding Molecular Simulation," Academic Press, San Diego（１９９６）．
• 上田 顕，「コンピュータシミュレーション」，朝倉書店（１９９０）．
• M. P. Allen, D. J. Tildesley，“ Computer Simulation of Liquids," Oxford Science, Oxford（１９８７）．
• J. P. Hansen, I. R. McDonald, “Theory of Simple Liquids," Academic Press, London（１９８６）．
• 岡崎 進, 吉井 範行, コンピュータ・シミュレーションの基礎（第2版）, 化学同人（２０１１）．
力学• H. Goldstein，「ゴールドスタイン新版古典力学（上），（下）」（瀬川富士，矢野 忠，江沢康生 訳），
吉岡書店（１９８３）．
数値積分• G. J. Martyna, M. E. Tuckerman, D. J. Tobias, M. L. Klein, Mol. Phys., ８７，１１１７（１９９６）．
Particle mesh Ewald法• T. Darden, D. York, L. Pedersen, J. Chem. Phys., ９８，１００８９（１９９３）．
FMM• L. Greengard, V. Rokhlin, J. Comput. Phys., ７３，３２５（１９８７）．
SHAKE• G. Ciccotti, J. P. Ryckaert, Comp. Phys. Rep., ４，３４５（１９８６）．
• P. Gonnet, J. Comput. Phys. ２２０, ７４０（２００７）．
MODYLAS• Y. Andoh, et al., J. Chem. Theory Comput. ９, ３２０１（ ２０１３）• Y. Andoh, et al., J.Chem. Phys. １４１, １６５１０１ (２０１４).