SENAC Vol. 48, No. 3（2015. 7）— 7 —［共同研究成果］

界面反応の第一原理シミュレーション
柳澤将 1,2, 稲垣耕司 2, 濵本雄治 2, 木﨑栄年 2,森川良忠 2
1
琉球大学理学部, 2 大阪大学大学院工学研究科

1. 緒言
本研究では半導体デバイスや有機デバイス界面、燃料電池電極界面など、応用上重要な界面での電子状態と反応過程を高精度かつ高効率に計算し、界面の性質を支配する要因を明らかにし、それによって、より良い界面を設計する為の指針を与えることを目指している。本
年度は第一原理電子状態計算プログラム STATE の高速化、特にレプリカ並列によって反応経路上の複数のレプリカを一つのジョブによって効率的に最適化し、エネルギー最小の反応経路を求める計算を可能にした。MPI プロセスでデータのやり取りを行う際に問題が生じていたが、計算機科学のエキスパートの助けを借りてこの問題を解決することができ、これによって、半導体表面エッチング過程や二酸化炭素還元触媒反応などの計算が効率的に行えるようになった。具体的には GaN 表面の水によるエッチング過程の研究、および、CO2 の銅表面上での解離反応過程、水素化反応過程の研究を行った。CO2 の銅表面上での解離反応過程については論文化するところまで進み、出版を行った。
また、有機機能性材料として重要なピセンの単結晶の電子状態について研究を行った。密度汎関数法により求めた固有関数・固有エネルギーを 0 次とし、多体摂動法の GW 近似に基づく自己エネルギーの計算を行った。GW 近似で得られる準粒子エネルギーは光電子分光の実験で観測される占有・非占有状態と物理的に対応する。GW 近似は、DFT 計算に比べて計算機負荷が膨大だが、我々は、計算対象のサイズに対するスケーリングで有利な GW space-time 法のプログラムを東北大学サイバーサイエンスセンターSX-9(以下 SX-9)上で整備し研究を進めた。本研究の計算シミュレーションは演算性能に加えて 2TB 以上の大きな主記憶が必要であり、SX-9 上で 32〜64 CPU を使用することによって初めて可能となった。

2. GW 近似による高精度なバンド計算の有機半導体結晶への適用
有機半導体を材料とする電子デバイスは、次世代の電子材料として国際的に研究開発がさかんである。最近の実験で、ルブレン、ペンタセン、ピセンなど一部の有機半導体の単結晶や薄膜で無機半導体に劣らない正孔移動度が報告され、従来考えられたホッピングによるキャリア伝導だけでなく、バンド伝導的なキャリア伝導機構の寄与の可能性が指摘されている。
著者らは最近、GW 近似による高精度なバンド計算を SX-9 上で、ルブレン、ピセンおよびフタロシアニンの単結晶について行い、一般にファン・デル・ワールス(vdW)力で緩く結合する分子間に、パイ電子雲が非局在化して広がることについての知見を得てきた[1-3]。本稿で

— 8 —

SENAC Vol. 48, No. 3（2015. 7）は、それらの研究例のうち、ピセンの単結晶に関する成果[2]を報告する。

2.1

計算方法

ピセン(分子式: C22H14)の単結晶の単位格子は 2 分子からなる単斜晶(図 1)で、格子定数は実験値に固定した。STATE プログラムコード[4]を用い、格子内の原子位置を最適化した。原子核のクーロンポテンシャルをノルム保存擬ポテンシャルで表し、価電子波動関数を平面波基底で展開した。電子の交換・相関ポテンシャルには PBE 汎関数を用い、分子間の vdW 力の記述のため、半経験的に vdW 力を PBE 汎関数に取り込んだ[5]。
こうして得られた固有関数・固有エネルギーを 0 次とし、多体摂動法の GW 近似に基づく自己エネルギーの計算を行った。GW 近似で得られる準粒子エネルギーは光電子分光の実験で観測される占有・非占有状態と物理的に対応する[6]。GW 近似は、DFT 計算に比べて計算機負荷が膨大だが、我々は、計算対象のサイズに対するスケーリングで有利な GW space-time 法[7]のプログラムを SX-9 上で整備し研究を進めた。本研究の計算シミュレーションは主に SX9 上で 32〜64 コアを使用して行われた。

