計算科学技術特論Ａ

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

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

MD計算の例

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

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

相互作用（力）•分子間相互作用（力）Vintermolecule (r N )   V (rij )
2体力の和
i j

Lennard-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、
高分子

座標変数

運動方程式

質点の位置
座標

ニュートンの第二法則

質量中心の位置

ニュートンの
第二法則

オイラー角
四元数

オイラー方程式

原子の位置座標

ニュートンの第二法則

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

デカルト座
標

一般化座標

拘束動力学

原子の位置座標
フレキシブ
ルモデル

分子内ポテンシャル

ニュートンの
第二法則

デカルト座
標

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

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

対象
タンパク質

固体
不均一系、界面

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

両親媒性分子水溶液 球状ミセル、棒状ミセル、ウイルス
ミセル、単分子膜、二 Ｌ膜、ＬＢ膜、ベシクル
分子膜、バイセル
溶液内化学反応、
化学反応
ガラス
電解反応、
自由エネルギー曲面
超臨界流体
分子系、混合系、クラス 金属、半導体
ターの生成
高分子

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

系・現象
折りたたみ、酵素、イオンチャ
ンネル、阻害剤、タンパク質複合体
脂質膜、ラフト、抗生物質、薬物透過、吸収カプシドタンパク質

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

電極、液体金属
対流、ずり、核生成、流体力
学

ＭＤ計算の歴史
発表年
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

Karplusポテンシャル名

開発者

GROMACS

Hess

AMBER

Kollman

TINKER

PonderCHARMM

Karplus

NAMD

Schulten

OPLS

Jorgensen

DESMONDShaw

GROMOS

Berendsen

LAMMPS

Plimpton

MODYLAS

Okazaki

GENESIS

SugitaMRABLE

IkeguchiMD計算の方向性• 大規模化（原子数増）
高並列コンピュータ＋高並列用ソフトウェア• 高精度化（アンサンブル増）
高並列コンピュータで並列同時実行＋汎用ソフトウェア• 長時間化（ステップ数増）
専用計算機＋専用ソフトウェア
ANTON by Wikipedia

大規模化・高精度化
ユーザー• 問題設定• 必要な計算資源はどの程度か• ソフトはあるか• プレ・ポスト処理できるか
初級の開発者• 機能拡張を自分で行う• ソフトはソース書き換え可能か（ ソース公開。技術的にも）上級の開発者• 問題規模• アルゴリズムの選択• データ構造• 非並列のしっかりしたプログラム作成• 並列化（MPI, thread…）• 開発期間内容
ＭＤ計算の要素
各要素の高速化、並列化
ＭＤ計算とは• MD計算の流れ、取り扱う原子数と周期境界条件、相互作用（力）、分子モデルと運動方程式• ＭＤ計算の例、歴史、現状
非並列のＭＤ計算• 短距離相互作用（bookkeeping, cell index法）• 長距離相互作用（クーロン相互作用） Ewald, PME法
• 効率化のための方法
拘束動力学（Δtを大きくする）
Multiple time step（Δtを大きくする）
On the flyによる最適化（収束計算をなくす）
高速化のために、キャッシュの有効利用、キャッシュミス
ＭＤ計算の並列化• 高速多重極展開法(FMM)
データ構造 メタデータ法、通信量の最小化、データ局在化による計算効率の向上
ベンチマークMODYLAS非並列の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 j
1
1
′
2 
n i
j 4 0 | ri  r j  Ln |（減衰の遅い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法
セルインデックス法

Q Q erfc( | ri  r j  Ln |)
1
′i j

2 n i j 4 0| ri  r j  Ln |実空間での相互作用

LJと一緒に計算する

VN(2) =

2
2
exp( | 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 ) =2
G
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
Dt = 2fs

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

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

時間発展演算子

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

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

時間発展演算子

eiLDt  eiLheavy Dt /2 ei  Lr Llight  Dt eiLheavy Dt /2
iLheavy Dt /2eiLlight Dt  t 2 i  Lr  t iLlight  t 2 n iLheavy Dt /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 ii
Ei = Tij  μi

μi =  Ei

拡張系のラグラジアン
粒子と双極子モーメントの自由度がある。
1 N
1 N
2
L  r , μ  =  mri   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の高並列になると効率が悪い• 領域分割：基本セルを複数の領域に分割⇒並列化が大変。
⇒通信は隣接領域間のみ。
⇒超並列でも効率がよい。
１
データ構造を大幅に変更
５
ほぼ最初から書き直し

２

３

４

６

７

８

１０ １１

１２

１３ １４ １５

１６９

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

Cell indexY・セグメントデータ
一塊の原子団のデータ
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キャッシュに載せておくことが重要

