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

課題番号 jh140029-NA15

並列フラグメント分子軌道計算プログラム OpenFMO のマルチプラットフォーム化
渡邊寿雄 （東京工業大学）

概要
京コンピュータなどの SMP クラスタ型超並列計算機に最適化された並列フラグメント分子軌道（FMO）プログラム OpenFMO をベースとして，GPGPU や Intel Xeon Phi などを搭載したアクセラレータ型スパコンへのマルチプラットフォーム化を進めた．マルチプラットフォーム化においては，JHPCN の枠組みを活用し，複数拠点の多様なアクセラレータ型スパコンを利用した．Hartree-Fock 計算ルーチンの GPGPU 化と OpenFMO への組み込みが行われ，FMO 計算での高速化を確認した．
またアクセラレータをより効率的に活用するために，分子軌道計算のカーネルコード（性能ボトルネックである 2 電子反撥積分）の簡素化を検討した．簡素化により高い SIMD 演算率が達成された一方で，計算量が増加したため，実行時間としては既存プログラムよりわずかに増加する結果となった．

1. 共同研究に関する情報

筑波大学・計算科学研究センター(1) 共同研究を実施した拠点名

梅田 宏明

東京工業大学 学術国際情報センター

産業技術総合研究所・ナノシステム研究部門

京都大学 学術情報メディアセンター

長嶋 雲兵(2) 共同研究分野

コンピュータシミュレーションによって超大規模数値計算系応用分野

原子や電子を露わに扱い化学反応を解明・予(3) 参加研究者の役割分担

測する計算化学は，2013 年ノーベル化学賞

東京工業大学・学術国際情報センター
渡邊 寿雄
青木 尊之

簡素化コード作成支援

2. 研究の目的と意義

九州大学 情報基盤研究開発センター

GPU 化/性能評価を受賞し，方法論やプログラムの開発から課題とりまとめ様々な応用まで幅広く研究・開発が行われて簡素化コードの作成/評価いる．その中でもフラグメント分子軌道(FMO)法は，たんぱく質などの巨大分子に対GPU 化の支援

京都大学・学術情報メディアセンター

する第一原理電子状態計算を可能にする計中島 浩Intel Phi 化の支援

算手法として注目されている．いくつかの

河合 直聡

Intel Phi 化/性能評価

FMO プログラムではスパコンへ向けた最適

九州大学・情報基盤研究開発センター化も行われており，稲富(九州大学)らにより

青柳 睦

簡素化コード作成支援

開発された FMO プログラム OpenFMO は，本田 宏明

Intel Phi 化/性能評価

京をはじめとした SMP クラスタ型並列計算

簡素化コード作成支援機で効率的な超並列実行が可能なプログラ

九州大学・大学院システム情報科学研究院
稲富 雄一

ムである．

OpenFMO のオリジナル

昨今，様々なアーキテクチャのスパコンがコードの改良登場し，特に省電力性能に優れたアクセラレ
1

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

ータ型スパコンはエクサへ向けて重要な役

ろう．

割を果たすことが期待されている．しかしな

本研究課題の目的は，並列フラグメント分がら，同じ計算化学分野のプログラムである子軌道計算プログラム OpenFMO のカーネ

分子動力学法がアクセラレータ化に素早くルコードの簡素化と，そのアクセラレータ化対応したのとは対照的に，分子軌道法のプロにより，実用的な FMO プログラムの開発をグラム群のアクセラレータ化はその複雑な行うことである．

カーネルコードが原因で遅れている．そのカ

3. 当拠点公募型共同研究として実施した意義

ーネルコードの複雑さを解消すれば，既存の本研究課題の目的の 1 つであるアクセラレ

GPU や Intel Xeon Phi などへのアクセラレータ化を行うには，多用なアクセラレータ型ータ化が容易になるのみならず，今後新たなスパコンを利用できる環境が不可欠であっアーキテクチャのアクセラレータが登場した．当拠点公募型共同研究においては，複数た際にも他分野のアプリケーションに後れ拠点を利用することで，各拠点が持つ多様なを取ることなく，アクセラレータ化が可能にアクセラレータ型スパコンの利用が可能でなる．