図 1. ピセンの単結晶の図。a,b,c は結晶軸を表し、b 軸方向に分子がヘリングボーン状に積層する様子を示す。実線は単位格子を表し、単位格子の長さは、a: 0.8480 nm, b: 0.6154 nm, c: 1.3515nm。格子ベクトルのなす角度は、α = γ = 90°、β = 90.46°。

2.2

計算結果・考察

PBE によるバンド構造(図 2 の点線)は先行研究[8]とよく一致し、最高占有軌道に由来するバンド(HOMO バンド)がΓ-X, Γ-Y 方向で 0.2-0.5 eV 程度の、比較的大きなバンド幅を与えた。

界面反応の第一原理シミュレーション— 9 —た。最低空軌道由来のバンド(LUMO バンド)は、それに比べて分散が小さい。
次に、GW 近似によるバンド計算の結果について述べる(図 2 の実線)。本研究では SX-9 を利用し大規模な系での計算実行を可能にしており、収束の達成も容易になって、以降の議論に十分な計算値の収束 (HOMO バンド幅: 0.01 eV 以内)に達している。
直接バンドギャップは PBE より 1.4 eV 増加して、3.8 eV 程度となり、DFT で一般に過小評価されるバンドギャップについて大幅な改善が見られる。バンドギャップの実験値の報告はないが、結晶の励起エネルギー(光学ギャップ)の実験値が 3.1-3.3 eV 程度で、バンドギャップとのエネルギー差 0.5-0.7 eV が励起子結合エネルギーに相当すると考えると、有機結晶の一般的な励起子の描像に近く、バンドギャップの計算値は妥当と考えられる。
HOMO バンドをΓ点付近で放物線にフィットすることによって、Γ-Y 方向の有効正孔質量[m h */m e ]は、PBE で 1.26、GW 近似で 1.06 と見積もられた。GW 近似による固体中の多体電子間相互作用の効果により、正孔キャリアは、より自由電子的になると説明される。このことは、ピセンで報告された数 cm2 V-1 s-1 程度の比較的高い正孔移動度[9,10]が、バンド伝導的な機構の寄与に由来する可能性を示唆すると考えられる。
上述のバンドの曲率の増加に対応し、特にΓ-Y 方向で、HOMO バンド(Hu)の幅が、GW 近似の自己エネルギー補正で 0.1 eV ほど増加するのが分かる(図 2)。これは、著者によるルブレン単結晶の先行研究[1]でも明らかになったように、対応する結晶中の方向(ここでは b 軸の、分子がスタックする方向)で、各分子の HOMO に由来する状態が分子間位置にまで非局在化することに対応している。実際、Hu バンドの Y 点での波動関数の空間分布を見ると、分子間で少なからず混成して広がる様子が分かる(図 3)。
最近の角度分解紫外光電子分光(ARPES)の実験によって、結晶中の 2 分子に由来する 2 枚のサブバンド(Hu, Hl)が分解されて観測され、高エネルギー側の Hu バンドのΓ-Y 方向の分散が0.18 eV に対して、低エネルギー側の Hl バンドのΓ-Y 方向の分散は 0.15 eV 以下、と報告された[11]。本研究の計算結果は、Hl の方が Hu よりも分散は小さく、実験結果を定性的に再現できているが、実験値よりも倍程度大きな値であり、定量的に一致していない。その不一致の原因として、室温付近の実験であるために ARPES のピークに熱によるゆらぎの効果が現れるのに加え、電子−振動結合の効果もバンド幅の大きさと同程度のオーダーで寄与し、その結果、絶対零度下の精密な電子状態計算の結果とずれてくることが挙げられる。
今後は、電子状態計算の精密化に加え、電子−振動結合の評価も含めたバンド構造の再現も可能にし、実験結果の解釈や、有機半導体中のキャリア伝導の支配要因の解明に寄与していきたいと考える。

