— 1 —SENAC Vol. 50, No. 3（2017. 7）［共同研究成果］

有機-金属界面での電子準位接続の精密決定に向けた
GW space-time コードの高速化
琉球大学理学部 柳澤 将

1. 緒言
有機分子からなる半導体は、材料の柔軟性や印刷工程に近い安価な製造工程、低消費電力などを利点とし、次世代エレクトロニクスの材料の候補として注目されている。
材料の効率・性能を本質的に支配すると考えられる基礎電子物性として、i) 結晶での注入電荷(電子・正孔) の移動度の支配要因[1]、ii) 有機層と金属電極層の接触する界面での電荷注入障壁の形成要因[2]、という問題が注目されてきた。
i)について、ルブレンなどの一部の有機分子単結晶・薄膜で、無機半導体で典型的な伝導機構（バンド伝導）の重要な寄与の可能性が指摘された[3, 4]。報告された角度分解光電子分光によるバンド構造や、そこから示唆される正孔伝導の性質について洞察を得るため、筆者らはルブレン単結晶において高精度なバンド計算法(GW 近似[5])を適用した。そこから、バンド伝導的な正孔移動の起源として、隣接分子間に非局在化する分子軌道の役割を指摘した[6]。他の有機結晶にも GW近似の適用を進め、結晶構造と電子状態との関係[7,8]や、結晶中の注入電荷に対する動的分極の効果の大きさについても知見を得てきた [9,10]。
一方、ii)の有機–金属界面での準位接続の問題でも、GW 近似の適用により、準位接続を支配する要因について洞察を得られると期待できる。特に、有機–金属界面では、結晶中の注入電荷が金属表面に近づく際に受ける鏡像ポテンシャルの効果の重要性が、GW 近似の適用により指摘された[11]。しかしながら、その後、GW 近似が有機–金属界面の準位接続に適用された例は、一部の典型的な場合を除くと少なく、長距離的な非局所ポテンシャルの計算を要する GW 近似の計算を、表面・界面系に適用するのが容易でないことを示唆していると言える。
本共同研究では、筆者がこれまで有機結晶系に適用を進めてきた GW space-time(GWST)法[12-14]を、一般的な有機–金属界面や表面の電子準位を決定する問題に適用できるよう、計算プログラムを高速化・高度化することを目的とした。本稿では、計算法やプログラムの概要紹介とともに、本研究で行ったプログラム上の工夫や、計算方法の高度化による成果について示し、今後の高速化・高度化の指針も述べる。

2. 理論、計算法

GW 近似では、固体内の注入電荷（電子・正孔）の周囲に反対符号の分極雲が動的に誘起されることによる、電子間クーロン相互作用の遮蔽を正しく記述する。このような取り扱いのため、GW 近似で計算されるエネルギー（準粒子準位）は、光電子・逆光電子分光の実験で測定される

— 2 —

SENAC Vol. 50, No. 3（2017. 7）

注入電荷の準位と物理的に対応する。
遮蔽された相互作用 W は、固体内の微視的な誘電関数と、(裸の)クーロン相互作用coulomb とによって次のように表される。

W   1 coulomb
ここで、これら物理量は、本質的には空間での任意の 2 点(r, r’)にある粒子間の相互作用を表すため、一般に W ( r, r’)のような非局所な量の計算を要する。固体のバンド計算で広く用いられる密度汎関数理論(DFT)の計算では一般に局所的な一体ポテンシャルを取り扱い、計算に必要な資源は N 2 から N 3 のオーダーで増えるが、GW 近似は N 4 のオーダーであり、DFT 計算よりも膨大になる。
一般に計算資源を要する GW 近似の計算プログラムでは、数値計算を実行可能にするため、プログラムごとで様々な工夫がされている。一般に周期的な物質系の電子状態計算では、波動関数は波数 k を良い量子数とするブロッホ波で、平面波基底で展開できるため、計算では高速フーリエ変換(FFT)が使われる。GW 近似も同様であるが、GWST 法では、畳み込み積分を行なわずに単純な積で計算が済むよう、計算量の変数を FFT でスイッチする。たとえば、遮蔽相互作用 W
は逆格子点(G,G’)と、虚振動数 iを変数として以下のように表される。