あり，本研究課題に必要な環境が揃っていた．

カーネルコードの複雑さを解消する一つ

また，各拠点の研究者の協力により，当初かのアイディアとして，電子状態の基底関数展ら予定していた GPU 化のみならず，Intel開を最も単純な関数(s 型 Gauss 関数)にて再Xeon Phi 化も実施可能になった．

展開する方法がある．このアイディアは非常

また，本研究課題のもう 1 つの目的であるに古く(1956 年 H. Preuss らの Gaussian

OpenFMO のカーネルコードの簡素化のメ

Lobe 法など)，アルゴリズムやカーネルコーリットは，これまでアクセラレータ化が遅れドの簡素化とそれに伴うプログラムの高速てきた分子軌道法のプログラム群のアクセ化が可能である．この方法に於いては，計算ラレータ化が容易になることに加えて，計算量の増大が起こるため，現在の主要プログラ化学者と計算機科学者によるプログラムのムでは採用されていない．この計算量の増大共同開発を容易にすることも大きなメリッを克服するには，簡素化したカーネルコード

トである．当拠点公募型共同研究においては，の高速化が必要であり，今後さらに SIMDこのような学際的な共同研究を推奨してお（Single Instruction / Multiple Data）長がり，本研究課題での共同研究の実施に適して増大するアクセラレータの活用が有効であいた．

ろう．分子軌道計算の性能のボトルネックはまた本研究課題で開発したプログラム群

2 電子反撥積分であるが，簡素化したカーネは公開予定であり，OpenFMO プログラムとルコードではそのループ長が揃うため，して直接利用されることに加えて，マルチプSIMD 演算の有効利用が可能になる．

ラットフォーム化されたカーネルコードがカーネルコードの複雑さを解消する他の計算化学分野の他のプログラムのマルチプ方法としては，基底関数展開に用いる関数をラットフォーム化にも大きな波及効果を及縮約された基底関数ではなく原始基底関数に基づいて取り扱う方法がある．これにより，カーネルコードの最内ループが容易にアンぼすことが期待される．

4. 前年度までに得られた研究成果の概要

ローリング可能となると共にループボディ新規課題のため，記入不要部にデータ依存性がなくなることによりSIMD 演算の効率的利用が可能になるであ

5. 今年度の研究成果の詳細
2

関数値

関数値

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

x 軸座標
図 1a p 型関数の s 型関数による再展開(ア)OpenFMO のカーネルコードの簡素化１x 軸座標
図 1b px 型関数の s 型関数による再展開における軌道中心位置と誤差: 軌道指数: 規格化定数【概要と目的】
カーネルコードが複雑である原因の 1 つは，

d : 軌道中心のずれ

HF 法や DFT 法におけるホットスポットである二電子反撥積分計算ルーチンがガウス型

px 型基底関数は共通の軌道指数αと中心位基底関数の４つ組の積分タイプ毎に必要と置のずれが d である２つの s 型基底関数の符なる点にある．その複雑さを解消するために，号が異なる線形結合で近似的に表される．軌
道指数αは中心位置のずれ d に依存し，中心最も単純な s 型基底関数による基底関数の再位置のずれ d が 0 に近づくと軌道指数αは無展開による OpenFMO のカーネルコードの限大になり，d が 0 の極限では px 型基底関数簡素化を行った．カーネルコードの簡素化にに一致，つまり再展開時の近似誤差がゼロにおいては混合精度計算の検討も不可欠であなる．一方で，d が増えるに従い，元々の px

り，それに伴い，下記 4 点の検討が必要であ型基底関数からのずれが大きくなる．

る．

一方で d が減少することで px 型基底関数を(1) 再展開時の近似誤差による精度低下

構成する 2 つの s 型基底関数の重なりは大き(2) 再展開後の基底関数の線形従属性によるくなり，結果として線形従属性による積分精度の低下が起きる．この積分精度の低下によ

二電子反撥積分の精度低下(3) 再展開後の基底関数の球面調和性の低下って SCF 計算が収束しない現象は以前より報(4) 再展開後の SIMD 演算活用による高速化告されている．そこで px 型基底関数を s 型基底関数で近似する際の一定の誤差を許容し，今年度は，(1)～(3)による精度低下について積分精度を確保する固定値 d を採用すること検討を行った．

