学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

jh150015-NA11

分子動力学計算ソフトウェア MODYLAS のメニーコアアーキテクチャ対応並列化に関する研究
安藤 嘉倫（名古屋大学 工学研究科 計算科学連携教育研究センター）
概要	  本研究では、汎用分子動力学計算ソフトウェア MODYLAS について、次世代のメニーコアアーキテクチャ (FX100 および Xeon Phi) の性能を発揮させるための並列化チューニングを行った。本年度は特にスレッド並列計算について、まず問題点抽出のためのコード性能評価を行った。その上で MODYLAS を用いた分子動力学計算における２つのホットスポット (粒子対での相互作用計算 p2p および高速多重極展開法での多極子展開係数から局所展開係数への変換 M2L) でのスレッド間ロードンバランスの解消を目的として、スレッド並列性能向上のための並列アルゴリズム開発およびコーディングを行い、顕著な性能向上を実現した。

1.

共同研究に関する情報	 料や高分子材料 (図 1) を設計する際に MD 計(1) 共同研究を実施した拠点名	 

算により得られる知見が不可欠になりつつあ名古屋大学情報基盤センター	 

る. しかしながら, 長距離原子間相互作用を東京大学情報基盤センター	 

含めた実用の研究において MD 計算であつか	 える原子数および計算時間は,「京」コンピュ(2) 共同研究分野	 

ータといった最新鋭のスーパーコンピュータ



超大規模数値計算系応用分野	 

を用いたとしても 1 千万原子系 (空間サイズ	 として 50 ナノメートル立方程度) に対する数(3) 研究者の役割分担	 

100 ナノ秒の計算が限界である[1]. より大規模・全体統括	 安藤	 

かつ長時間な MD 計算を行うことで, 学問上・MPI 並列化性能統括	 荻野	 のブレークスルーだけでなく, より高精度な・分子動力学計算のアルゴリズム開発	 	 

材料設計が可能になると期待される.	  安藤、吉井、藤本、遠藤、篠田、岡崎	 ・並列化コーディング	 (MD 計算全般)	 	 	  安藤,	 藤本,	 遠藤	 ・プログラム性能評価	 大島、片桐	 ・Xeon	 Phi 向け並列化コーディング	 大島	 ・自動性能チューニング技術提供	 片桐	 ・並列化コーディングおよび並列アルゴリズム	 	  開発	 小村、鈴木	 
2.

研究の目的と意義	 研究の目的	  分子動力学 (MD) 計算は, 化学,物理, 生物, およびウィルス学といった様々な学問分野において実験とならぶ解析ツール

図 1	  MODYLAS を用いた計算対象の例。燃として広く普及している. 加えて工業分野に料電池に用いられるナフィオン高分子分離膜

おいても分子の特性を活かしたナノ機能性材(全原子モデル)。	 
1

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月	  本研究では「京」コンピュータおよび FX10に対応できると期待される. さらに大規模な上で稼働実績のある汎用分子動力学ソフトウ

ハイブリッド並列においてプロセス数とスレェア MODYLAS[2,3]に対して, 次世代のメニッド数の最適化は様々なチューニングパラメ

ーコア計算機, さらには将来のエクサスケーータの関与する複雑なプロセスであり, 本研

ルマシンの性能を発揮させるための並列チュ究成果をベースとした今後の研究において,ーニングを行った. 必要に応じてハードウェ

このコストを自動性能チューニング技術によアに適合したアルゴリズム開発についても行って削減することを図る. 近い将来、次世代

い, 将来的にエクサスケールマシンとスーパーコンピュータを高効率に動作させ,MODYLAS によって長距離原子間相互作用数億~10 億原子系での MD 計算を数 100 ナを含めた数億から 10 億原子系での実用的なノ秒~数μ秒のオーダーで行うことにより,MD 計算を可能にすることを目標とする.

上記学問分野におけるより現実的な MD シミ

研 究 の 意 義 	  分 子 動 力 学 ソ フ ト ウ ェ アュレーションにもとづく新規な科学的発見,

MODYLAS は「京」コンピュータの全ノードおよび工業分野における実用的な MD シミュ