W G, G; i   1 G, G; i coulomb G, G
粒子の準粒子エネルギー演算子は、  r, r; i  iG r, r; i W r, r; iのように、実空間点(r,r’)と虚時間 iに依存する。このようにして計算は簡素化され、計算速度は一般に FFT の速さに依存する。ただし、実空間グリッド上のデータが必要なため、一般にメモリやディスクの要求が多くなっている。
現状では多量のデータのディスク I/O による遅延を避けるため、計算データを極力メモリにロードして計算を行なっており、多くのメモリを使用できる SX-ACE 上で主に使用している。元々のプログラムは英・独の研究グループで開発され[12-14]、並列計算は不可であったが、筆者がフラット MPI を実装している。

3. 結果・考察
3-1. プログラムの効率上の工夫
GWST コードによる大規模な実行を見据え、SX-ACE 上や PC クラスター上での動作の検討・効率化を進めた。東北大学サイバーサイエンスセンターの技術専門職員の助力を得て、 LX406Re-2 と SX-ACE での以下のテストをもとにコードの改善を進めた。

有機 - 金属界面での電子準位接続の精密決定に向けた
GW space-time コードの高速化— 3 —テスト A(小規模: シリコン単結晶) :4〜24 プロセス並列、メモリ使用: 30 GB 程度、演算時間数分
テスト B(中規模: ナフタレン単結晶) :144 プロセス並列、メモリ使用: 400 GB 程度、演算時間 30-50 分
テストの結果、テスト A は LX 406Re-2 上での実行が有利である一方、将来的に目指す大きな系のシミュレーションには SX-ACE の方が有利であることが分かった。コード改良前の段階でSX-ACE での平均ベクトル化率は A: 83.7 %、B: 95.9 %となった。小規模な A の方がかなり低い
ベクトル化率であるのは、ファイルからのデータ読み込み・ロードなど、比較的ベクトル化しにくいプロセスが実行の少なからぬ部分を占めるためで、数値演算の割合が増すテスト B でベクトル化率が上がる。それでもベクトル化の余地がまだあり、ベクトル化率が低いパートを改善した。
ファイルから、一電子軌道などの数値を読み込みストアする部分で、インデックス配列を用いて配列にストアしていた。技術専門職員によって、一旦、一時配列に一括して読み込んだ後にインデックス配列を用いたストアに一時配列から読んでストアするよう改変された。結果、平均ベクトル化率は若干下がったものの、SX-ACE 上でのテスト A で 5%程度の全演算時間短縮、テスト B で 250%の全演算時間短縮が見られ大幅な効率化につながった。しかし LX 406Re-2 では効率化せず、テスト B で 10%程度の全演算時間増につながった。演算の重複化が原因と考えられ、LX 406Re-2 では改良前の方が効率的だが、LX 406Re-2 上の実行で想定されるデータ量に応じ、SX-ACE と別のコードを実装する必要があると考えられる。
それ以外にも、do-loop 内の依存関係でベクトル化が不可であったパートについて、コードの見直し・改訂を行ない、テスト A で 90 %、テスト B で 98 %までベクトル化率が上がり、実行時間もさらに 10%前後の短縮が見られた。
GWST コードでは高速フーリエ変換が多用されるが、特にグリッドサイズの異なる二種類のフーリエ変換(r と G および実格子点 R と k)が必要であり、適切な FFT コードの選択が重要である。
一般に r と G のフーリエ変換は数万〜十万程度の次元に対し、もう一方は数 10〜100 程度の次元となる。前者の FFT は、特に ASL ライブラリを使用すると高ベクトル化率での高速処理が可能だが、後者の FFT ではベクトル化のロードが負荷となり、ASL の使用だと効率が悪い。これまで、小グリッドサイズの FFT には FFTW を使用したが、技術専門職員による解析の結果、多グリッドの FFT での高ベクトル化率と、比較的グリッドが少ない場合でも性能を有する FFTE[15]を使用するのが良さそうであることが分かった。GWST コードでは R と k の FFT のルーチンが数万回以上コールされるが、FFTE はソース自体入手が可能で、ソースからコンパイルしインライン展開を行なうことによるルーチンのロードの負荷もある程度軽減できると見込まれる。テスト A、B で試した限りでは大幅な高速化は見られなかったが、将来的に、テスト B よりもさらに— 4 —

SENAC Vol. 50, No. 3（2017. 7）