で，SCF 計算が収束することを確認した．
(3)の再展開後の基底関数の球面調和性の【結果】低下については，再展開時の s 型基底関数のs 型基底関数で px 型基底関数を再展開する張る空間の等方性に着目し，3 つの p 型基底場合の式を以下に示した．

関数を 6 つの s 型基底関数にて再展開する際に球面調和性の低下が最も大きいことを確認した．球面調和性の確保には，再展開時のs 型基底関数が張る空間の等方性が不可欠な3

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

ため，s 型基底関数を追加して展開することるカーネルコードの更なる簡略化により，で球面調和性の担保を行うことを検討中でより効率的な SIMD 演算の実現を目的としある．

た．
【効率的な SIMD 演算に向けたコード変更】

OpenFMO のカーネルコードの簡素化２(イ)

一般に，コンパイラの自動並列化機能を【概要と目的】利用したコード中の多重ループ箇所の SIMD

今後開発されるメニーコアプロセッサに化オブジェクトコードの生成には，最内ル於いては，搭載される SIMD 演算器の SIMDープが容易にアンローリング可能であるこ長が増大すると予想されるため，その特性とならびにループボディ部にデータ依存性を活かして効率的に計算するためにはが少ないことが必要条件となる．しかしな

SIMD 演算の割合を高めることが必要であがら，一般的な G 行列計算コードは 8 重ルる．しかし，現在の分子軌道法計算プログープで構成され，下位 4 重ループ長が上位ラムでは，そのボトルネックとなる二電子4 重ループの変数に依存して 1～6 程度のFock 行列(G 行列) 計算において効率的な範囲で変化し，それが実行時に決まる．まSIMD 演算が難しいことが知られており，現たループボディ部では条件分岐を含むデー状では効率的なメニーコアプロセッサ利用タ依存関係のある複雑な計算が行われていが困難である．

る．

そこで本研究では前項にて示した通り，sこれらの問題に対し本研究では，これま型基底関数による基底関数の再展開によるでの縮約基底ではなく原始基底に基づく方カーネルコードの簡素化によってそのルー法を試み SIMD 演算の効率的利用を図る．こプ長を揃えつつ，SIMD 演算器の WAY 数にの際，下位縮約ループをアンローリングし対して十分なデータ並列性を得ることで，上位シェルループへの組み込みを行うといSIMD 演算の有効利用を目指している．

うアプローチをとる．さらに，ループボデ

しかしながら，効率的な SIMD 演算を可能

ィのコードに対し配列表記や専用ディレクにするには，それだけでは不十分であり，ティブを利用してコード変更を行うことで，条件分岐がないことやメモリ中の連続領域より明確にコンパイラへ SIMD 演算箇所をへの被演算データの配置など更なる SIMD 演伝える工夫をする．

算の阻害要因の排除が必要不可欠である．
そのため，本研究では，基底関数展開に【性能評価】

用いる関数を縮約された基底関数ではなく本研究では，九大の稲富らにより開発さ原始基底関数に基づいて取り扱うことによれているフラグメント分子軌道法プログラ

Original

Loop Coalescing

50

Loop Coalescing + Array Notation

40
クロックサイクル数

SIMD Ratio [%]100806040200GlySerGlu

Phe

Trp

図 2 各アミノ酸における SIMD 演算率30

Original
Loop Coalescing
Loop Coalescing + Array Notation

20100GlySerGlu

Phe

Trp

図 3 各アミノ酸でのクロックサイクル数4学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月

ム OpenFMO のカーネルコードを使用した．

は確実である．今後メモリ上のデータ配置

既存の RHF 法プログラムに対する評価の方法の考慮など SIMD 化コードの最適化も同ため，アミノ酸５種(Gly，Ser，Glu，Phe，時に行うことにより，更なる高速化を検討

Trp) に対し上記のプログラムを用いた 1予定である．

回の繰り返しのみの RHF 法テスト計算を行(ウ)った．計算には Intel Xeon E5-2620(6 core)【概要と実装】