規 模 で の 並 列 化 計 算 に 対 応 し た

レーションによる材料設計に貢献する.

MPI/OpenMP ハイブリッド並列化チューニングがすでに施されており, 次の目標は数 10

3. 当拠点公募型共同研究として実施した意義	  	 PFLOS (ポスト T2K) ないしエクサスケール	  MD 計算ソフトウェアがハードウェアの日マシンでの高効率な並列化計算への対応にあ進月歩の進化に対応するためには, 最新のアる. 一方スーパーコンピュータの進化はノーーキテクチャの基本性能および対応コーデド数の増加が頭打ちになり, ノード当たりのィングに詳しいコンピュータ・サイエンスのコア数および SIMD ベクトル長を増加させる研究者と, MD シミュレーションを用い実際ことでシステム全体の演算性能を向上させるの研究を実施している者とが, 互いに協力し方向にある. これら次世代のマシンでは, コた学際領域分野の研究が不可欠である. 本研

ア数および SIMD 長の増加によって現行のプ

究では両分野の専門家が共同研究者としてログラムの各所でスレッド並列性および参画している. 本課題に参加するコンピューSIMD 並列性が確保できなくなる。 本研究でタ・サイエンスの研究者はハイパフォーマンのスレッド並列性確保のためのチューニング,スコンピューティング(HPC)全般を専門と計算負荷の均等化によるスレッド並列化効率し, 特に自動性能チューニング(AT)に詳しいの向上, および階層キャッシュの最適化につ研究者（片桐）、最新のメニーコア技術に詳いての研究は, エクサスケールマシン上でしい研究者(大島),スレッド並列(小村,鈴木)

MODYLAS による高効率なスレッド並列化および MPI 並列化技術に詳しい研究者(荻MD 計算を実現するために不可欠である。

野) から構成される. 複合・階層的な最新の	  本研究では、FX100 をベースに 32 コアにメニーコア型の計算機システムを使いこな対応したスレッド並列化コーディングを実現す上で不可欠な人員構成である。一方, MDした上で, さらなるコア数増加に対応できるシミュレーションの研究者は, MD 計算ソフ

よう Xeon Phi をベースにした 50 コア以上のトウェア MODYLAS の基本設計に最初からスレッド並列化を見越した研究を行った。研

関与しかつ MD 計算の基本原理に詳しい研究成果は顕著であり, 今後研究をより深化さ究者(安藤,吉井,岡崎,篠田)および現行のせることでポスト T2K およびエクサスケー

MODYLAS 並列化の内容に詳しい研究者(安ルマシンでのメニーコア並列環境へスムーズ

藤,藤本,遠藤)から構成される. 両者の協業に
2学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

より, プログラミングレベルでの並列性能チの energy_direct 疑似コードを示す。この粒度ューニングだけでなく, 新規並列アルゴリズ不足を解消するため 2 つの観点からスレッドム開発・実装による性能向上を迅速に達成す並列の粒度向上および負荷均等化を試みた.	 

ることができた.	 

4. 前年度までに得られた研究成果の概要	 
本年度よりの開始課題のため該当せず。	 	 

5. 今年度の研究成果の詳細	 	  当初 4 月 1 日より利用可能予定であった名大情報基盤センターの FX100 の運用が 9 月に

図 2 高速多重極展開法 (FMM)の主要演算.

ずれ込んだため、スレッド並列の事前評価をホットスポットは近傍粒子との粒子対計算

まず東大情報基盤センターの FX10 上, およ

p2p [energy_direct] および多極子から局所び Xeon Phi 上で行った。

展開係数への変換 M2L [energy_fmm_2 およ	  MODYLAS では長距離静電相互作用を高び_3]. 演算はほか点電荷から多極子展開係速多重極展開法 (FMM, 図 2)を使って計算数への変換 p2M [calc_fmm], 多極子のマーしている。MD 計算においてのホットスポッ

ジ M2M [merge_fmm], 局所展開中心の移動

トは以下の二箇所,