大規模な系を対象とし、よりグリッドサイズが大きくなると利点は増すと期待される。SX-ACEの自動並列化機能が FFTE で有効なようであり、ハイブリッド並列の導入による効率化も期待される。
並列化について、GWST コードの現状はフラット MPI 並列のみである。現状、100 並列度程度まではまずまず効率的だが、それを越えると効率がかなり低下する。テスト B 以上の規模の計算でもっとも時間を要するのは微視的誘電行列 GG ’の逆行列化で、現在 ScaLAPACK を使用している。ScaLAPACK による複素行列の LU 分解・逆行列化は、ある程度の並列度・ノード数で効率的だが、現状ではコードの序盤のブロック分割から ScaLAPACK での 2 次元サイクリック分割に合うようにノード間で冗長な通信処理を追加しているため、コードの見直しを行って並列化効率を上げる必要がある。対策として、SX の自動並列化機能を用い、OpenMP/MPI のハイブリッド並列を実装する作業を進め、通信の軽減による効率化を目指している。

3-2. 計算精度・対象の高度化
計算アルゴリズム上の改善に加え、物性予測の精度の向上や、計算対象をより高度化する工夫も行った。上述の効率化の達成により、様々な有機半導体結晶の GW 近似による計算が可能になった。
通常、DFT 計算による解をゼロ次として準粒子エネルギー補正を一度行う(G0W0)。しかし、DFT 計算によるゼロ次エネルギーでのエラーが大きい場合、(n−1)回目の準粒子補正のエネルギー(n ≥ 2)を、n 回目の準粒子補正におけるゼロ次エネルギーとして用いることを繰り返す自己無撞着(self-consistent)な計算を行うことで、精度を向上させることができる(evGW)。このような計算は元の GWST コードにすでに組み込まれていたが、これを現状のコード向けに改良し使用できるようにした。結果、より実験値に近いバンドギャップの予測を達成した(表 1)[10]。
また、現状のプログラムでは表面(2 次元周期系)などの、非周期性を持った系への適用が計算のコスト上、容易ではない。そのため、よく規定された真空準位を GW 近似の計算で定義しイオン化ポテンシャル (IP) や電子親和力 (EA) を決定するのが難しい。しかし、IP と EA の差(バンドギャップ)が 3 次元周期系(バルク)の性質で本質的に決まり、表面の分子配向に依存する静電的効果で IP・EA が真空準位に対してシフトすると仮定すれば、GW 近似による 3 次元バルクのバンドギャップと DFT 計算による表面系(周期的スラブモデル)での静電効果の見積もりを合わせて、表面の配向に依存した IP・EA を見積もることができると考えられる。
図 1 に、そうして見積もった IP・EA の例(ペンタセン結晶)を示す。ここでは、GW 近似で記述される動的分極効果により、結晶・バルク相で、注入された正孔・電子が周囲の分極雲により安定化される(図 1 の Ep に相当。大きさは、正孔と電子で共通と仮定)。
一方、表面配向に依存する静電効果により、真空準位に対して同じ方向にシフトする(図 1 のWes)。Wes は、孤立分子と表面のスラブモデルのエネルギーギャップ中心の差として見積もった。

有機 - 金属界面での電子準位接続の精密決定に向けた
GW space-time コードの高速化— 5 —Wes は正孔・電子で共通で、表面分子の電気四重極モーメントによって決まると議論されている[16,17]。
このような見積もりにより、実験的に確認された、異なる表面の分子配向を反映した IP・EAの変化を計算で再現することができ、IP・EA の数値も実験値[17]とまずまず一致する[18]。相違の原因として、Ep が正孔・電子とで共通と仮定したことなどが考えられる。将来的には、表面のスラブモデル自体の GW 近似の計算によって IP・EA の決定を目指し、表面での電荷注入準位を支配する物理的要因の解明を目指す。

4. 結言
本共同研究では、次世代電子材料として注目される有機エレクトロニクス材料との関連で、有機–金属界面や有機結晶表面などのより複雑な物質系での電子準位の理論的決定を目標に GW
space-time コードの効率化・高速化を進めた。一般に、非局所な物理量を計算するため、DFT に基づく従来の第一原理計算シミュレーションよりも大量の数値データを取り扱う必要がある。現在、数百程度の並列度で SX-ACE 上での実行を主としており、実際、中程度以上の規模の計算では PC クラスターよりも SX-ACE で高速に計算が実行される。本共同研究でもベクトル化率向上
をはじめ、SX-ACE 上での効率向上を達成した。当面は、非周期性を持った系への適用性向上のため、SX-ACE 上での 1000 並列度以上での効率化を目指す。誘電行列の逆行列化を、ノード内のスレッド並列・ノード間の MPI 並列のハイブリッド並列によって並列効率を高めることを目標に進める。より将来的には、ベクトル機実行も含め、メニーコア計算機向けの高効率化も考えていく。

