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

jh180069-NAH

GW space-time コードの大規模な有機-金属界面への適用に向けた高効率化
柳澤将（琉球大学）

概要 GW space-time コードによる大規模な表面・界面系の効率的な計算を可能にする
ため、今年度は、主に複素誘電行列の、LU 分解による逆行列化計算のノード内タスク並列化を検討した。結果、スレッド数の増加に伴う並列化の効果が認められ、行列サイズの増加に応じてその効果が大きくなることも認められた。表面などの二次元系の電子状態に向けた計算法についても、非周期的な方向でクーロンポテンシャルを打ち切る方法の実装を進めており、表面・界面系への適用が可能になってきている。これまでの共同研究によるプログラムの改善により、表面・界面系への適用も含め、以前よりも適用可能な問題・サイズが増え、応用計算の成果も出てきている。

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

模な計算資源を活用し有機電子材料系に適用し(1) 共同研究を実施した拠点名た例や、表面・界面系に適用した例は少ない。本

東北大学サイバーサイエンスセンター

研究計画の達成により、その状況を改善し、大規模な GW 近似の計算による国内の大型計算機資(2) 共同研究分野
n

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

o

超大規模データ処理系応用分野

GW 近似は物質中のグリーン関数に基づく方o

超大容量ネットワーク技術分野

法で、一体有効ポテンシャルではなく、電子間相

o

超大規模情報システム関連研究分野

互作用をより精度良く考慮した電子準位（純粒源の使用を促すことにつなげたいと考える。

(3) 参加研究者の役割分担

子エネルギー）を計算し、半導体のバンドギャッ・柳澤 将（琉球大学・代表）プログラム開発・テプなどの予測の信頼性が高い。しかし、一般に、ストおよび、実物質系の大規模計算の実行。

空間中の任意の二点間の非局所的な多体電子間・江川 隆輔（東北大学・副代表）プログラム効率相互作用の記述が必要で、バンド計算で頻出の化のための実装および助言。

密度汎関数理論（DFT）の計算ほどプログラムコードの効率の良い並列化が容易ではなく、計算2. 研究の目的と意義に時間がかかる。本研究計画の達成で、GW 近似

本研究では、有機電子材料と関連のある物質をはじめとする、DFT 計算よりも高度な第一原系の基礎的な電子物性の理解に資するため、量理計算法の大規模計算に向けた実装を促すこと

子力学や電磁気学などをもとに物質の電子物性にもつながると考えられる。

を予測・再現する第一原理電子状態計算のプロ

また、有機半導体の電子物性の実験研究は日グラムの整備・開発を進める。特に、有機分子膜

本では歴史と実績があり、現在でも分野を先導と金属電極との接触する領域（界面）でのキャリする最先端の成果が得られている。その一方、実ア（電子・正孔）注入障壁の正確な予測に向けた

験的観測の難しい電子状態に関する理論的解釈第一原理計算（GW 近似）を可能にするため、大でも先端の成果を得て、総合的に有機半導体の規模な並列計算を効率的に実行することを目的

研究分野を国際的にリードする体制の構築にもに、MPI 並列と OpenMP とのハイブリッド並列つなげることが必要である。本研究では、コードの実装を行ってプログラムの効率向上を進める。

開発の進展とともに、大規模な計算機資源の利GW 近似の計算は最近、国内外の物性の理論用により、実験事実の理論計算による洞察・解釈研究で広く実行される一方、日本国内では、大規の深化や、新規の成果も目指す。
1

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

このような目的・動機の下、昨年度から GWている。また、物性の計算科学コミュニティにも

近 似 の 計 算 プ ロ グ ラ ム コ ー ド GW space-time

成果を広めることで、同センターでの開発プロ（GWST）をもとに研究を開始した。昨年度は、グラムの使用を広めることができる点も意義が複素誘電行列の逆行列化に使用する ScaLapackあると言える。のライブラリの使用をやめ、Lapack のライブラ

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

リを使用してノードごとに独立して逆行列化を行うようコードを改良し、一定の成果を得た。