L2L [energy_fmm_2 および_3], および局所	  ・Lennard-Jones 相互作用および静電相互展開係数をもちいた点電荷上の電場の計算	  	  作用の粒子対計算 p2p [energy_direct]

L2p [energy_fmm の一部] から構成される。	 	  ・多極子から局所展開係数への変換 M2L	 

[energy_fmm_2 および energy_fmm_3]

である。ここで[...] 内はサブルーチン名。この二箇所を本課題でのスレッド並列最適化の対象とした。
5.1 FX10 および FX100 での最適化	  まず energy_direct について, 基本的性能を調べるため FX10 での 16 スレッド実行時の演

図 3	  FX10 での energy_direct 16 スレッド実行

算性能およびスレッド間インバランス性能時のスレッド番号別演算性能。改良その 1 にの測定を行った。そのデータを元に問題点をより 30%のインバランスが 10%に削減された。

抽出しスレッド並列性能改善に取り組んだ。	  図 3 にあるように, energy_direct のオリジ
ナルコードでは 16 スレッド実行時にスレッド間に最大 30%のインバランスが生じていた。その原因として図 5 左にあるように, オリジナルコードでは対相互作用の自原子(iatom) に関するループが iatom の属する小領域 (サブセル) 内の原子数 (平均 40 程度)

図 4	  FX100 での energy_direct 32 スレッド実しかなく, 16 スレッド実行時に粒度不足に陥行時のスレッド番号別演算性能。改良その 2っていた。参考のため図 6 上にはオリジナルによりインバランスは 1%以下に削減された.
3学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月	  一つ目の観点は iatom についてのループの伸張である. 図 5 右にその概念図を, 図 6 下に具体的なコードを示す. 該当 do ループにて参照する小領域を z 軸方向に複数のサブセルにまたがるよう広げることで iatom に関するループ長を伸張した。コード改良の結果, 図 3にあるように 16 個のスレッド間でのロードインバランスが 10% に縮小され, かつループ長伸張の効果によりスレッド単体性能についても 25%程度向上した。この方法では相手

図 7	  energy_direct の疑似コード。改良コード

jatom 粒子についての参照範囲も拡大するた

その 2。

めキャッシュミス率の増大が心配されたものの, 現状深刻な問題とはなっていない。	  二つ目の観点からの改良では, スレッド並列対象箇所を iatom のループではなく iatom の属するサブセルに関するループに移動した.
その具体的なコードを図 7 に示す. この方法
では, 分割対象ループに nowait を指定しているためスレッド数の上限は [MPI プロセス内のサブセル数 Nsub]25 になる. Nsub=23 のとき200, 43 のとき 1600 であり, メニーコアに対し
図 5	  energy_direct コード改良の概念図. 中心て十分な並列度が確保される。図 4 にあるよの緑および赤い立方体が iatom の所属セル,うに, 32 コア実行時 (Nsub=43) でのスレッド端の赤および青い立方体が jatom の所属セル.

間ロードインバランスが 1%以下と劇的に縮小され, かつスレッドあたりの計算粒度が約スレッド数倍になった効果などによりスレッド単体性能についても向上した。その結果
Nsub=43 のとき 1.5 倍 (16core) および 1.7 倍(32core) の高速化が実現された. 我々はさらに性能を向上すべく i-square 単位, i-column 単位のスレッド並列化についても検討しており,問題規模ごとのこれら方法の最適解を探っている最中である.	  もう一方のホットスポット energy_fmm_2
についてはプロセス数に応じて M2L 演算をFMM 階層ごと異なる相互作用計算対象セル数 (すなわち, do ループ長さ) で行うという特徴がある。しかしながら、オリジナルコー

図 6	  energy_direct の疑似コード。オリジナル

ドではスレッド並列のチャンクサイズを階層コード(上), および改良コードその 1 (下)によらず一定かつ 8 スレッド実行に最適化した値としていた。そのため並列対象 do ループ4

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

粒度の大きい下層部についてはインバランスるユーザー定義関数 algndr (ルジャンドル倍が小さいが, 粒度の小さい上層部ではインバ多項式) での整数演算の効率化, および if 文ランスが大きくなる傾向があった。階層ごとの削減についても検討している。