が 2 CPU の PC を用い，2 ならびに 4 WAY の筑波大学の梅田らを中心として，HF 計算コ倍精度浮動小数点演算が可能である．

ードの GPGPU 化を行っている．GPGPU 化を行Intel Compiler14.0.0 の -O3 最 適 化 で のった Fock 行列計算コードは OpenFMO におい

PAPI ライブラリによるハードウェアカウてフラグメント(及びフラグメントペア)のンタの取得結果から，G 行列計算における計算で利用される HF 計算コードをスケルト

SIMD 演算率は 5%以下であった．

ンコードとして切り出したものである．C 言

次に多くの積分タイプの中から(ss|ss)積語で書かれたスケルトンコードは簡潔であ分 に 関 わ る G 行 列 計 算 に 注 目 し ， Loop

り機能追加や改良も容易となっている．またCoalescing＋ループボディ SIMD 化した場合このコードへの改良点を OpenFMO プログラムの SIMD 演算率と計算に要したクロックサイに導入することも比較的容易であり，これにクル数をそれぞれ図 2，図 3 に示した．図 2より GPGPU 化した Fock 行列計算コードを FMOが示す通り，ループボディ部の変更に従っ計算に適用することが可能となる．て SIMD 演算率は全分子で 96%超を達成する

多くの計算コアを持つ GPU では行列への加ことができた．しかしながら，クロックサ算にコストの高い排他制御が必要となるたイ ク ル 数 を 見 る と ， オ リ ジ ナ ル < Loop

め，これを避けるようなアルゴリズムを開発

Coalescing < Loop Coalescing＋SIMD 化 とした．また GPU スレッドの SIMD 動作が効果

実行時間は最大で 1.6 倍に増加している．

的に実行可能となるような最適化や，CPU と

これは Loop Coalescing＋SIMD 化により浮GPU の同時実行についても実装している．さ

動小数点演算以外の部分によるオーバーヘ

らに MPI による並列化と合わせることで複数ッドのためである．

GPU を利用した高速化も実現した．

今回用いた手法で達成した高い SIMD 演算

組み込みに利用したバージョンの OpenFMO率が示す通り，現在ある 8 ウェイの SIMD 演

プログラムは動的プロセス生成を利用する算器を持つプロセッサだけでなく，今後開MIMD プログラムとなっており，マスタ MPI発が予想される 16 ウェイや 32 ウェイのプロセスとメモリサーバ MPI プログラム，そSIMD 演算器を持つプロセッサなど，より長して多数のワーカ MPI プログラムからなって

い SIMD 長を持つプロセッサを利用した場合いる．我々のコードを実装するのは，このうでも効率良く SIMD 演算器を活用できること
表1

OpenFMO のアクセラレータ化（GPU 化）

GPGPU 化 Fock 行列計算コードを用いた FMO 計算の実行時間(秒)Algorithm
CPU (direct)
CPU (in-core)
GPU+CPU (direct)
GPU+CPU (in-core)

SCC
219
212
196
195(Ala)10
DimerSCF
165
112
10685

5

Total
392
333
320
299

SCC
1,418
1,518
1,342
1,391

CrambinDimerSCF
1,606
1,446
1,269
1,229

Total
3,115
3,069
2,746
2,749

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

ちでワーカプログラムだけである．GPGPU 化

で表中の実行時間は Fock 行列計算が含まれコードの組み込みは 1)GPU の初期化，2)分子る二つの計算手順(SCC 計算部分，dimerSCFデータの転送，3)Fock 行列計算ルーチンの差計算部分)，および FMO 計算全体に対するマし替え，4)分子データの消去，5)GPU の終了

スタプロセスにおける経過時間である．我々処理の各コードを OpenFMO ワーカコードに挿が実装した GPGPU 化 Fock 行列計算コードは入することで実装できる．

CPU による direct 計算より速いだけでなく，in-core 計算手法を用いた場合よりも高速に【ベンチマーク】動作していることがわかる．