静電相互作用• 電気を帯びた系：水、電解質溶液、タンパク質、脂質な
どの両親媒性分子、イオン液体、、、• 周期境界条件 カットオフ
Ewald法 O(N2)
Particle Mesh Ewald法 O(NlogN)
Darden,et al. (1993),Essman, et al. (1995)

Multilevel summation法 O(N)D. J. Hardy, et al. (2015)

高速多重極展開法（FMM） O(N)Greengard, Rokhlin (1987)

レ
多

極

子

モ

ー

メ

ン

ト
局

M2M

所

展

開

係

ル

３

数

G

L2L

第二隣接イメージセル

M2L

ベ

第三隣接イメージセル

FMMの概要

レ

ベ

ル

４

F

M2L

D
D

P2M

L2P

D

B

レ

C

D

B
A

B

互

作

用

の

計

算

方

法

ル

４

B

E

子

B
i

．

D
C
B
A

基本セル

相

原

ベ

第一隣接イメージセル

P2M

原

子

と

サ

ブ

セ

ル

と

の

相

互

作

用

．

GreengardらによるFMMの定式化
Pr
Q

球面調和関数 n  m ! m
Pn  cos   eim n  m !

Ynm   ,   =  , ,  ×nYnm   ,  1n m=  Yn   ,   r n1
r  n
=0 m= nnYnm   ,  qin m  P  = =   qi  Yn   ,  r  n =0 m = n
r n 1多極子モーメント

M nm =  qi inYn m  i , i i

 r,, FMMにおける変換• P2M
M nm =  qi inYn m   i , i 
i

• M2M
j

M jk =

M kj nmi k  m  k  m Anm Akj nm  nYn m   ,  

n

 

Anm =

Akj

n = 0 m = n 1  n n  m  n  m • M2L
M nmi k  m  k  m Anm Akj Ynmjk   ,  

n



Lkj =



 1 n Anmjk  n  j 1

n =0 m = n

• L2L


Ljk =n 

Lmn i m  m k  k Anmjk Akj Ynmj k   ,    n  j 1 n  j Anm

n = 0 m = n

• L2P
i =

qi



j

LkjY jk  i , i ri j


4 0
j = 0 k = j

nmax

j





fi = qi   Lkji ri nYjk i ,i 
j =0 k = j

 i ri jY jk  i , i   ri j 1Y jk11  i , i 

Solid harmonics化
FMMにおける変換

Legendre陪関数Pm  cos   =   sin  

•

dm

m

 d cos  m

P2M

M m =  qi Rm  ri P  cos  

 m

i


1
d
m
2
=    sin  
cos

1

 m
2 !
 d cos  



•

M2M



Singular solid harmonics

  m    m ! m

r

 1

M  =   M m R

•

M2L

nmax



L =   M  S   m  mrPm  cos  eim r , ,   =


m
!



S m  r ,  ,   =  1 m

 =0  = Regular solid harmonics
Rm

nmax

P  cos  e im = 0  = 

•

L2L

nmax



L =   L Rmm
•

 =   = 

L2P

=

1

nmax


m


LM
4

m


0  = 0 m = 1=   S  m  rS  Rm  rR 
rS  rR = 0 m = 

 1 m1 m1 
 2 R1  R1 nmax 
m 1
m1
m1 
fi = qi   L  R1  R1
 2

 =0 m=


Rm1













実部、虚部の対称性の利用（実数化）• Regular, singular solid harmonicsや多極子モーメント、
局所展開係数の実部、虚部の間にある対称性を利用してこれらの要素数を減らし、変換の演算量の削減を
図る。
nmax



L =   L Rmm


 =   = 

実部、虚部の間にある対称性

Rnm =  1 m Rn*mmの非負部のみ取り扱えばよいn  m  n0mn

Snm Lmn M nm も同様の式が成り立つ
要素数 nmax  1 2 =25 (nmax=4)

 nmax  1  nmax  2 
2=15

実部、虚部の対称性の利用（実数化）
•

M2M

nmax



m

M  =   M

m





R

演算量＝要素数2
演算量の減少 (nmax=4)

 = 0  = 

•

M2L

nmax

2 2





m


L =  M S•

  nmax  1   m   


152
1
=
252 2.78

 = 0  = 

L2L

nmax



  nmax  1  nmax  2  2L =   L Rmm


 =   = 

実数化

nmax

2



L =   L Rmm


m

Rm =  1  Rm =   = nmax

 0  m    m 1    m =   L R    L R    L R  
 = 
 =1
 = 


Rm =  1 m 1

Rm

nmax

 0  m    m       m =   L R    L R    L R  
 = 
 =1
 =1