固体・結晶系の準粒子エネルギーを正確に予

本年度はその成果に続き、ライブラリの逆行測できる GW 近似の計算プログラム GWST を、列化コードをノード内でタスク並列実行するこ

有機半導体系の表面・界面の問題にも適用できとから検討を開始し、それ以外の箇所のノード

る よ う 、 MPI 並 列 と ノ ー ド 内 タ ス ク 並 列

内タスク並列化についても検証・実装を進めて（OpenMP）とのハイブリッド並列の導入によっいる。本報告書では、逆行列化コードとしてて、並列化効率の向上を試みた。プログラム中のLapack のライブラリの LU 分解を用いた逆行列複素誘電行列の逆行列化計算において、化コードでの、ノード内タスク並列のテスト結

ScaLapack のライブラリの使用をやめ、並列度が果を報告する。

増すにつれて通信時間の増大で発生するボトルまた、GWST プログラムによる比較的規模のネックの解消を進めた。計算データは、ノード・大きい応用として、GWST コードの適用目標でプロセス間では独立な変数についてノード間にある表面系について必要な実装・テストを開始

分散させ、さらにノード内でタスク並列化を行したので、その成果も報告する。

うようにした。結果、現在テスト可能な 224 程度の並列度で、計算時間の短縮と並列化効率の改3.

当拠点公募型共同研究として実施した意義善が見られた [29 年度 JHPCN 最終報告書および2016 年度まで、東北大学サイバーサイエンス

29 年度 JHPCN 最終報告会]。

センター（東北大学 CC）と共同研究を進め、GW

表面の周期的スラブの GW 近似による取り扱近似プログラムの整備・開発や、それを用いた有いは、まだ実行できてないが、表面での有機分子機半導体の電子物性計算で少なからず成果を収の配向に依存する静電ポテンシャルの効果と、めてきた。特に、通常の PC クラスタよりも 1 コ

有機半導体の結晶内部での注入電荷への誘起分アあたりのメモリの多い SX-ACE を使用するこ極の効果とを切り分けて扱う工夫を行なった。とで、比較的大きなサイズの問題にもプログラ

そのような計算によって、有機半導体の表面のムを適用できる利点がある。また、平面波基底関構造に依存した電荷注入準位の正確な決定が可数に基づく第一原理計算コードであるため、ベ能となることが分かった[S. Yanagisawa, AIP Conf.

クトル化で得られる利点も大きい。

Procs. 1906, 030014-1-4 (2017); K. Yamada, S.

昨年度より、東北大学 CC からプログラム開Yanagisawa et al., Phys. Rev. B 97, 245206-1-8(2018)]。

発の支援を受けながら、効率のよい大規模な並列計算を可能にし、計算対象をバルク結晶だけでなく表面・界面などの並進対称性のない系に5. 今年度の研究成果の詳細まで拡大するため、計算機資源の使用料金で大A. OpenMP 導入による、MPI 並列とのハイブリきなメリットのある貴拠点公募型共同研究としッド並列の導入: LU 分解による複素行列の逆て進めている。特に、今年度からは東北大学 CC

行列化計算のタスク並列の関係者との定期的なディスカッションも開始

昨年度の、逆行列化計算の部分の MPI 並列化しており、より綿密な共同研究の体制を確立しでの改善に続き、ハイブリッド並列の実装に向2

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

け、逆行列化計算のノード内タスク並列の検討

表 1-2. SX-ACE 上 で の Lapack の libraryを行った。今年度前期は、29 年度最終成果報告(ZHETRF/ZHETRI: Mathkeisan) の計測結果で、会でも報告した、8000 から 20000 次元程度の複14965 次元の複素行列の場合。それ以外の条件は、素行列の逆行列化について、タスク並列化のテ

表 1-1 と同じ。

ストを進めた。
これまでの MPI 並列化コードで用いた Lapackのインターフェースを変えずとも、新たにコンパイルオプションを追加するだけで Mathkeisanによるタスク並列化が起動できるので、その性能を検討した。