OpenFMO に実装した GPGPU 化 Fock 行列計算In-core 計算手法よりも高速ではあるが，コードの性能評価は筑波大学の HA-PACS ベーGPU による高速化は計算全体で 1 割程度でしスクラスタを用いて行った．HA-PACS ベースかない．今回の実装では Fock 行列計算部分クラスタの計算ノードには 2 台の 8 コアのみを GPGPU 化しているため，FMO 計算におIntel E5 CPU(Sandy Bridge-EP, 2.6GHz)と 4けるもう一つのボトルネックである環境ポ

台の Fermi 世代の GPGPU(NVIDIA M2090 GPU)，テンシャル計算部分が高速化されていないおよび 128GB のメモリが搭載されており，そことが原因である．現在，Fock 行列計算と同れらが InfiniBand により接続されている．

様に二電子積分を利用する 4 中心項部分につ複数 GPU を活用するためノードごとに 4MPI

いての GPGPU 化を検討中であり，これによりプロセスを起動し，それぞれのプロセスがさらに大きな高速化が期待できる．

OpenMP 並列で 4 CPU コアと 1 台の GPU を利用することとした．コンパイルや実行には(エ)OpenFMO のアクセラレータ化（Intel Xeon

Intel コ ン パ イ ラ 14.0, CUDA 5.0.35,Phi 化）

mvapich2 1.8.1 を そ れ ぞ れ 利 用 し た ．

OpenFMO の Intel Phi 化については，京都

OpenFMO では動的プロセス生成や片側通信を大学 学術情報メディアセンターが JHPCN 採利用するため，MPI 処理系である mvapich2択課題向けに募集し採択された大規模計算については適切に環境変数を設定している．

プログラムの高度化・高性能化支援制度を利用して実施した．