nmax




m 0
m
m2
m
=   1 LR    L Rm   1 Rm  L Rm   1 Rm 
 = 
 =1

nmax



m 1 0m 1m 1m
 i   1 L R    L Rm   1 Rm  L Rm   1 Rm 
 = 
 =1


 

 













直交座標の漸化式の利用• 従来は、P2M,L2Pにおいて、原子座標として直交座標を球面座標に変換し、球面調和関数を計算していた。
• ここでは、直交座標から直接solid harmonicsを求める漸化式を導入し、座標変換を消去する。
rでスケールされたLegendre陪関数Pm  r , z  = r   r sin   m Pm  cos  Regular solid harmonicsPm r , z  x  iy m

= r Pm cos  e

Rm  x, y, z  =im1Pm  r , z  x  iy m   m !

Singular solid harmonics
   m !  m
m
S m  x, y , z  =P  r , z  x  iy 
2  1rm≧0についての漸化式 1 m  2m  1 !!
Pm  r , z  =  2  1  zPm1  r , z  2  1
  m 1 2 m
zPm1  r , z  r P  2  r , z 
m
m

 = m  = m  1
  m  2K. Nitadori, Particle mesh multipole method: An efficient solver for gravitational/electrostatic forces based on multipole method and fast
convolution over a uniform mesh.

多重ループの１重ループ化• M2M,M2L,L2Lにおいて、solid harmonicsの添え字についての多重ループが存在する。
• １重ループ化しループ長を伸ばすことにより、SIMDなどのマシン性
能を引き出す。
２重ループ
nmax



Lm =   M  S  m    nmax  1  nmax  2  2 = 0  = 

最内側ループ長 １５ → ２２５
１重ループ

2

YN7

Benchmark test
•
•
•
•
•

•

水系（19530原子）1000MDステップ実行
周期境界条件
NVEアンサンブル
実行環境
京コンピュータ
8ノード（8MPI, 1Thread)
MODYLAS使用
FMMの展開次数 nmax=4
球面調和関数版
改良版

M2L:
60%
non FMM
35%

FMM
65%

L2P:
3%
M2M:
0%
comm_global:
0%

P2M:
1%

L2L:
0%
comm_local:
1%

FMM_Ewald:0%

球面調和関数使用時、FMMは全体の６５％。
M2Lだけで６０％を占める。
P2MやL2Pで４％ほど。
FMMの中ではM2L,P2M,L2Pが９９％を占めている。

ス ラ イ ド
Y N 7

4 8
計算条件、全体に対するFMMの割合を⽰す。
YOSHII Noriyuki, 2017/05/11

対称性を利用した実数化（ループ2重） → ループの1重化 → 漸化式利用による座標変換の削除
時間
(秒)
20
18
16
14
12
10
8
6
4
2
0

時間
(秒)
0.4
0.35
0.3
0.25
0.2
0.15
0.1
0.05
0

P2M

時間
(秒)

対称性
座標変換あり漸化式
座標変換なし
３．５倍

M2M

900
800
700
600
500
400
300
200
100
0

時間
(秒)
0.7
0.6
0.5
0.4
0.3
0.2
0.1
0

M2L
対称性 ２．４倍（理論値 ２．８倍）

L2P

時間
(秒)

対称性
座標変換あり

60
50
40

漸化式
座標変換なし
６．９倍
３．５倍

30
20
10
0

L2L

時間
(秒)
1000
800
600
400
200
0

FMM
２．４倍

1600

時間 1400(秒) 1200
1000
800
600
400
200
0

削減部分
38%

L2L:
0%

M2L:
25%

L2P:
1%

FMM
27%
FMM_comm_local:
1%
FMM_comm_global:
0%

non FMM35%
FMM_Ewald:
0%

M2M:
0%

P2M:
0%

YN3

計算精度
• 変更前後における計算結果の比較(1step)
MODYLASの熱力学量に関する標準出力,
M2L(single loop)
M2L(double loop)
original MODYLAS
M2L(single loop)
M2L(double loop)
original MODYLAS

Hamiltonian
potential-E-3.901749600398E-16-4.449136174563E-16-3.901749600398E-16-4.449136174563E-16-3.901749600398E-16-4.449136174563E-16
kinetic-E
temperature
5.473865741657E-17
2.030212809202E+02
5.473865741657E-17
2.030212809202E+02
5.473865741657E-17
2.030212809202E+02• 12ケタ一致。
• 正しい静電相互作用計算を行えている。

ス ラ イ ド
Y N 3

5 1
計算が正しく⾏えていることを⽰す。
YOSHII Noriyuki, 2017/05/11