ZHETRF (秒)ZHETRI (秒)no OpenMP

76.1 (1.00)

150.6 (1.00)

OpenMP (2)42.7 (1.78)

134.2 (1.12)

OpenMP (4)25.6 (2.97)

79.8 (1.89)

表

GWST コードでは計算される倍精度複素誘電1-3. SX-ACE 上 で の Lapack の library(ZHETRF/ZHETRI: Mathkeisan) の計測結果で、行列がエルミート行列なので、LU 分解による逆20121 次元の複素行列の場合。それ以外の条件は、行列化計算には、メモリ節約を主な理由に、上三表 1-1 と同じ。

角・下三角行列を扱う Lapack のライブラリのコード（ZHPTRF および ZHPTRI）を使用してきた。
しかし、これらのコードをタスク並列化しても、並列化の効果が非常に小さいことが確認された。
その主な理由は、下位で並列化のためにコールされる ParBLAS の処理よりも、LU 分解・逆行列

ZHETRF (秒)ZHETRI (秒)no OpenMP

181.3 (1.00)

361.1 (1.00)

OpenMP (2)99.2 (1.83)

234.5 (1.54)

OpenMP (4)57.5 (3.15)

194.9 (1.85)

としても問題はないことが分かり、結果、それほ化のコード自身の処理のコストの方が非常に大ど使用メモリは増加しない。表 1-1 から 3 に、測きいため、と考えている。

定されたタスク並列化の効果を示す。経過時間

そこで、LU 分解と逆行列化のコードを、それの計測には、簡易性能解析機能 ftrace を使用し、ぞれ、二次元配列をより高速に処理可能な行列の次元は 8279、14965、および 20121 の三種

ZHETRF・ZHETRI に変更した。二次元配列を扱類を検討した（対象物質はナフタレン: C10H8）。

うために書き換えるコードの箇所はそれほど多

まず OpenMP なしでの単体性能で、これまでくなかったが、使用メモリの増大が懸念された。の ZHPTRF・ZHPTRI のコードと比べて 1.5 倍以コードの各箇所の計算精度や計算量の有効数字上の経過時間の軽減があることを確かめた。

を考慮した結果、一部の一時データの精度を落OpenMP によるタスク並列化は、SX-ACE の1-1. SX-ACE 上 で の Lapack の library

最大値の 4 スレッドまで検討したところ、スレ(ZHETRF/ZHETRI: Mathkeisan)の計測結果。タスッド数の増加に伴う効率化が認められた。特に、ク並列なし、またはタスク並列（2 または 4 スレ

LU 分解（ZHETRF）については、行列の次元のッド）による 8279 次元の複素誘電行列の LU 分解増加とともに、並列化の効率が増す様子が伺え（ZHETRF）
・逆行列化（ZHETRI）での経過時間。

る。たとえば、4 スレッドのタスク並列化によっ

カッコ内は、それぞれの経過時間でタスク並列化て、8279 次元ではタスク並列化しない場合の 2.5なしの経過時間を割ったもの。

倍の高速化だったのが、20121 次元では 3.2 倍の表

ZHETRF (秒)ZHETRI (秒)高速化となった。次元の増加とともに、演算処理

no OpenMP

22.0 (1.00)

42.8 (1.00)

時間の増加により、それ以外のオーバーヘッド

OpenMP (2)13.3 (1.65)

24.4 (1.75)の占める割合が減っていると考えられる。さらOpenMP (4)8.8 (2.50)

21.9 (1.95)に次元が増すと、理想的な 4 倍の高速化に近づ3

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

くと期待される。ともその原因を議論し、今後のさらなるタスク

しかしながら、それに続く逆行列化（ZHETRI）並列化の導入に生かしたいと考える。

では、LU 分解のように行列サイズの増加にともB. 非周期系でのクーロン相互作用カットオフのなう並列化の効率化が、必ずしも認められないのが分かる。基本的にいずれの行列サイズでも実装: 表面系や、孤立分子系への適用