性能評価としてグリシンの 10 量体(112 原子, 5 フラグメント)及び Crambin タンパクOpenFMO の最適化は，まず Xeon E5-2697 v3(642 原 子 , 20 フ ラ グ メ ン ト ) の(HSW 2.6GHz 14 コア）の２ソケットシステム

FMO-HF/6-31G(d)計算を行った．それぞれ計上で最適化し，次に同システム上に実装した

算には HA-PACS ベースクラスタの 2 ノード

Xeon Phi 7120A (1.238 GHz 61 コア) 上での(8MPI プロセス)および 8 ノード(32MPI プロ単独実行（native）向けに最適化した．

セス)を用いている．HA-PACS クラスタではノHSW 向けの最適化には，インテルの最新コードあたり 128GB と大きなメモリ空間を利用ンパイラ(15.0.2,164)を用いて OpenMP で並可能なため，FMO 計算のフラグメントサイズ列化し，HAVX2 命令を使用するようにコンパ程度であれば二電子積分をメモリに保存しイルした．Phi 向けには，同コンパイラにてておいて何度も利用する in-core 手法が可能

OpenMP で並列化し，512 ビットの SIMD 命令になる．そこで性能評価では CPU や GPU によを使用するようにコンパイルした．

る direct 計算だけでなく，この in-core 計

最も時間を要する関数のループをベクト算手法も比較の対象とした．表 1 にはそれぞル化できるように書き換え，内部関数でベクれの計算手法による実行時間を示した．ここ

トル化を行った結果，HSW 上で 3.3%、Phi 上
6

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

7. 研究成果リスト

で 20.1%の性能向上が得られた．

(1) 学術論文

Phi 上での実行では，HSW と同様に誤差関数テーブルがメモリ上に乗るようにできな

1. 梅田宏明, 塙敏博, 庄司光男, 朴泰祐,かったため，HSW の 11 分の 1 程度の性能とな

重田育照, "GPGPU クラスタ上での FMO 計

った．テーブルサイズは展開項数と要求精度

算 の 性 能 評 価 ", Journal of Computerに依存し，計算量を増やす代わりにテーブル

Chemistry, Japan, 13(6), 323-324(2014)(2) 国際会議プロシーディングス

サイズを小さくすることが可能なため，メモなしリ上に乗るように最適化する必要がある．

(3) 国際会議発表

HSW と Phi においては同じ条件での比較ができなかったため，バッファーメモリサイズ

1. H. Umeda, "Fock matrix preparation

を同じ 32MB にして比較した場合でも Phi 上with GPGPU for fragment molecular

での性能は HSW の半分弱であった．Knights

orbitalCorner には誤差関数の演算をベクトル化す

connection workshop @SC14, 招待講演,る機能があるため，その活用による更なる高2014.11.19, ニューオリンズ(アメリカ)(FMO)

calculation",HPC(4) 国内会議発表

速化を今後検討している．

1. 渡邊 寿雄，稲富 雄一，梅田 宏明，本田

6. 今年度の進捗状況と今後の展望

宏明，青柳 睦，長嶋 雲兵「ガウス型基カーネルコードの簡素化に関しては，当初底関数の s 型関数による再展開についての実施計画にあった s 型基底関数による基底の考察」
第 8 回分子科学総合討論会 2014，関数の再展開による簡素化に加えて，基底関

ポスター発表 4P128, 2014.9.21-24，広島

数の縮約を解いた原始基底関数の露わな利

大学 東広島キャンパス用による簡素化にも取り組んだ．現時点では，簡素化されたカーネルコードにおける演算

2. 梅田 宏明，塙 敏博，庄司 光男，朴 泰祐，重田 育照「GPGPU 化 Fock 行列計算ル数の増加に SIMD 化による高速化が追い付いておらず，実行時間の増加を招いているため，更なる最適化の検討が必要である．

ーチンの OpenFMO への組み込み」第 8 回
分子科学総合討論会 2014，ポスター発表
4P119, 2014.9.21-24，広島大学 東広島キ

カーネルコードの GPU 化と，GPU 化されたャンパス

カーネルコードの取り込みによる OpenFMO プ

3. 渡邊 寿雄「ガウス型基底関数の s 型関数ログラムの GPU 化については計画通りに達成による再展開についての考察」第 2 回できた．OpenFMO プログラムの GPU 化につい

CUTE シンポジウム・コンピュータ化ては，カーネルコード以外の環境ポテンシャ

学，招待講演, 2014.10.30-31，三重大学

ル計算部分の GPU 化による高速化を検討中で地域イノベーション研究開発拠点あり，さらに大きな高速化が期待できる．

4. 梅田 宏明「FMO 計算のための FockOpenFMO の Intel Phi 化については，既存行列生成ルーチンの GPGPU 化」第 2

コードを Phi 上での単独実行（native）向け

回 CUTE シンポジウム・コンピュータに最適化を行った．しかしながら，性能的に化学，招待講演, 2014.10.30-31，三重大は満足いくものではなく，今後は誤差関数テ

学 地域イノベーション研究開発拠点ーブルサイズの最適化や，誤差関数の演算を

5. 今村憲亮，本田宏明，稲富雄一，南里豪

ベクトル化するなど更なる高速化を試みる

志「二電子 Fock 行列の原始基底計算アル必要がある．

ゴリズムの提案」，日本コンピュータ化学7学際大規模情報基盤共同利用・共同研究拠点 平成 26 年度共同研究 最終報告書 2015 年 5 月会 2015 春季年会，ポスター発表 1P14,2015.5.28，東京工業大学 大岡山キャンパス
6. 梅 田 宏 明 , "Development of FMO
program

forrecentHPC

systems:

K-computer and GPGPU cluster", ポス
ト京向けアクセラレータについての勉強会, 招待講演, 2014.4.17, 岡崎
7. 梅田宏明, "フラグメント分子軌道法における Fock 行列計算の GPGPU 化", GTC
Japan 2014, 招待講演, 2014.7.16, 東京
8. 梅 田 宏 明 , 塙 敏 博 , 庄 司 光 男 , 朴 泰 祐 ,"GPGPU 化した Fock 行列計算ルーチンの OpenFMO への組み込み",日本コンピ
ュータ化学会 2014 春季年会, ポスター発表, 2014.5.30, 東京
9. 梅 田 宏 明 , 塙 敏 博 , 庄 司 光 男 , 朴 泰 祐 ,"GPGPU クラスタ上での FMO 計算の性能評価", 日本コンピュータ化学会 2014
秋季年会, ポスター発表, 2014.10.19, 郡山(5) その他（特許，プレス発表，著書等）
なし

8