FMM• FMMを、Greengard & Rokhlinに従う球面調和関数を基底としたものから、Solid harmonicsへと変更した。
• 実部と虚部の対称性を利用による演算量の削減を図った。
• これに伴いP2M、M2M、M2L、L2L、L2P変換の演算量が減少。 特にM2Lが2.4倍速くなった。
• 直交座標から直接solid harmonicsを評価できる
漸化式を利用することにより、座標変換を消去した。これによりP2M,L2Pが3.5倍速くなった。
• FMM全体で2.4倍高速化することに成功した。

ＭＤ calculation
Empty capsid
ＶＰ１、ＶＰ２，ＶＰ３，ＶＰ４
６０ eachSphingosine
６０
Ｈ ２Ｏ
１、８８４、２１８
Ｎａ＋
５，３１０
Ｋ＋
１３６
Ｃｌ－
５、２０６
６、４８０、２３６ atoms

ＮＰＴ ensemble
Ｐ＝１ａｔｍ、Ｔ＝３１０．１５ Ｋ
ＳＨＡＫＥ／ＲＡＴＴＬＥ／ＲＯＬＬ
ＣＨＡＲＭM２２ ｗｉｔｈ ＣＭＡＰ
ＬＪ １２ Å ｃｕｔ ｏｆｆ
Coulomb ＦＭＭ

ＰＤＢ １ＨＸＳ
ｐＨ ７．４
positive charges
negative charges
Dｔ = ０．５ ｆｓ／２ ｆｓ／４ ｆｓ
２００ ｎｓ

４、８６０
５，１００

Negative pressure

instability

Local pressure from virial theoremNoutside

N

N kT
1

rij FijP =V
3V i j iinside( ri V )Z. S. Basinski, M. S. Duesbery, R. Taylor,Can. J. Phys. 49, 2160(1971)- 34 atm
c.f. deinition by Harajima

Definition of local pressure by Irving-Kirkwood
PN ( R ) = k BT ( R )  
m

PT ( R ) = PN ( R ) 

j

k

l

R d PN ( R )2dR(e R  f)(e R  rjljk )1
m jljk ( R )8mR
c jljk r jljk
  P = 0 ： mechanical balanceOrigin of the negative pressure
Local electrostatic pressure fieldby the capsid

圧力プロフィール
Pressure profile

Local pressure from virial theorem 240 excess negative charges at pH=7.4
N

P =

N

N kT
1

rij FijV
3V i j i
( ri V )1 N
Q ( ri ) =rij Eij3V j i
P( ri ) = qiQ ( ri )

Exchange of water
Cumulative number of transfer: from outside to inside: from inside to outside

Stroboscopic picture of a watermolecule crossing the capsid.

1. Water exchange is in thermal
equilibrium and is very rapid.

→ Chemical potential of water molecules is the same
between the inside and outside of the capsid.

2. The capsid doesn’t allow the ions to pass through.
→ The capsid works as a semipermeable membrane.

プロジェクトに育てられたMODYLAS•

NAREGIプロジェクト
自前で汎用プログラムを作ろう！（2004年秋）（自前コードがないと、開発した手法のテスト、最先端の手法の導入といった方法論に関する研究ができない。）
開発方針：多少遅くてもよい、高精度の計算ができなければならない。
まずは自分たちの研究用として利用。
数ノードで並列化MPI＋SMP（富士通SR8000）、数万原子系に対応。
⇒ループ分割、簡単。
スパコンを占有した研究のための実証計算。
分子研SR11000（2005年正月）⇒両親媒性分子集合体の熱力学的安定性についての
自由エネルギー計算。
ウイルス系（1000万原子）への挑戦⇒高速多重極展開法（FMM)実装（2007年春）。

プロジェクトに育てられたMODYLAS• Nano Grand Challenge プロジェクト
領域分割化。
並列化を考慮した変数の取り直し。プログラムの書き直し。
高並列化へと歩み出す。
開発方針：領域分割化し、高並列に対応。
• 京プロジェクト
変数のコピーレス化、変数を通信用、演算用に分けない。
TOFUによる隣接通信の活用。
SIMDの活用、オンキャッシュ化
開発方針：1000万原子系を1MDステップ10ms以下に。（実際には5msにまで高速化）
ソフト公開、ソースコード公開。
ミニアプリ化。
• ポスト京プロジェクトの要求、マシンの仕様に合わせて最適化してきた。
将来のユーザーのニーズ、マシンのアーキテクチャーを見据えながら、時代に合ったソフトウェア開発を進めなければならない。

今日のまとめ
ＭＤ計算の要素
各要素の高速化、並列化
ＭＤ計算とは• 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. １４１, １６５１０１ (２０１４).