スレッド数の増加にともない並列化が効率化さ

一般に固体・結晶系の電子状態・バンドエネれるのは認められるが、たとえば 14965 次元でルギーの第一原理計算では、周期的ポテンシャスレッド数が 2 の場合に、タスク並列なしと比ル下の電子波（ブロッホ波）を効率よく計算する

較すると高速化の効果はわずか 1.1 にしかならために、周期的境界条件の下で平面波基底関数ず、これは他の行列次元の場合よりも低い。

を用いる。

2. SX-ACE 上 で の Lapack の library

表面などの二次元周期系では、表面垂直方向(ZHETRF/ZHETRI: Mathkeisan)の計測結果で、に周期性がないため、周期的繰り返し単位（ユニ22997 次元の複素行列の場合（パーフルオロペンットセル）を、その方向に沿って大きめに取り、

タセン: C22F14）。

隣接セル同士が相互作用しないように工夫する

表（スーパーセルの方法；図 1 を参照）
。スーパーZHETRF (秒)ZHETRI (秒)no OpenMP

269.3 (1.00)

540.2 (1.00)

金属などの表面では電子間相互作用が短距離的OpenMP (2)144.7 (1.86)

324.1 (1.67)

で小さいため、有効である。しかし、一般に固体

OpenMP (4)82.4 (3.27)

290.0 (1.86)

中の多体電子間のクーロン相互作用は、誘起分セルの方法は、自由電子の振る舞いに近い典型極で遮蔽されても長距離的である。人為的なセ
これまで述べたタスク並列化での特徴が、計ル間の相互作用を避けるには、100 Å以上のセ算対象の物質や原子種に依存する傾向であるかル長が必要と見積もられるが、そのようなサイどうかを確認するため、本年度後期には、ナフタ

ズのスーパーセルの計算は、今日の計算機環境

レン単結晶以外の有機分子電子材料もいくつかでも実行は容易ではない。

検討した。その例の一つ（パーフルオロペンタセ

その状況の改善のため、2006 年に、本来相互ン: C22F14）の結果を表 2 に示す。ここでは行列作用すべきでない非周期的な方向のクーロン相の次元が 22997 で、表 1-3 よりもやや大きな次元互作用を、適当な距離（カットオフ）で打ち切るだが、全体的な傾向は表 1-3 に近い傾向にあり、方法が提案された [C. A. Rozzi et al., Phys. Rev. B

次元数がやや大きくなったことでスレッド数の73, 205119 (2006); S. Ismail-Beigi, ibid. 73, 233103

増加に対する効率化が全体的に少し増している。

(2006)]。この方法では、たとえば二次元シート物

これ以外の物質についても、並列の効率化の特質の場合、面に垂直な方向のクーロン相互作用徴は大体同じ傾向にあるので、本研究で見出さが、シートからの距離|z|について、カットオフ zc

れた特徴は、物質によらず、誘電行列のサイズでを超えるとゼロになる（図 1）。本年度の前期・概ね決まっているようである。

後期では、後述の二次元系のクーロン相互作用

LU 分解と逆行列化のうち、逆行列化のコードカットオフの表式の組み込みと、その動作テスで並列化の効果が比較的小さい上、その効果もトを、典型的な二次元シート物質について行な

LU 分解に比べて安定しない理由は、現時点ではった。

不明であるが、より大きなサイズのアプリケーションに向け、さらに大きな次元での実行・テス
トも試していくとともに、東北大学 CC の関係者
4

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月i) ガ ン マ 点 (k= 0) に お け る 複 素 誘 電 行 列

eGG’(k=0)のe1,1 要素 (“head”)と、e1,2..N およびe2..N,1 要素（"wing”）を分離し、残りの(N - 1) ×(N - 1)の行列（"body”）だけで逆行列化計算を行う。
ii) eGG’-1(k=0)とクーロンポテンシャルの積から
遮蔽ポテンシャル W GG’ (k=0)を得る。”body”部分の計算はそのまま行うが、”wing”部分は数値的に無視しても問題ないことが分かって
図 1. 二次元系でのクーロン相互作用打ち切りのおり、ゼロと近似する。”head”部分は、ここで概念図。表面からの距離|z|が zc を超えると、クーは計算しない。
iii) 次に W GG’を逆フーリエ変換し、W rr’を得る。