謝辞
本共同研究を進めるにあたり、東北大学サイバーサイエンスセンター共同利用支援係技術専門職員の山下 毅氏には、大変お世話になりました。この場を借りて御礼申し上げます。GWST コードの実行および整備・開発について、ドイツ・ユーリッヒ固体物理研究センターの Arno
Schindlmayr 研究員(現ドイツ・パーダーボルン大学 教授)にも感謝申し上げます。

参考文献[1] V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey, and J.-L. Brédas, Chem.
Rev. 107, 926 (2007).
[2] H. Ishii, K. Sugiyama, E. Ito, and K. Seki, Adv. Mater. 11, 605 (1999).
[3] J. Takeya, M. Yamagishi, Y. Tominari, R. Hirahara, Y. Nakazawa, T. Nishikawa, T. Kawase,T. Shimoda, and S. Ogawa, Appl. Phys. Lett. 90, 102120 (2007).
[4] S. Machida, Y. Nakayama, S. Duhm, Q. Xin, A. Funakoshi, N. Ogawa, S. Kera, N. Ueno,

— 6 —

SENAC Vol. 50, No. 3（2017. 7）

and H. Ishii, Phys. Rev. Lett. 104, 156401 (2010).
[5] L. Hedin, Phys. Rev. 139, A796 (1965).
[6] S. Yanagisawa, Y. Morikawa, and A. Schindlmayr, Phys. Rev. B 88, 115438 (2013).
[7] S. Yanagisawa, Y. Morikawa, and A. Schindlmayr, Jpn. J. Appl. Phys. 53, 05FY02 (2014).
[8] S. Yanagisawa, K. Yamauchi, T. Inaoka, T. Oguchi, and I. Hamada, Phys. Rev. B 90,245141 (2014).
[9] S. Yanagisawa, S. Hatada, Y. Morikawa, J. Chin. Chem. Soc. 63, 513 (2016).
[10] S. Yanagisawa and I. Hamada, J. Appl. Phys. 121, 045501 (2017).
[11] J. B. Neaton M. S. Hybertsen, and S. G. Louie, Phys. Rev. Lett. 97, 216405 (2006).
[12] M. M. Rieger, L. Steinbeck, I. D. White, H. N. Rojas, and R. W. Godby, Comput. Phys.
Commun. 117, 211 (1999).
[13] L. Steinbeck, A. Rubio, L. Reining, M. Torrent, I. D. White, and R. W. Godby, Comput.
Phys. Commun. 125, 105 (2000).
[14] C. Freysoldt, P. Eggert, P. Rinke, A. Schindlmayr, R. W. Godby, and M. Scheffler, Comput.
Phys. Commun. 176, 1 (2007).
[15] http://www.ffte.jp/[16] B. J. Topham and Z. G. Soos, Phys. Rev. B 84, 165405 (2011).
[17] H. Yoshida, K. Yamada, J. Tsutsumi, and N. Sato, Phys. Rev. B 92, 075145 (2015).
[18] S. Yanagisawa, AIP Conference Proceedings, in press.

— 7 —

有機 - 金属界面での電子準位接続の精密決定に向けた
GW space-time コードの高速化

表 1. 典型的なオリゴアセン結晶でのバンドギャップ。DFT 計算でのゼロ次エネルギーに一度だけ準粒子補正した場合(G0W0)と、自己無撞着に繰り返した場合(evGW)。単位は eV。
ナフタレン

アントラセン

テトラセン

G0W0

4.72

3.37

2.40evGW

5.70

4.17

3.09Exp.

5.0〜5.5

3.9〜4.2

2.9〜3.4

図 1 (a): ペンタセン薄膜相の(001)面における電荷注入準位の決定。イオン化ポテンシャル・電子親和力は、それぞれ 4.70 eV・1.99 eV。実験値[17]は、4.90 eV・2.35 eV。