チャンクサイズの最適化を行うことでこのインバランスを大幅に解消することができた(5.2 節で後述)。	  ほか、FX10 での基本プロファイル測定デー
タより明らかになった問題点として、スレッド間のバリア同期の回数の多さとそれに伴うバリア同期に要する時間の割合が大きい点があった。例えば図 8 にあるように, 主要ルーチンの p2p 演算 energy_direct で 12％, L2p 演

図8FX10 上での基本プロファイル測定結果。

算 energy_fmm で 56%, p2L 演算 calc_fmm でインプット pyp222 はタンパク質, イオン, お55%のバリア同期があり, プログラム全体でよび水を含む約 15 万原子系。8 プロセス 16も 12%と非常に大きかった。バリア同期が多スレッド実行。赤色が性能上改善の余地のあい原因として, これらサブルーチンにおけるる箇所。

リダクション指示句を使った変数 (例えばポテンシャルエネルギー) の足し込みがあると考え, 対象変数をスレッド数の次元を持つ配列として変数の足し込みをループ外で行うよう書き換えた。さらに, 依存性を吟味した上で該当 do ループ末に nowait 指示句を入れる,
特に energy_direct についてはスレッド並列対

図 9

FX100 上での基本プロファイル測定結

象の do ループ位置を変える (図 7), などの対果（改良コードその 2）。インプットは図 10処を施した。その結果、図 9 にあるようにと同じ。8 プロセス 16 スレッド実行 (Nsub=43)。

energy_direct および energy_fmm におけるバリ青色がバリア同期削減に成功した箇所。

ア同期を大幅に減少させることができた。しかしながらプログラム全体のバリア同期の割5.2 Xeon Phi での最適化合は依然 9%と高く, さらなる性能向上には	  図 10 には energy_direct, energy_fmm_2, お

ここで行ったと同様の改良をプログラムの該よび energy_fmm_3 へ上記改良を加えた上で当箇所に広く適用する必要がある。の Xeon Phi 上 60, 120, 180 および 240 スレッ	  また詳細プロファイル測定データから, 主ド実行での性能測定結果を示す。さらに図 11

要ルーチン以外で L1 キャッシュの dm ミスヒには例として 240 スレッド実行時のスレッドット率が顕著なサブルーチンが残されている

番号別経過時間を示す.

ことがわかった。これらルーチンについてデ	  energy_direct オリジナルコードについてスータの事前ソートによる不連続性の解消, おレッド数の増加にともない経過時間が増加しよびデータ量のインバランスの調整による改た(図 10). これは図 11 上にあるように do ル

善を計画している。さらに MIPS 値測定からープ回転数の枯渇によるスレッド間負荷インサブルーチン energy_fmm, calc_fmm などで整バランスが発生したためである。今回の改良

数演算が目立っており, そこで使用されていその 1 を施すことによって 120 スレッドに至5学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月るまでの高速化が達成されたが, 180 スレッド以上では飽和した (図 10)。また FX100 上では改良その 2 が性能的に優位であったが,Xeon Phi 上では反対の傾向となった。予想に反する結果のためその原因については現在調査を進めている。	  一方 energy_fmm_2 および_3 について, オ
リジナルコードでは最下層および二階層目については均一なスレッド並列性が確保される
一方, より上層においてはチャンクサイズの不備によりインバランスが生じていた(データ未掲載)。今回チャンクサイズの最適化を行った結果, 図 11 下にあるように 3 階層目以降においても均一な並列性を確保できるようになった。ただし 240 スレッドすべてを活用で	 きておらず, チャンクサイズの計算式には改図 10 Xeon Phi で の energy_direct (p2p),善点がのこる。今後 Xeon Phi のコア数がさらenergy_fmm_2, energy_fmm_3 (M2L+L2L) のに増えた場合, energy_fmm_2 および_3 ともに測定値. energy_fmm については階層ごとの合並列性を保つ上で該当 do ループ長のさらな計値. スレッド数 60, 120, 180 および 240.