ロン相互作用はゼロとなる。

誘電関数の”head”部分は、遮蔽されないクー周期的ポテンシャルの下でのブロッホ波の計ロンポテンシャルからの寄与を与えるので、算では、実空間での表式よりも、逆格子空間の表その解析的表式を足し合わせるだけで、遮蔽式で計算する方が数値計算の収束に有利である。

相互作用の計算が完了する。

二次元系の、カットオフ zc のクーロン相互作用

このようにして、コード中での逆空間表示のの表式は、フーリエ変換により波数 k の関数とクーロンポテンシャルの計算における数値的発して、逆格子基底 G の下で次の表式となる。

散を回避することができた。

𝜐" (𝐤 + 𝐆) =4𝜋.1 − 𝑒 2(345)6789 cos[(𝑘 + 𝐺)8 𝑧" ]B|𝐤 + 𝐆|ここで(k + G)xy は、逆格子ベクトル k + G のxy 平面への正射影のノルム、(k + G)z は z 成分である。また、カットオフ zc をスーパーセルの面垂直方向の長さ Lz の半分と設定することで表式が簡略化されているが、一般性は失われていない。
この表式の数値的扱いで注意すべき点は、表式の分母に現れるベクトルのノルムが数値的に図 2. GaN 二次元シートにおける、静的な巨視的ゼロになり発散しうることである。そこで数値誘電関数の波数 q への依存性。

的発散を数値テクニックで回避し、クーロン相互作用カットオフの計算を実装する必要がある

この実装のテストとして、二次元系のテストが、GWST コードではその実装が行われていな計算を行った。GaN の二次元シートにおいて、計かった。

算された静的な巨視的誘電関数(e00-1(q, w=0))を本研究では、GWST コードの次の特長[C.

図 2 に示す。また、ガンマ点 (k= 0) での直接バ

Freysoldt et al., Comput. Phys. Commun. 176, 1

ンドギャップは 5.4 eV となった。これらの結果(2007)]を生かし、以下のようにクーロン相互作は、以前、他の GW 近似計算コードで得られた結用カットオフの実装を試みた。

果[S. Ismail-Beigi, Phys. Rev. B 73, 233103 (2006)]
5

学際大規模情報基盤共同利用・共同研究拠点 平成 30 年度共同研究 最終報告書 2019 年 5 月と数値的によく一致している。この結果を見る

バンドギャップを正確に再現できたため、その限り、実装が成功していることが示唆される。

成果を盛り込み、Book chapter として出版するこ

一方、新たなテスト対象として、半導体材料とができた（「7. 研究成果リスト(5) その他」
、として注目されるリンの二次元シート（フォスを参照）。

フォレン）を検討した。本研究でのクーロンカットオフの導入の仕方に近い方法 [F. A. Rasmussen

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

et al., Phys. Rev. B 94, 155406 (2016)] でのバンド

今年度は、GWST プログラムコードのハイブギャップの計算値（1.8 eV）が最近報告された。

リッド並列化に向けた取り組みとして、主に、本研究でも 1.9 eV と計算され、まずまず一致しSX-ACE（1 ノード内 4 コア）上での LU 分解にていると言える。

よる逆行列化計算の、ノード内タスク並列化のしかし、本研究での場合、2 次元ブリルアンゾ検討を行った。その目的に向けて、これまでよりーン内の k 点サンプル数を増やして電子間相互も高性能なライブラリのルーチンを使用し、ノ作用のおよぶ範囲を大きくすると、バンドギャード・コア単体での性能向上と、SX-ACE の最大ップが数値的には無視できないほど減少する傾スレッド数まで並列度を上げることによる並列向が見られた。もしクーロンカットオフが正し化の両方を達成できた。