— 10 —

SENAC Vol. 48, No. 3（2015. 7）

図 2. ピセン単結晶のバンド計算の結果。Γ-X、Γ-Y、Γ-Z 方向がそれぞれ結晶 a, b, c 軸に対応。バンド図は、Γ点の HOMO バンドエネルギーを基準とするエネルギー分散を表し、PBE によるバンド図を点線で、GW 近似によるバンド図を実線で示した。Hu (Hl) は、単位格子中の各分子の最高占有軌道(HOMO)由来の 2 つのバンドのうち、エネルギーが高い(低い)バンドを表す。

図 3. ピセン単結晶での HOMO(Hu)の波動関数の等高面図。 等高面は大きさ 2.21 Å−3/2 に対応し、色の違いは波動関数の符号の違いを表す。単位格子を実線で示し、右図では、左図の太い点線による等高面の断面を表す。

— 11 —

界面反応の第一原理シミュレーション

3.

触媒援用加工法におけるエッチング化学反応解析

本章では触媒援用加工法(CAtalyst Reffered Etching method; CARE)におけるエッチング化学反応解析に第一原理量子シミュレーションを適用した例について述べる。CARE 法は、Pt等の触媒をエッチング液中で加工対象表面に接触・搖動させ、その部分でのみエッチング液中の化学種を活性化して表面原子をエッチングすることにより、原子レベルで平滑な表面を得る方法(図 4）[12]であり、現在、SiC、GaN などのワイドバンドギャップ半導体に対しての加工技術開発が精力的に進められている。本研究の目的は第一原理量子シミュレーションを用いてエッチング反応を解析し、平坦化や触媒機能などのメカニズムを明らかにすることにより技術開発の方向性に示唆を与えることにある。ここでは LED やパワー制御半導体デバイスで重要な材料である GaN 表面の水による加工で平坦化される機構を明らかにした例について述べる。
計算には我々が開発した第一原理分子動力学シミュ押し付けて搖動

レーションコード STATE を用いた。STATE は密度汎関数法に基づいており、波動関数を平面波展開法、原子核ポテンシャルをウルトラソフト擬ポテンシャル、交換相関項を GGA-PBE の方法でそれぞれ扱っている。反応の活性

Catalyst
Solution

化エネルギーは CI-NEB 法によって求めた。この手法では、2 つの準安定原子構造の間をつなぐ反応経路に沿っSampleて複数の原子構造を仮定し、反応経路方向以外の自由度について構造最適化しつつ反応座標を収束させていくことで最小エネルギー障壁の反応経路を探索する。また障壁近くの構造でのみ反応座標方向に極大エネルギーの構造を探索することで遷移状態を求める。複数の本研図 4. CARE 加工法の概念図
液中で触媒とサンプルを接触させ、そこでのみエッチング反応を起こさせることにより加工対象物の表面を平坦化する。

究では、1 つの MPI プロセスでこれらの複数の構造を同時に扱って反応経路最適化を行うプログラムを開発し、これを利用して解析を行った。ここでは、水分子の GaN 表面への解離吸着の反応性が表面のステップ部分とステップから原子が欠損したキンク部でどのように異なるかで平坦化メカニズムを調べた結果について示す。モデルは計算の簡単化のため現実のウルツ鉱構造でなく閃亜鉛鉱構造を用いている。 (111)面のステップ-テラス構造とステップのキンク部を模した構造で生じる表面への水分子の解離吸着反応をそれぞれ NEB 計算により解析した(図 5、6)。活性化障壁はステップ部、キンク様部でそれぞれ 1.18eV、0.81eV となっておりキンク様部での反応が容易である結果となった。この差は、図からわかるようにキンク様モデルの場合は N 原子が回転することにより Ga を束縛する作用が弱められているためであると考えられる。この解析は水分子 1 つが解離吸着するエッチングのごく初期の過程だけであるが、水分子が次々に解離吸着して Ga をエッチングすると予想される反応過程において、

— 12 —

SENAC Vol. 48, No. 3（2015. 7）

キンク様部分での反応では変形の自由度が大きく、原子構造の構造緩和の程度が大きいため障壁が下がることが容易に推察できる。この結果から、表面平坦化のメカニズムが、エッチングが表面ステップ構造のキンク部で起こりやすいためと理解される。

(a)Initial model(b) Meta-stable(c) Final

図 5. GaN ステップサイトへの水分子の解離吸着過程。球は原子を表し、色により元素が区別される。紫：Ga、灰：N、赤：O、水色：水素、黄：吸着する水分子の O。表面を上から見ている。奥にある原子ほど薄く描かれている。Ga-N のバックボンドが切れて Gaに水分子が吸着した(b)後、水分子の H が N を終端する。活性化障壁は 1.18eV であった。

(a)Initial model(b) Meta-stable(c) Final

図 6. GaN キンク様サイトへの水分子の解離吸着過程。図 2 と同様の表示である。ステップ部での吸着と同様に Ga-N のバックボンドが切れて Ga に水分子が吸着した(b)後、水分子の H が N を終端する。このモデルでは Ga の動きに伴って結合している N が回転していることが見て取られる。活性化障壁は 0.81eV であった。

今回の解析では、触媒の作用や反応点周辺に存在する他の水の分子は取り扱っておらず、今後の課題となっている。特に、水は水素結合を介してネットワークを形成しており、これ界面反応の第一原理シミュレーション— 13 —が反応において重要な役割を果たしていることは言うまでもない。周囲の水を含めた系全体での反応障壁を解析するためには、時々刻々変化するネットワーク構造に対応した反応障壁を熱力学的に平均化する必要があり、非常に多量の計算資源が要求される。今後、ブルームーン法やメタダイナミックス法などを用いた自由エネルギー障壁計算を効率化するプログラム開発を実施する予定である。

謝辞
本研究は、東北大学サイバーサイエンスセンターのスーパーコンピュータ SX-9 を利用することで実現することができた。また、研究にあたっては同センター関係各位に有益なご指導とご協力をいただいた。本研究の一部は JSPS 科研費 23226004 の助成を受けた。

参考文献[1] S. Yanagisawa, Y. Morikawa and A. Schindlmayr, Phys. Rev. B 88, 115438 (2013).
[2] S. Yanagisawa, Y. Morikawa and A. Schindlmayr, Jpn. J. Appl. Phys. 53, 05FY02 (2014).
[3] S. Yanagisawa, K. Yamauchi, T. Inaoka, T. Oguchi and I. Hamada, Phys. Rev. B 90, 245141(2014).
[4] Y. Morikawa, H. Ishii and K. Seki, Phys. Rev. B 69, 041403 (2004).
[5] S. Grimme, J. Comput. Chem. 27, 1787 (2006).
[6] L. Hedin, Phys. Rev. 139, A796 (1965).
[7] M. M. Rieger et al., Comput. Phys. Commun. 117, 211 (1999); C. Freysoldt et al., Comput. Phys.
Commun. 176, 1 (2007).
[8] T. Kosugi, T. Miyake, S. Ishibashi, R. Arita, and H. Aoki, J. Phys. Soc. Jpn. 78, 113704 (2009).
[9] N. Kawasaki, Y. Kubozono, H. Okamoto, A. Fujiwara, and M. Yamaji, Appl. Phys. Lett. 94,043310 (2009).
[10] H. Okamoto, N. Kawasaki, Y. Kaji, Y. Kubozono, A. Fujiwara, and M. Yamaji, J. Am. Chem. Soc.
130, 10470 (2008).
[11] Q. Xin, S. Duhm, F. Bussolotti, K. Akaike, Y. Kubozono, H. Aoki, T. Kosugi, S. Kera, and N.
Ueno, Phys. Rev. Lett. 108, 226401 (2012).
[12] 応用物理:原英之,佐野泰久,有馬健太,山内和人,触媒基準エッチング法, 77, 2, 168-171,(2008).