る伸張は不可欠である。しかしながら 1 階層ごと M2L を行う現在の実装では, 相互作用計算対象セルのループ長は数千程度に限られている。メニーコアをより有効に活用する改良の候補として, ループ融合により階層を跨いだループに書き換えるなどを考えている。

6.

今年度の進捗状況と今後の展望	 	  本年度は, 基礎評価とメニーコア対応コード作成を目的として 4 期に分け研究を進める計画であった。FX100 が 4 月 1 日より使用できないとのアクシデントはあったが, 全般に順調に当初計画を進めることができた。今年

図 11 	  Xeon Phi で の energy_direct ( 上 ),度与えられた計算機資源については,FX10

energy_fmm_2 および energy_fmm_3 (下) 測を追加分を含め 100%, FX100 についてもほぼ定値。240 スレッド実行。_3(4), _3(3), …とあ100%消費した。

る( )内は FMM の階層番号。0 が最下層。	  第一期 (4-6 月) および第二期 (7-9 月) において, スレッド並列の事前評価を FX10 お

二箇所 (energy_direct, energy_fmm_2 およびよび Xeon Phi 上で行い, 図 3, 10 および 11 に_3) についてのコード改良およびコードの高あるようなプログラムの傾向を把握した。	 

度化を行った。energy_direct については, まそののち, 第二期において, ホットスポット

ず図 5 に示した方法でのスレッド並列対象
6

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月

iatom ループ長の伸張により 16 スレッド実行

び Xeon Phi のワイド SIMD 幅を活用できるよ時のインバランス解消および演算性能の向うな性能チューニングを開始する。ワイド上を得た。energy_fmm_2 および_3 についてSIMD 幅への対応は, メニーコアへの対応とはチャンクサイズを階層ごと最適化するこ

ともに次世代スーパーコンピュータの性能とでインバランスを低減させた。

を引き出すために不可欠である。それととも	  後半の第三期および第四期では, FX100 おに, 改良されたコードを用いた FX100 上でのよび Xeon Phi 上においてさらなるスレッド大規模性能テスト, および実際のサイエンス並列性能向上に取り組んだ.

Energy_direct

研究への応用を行う。ウィルスや高分子などについて図 7 にあるセル単位スレッド並列化を題材に, 新規サイエンス開拓を目的としたを実装することで, 32 スレッド実行時のインFX100 を全ノード規模利用した原子数 1000

バランスを 1% 以下に抑えることに成功し万 ~ 1 億オーダーでの大規模 MD 計算を計た. 結果として energy_direct オリジナルコー画している。	 

ドにくらべ 1.5 倍 (16core) および 1.7 倍(32core) の高速化が実現された。さらに Xeon

7.

Phi において 120 スレッドまでの高速化が実(1) 学術論文	 現された（図 10）。本来優位のはずの改良そ該当無しの２の性能が芳しくない, チャンクサイズの研究成果リスト	 

(2) 国際会議プロシーディングス	 

計算式に改善の余地がある, などの問題点は該当無し

残るものの 100 スレッド以上での並列性を確(3) 国際会議発表	 

保したことは本研究の大きな成果である。な	  	  該当無し

お energy_direct については, 改良その 1 およ(4) 国内会議発表	 

び 2 の他にも並列化方法を複数提案し最適解・安藤嘉倫, 吉井範行, 藤本和士, 小嶋秀和,を探っている最中である。

山田篤志, 岩橋建輔, 水谷文保, 岡崎進, 高並	  採択済みの平成 28 年度 JHPCN 課題では、

列対応汎用分子動力学シミュレーションソフ

引き続き Xeon Phi での性能向上, および５章

ト MODYLAS による大規模分子動力学計算,に挙げた問題点 (ホットスポット以外でのスHPCS2015, 東京 (2015).

レッド間のバリア同期, キャッシュミス率,・安藤嘉倫, 吉井範行, 岡崎進, 汎用 MD ソ

整数演算の最適化など) の解決に取り組む。

フトウェア MODYLAS の異方的な MPI プロ

バリア同期削減については今回行ったと同セス分割および基本セル分割への拡張, 第 29

じ改良をプログラムの該当箇所に広く適用

回分子シミュレーション討論会, 新潟することでプログラム全体のさらなる性能(2015).