く導入されていると、上述の 2016 年の論文でもしかしながら、LU 分解の箇所に比べて、それ見られるように、k 点サンプル数を増やしてもバに続く逆行列化のタスク並列化では、スレッドンドギャップはほぼ全く変わらないのが正しい。

数の増加に対して並列化の効率が安定に増加しその問題の理由を検討した結果、複素誘電行列ない場合がある、などの問題が見られた。この問

eGG’ (k)のブリルアンゾーン内の波数 k への依存

題については、東北大学 CC の関係者との議論や、性、特に、ガンマ点 (k= 0) 近傍での振る舞いは場合によっては計算機ベンダーの協力も得なが正しく記述できているが、その逆行列にクーロら、その原因と解決策を考えていきたい。

ンポテンシャルをかけて遮蔽ポテンシャル W GG’

今後は、東北大学 CC の関係者とのディスカを構築する際には、k → 0 の極限でのガンマ点近ッションの機会を得ながら、計算プログラムの傍の数値的扱いを当初予期していたよりも、注他の箇所の効率の解析も依頼し、共同研究関係意深く行う必要があることがわかった。それでをさらに深めることでプログラムコードの効率も、そのための改善策はほぼ分かっているため、化を進めていく。

来年度の継続課題を開始後間もなくには、より本研究計画での GW 近似の計算による達成目

信頼性が向上したクーロンカットオフの取り扱標の一つとして、表面系への適用を可能にするいができることが期待できる。

取り組みも行った。ここでは、従来提案された、クーロン相互作用カットオフの実装を GWST コC. その他の実績ードに初めて導入した。その際には、数値的な発昨年度からのプログラムコードの改良により、散を回避するための工夫も行った。この作業にこれまでよりも計算が速くなって、比較的大きより、表面系についても、信頼のできる取り扱いな単位胞の有機分子結晶でも、計算結果の収束が可能になりつつある。一部、数値的な修正が必のための数値パラメータのチェックに要する時要な点も改善策はわかってきており、来年度に間を、以前よりもかなり短縮できることが分かは、実際に表面や界面での電荷注入準位のシミった。このことを生かして、以前出版した論文で、ュレーションが可能になると予期される。

計算時間の制約により詳細に検討できなかった計算プログラムが以前より速くなったことに有機半導体結晶での GW 近似によるバンド構造・より、表面・界面のシミュレーションも含め、新6

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

たな応用計算が射程に入ってきており、そのような成果も昨年度から引き続き得られている（「7. 研究成果リスト(5) その他」を参照）。継続課題となる 2019 年度も、共同研究利用のために提供される大規模な計算機資源を活用し、多様な成果の創出を図りたい。
また、「5. 今年度の研究成果の詳細」の項目

A, B それぞれについて、ご協力をいただいた
東北大学サイバーサイエンスセンターの山下毅氏、および、ドイツ・バーダーボルン大学のArno Schindlmayer 教授に、この場を借りて感謝申し上げます。
7.

研究成果リスト(1) 学術論文(2) 国際会議プロシーディングス(3) 国際会議発表
S. Yanagisawa, “Enhancement of the GW space-time
program code for accurate prediction of theelectronic properties at surfaces and interfaces inorganic electronics materials”, 29th Workshop on
Sustained Simulation Performance, Cyberscience
Center, Tohoku University, March 19-20, 2019
(Invited Talk).
(4) 国内会議発表
柳澤将,「有機半導体の結晶構造と電子状態: 分子間ファン・デル・ワールス力の記述と、多体摂動論によるアプローチ」, 第二回琉球大学計算科学シンポジウム, 2018 年 10 月 19〜21 日,琉球大学工学部 2 号館．
(5) その他（特許，プレス発表，著書等）
S. Yanagisawa and I. Hamada, A Chapter inTheoretical Chemistry for Advanced Nanomaterials-Functional Analysis

byComputationandExperiment (Springer Nature, 2019), "Nano-scale
first-principles electronic structure simulations ofmaterials relevant to organic electronics"

7