向上が期待できる。さらに開発した複数のコ・安藤嘉倫, 鈴木惣一朗, 大島聡史, オーガナードについて, 自動性能チューニング技術をイズドセッション「分子動力学計算ソフトウ用い, 問題規模(原子数), プロセス数, コア

ェア MODYLAS のメニーコアアーキテクチ数といった実行条件ごと最適なコードを探ャ対応並列化に関する研究」, HPCS2016, 仙索する。また FX100 で得られた性能測定デー

台 (2016).

タをもとにエクサスケールマシンでの性能(5) その他（特許，プレス発表，著書等）	 

予測をするための数値モデルの構築も計画

安藤嘉倫,“汎用分子動力学計算ソフトウェしている。	 

ア MODYLAS”, 日本機械学会 計算力学部門	  さらに次年度 JHPCN 課題では FX100 およ

CMD Newsletter, No.54, 28-30 (2015).
7

学際大規模情報基盤共同利用・共同研究拠点 平成 27 年度共同研究 最終報告書 2016 年 5 月	 
参考文献[1] Y.Andoh, N.Yoshii, et al., “All-atom molecular
dynamics calculation study of entire poliovirus empty
capsides in solution”, J. Chem. Phys., 141, 165101(2014).
[2] Y.Andoh, N.Yoshii, et al., “MODYLAS: A highly
parallelized

general-purpose

moleculardynamics

simulation program for large-scale systems withlong-range forces calculated by fast multipole method(FMM) and highly scalable fine-grained new parallel
processing algorithms”, J. Chem. Theory Compt., 9,3201-3209 (2013).
[3] N. Yoshii, Y.Andoh, et al., “MODYLAS: A highly
parallelized

general-purpose

moleculardynamics

simulation program”, Int. J. Quantum Chem., 115,342-348 (2015).	 	 

࣭┠Ᏻ࡜ࡋ࡚ 8 ࣮࣌ࢪ⛬ᗘ࡛ࡲ࡜ࡵ࡚ୗ
ࡉ࠸࣭ᴫせ㸪࠾ࡼࡧ㸪1.㹼5.ࡢぢฟࡋ㡯┠࡟ᚑࡗ࡚グ㏙ࡋ࡚ୗࡉ࠸㸬⾜㛫ࡸᩥᏐࢧ࢖ࢬ➼ࡢㄪᩚࡣ⮬⏤࡛ࡍ㸬࣭ᚲせ࡛࠶ࢀࡤ㸪ཧ⪃ᩥ⊩ࡣヱᙜࡍࡿ㡯ࡢ୰࡟グ㍕ࡋ࡚ୗࡉ࠸㸬࣭୰㛫ሗ࿌᭩ࡢෆᐜ࠿ࡽࡢ㐍ᤖࡀศ࠿ࡿࡼ࠺㸪᭱⤊ሗ࿌᭩ࢆసᡂࡋ࡚ୗࡉᙜሗ࿌᭩ࡣ㸪2014࠸㸬ᖺ 5 ᭶୰࡟ Web ୖ୍࡛⯡࡬බ㛤㛤ጞࡋࡲࡍ㸬බ㛤᪥࡟ࡘ࠸࡚ࡣ㸪ู㏵ࡈ㐃⤡ࡋࡲࡍ㸬Ⓨ⾲ෆᐜ࡟≉チฟ㢪ணᐃࡢⓎ᫂ࡀྵࡲࢀࡿሙྜ࡟ࡣ㸪Web බ㛤᪥ࡲ࡛࡟ᡤᒓᶵ㛵ෆ࡛ࡢᡭ⥆ཬࡧ≉チฟ㢪ࡀ᏶஢࡛ࡁࡿࡼ࠺㸪ྍ⬟࡞㝈ࡾవ⿱ࢆࡶࡗ࡚Ⓨ᫂ᒆฟࡉࢀࡿࡼ࠺࠾㢪࠸࠸ࡓࡋࡲࡍ㸬

目
安として
1~
2
節は
2
ペー
ジ
以内，
全体
で
6
ペー
ジ
程

8