Comments
Description
Transcript
バンド構造を用いた材料開発（実践編）
応用物理学会結晶工学分科会結晶工学スクールテキスト第９版(2008 年) バンド構造を用いた材料開発（実践編） 東京工業大学応用セラミックス研究所 神谷利夫 1. はじめに 電子状態理論（バンド理論）の基礎を学んだわけですが、それでは、この知識を実 際の材料に応用するにはどうしたらよいでしょうか。私たちが電子構造を知るために は量子理論に基づいた膨大な計算が必要となり、コンピュータで数値計算することに なります。このようなソフトウェアを総称して「量子計算ソフト」と呼ぶことにしま すが、現在では第一原理計算が使われる例が多くなっています。最近の第一原理計算 の精度は非常に高く、ある程度の「材料設計」に使えるレベルに到達しています。多 くの量子計算ソフトが入手できるようになっていますが、それぞれに特徴があり、対 象とする材料と物性によって使い分ける必要があります。ここでは、半経験法分子軌 道法ソフトである MOPAC7 1 、第一原理密度汎関数法バンド計算ソフトである WIEN2k 2、VASP 3、CASTEP 4 を使った結果を例に、実際の材料を理解する、ある いは研究する上でこれらをどのように利用していけるかを説明していきます。 なお、ここで掲載するデータや図は、精査して論文に掲載されたものも、本原稿を 書くために突貫で計算してみた、ほとんど検証していないものまで混在しています。 また、私の不勉強ゆえに不正確な記述もあることは否定できません。書いてあること を鵜呑みにするのでなく、あくまでも、量子計算の可能性を述べるための参考データ と考えていただければ幸いです。また、各計算結果の後に (VASP, PAW, PBE) などと 書いてありますが、これは用いたソフトェアと計算モデルを記述しています。この場 合は、VASP を用い、PAW 法で PBE 汎関数を用いたという意味です。 本稿は、1~8 章までは長い前置きになっています。 「量子計算で何ができるのか」だ けに興味がある人は、途中を飛ばして 9 章から読んでみてください。脚注には知って いると役に立つかもしれないことを思いつくままに書きましたが、最初に本文を読む 時には無視して読み進めていただければと思います。 2. 量子計算で何ができるか？ 1) 原子・電子の状態 (電子準位、波動関数、電子密度、電位分布など) を視覚化で きる ・講義で学んだ理論を直観的に理解できる ・自分が組み立てた理論・モデルが正しいかどうかを検証するツールになる 2) 精度の高い全エネルギー計算を行える ・絶対零度の基底状態 5で安定な構造がわかる ・生成エネルギーが計算できる 合成できるはずのない構造をスクリーニングできる 標準状態で安定でない構造を安定化させる条件についてヒントをくれる ・力学物性 (弾性テンソル、格子振動分散など) の計算を行える 赤外分光、Raman 分光スペクトルなどの同定 3) 電子構造を定量的に計算できる (バンド構造) 1 2 3 4 5 Winmostar, http://winmostar.com/ http://www.wien2k.at/ http://cms.mpi.univie.ac.at/vasp/ Materials Studio 4.1, http://www.comtec.daikin.co.jp/SC/prd/material/mol-orbital/index.html 一番全エネルギーが低くなる原子配位、電子構造の状態 ・直接遷移か間接遷移か、許容遷移か禁制遷移か ・光学スペクトル (誘電関数 = 屈折率、吸収係数) ・キャリア輸送特性 (有効質量、フェルミエネルギー、状態密度) 散乱モデルを別途考慮すると、電子伝導度、電子比熱、熱電係数などの 計算もできる ・化学結合性状に関する定性的な知見 共有結合性、投影(部分)状態密度 ・フェルミ面の形状 ・スピン分極した系: スピン配置、スピン配列、自発分極など 4) 電子－格子結合系の物性: 誘電率 (Bery 位相)、圧電定数など 「バンド計算法」が与えてくれる情報は、バンド構造だけではありません。あくま でも結晶について効率的で精度が高い計算ができるがゆえに「バンド理論」が使われ ているだけであり、バンド構造はその一部に過ぎません。 3. 分子軌道法とバンド計算 1) 分子軌道法 その名の通り、分子の電子構造を量子論にのっとって計算する方法です。構造 モデルは必ず端面をもちます。電子構造は、縦軸にエネルギーをとって占有状態、 非占有状態、場合によってはスピン状態を表示した「エネルギー準位図 (energy level diagram)」で表わされます。横軸にエネルギー、縦軸にそのエネルギーをも つ電子の数をプロットした「状態密度 (density of state: DOS)」も使われますが、 分子軌道法の場合には、両者はほぼ同一です。状態密度を各原子や原子の波動関 数へ分解したものを「部分状態密度 (partial DOS: PDOS)」と呼び、化学結合がど の軌道から作られているかを図的に理解するためのツールになります。 基底関数6には原子の波動関数 (原子軌道, atomic orbitals: AO) が使われること が多く、これを LCAO 法(Linear Combination of AO, 原子基底の一次結合)といい ます。そのため、原理的に原子ごとに分解した PDOS を見積もるのに適していま す。実際の計算では、解析的に積分を実行できる Gauss 関数で AO を表現した GTO (Gaussian-Type Orbitals)が使われることがあります (Gaussian03 7) が、数値 解析で求めた原子軌道も基底として使われています (DV-Xα8)。 理論モデルは Hartree-Fock (HF) 近似が使われることが多い (Gaussian03) ので すが、電子相関を考慮していないため、配置間相互作用 (Configuration Interaction: CI) などと組み合わされることも多い。密度汎関数法も使われます (DV-Xα法, Gaussian)。 2) バンド計算 周期的な構造をもつ結晶を扱います。分子と異なり、周期条件のおかげで構造 モデルは端面を持ちません。Bloch の定理を利用して計算対象を単位格子内の原 子だけに減らしていますが、その代償として、バンド計算に固有のパラメータ「波 数ベクトル k」が現れます。エネルギー固有値は k の関数 E(k)になり、電子構造 は、横軸を k、縦軸を E(k)とする「バンド構造」となります。バンド構造を使う 6 7 8 量子計算では、複雑な波動関数を原子軌道や平面波を足し合わせて(線形結合を作って)近似しま す。足し合わせる基本になる関数を「基底関数 (basis functions)」と呼びます。 http://www.gaussian.com/ http://www.dvxa.org/ ことにより電子構造を直観的に理解しやすくなりますが、多くの学生からはわけ がわからないと嫌われ、「スパゲッティ」と呼ばれることもあります。バンド構 造は対称性の高い特定の方向だけでしか表示できず、結晶の全体の電子構造を理 解する場合には見落としが出る可能性があります。バンド計算における DOS は すべての k に関する積分9となるため、全電子状態の情報を判断することができ ます。分子軌道法と同様に、PDOS によって化学結合を作っている波動関数の知 見が得られます。 ただし、結晶の場合、特にエネルギーの浅い伝導帯では、格子間位置に拡がっ た波動関数が無視できなくなり、これを原子の PDOS に帰属させることが困難で す (あるいは帰属させるべきではない)。これを「Interstitial (隙間の の意味)」と 呼ぶことがあります (WIEN2k など)。 この事情から、金属や半導体では原子基底だけでは波動関数の精度が出せませ ん (特に伝導帯など浅い準位)。そのため、エネルギーの浅い準位（非占有準位 と浅い占有準位をあわせて「valence states (価電子状態)」と呼ばれる）に平面波 を基底関数として用いる計算法が主流であり、現在では擬ポテンシャル平面波 (PseudoPotential Plane Wave Method: PP-PW) 法、(線形化した)補強された平面波 ( (Linearlized) Augmented Plane Wave: L/APW) 法、PAW (Projected Augmented Wave) 法などが主に使われています。ただし、原子基底を拡張し、「拡がった基 底 (diffused functions)」や主量子数の大きい AO などを加えることで、LCAO ベ ースでの精密計算も行えます。 後で詳述しますが、平面波を基底に用いると、波動関数がどの原子に帰属して いるかを判断するのに大きな任意性がでてきます。つまり、平面波基底での PDOS には大きな任意性がでるという問題があります。また、単位格子が大きく なると、必要な平面波の数、ひいては固有値問題を解く際の行列が大きくなり、 計算が極端に膨大になります。これらの問題を改善するため、平面波基底の場合 も、Wannier 関数をつくることで原子に局在した基底に作りなおすなどの工夫も 行われています (Order-N 法など)。 密度汎関数法が主流です (WIEN2k, VASP, CASTEP など) が、HF 法を使えるソ フトウェアもあります (CRYSTAL)。 4. 第一原理法、半経験法、経験法 1) 第一原理法 (first-principles method)10 量子論の基礎方程式をなるべく忠実にコンピュータを使って解く方法。基本的 に入力は原子の種類と座標の「構造パラメータ」のみです。結晶の場合には、単 位格子の格子定数 (軸長 a, b, c, 軸角 α, β, γ) と、単位格子内にある原子の種類 と座標を入力しますが11、空間群を指定することで、最小の原子座標だけの入力 で済ませられる場合もあります12。ただし、第一原理計算といっても、構造パラ メータ以外に全くパラメータを入力する必要がないわけではありません。計算の 9 正確には「第一ブリルアンゾーン内の積分」 ”最初から”という意味のラテン語から派生した”ab initio 法”とも呼ばれます。明確に区別されて 使われてはいませんが、前者はバンド計算、後者は分子軌道法で使われることが多いようです。 11 分子軌道法の場合は、座標をデカルト座標で入力する”XYZ 形式”と、原子間距離と角度で分子 の形状を規定する”Z 行列形式”が広く使われています。 12 空間群を用いた結晶構造の理解が必要になりますので、結晶学や群論の教科書を参考にしてく ださい。Bible は『International Tables for Crystallography, Vol. A-D』です。 10 モデルの種類、精度を決めるパラメータ群を適切に入力する必要があります。た だし、これらのパラメータを多少変えたとしても計算結果に与えることは無く （というよりも、数値を変えても計算結果が変わらない程度に高い精度の計算を していることが、正しい結果を得るための必要条件）、本質的な計算結果に影響 を与えるのは、構造パラメータだけです。 HF 法 (CRYSTAL06), FL/APW 法 13 (WIEN2k), 擬ポテンシャル/平面波法 (VASP, PWscf 14, CASTEP), PAW 法 (VASP), Gaussian03, GAMESS などがありま す。DV-Xα法は任意性のあるパラメータαを含むために第一原理ではないという 意見もありますが、本稿ではこの範疇とします。 2) 半経験法 (semi-emprical method) 物理方程式は Schrödinger 方程式を用いるものの、計算に時間がかかる項－エ ネルギー積分や共鳴積分－を経験パラメータに置き換えることで計算時間を短 くした方法。構造パラメータの他に、例えばエネルギー積分、共鳴積分といった、 実験などによって決めた経験パラメータを入力する必要があります。計算結果は 経験パラメータによって大きく変わる（逆に言えば、都合の良い結果を出せるよ うに変えることができる）ため、これらの選択が計算の信頼性を決定します。実 験結果に合うようにパラメータを再設定しなおす必要がでることも多く、この時 間の方が第一原理計算よりも長くかかることも珍しくありません。 分子軌道法で使われる（拡張）Hückel 法, CNDO 法 (MOPAC)、バンド計算へ の拡張である Tight Binding 法などがあります。 3) 経験法 (empirical method) 量子論的な効果を明示的に考慮せず、原子間の相互作用ポテンシャルの形を実 験から決め、古典力学の運動方程式(Newton の運動方程式)を用いて材料中の原子 の運動をシミュレーションする方法。 格子力学(lattice dynamics: LD) 法、分子動力学 (molecular dynamics: MD) 法、 モンテカルロ(Monte Carlo: MC) 法などの多くがこれを採用しています。これら には量子計算法もあり、区別する場合には「古典 MD 法」などと呼ばれます。 5. Schrödinger 方程式と密度汎関数理論 実際に量子計算を行って結果を解釈する際に問題を起こすことは少ないのですが、 Schrödinger の波動力学にでてくる一電子 HF 方程式と密度汎関数 15 理論 (Density Functional Theory: DFT)にでてくる Kohn-Sham (KS)方程式を混同されることが多いよ うなので注意を書いておきます。これは、それぞれの理論の基礎方程式を変数分離し て一電子化した方程式が、下記のように形式上似ているために起こるのでしょう16。 ⎧ 1 2 ⎫ (1) HF 方程式: ⎨− ∇ l + Vext (r ) + Ve −e (rl ) + V Xl (rl )⎬φl (rl ) = ε l φl (rl ) ⎩ 2 ⎭ ⎧ 1 ⎫ KS 方程式: ⎨− ∇ 2 + Vext (ρ (r )) + Ve −e (ρ (r )) + VXC (ρ (r ))⎬φ (r ) = εφ (r ) (2) ⎩ 2 ⎭ ここで Vext は原子核と電子の相互作用で、単純な Coulomb ポテンシャルです。Ve-e は 電子－電子間の Coulomb ポテンシャルです。VX と VXC が違って見えますが、C は電 13 14 15 16 F は Full potential 法、L は Linearlized (線形化法) の意味 http://www.pwscf.org/ 関数 g (x) が関数 f (x) の関数 g( f (x) ) であるとき、g を f の「汎関数」といいます。 原子単位を使っているため、プランク定数や電子質量などの記号が消えています。 子相関 (electron correlation) の意味で、HF 近似では考慮されていないために入ってい ないだけです。ですから、両者はこの段階の解析では同一になり、両者を同じものと 思っても問題は起こりません。 しかし、Schrödinger 方程式が、各電子の座標と運動量からつくったハミルトニアン 17 に物理量の交換関係18を適応して導出することからも明らかなように、Schrödinger 方程式というのは各電子の座標を変数としています。 DFT はまったく異なる概念を基礎とした、別の形式の量子理論です。 「ある電子系 において電子密度ρ(r)が決まると、基底状態の外部ポテンシャル V(r)、系の全エネル ギーE は一意的に決まる」という Kohn-Hohenberg の定理に基づいています。つまり、 基本的な概念として波動関数は現れず、電子密度が基本的な量となります。電子密度 は空間座標だけの関数なので、KS 方程式にでてくる「波動関数のようなφ(r)」の r も 空間座標です。 HF 近似と DFT には他にも、似ているけれども違うことがいくつかあります。HF 近似で得られるエネルギー固有値 ε i はイオン化ポテンシャル ε i= E N − E N −1 (EN は電子 が N 個あるときの全エネルギー) や電子親和力 ε i= E N +1 − E N に対応します19が、DFT では ε i= ∂E / ∂ni (E は全エネルギー、ni は i 番目の軌道を占有している電子数で、非整 数として扱う)と、化学ポテンシャルのような量に対応します20。そのため、DFT で得 られるエネルギー固有値は、光電子分光で測定される束縛エネルギーよりも小さくな る21ため、遷移状態法などによってイオン化ポテンシャルを求める必要があります。 6. バンドギャップができるいろいろな成因 これまで、 「(ほぼ)自由な電子理論 (Nearly Free Electron Theory: NFE)」で、ブリルア ンゾーン (Brillouin Zone: BZ) 境界では電子の進行波と反射波が干渉してバンドギャ ップを形成するということを学びました。しかしながら、私たちが普段問題にする「バ ンドギャップ」は、その下のエネルギー準位が電子で完全に占有されていて、その上 は電子が全く存在しない「バンドギャップ」のことです。この講義では、後者の狭義 のものだけをバンドギャップと呼びますが、その多くは NFE ででてくるバンドギャ ップとは異なる成因で形成されています。 1) 共有結合（等極結合）Si など 結合、反結合軌道のエネルギー分裂によってバンドギャップが形成される。 2) イオン結合（異極結合）酸化物など 22 イオンのエネルギー準位の違いによってバンドギャップが形成される。 原子の 2p 軌道と 3s 軌道には大きなエネルギーギャップ (energy gap) があるこ とと、本質的に同じ。 3) BZ 境界での干渉： Bragg 反射 17 全エネルギーを粒子の座標と運動量であらわしたものと理解しても大丈夫です。 座標と運動量についての xPx - Pxx = ih のような関係。単なる数値同士ではゼロであるが、特定 の物理量同士ではプランク定数のオーダーでゼロではない、というのが量子論の本質。このよ うな物理量は「互いに共役な物理量」と呼ばれ、ラグランジアンから見つけることができる。 19 Koopmans の定理: T. Koopmans, Physica 1 (1933) 104. 20 Janak の定理: J.F. Janak, Phys. Rev. B 18 (1978) 7165. 21 例えば W.R.L. Lambrecht et al., Phys. Rev. B 50 (1994) 14155. 22 固体中のイオンのエネルギー準位を正しく理解するには、 「マーデルングポテンシャル」を考 慮する必要があります。 18 NFE ででてくる広義のバンドギャップ。典型的な例は「パイエルス転移 23」。 4) 強電子相関系材料 上記の 1) ~ 3)によるバンドギャップは、HF 法や DFT などの第一原理バンド計 算でも再現できます24。しかし、これらの計算では現実の材料で現れるバンドギ ャップを再現できないことがあります。典型的な原因に、電子相関を十分取り入 れられていないことが挙げられ、特に d 電子、f 電子を含み、局在性の強い電子 状態を作りやすい遷移金属酸化物で問題になります25。この場合、LDA(GGA)+U などのいわゆる「beyond DFT」を使うことが必須になります。 7. 基本知識: ブラベー格子(Bravais lattice)と基本格子(Primitive lattice) GaAs の結晶構造は「閃亜鉛構造」と呼ばれます (図 1) が、Ga と As を同じ原子に 置き換えると Si の結晶構造である「ダイヤモンド構造」になります。閃亜鉛構造は 8 つの原子を含む立方体の単位格子で表されます（ 「立方晶に属する」といいます）が、 図 1(A)からわかるように、Ga だけでなる面心立方構造 (Face centered cubic: FCC)と As だけでなる FCC 構造が(1/4, 1/4, 1/4)ずれて重なった構造です。しかし、面心立方格 子などは「ブラベー格子」といわれるもので、単位格子に含まれる原子の数を増やし て見かけの対称性が高い形の単位格子をとったものです (図 1(B))。実際には、結晶の 3 次元周期性を表現するのにもっと小さい、原子数の少ない格子を取ることができ、 そのような最小の格子を「基本格子」といいます（図 1(C)）。実際の量子計算では、 扱う原子数を減らして計算時間を短くするため、ブラベー格子ではなく「基本格子」 を使います。計算結果も基本格子ででることが多いので注意が必要です。 文献などの格子定数は、通常はブラベー格子で表され、Si の場合は 296 K において aC = 0.5431 nm です(C は立方晶、cubic の意味)。一方、基本格子である菱面体格子 (Rhombohedral lattice) では aR = 0.3840 nm, αR = 60o になります(R は菱面体の意味)。 (A) Si,As (B) Si,Ga Si,Ga (C) Si,As cC bC aC Si,Ga cR bR aR Si,Asf 図 1 Si, GaAs の単位格子。(A) 2 つの面心心立方構造が重なってできています。 (B) ブラベー格 子である面心立方格子。(C) 基本格子である菱面体格子。ブラベー格子では 8 つの Si 原子を含み ますが、基本格子では 2 つに減っています26。基本格子の体積はブラベー格子の体積の 1/4。 23 フェルミ面が BZ 境界ではなく BZ 内部にある場合、結晶の周期を大きくしてフェルミ面に BZ 境界が来るようにできることがあります。そうすると、バンドギャップが開き、その分、電子 系のエネルギーが低くなります。このエネルギー利得が結晶の周期性を変えることによる弾性 エネルギーの損より大きいと、この周期の構造が安定化します。代表的な例がポリアセチレン。 24 後述のバンドギャップ問題があるので、定量的に再現できるという意味ではありません。 25 バンドギャップを作る原子軌道の組み合わせにより「ハバードギャップ」 、 「電荷移動ギャップ」 などと区別されます。 26 単位格子の角にある原子は 4 つの単位格子、面上にある原子は 2 つの単位格子で共有されてい ることを考慮して、それぞれ 1/4 と 1/2 の重みをかけて原子の数を数えてみましょう。 8. 結晶構造の入力 図 2 には、WIEN2k を使って Si の電子構造を計算する際の、結晶構造の入力画面を 示しています。WIEN2k には、w2web という Web プログラムが付属しており、Web ブラウザを使うことで、Graphical User Interface (GUI) でのプログラムの制御、結果の グラフ表示ができます。これからわかるように、構造モデルとして入力しているのは 空間群(Spacegroup) Fd-3m No. 227 格子定数(Lattice parametes) a, b, c, α, β, γ 非等価な原子の種類と座標 (Inequivalent Atoms) Si, 0.125, 0.125, 0.125 だけです。構造パラメータではない RMT 27 や、計算精度に影響を与えるパラメータ 28 も入力しないといけませんが、これらは脚注を参考にしてください。 図2 WIEN2k の GUI (w2web)によって結晶構造モデルを入力している様子。 9. 実例を通して理解してみましょう 1) 全エネルギー: 生成熱など ほとんどの量子計算ソフトウェアでは、入力した構造モデルの全エネルギーを出力 してくれます。WIEN2k で PBE96 という密度汎関数 29を使うと、Si の全エネルギー 27 「Muffin-Tin (MT)球半径」と呼ばれる。L/APW 法では、ある半径(RMT)内の波動関数を原子の 波動関数とし、その外側を平面波で補強します。MT 球表面で両者の波動関数が連続的につな がる条件からエネルギー固有値が求まります。RMT を小さくとりすぎると原子の波動関数が MT 球内に収まりきらなくなり、逆に大きくとりすぎると、本来は平面波で表現しないといけ ない状態まで原子の波動関数だけで表現することになるので計算精度が悪くなります。そのた め、L/APW 法では、RMT を適切に選ぶことが重要です。w2web では、これを自動的に決める 機能ももっていますが、最終的には計算をしている本人が判断しなければいけません。 28 平面波の数も重要な計算パラメータです。打ち切りエネルギー (cut off energy) Ecut = h 2 k k2 max / 2me （波数ベクトル k を持つ平面波のエネルギーは h 2 k 2 / 2me ですね）で代表され、 29 論文等にも明記する必要があります。WIEN2k では Ecut の替わりに、最小の MT 球半径 R との 積 Rkmax を原子単位で入力するようになっています。これも論文等では Ecut に直しましょう。 DFT は(2)式の VXC (交換相関相互作用)などの具体的な形を決める方法を与えません。そのため、 前提や近似により多くのモデルがあります。DFT で単に「密度汎関数」というときは VXC のモ デルのことを指します (VXC も電子密度の汎関数です)。大別して、各空間座標の密度だけで VXC を表す局所密度近似 (Local density approximation: LDA)と、電子密度勾配を取り入れて非局所化 した一般化密度勾配近似 (Generalized gradiant approximation: GGA)があります。GGA にはさら に、PW91, PBE96 などが提案されています (名前は、提案者のイニシャルと報告年をとってい は-1160.141 Ry 30 が得られます。これは、基本格子に含まれる Si 原子 2 つ分のエネ ルギーです。1 原子あたり 8 keV になりますから、普段私たちが研究で接しているエ ネルギーから比べると数桁大きいことがわかりますね。これは、WIEN2k では L/APW 法を用いており、1s 軌道からのすべての電子軌道までエネルギーを加算しているため です31。私たちが普段測定するような物性量、例えば生成エネルギーなどと対応させ るには、化学反応式の両辺にでてくる物質の全エネルギーを計算し、そのエネルギー 差を求める必要があります。 例 1-1: Na の昇華熱 (Na (結晶) => Na (原子)のエネルギー差) (VASP,BPAW, PBE) ※ Na(結晶)の全エネルギー : E = -2.6203 eV/cell ※ Na(原子)の全エネルギー : E = -0.0007 eV/atom ※ Na(結晶) => Na(原子)のエネルギー差: ΔE = 1.3094 eV = 126 kJ/mol ※ RT = 2.49 kJ/mol (300 K)を足してエンタルピーにする: ΔH = 128 kJ/mol ※ 文献値: 108 kJ/mol 32 例 1-2: NaCl の生成エネルギー ※ NaCl(結晶) の全エネルギー: E = -27.2610 eV/cell (4NaCl) ※ Na(結晶)の全エネルギー : E = -2.6203 eV/cell (2Na) 33 ※ Cl2(分子)の全エネルギー : E = -3.5504 eV/cell (2Cl) ※ Na(原子)の全エネルギー : E = -0.0007 eV/atom ※ Cl(原子)の全エネルギー : E = -0.0183 eV/atom 34 ※ 凝集エネルギー NaCl(結晶) => Na(原子) + Cl(原子): -6.7962 eV/NaCl = -655.7 kJ/mol 文献値 641 kJ/mol ※ 生成エネルギー NaCl(結晶) => Na(結晶) + 1/2 Cl2(分子): -3.7301 eV/Na = 359.9 kJ/mol 文献値 411 kJ/mol 35 例 1-3: Si の凝集エネルギー (結合解離エネルギー) 上の例と同様、Si(結晶)と Si(原子)の SCF 計算を行い、全エネルギー差をとります。 しかし、WIEN2k, VASP の両方で Si(原子)の計算が収束しなかったので、ごまかしま す（本当はいけないやり方）。VASP で PAW 法, PBE を使うと、原子の全エネルギー がほぼ 0 eV になる (例 1-2 参考) ので、Si(原子)のエネルギーを 0 と置き換えます。 30 31 32 33 34 35 ます)。 Ry は Rydberg の原子単位で、エネルギーについては 1 Ry = 13.6 eV (水素原子の 1s 軌道エネル ギー)。他の量子計算プログラムでも出力、入力が原子単位になっているものが多く、また、 Hartree の原子単位も使われているので注意が必要です (1 Hartree = 27.2 eV)。 エネルギーの浅い価電子状態しか考慮しない擬ポテンシャル法では、全エネルギーの計算値は かなり小さくなります。 例えば VASP (PAW, PBE)で Si の計算をすると-10.844 eV が得られます。 『バーロー物理化学』より。 分子や原子の全エネルギーは、大きい単位格子内に一つの原子・分子を入れた構造をつくり、k 点をΓ点だけに限定した計算から求めます。 生成エネルギー: 化合物を構成元素の単体から作るのに必要なエネルギー 凝集エネルギー: 0K の凝集体を構成する原子に分けてばらばらにするのに必要なエネルギー 格子エネルギー: 0K の結晶を構成要素(原子・イオン・分子)に分けてばらばらにするのに必要な エネルギー。 凝集エネルギーより精度が非常に悪い。Cl2 の結合解離エネルギーのずれが大きいので、Cl2 分 子の計算モデルが適当でないようです (これ以上はまじめに計算していません)。 ※ Si(結晶)の全エネルギー: E = -43.3748 eV = -5.4219 eV/atom = 524 kJ/mol Si(原子)の全エネルギーを 0 と置き換えたので、これが凝集エネルギー 文献値 446kJ/mol 36 ※ 結合解離エネルギーは、これを結合数で割ればよい。Si 原子１つあたり 4 つ の Si-Si 結合をもつが、結合は２つの原子で共有されている。つまり、原子１ つあたりには２つ分の結合エネルギーがあるので、2 で割ればよい。 Si-Si の結合解離エネルギー: E = 262 kJ/mol 文献値 224 kJ/mol 37 2) 単純な結晶の安定構造（構造緩和計算）と体積弾性率 さて、全エネルギーが求まるということは、それが最小になるような格子定数や原 子座標を探すことで再安定な構造を計算できることを意味しています。このような計 算を構造緩和計算と呼びます38。ただしここで注意が必要なのは、DFT で求められる 電子状態は基底状態であること、原則として電子分布は絶対零度の Fermi 分布を仮定 しますし、原子核の運動は考慮していないので、得られた結果は、「基底状態、絶対 零度」での構造になります。 -551.424 -1160.137 -1160.139 -551.425 Exp.(RT) aC = 0.5431 nm VR = 270.5 a.u.3 (primitive cell) -1160.140 Opt. aC = 0.5472 nm VR = 276.67 a.u.3 -1160.141 -1160.142 図 3 -551.427 -551.428 Opt. ac = 0.4254 nm VR = 130.0164 a.u.3 -551.429 -551.430 -1160.143 -1160.144 Exp.(RT) ac = 0.422 nm VR = 126.9 a.u. 3 (primitive cell) -551.426 Energy / Ry Energy / Ry -1160.138 260 270 280 290 Volume / a.u. 3 300 -551.431 120 124 128 132 Volume / a.u. 3 136 140 単位格子体積を変えながら計算した Si (左)と MgO (右) の全エネルギー変化 (WIEN2k, PBE96) Si の例を見てみましょう。図 3(左)のように、V0 = 276.69 a.u.3 39に極小を持つエネル ギー曲線が得られます。横軸の体積は基本格子である菱面体格子のものですので、計 算しやすいように、4 倍してブラベー格子である立方晶になおしてから格子定数を求 めてみましょう。a = (4V0)1/3 = 0.5472 nm が得られます。文献値は 296 K で 0.5431 nm ですから、計算誤差は 0.8%にすぎません。ただ、実験値は 296K での値ですので、0 K での外挿値 (ac = 0.5426 nm) と比較するのが自然ですが、いずれにしても大きな違い 2 はありません。また、全エネルギーE と体積 V は、 E = E min + B0 (V /V0 ) の関係があり ますから、極小点の曲率から体積弾性率が B0 = 87.6 GPa と求まります。実験値は 298K で 97.88 GPa ということですから、11% 程度の誤差になります。 36 37 38 39 『理化学辞典』より。 化学便覧基礎編 II(第 3 版)より: (CH3)6Si2 で測定した値。 収束の安定性確保や計算時間の関係で、格子定数を固定して緩和計算がされることも多い。こ れと区別して、すべての格子定数まで最適化する計算を Variable-cell relaxation、Full relaxation と呼ぶことがあります。 ここでも原子単位が使われています。長さの単位は Bohr 半径 a0 = 0.05292 nm で、単位とし て ’bohr’ も使われます。 MgO の場合は、V0 = 130.02a.u.3 に極小を持ちます。MgO も FCC 構造を持ち、これ は菱面体格子の体積なので、Si の場合と同様の手順で a = 0.4254 nm が得られます。 文献値は室温で 0.422 nm ですから、この場合も計算誤差は 0.8%にすぎません。体積 弾性率は B0 = 157 GPa であり、室温の実験値 162 GPa と 3%程度の誤差で一致します。 ここでもうひとつ、量子計算の特徴と注意を書いておきます。図 3 をみてわかるよ うに、計算で得られるエネルギーは 1000 Ry 程度 (重原子を含む場合はさらに桁で大 きくなる)ですが、力学物性に効くのはせいぜい mRy 以下です。安定構造や物性を正 確に計算するには、量子計算には 8~10 桁以上の計算精度が必要といわれる理由です。 3) 一般的な構造緩和計算: 安定構造と局所構造 Si や MgO の結晶構造は対称性が高く、自由パラメータが格子定数の a だけしかあ りません。そのため、a を横軸にとった図 3 のようなグラフだけで最安定な構造を決 められますが、一般の結晶では自由パラメータが多くなり、グラフからは最安定な構 造を決められなくなります。その場合は、繰り返し逐次計算による最適化40によって 最小エネルギーを実現する自由パラメータを探し、最安定な構造を決めます。 表 1 に代表的な結果をまとめていますが、ほとんどの場合に 1~3%以内の精度で構 造パラメータが決まっていることが判ります。現在の DFT では、基底状態の電子構 造と結晶構造については、この程度の精度が得られます。逆に、これよりも実測値と の一致が悪い場合は、構造モデルが間違っている、計算条件 (Ecut、k 点数41、RMT な ど)が不適切、密度汎関数が適切でない、スピン分極を考慮しないといけない、など をまず疑います。これらが問題ない場合、例えば強相関電子系であるなど、DFT を適 用することの妥当性について再検討する必要があります42。 表 1 種々の結晶の構造パラメータ。カッコ内は緩和計算で得られた値 (VASP, PAW, PBE)。括弧 外が文献値。z(M)などは、M 原子の座標のうち自由パラメータであるもの。軸長の単位は Å。 Al (FCC) a = 4.04975 (4.0462) Ca (FCC) a = 5.5884 (5.51942 ) Mg (HCP) a = 3.2094 (3.1869) c = 5.2103 (5.19778) Na (BCC) a = 4.235 (4.20437) Si a = 5.41985 (5.46631) GaAs a = 5.65359 (5.7605) GaN (wurzite) a = 3.186 (3.24541) c = 5.176 (5.28965) z(N) = 0.375 (0.375783) NaCl a = 5.62 (5.65062) MgO a = 4.2109 (4.23617) CaO a = 4.8112 (4.83784) ZnO a = 3.2427 (3.25452) c = 5.1948 (5.21411) z(O) = 0.3826 (0.3816) In2O3 a = 10.117 (10.0316) SnO2 a = 4.738 (4.71537) c = 3.1865 (3.18356) β-Ga2O3 a = 12.23 (12.026) b = 3.04 (2.9927) c = 5.8 (5.7185) β = 103.7 (103.86) InGaO3(ZnO)1 a = 3.299 (3.29491) b = 5.714 (5.70415) c = 26.101 (25.4037) a = 11.989 (12.0284, 11.997, 11.9884) 12CaO·7Al2O3 (C12A7) α = 90 (α=89.9895, β=89.9334, γ=89.9619) 40 41 42 共役勾配(Conjugate gradient: CG)法、最速降下(Steepest descent: SD) 法などの最適化アルゴリズ ムや、熱運動を抑えた分子動力学(dumped MD)法などが使われます。 バンド計算では第一 BZ 内の異なる k ベクトルについて波動関数を計算し、全電子密度などを 求めます。この際に使う k 点の数を「k 点数」と呼ぶことがあり、精度に重大な影響を与えます。 実験の解析やモデルの方が間違えていることを見つけることもあります。DFT の結果と実験に 大きな違いが出た場合、新しい知見を得るチャンスかもしれません。 表 1 のうち、C12A7 については、文献にある結晶構造は立方晶ですが、緩和計算の 結果は三斜晶43になっています。これは、前者が X 線回折による構造解析で得られた もので、平均構造として立方晶となっているためです。実際の C12A7 では、12 個あ る等価な原子位置のうち、2 つを O2-イオンが統計的に44占有しています。そのため、 これらの O2-イオンの位置を決めて構造緩和計算をすると、対称性の低い正方晶や三 斜晶45になります。このように、固溶体46や乱雑構造をもつ材料の、実験的に知るこ とが難しい局所構造を知ることができるのも、構造緩和計算の大きな利点です。 ただし、多形47や相転移が近傍にあるような微妙な構造では、計算の結果には注意 が必要です48。たとえば鉄の場合には、LSDA49で計算すると FCC、非磁性構造が安定 になるとの結論が出ますが、実際には体心立方構造で強磁性です。これは GGA50を使 うことで正しい結果がでます51。また、ルチル構造の TiO2 においても、GGA の一種 である PBE を使うと安定格子体積が大きくなり、強誘電相が安定になってしまうと いう問題があります52。鉄の場合には GGA、ルチルの場合は LDA で正しい結果が得 られたわけです。必ずしも L(S)DA よりも GGA の方が正しいとは限りませんので、 注意が必要です。 4) 原子間ポテンシャル: MD 法との組み合わせ -4.7835 Total energy / 10-15 J/unit cell Total energy / 10-15J/unit cell -4.772 -4.774 -4.776 -4.778 -4.780 -4.782 -4.784 -4.786 -4.7845 -4.7850 -4.7855 0 0.2 0.4 0.6 Lattice parameter / nm 0.8 1.0 a=0.4nm -4.7860 -4.7865 a=0.39nm a=0.44nm -4.7870 a=0.421nm -4.7875 -4.788 図 4 -4.7840 0 0.02 0.04 0.06 0.08 0.10 0.12 0.14 Ionic displacement / fractional coordinate MgO について第一原理計算で求めたエネルギー曲線 (シンボル) と、それから求めたポテ ンシャル関数から求めたエネルギー曲線 (実線) 53 経験法の問題である「計算結果が経験パラメータの選択に大きく依存する」に対し 軸長 a, b, c がすべて異なり、軸角α, β, γ も 90o でないブラベー格子をもつ結晶。 「ランダムに」を格好よく表現しているだけです。 45 軸長が a = b, 軸角がα = β = γ = 90o のブラベー格子をもつ結晶。 46 半導体工学では「混晶」、金属工学では「不規則格子」と呼ばれる。 47 polymorph。同じ化学組成で結晶構造が異なるもの (ルチル構造とアナターゼ構造の TiO2 など)。 積層順序だけが違うものを特に「多型 (polytype)」という (閃亜鉛鉱構造とウルツ鉱構造の GaN など)。 48 計算精度にも注意が必要です。Ecut や k 点数を増やし、判断しようとしているエネルギー差よ りもこれらの効果が小さくなることを確認する必要があります。 49 Local Spin Density Approximation。DFT のもっとも簡単な形式の一つである局所密度近似 (LDA) をスピン分極計算に適用し、スピン密度の局所密度汎関数を VXC として計算する方法。 50 電子密度の勾配を考慮することで、L(S)DA の欠点を改良した方法。 51 C. S. Wang, B. M. Klein, H. Krakauer, Phys. Rev. Lett. 54 (1985) 1852; T. C. Leung, C. T. Chan, and B. N. Harmon, Phys. Rev. B 44 (1991) 2923. 52 B. Montanari and N.M. Harrison, Chem. Phys. Lett. 364 (2002) 528. 53 T. Kamiya, Jpn. J. Appl. Phys. 35 (1996) 3688 43 44 て、第一原理計算で原子配置の違う複数の構造で全エネルギーを計算し、MD 法など に必要な原子間ポテンシャルを決める方法があります。図 4 は、MgO についてイオ ン間ポテンシャルを求めたものです (CRYSTAL88, HF)。 5) アモルファス物質の構造を決める: 古典 MD 法でモデルをつくる 4) は、第一原理計算の結果を古典 MD 法に応用する例でしたが、逆の利用の仕方 もあります。結晶では有限個のパラメータで原子の位置を決められ、これらを回折法 などの構造解析で決定できます。しかしながら、アモルファス物質の場合は、電子密 度の動径分布関数のように、平均的な構造情報しか実験では得られません。 そこで「局所構造を決められる」という第一原理計算の利点が活きるわけですが、 構造緩和計算では、初期構造モデルの近くで全エネルギーが極小値を持つ構造 (local minimum) に収束してしまいます。そのため、なるべく実際に近い初期構造モデルを 作らないといけません。有限温度の動力学シミュレーションができる MD 法で、モデ ル構造を一旦融解させ、再固化させてアモルファス構造を作るのが直接的で単純な方 法です。しかし、第一原理 MD 法は計算が遅いため、時系列を追って溶融－再固化を シミュレーションすることはできません。そこでまず、同じ構成原子をもち構造が近 い結晶の構造や、アモルファス材料の実験データに合うように、大雑把に原子間ポテ ンシャルを決めます。これを用いて古典 MD 法で溶融－再固化シミュレーションを行 い、アモルファスの初期構造モデルを作ります。この構造から第一原理構造緩和計算 をすることで、アモルファス材料の安定構造や電子構造の計算ができます (図 5)54。 In Ga Zn In Ga Zn O O 図 5 MD 法と VASP (PAW, LDA)で求めたアモルファス InGaZnO4 の原子配列 (左)。84 原子を持つ 大きな単位格子55について計算した。右図は、最低非占有準位の波動関数 (正確には|ψ|2)。 6) バンド構造の読み方: フェルミエネルギー、バンドギャップ、間接遷移 ここでバンド構造を見てみましょう。でも、「わかりにくい」と評判なので、先に 「バンド構造図の読み方」をまとめておきます。図 6 左に WIEN2k で計算した Si の バンド構造図を示しますので、比べながら読んでみてください。 1) 横軸は波数ベクトル k を表す。ただし k をベクトル成分で表記してもわかりにく いので、対称性が高い k 点にはそれぞれ記号が定義されています56。実は、記号の 54 55 56 K. Nomura et al., Phys. Rev. B 75 (2007) 035212. 材料シミュレーションでは、基本となる単位格子を複数個つなげて大きい格子の基本構造をつ くることが多い。これを super cell と呼びますが、結晶学ででてくる本来の超格子 (superlattice) や、薄膜などででてくる人工超格子 (artificial superlattice)とは英語語源が異なります。 例えば、bilbao crystallographic server (http://www.cryst.ehu.es/) を見てください。この Web ページ 2) 3) 4) 5) 6) 7) 8) 意味がわからなくても重要な情報 (バンドギャップ、間接遷移など) は得られます。 縦軸は電子のエネルギーを表し、下へ行くほど電子のエネルギーが深く、安定に なります。 エネルギー固有値は波数ベクトル k の連続関数 E(k)になるため、複数の複雑な曲 線で描かれています。 波数ベクトル k をもつ電子は、飛び飛びのエネルギーしかとることができません57。 これは原子の中の電子と同じですね。ただ、原子の中の電子は、たとえば 1s 軌道 と 2s 軌道の中間のエネルギーを持つことはできませんが、結晶では k によってエ ネルギー固有値が変わるため、 「中間のエネルギー」も連続的にとることができる ようになります。エネルギー固有値は k の連続関数 E(k)になり、これを「エネル ギーバンド」と呼びます。 スピン分極を考慮しない場合は、一本58のエネルギーバンドには、単位格子あたり up スピンと down スピンの 2 つの電子が入ることができます。スピン分極したバ ンドの場合は、それぞれ 1 つずつです。 大雑把には、逆空間における k を、電子が進む方向とみなせます59 特に明記されていない場合、縦軸のエネルギーの原点はフェルミエネルギー (Fermi energy: EF)とされています。つまり、0 K において、0 eV より下の軌道は電 子が詰まっている被占有準位で、それより上の軌道は原子が無い非占有準位です60。 金属であるか、絶縁体（半導体）であるかは、EF がバンド構造図の曲線 E(k)と交 差する場所があるかどうか、あるいは後述の状態密度が EF で有限の値であるかど うかで判断できます。 上でまとめたように、図 6 のようなバンド構造図の横軸は、逆空間（第一 BZ）内 の波数ベクトル k を表しています。ただし、3 次元に拡がる波数ベクトルのすべてを − から”KVEC”を選び、空間群の番号(Si の空間群は Fd 3 m で、空間群番号は 227 です)を入力する と、図 6(右)のような第一 BZ の絵と、対応する点の記号、座標が表示されます。これでもまだ わかりにくいと思いますが、WIEN2k の場合は、XCrySDen (http://www.xcrysden.org/)をインスト ールすることで、バンド構造を計算する k 点を GUI を使って選ぶことができます。 57 「結晶の固有状態としては」そのエネルギーをとりえない という意味。電子は結晶中でも原子 による散乱を受け、進行波成分は減衰します。ところが、E(k)のエネルギーをもつ場合に限り、 3 次元周期構造による干渉が散乱波を打ち消し、進行波だけが残ります (薄膜による光の干渉と 同じ。表面反射があっても、干渉条件を満たす場合は透過率が 100%になります)。そのおかげ で、全く散乱を受けずに進行することができ、これが「エネルギーバンド」になります。これ 以外のエネルギーでも波動関数は存在しますが、波数ベクトルが虚数になります。そのような 電子の進行波成分は、原子による散乱によって波数の虚部に反比例した減衰長で減衰します。 厚さが減衰長程度の薄膜の場合 (例えばトンネル磁気抵抗素子) にはこのような減衰状態も重 要になります。 58 第一 BZ 全体で「一本」であり、バンド構造図の一区画についての「一本」ではありません。 59 k は運動量のようにふるまいますが、運動量保存則を満たさない(ブラッグ条件 k’ = k+G )こと hkl から理解できるように、厳密な意味での運動量ではありません。 60 金属の場合は EF がバンドを横切るので EF は一意的に決まります。絶縁体や真性半導体ではバ ンドギャップ中のどのエネルギーもこの条件を満たすので一意的に決まりません。半導体物理 で使われる EF は化学ポテンシャルの意味があり、Fermi 分布関数と電荷中性条件から決まり、 真性半導体ではほぼバンドギャップの中央に来ます。バンド計算では化学ポテンシャルとして EF を扱うことが少ないため、単純に、 「絶対零度での電子配置における最高の電子のエネルギー」 をフェルミエネルギーと呼ぶことも多くあります。そのため、VBM をエネルギー原点にとった バンド構造図が良く使われます。 表示するのは無理がありますから、結晶格子の中でも対称性の高い方向についてだけ 抜き出して描きます。横軸の記号、W, L, Γ, X, K は逆空間における対称性の高い点（実 空間における「方位」に対応します）、Λ, Δは対称性の高い逆空間の方位を表してい る (図 6 右) のですが、上記脚注のように、データベースなどで調べられますから、 この記号を覚える必要はありません。ただ、Γ点は逆空間の原点、つまり、k = (0, 0, 0)、 X 点は x 軸方向の BZ 境界 (単純格子の場合は (1/2, 0, 0) 61) を表しているということ くらいは覚えておきましょう (脚注62を読むと重要性がわかります)。 電子のエネルギー (C) (A) W 図6 ↑0.53 eV (B) L Λ Γ Δ L WK 波数ベクトル W L Λ Γ Δ L WK 波数ベクトル Si のバンド構造。(左) DFT で計算した結果。(中) 非占有状態を平行移動させてバンドギャ ップ値を実測値にあわせたもの。(右) 第一 BZ と対称性の高い k 点の位置と記号。 図 6(左)をエネルギー軸の方向から見てみましょう。}(A) と }(C) で示した領域は どのエネルギーでも E(k)曲線にかぶりますが、}(B) の領域では、そのような E(k)曲 線がありません。つまり、どのような k をもつ電子でも、}(B) の範囲のエネルギー をとることはできないことを意味しています。このようなエネルギー範囲を禁制帯 (forbidden gap) やバンドギャップ (bandgap) とよびます。 この図から、シリコンはΓ点（k=(0,0,0)）に価電子帯上端（Valence Band Maximum: VBM）をもち、Γ-X 軸の途中に伝導帯下端（Conduction Band Minimum: CBM）を持つ ことがわかります。バンド構造図が便利なのは、k 選択則63を理解しやすいことです。 この場合は VBM を形成している電子準位の k 点と CBM を形成している k 点が異な りますので、最小のバンドギャップでは光学遷移は非常に弱く64、 「間接遷移型半導体」 65 であることがわかります 。間接遷移が起こるためには、運動量保存則をみたすため、 61 62 63 64 65 一般に k 点の座標は、逆格子ベクトル a*, b*, c*を単位とした内部座標で表されます。結晶学や 回折学では逆格子ベクトルは a* = b × c / V のように定義されます (V は単位格子体積) が、物性理 論やバンド理論では a * = 2π b × c / V と、係数が 2π だけ違うことに注意してください。この余計 な係数のおかげで BZ 境界の k を 1/2 のような単純な数値で表現できるようになります。 LCAO での結晶軌道は Bloch の定理を使って e ik ⋅ri φ (r − ri ) (φ (r - ri)は原子の波動関数、ri は i 番目の単位格子の座標)で表されますが、Γ点では L + φ (r − a ) + φ (r ) + φ (r + a ) + L と簡単にな ります(a は格子ベクトル)。つまり、それぞれの単位格子内の原子の波動関数を単純に足したも ので、φが s 軌道であれば結合軌道になります。X 点は(1/2, 0, 0)ですから、a 軸方向の波動関数 だけを抜き出すと、L − φ (r − a ) + φ (r ) − φ (r + a ) + L になり、交互に係数が+1,-1 と変わること が判ります。これは分子軌道でいうところの反結合軌道です。 光による電子遷移 (双極子遷移) は同じ k 点の間でだけ起こるという遷移の選択則。 一次摂動 (双極子遷移) での遷移確率はゼロですが、高次の摂動やフォノンとの相互作用によ って遷移が起こります。ただし遷移確率は非常に小さくなります。 実験的に区別する場合には、吸収スペクトルα(E)が α (E ) = (E ± hω ph − E g )2 / ω 2 (間接遷移)であ るか α (E ) ∝ (E − E g )1 / 2 / ω 2 (直接遷移)であるかを利用することができます。 ∑ 波数ベクトルの差に等しいフォノンを吸収あるいは放出する必要があり66、光学吸収 係数が非常に小さくなります。また、発光ダイオード(Light-emitting diode: LED)では、 伝導帯に電子を、価電子帯に正孔を注入し、これらが再結合する際に出てくるエネル ギーを光として放出させますが、間接遷移型の場合はほとんど発光しない67ので LED の発光層としては使えません。バンド構造によって、このようなことまで読み取れま す。図 7 に WIEN2k (PBE96)で計算した典型的なバンド構造をまとめますので、それ ぞれ、金属であるか半導体であるか、直接遷移型であるか間接遷移型であるか、考え てみてください。 6 (B) (A) (C) (D) エネルギー / eV 4 2 0 -2 -4 -6 Γ 図7 M K Γ AW L Γ X WK W L Γ X WK W L Γ X WK (A) Mg, (B) Si, (C) GaAs, (D) TiN のバンド構造。 7) バンド構造の読み方: バンドギャップ問題、real state と virtual state さて、図 6 左から読み取れる最小のバンドギャップ（間接バンドギャップ）は 0.6 eV でしたが、これは、300 K での実験値 1.1242 eV よりずいぶん小さいですね。絶対零 度への外挿値ではさらに大きく 1.17 eV になりますから、温度のせいではありません。 これが DFT の最大の問題と言われる、いわゆる「DFT のバンドギャップ問題 (Bandgap problem)」です。DFT ではバンドギャップをかなり過小評価しますので、バンドギャ ップの定量的な議論には使えません。逆に HF 近似ではバンドギャップを過大評価す る「バンドギャップ問題」があり、第一原理計算でバンドギャップを正確に計算する のは現在でも難しい問題です。 これは、非占有準位の意味を考えてみると理解できます。被占有準位というのは、 現実に電子がいる状態そのものです。ところが非占有状態というのは、量子理論の方 程式を数値的に解いた結果として現れた、数学的な解にすぎません。計算していると きには存在を考慮していない電子の状態ですから、その結果に物理的意味はありませ ん。そのため、被占有軌道を real state、非占有軌道を「仮想軌道 (virtual state)」と呼 ぶことがあります。ですから、物理的意味の無い仮想軌道とのエネルギー差であるバ ンドギャップの計算値が実測値に合わないからといって、DFT や HF 近似の精度が悪 いということではありません。バンドギャップと非占有軌道は DFT の前提外である 66 67 前の注釈で「k は運動量保存則を満たさないので厳密な意味の運動量ではない」と書きました が、Bragg 条件を満たす場合などを除くと、通常の運動量と同じ保存則が出てきます。 注入された電子や正孔が両方とも高いエネルギーをもって同じ k 点の状態にあれば、k 選択則 を満たした許容遷移で発光再結合しそうです。しかし、高エネルギー電子(ホットエレクトロン) は fs オーダーの短い時間でフォノンなどとの非弾性散乱によって一番エネルギーの低い状態に 緩和します。一方、許容遷移でも再結合時間は ps 程度になりますから、ほとんどの場合、再結 合は CBM と VBM まで緩和した後で起こります。 ということです68。 ただし、「仮想軌道」とはいっても、経験的に、非占有準位の構造が励起状態の電 子構造をかなり良く反映していることが確認されています。そのため、これを「伝導 帯」の電子構造として解釈するということが一般的に認知されています。「バンドギ ャップ問題」に対しては、実験値のバンドギャップを再現するように、非占有準位全 部を適当なエネルギーΔsc だけ平行移動してしまう、という対応がされます 69 (図 6 右)。 8) 有効質量 先に「渡した値の興味のあるバンドギャップは(ほぼ)自由電子(NFE)近似のものとは 違う」と書きましたが、これはあくまでもバンドギャップだけの話です。NFE は、よ り複雑な材料のバンド構造を理解するにも非常に有用ですが、それは、バンドギャッ プではなく、伝導帯や価電子帯の形を理解したり、それから「有効質量近似」という アイデアを提供してくれることにあります。自由電子バンドギャップは図 7 のぞれぞ れのバンド構造について、CBM と VBM 付近の E(k)は、2 次関数的になっていること に気づいたでしょうか？ LCAO に基づいてバンド理論を展開することも可能ですが 70 、その場合、一次元結晶でのエネルギー分散 E(k)は (3) E (k ) = ε 1 − 2 h12 cos(ka ) ~ ε 1 − 2 h12 + h12 a 2 k 2 で与えられます。これを自由電子の E (k ) = E0 + (h 2 / 2me )k 2 と比較すると、あたかも電 子が質量 m* = h 2 / 2 h12 a 2 を持っているとみなすことができます (図 8)。事実、半導体 では、結晶中の電子を「誘電率εr の真空中を運動する質量 m*、電荷-e の電子」、正孔 を「(同)電荷+e の電子」として扱っても良い近似となることが知られています。その ため、m*を有効質量と呼び、このような近似を「有効質量近似」と呼びます。 有効質量が小さいほど電子が動きやすいわけですが、LCAO 理論からは、「電子の 動きやすさと波動関数の直観的な関係」も得られます。(3)式から、有効質量を小さく するには共鳴積分|h12|が小さい材料を設計すればいいことがわかります。|h12|は隣り合 った原子の軌道関数とハミルトニアンの積分ですから、これを大きくするには、波動 関数同士の重なりが大きくなるような結晶構造をみつけるか、波動関数の拡がりの大 きい原子（一般には周期表の下の元素）に置き換えることを考えれば良いわけです。 図 8 には、SnO2 を例に有効質量の解析を示しています。SnO2 は TiO2 と同じルチル 構造を持ち、正方晶に属し、Γ-X とΓ-Z 方向で異なる有効質量をもちます。図 8 右は Γ-Z 方向の CBM バンドについて、E(k)の 2 階微分係数から計算した m*(E)を重ねてプ ロットしています。これから、Z 方向、つまり c 軸方向の有効質量 me//c はΓ点で 0.182me と得られます。文献値は 0.234me なので少し過小評価しています。原因の一つは、イ オン性結晶の有効質量が、電子の電荷と光学フォノンが相互作用して重くなることと 68 非占有軌道の計算精度を上げるには、ここに電子が入った場合、つまり、励起状態を考えて計 算しないといけません。CI, 時間依存 DFT (time dependent DFT: TDDFT) 法, GW 近似などが必要 とされますが、計算時間が膨大になることや、TDDFT のように適用可能範囲がまだわかってい ないものもあり、まだ実験研究者が手軽なツールとして使うほどは簡単ではないようです。 69 物理学者というのは、普段は厳密に数式を扱っているように思えますが、必要なときには極端 な単純化や、非常に荒っぽい手段も平気で使います。ただし、このΔsc に「Scissors operator」と いう格好良い名前をつけて、高尚なイメージを与えることも忘れません。 70 化学をベースにしている人や化合物半導体を扱っている人の場合は、自由電子理論よりもむし ろわかりやすいのではないでしょうか 考えられます71。そのため、バンド構造72から求めた有効質量を「バンド有効質量」 と呼んで区別することがあります。また、k がΓ-Z 方向へ変化するにつれ、有効質量 が増加して発散し、その後で負の値になっています。これは E(k)の変曲点で曲率がゼ ロになるためで、これを境に m*が正の側を「電子バンド」、負の側を「正孔バンド」 と呼ぶことがあります。 2 LCAO バンド 1 5 -5 -10 m* -1 0.5 4 0 3 2 m* 1 -15 0 -0.5 図8 -0.3 -0.1 2π/a 0.1 0.3 0.5 Γ Z A M Γ X R Z -2Γ0 Energy / eV 1.5 1 0 m* / me Energy / eV Energy / |h12| Energy / 2|h12 | 自由電子 6 自由電子 5 2 7 Z 0 (左) LCAO 法により求めたエネルギーバンドと自由電子近似したバンド。h12 は共鳴積分で、 バンド幅は 4|h12|になります。(中) SnO2 のバンド構造 (VASP, PAW, PBE)。(右) SnO2 のΓ-Z の CBM エネルギーバンド(●)、自由電子近似バンド(実線)、電子の有効質量(◇)。 有効質量を見積もる簡単な方法もあります。一つのエネルギーバンドの最低エネル ギーと最大エネルギーの差を「バンド幅73 (W)」と呼びますが、(3)式から W = 4|h12| になることがわかります。単純なエネルギーバンドであればこの値はバンド構造図か らすぐにわかりますので、 m* = 2h 2 / Wa 2 から計算できます。SnO2 の例では、バンド 幅 W = 0.545 eV、格子定数 c = 0.3227 nm ですから、m//* c = 2h 2 / Wc 2 = 0.268me が得られ ます。E(k)から計算した値より結構過大評価しますので、大事なデータの場合はちゃ んと E(k)曲線から求めましょう。 9) 状態密度 図 9 左は状態密度74 D(E) と呼ばれるもので、Si について計算したものです(図 6 左 図と対応させて見てください)。横軸に電子のエネルギー、縦軸に、その付近のエネ ルギーを持つ電子が何個存在するか75をプロットしたものです。縦軸の単位は [状態 数/eV] なのでちょっとわかりにくいので、右図をみてください。これはエネルギー –15 eV から E までの準位を占有している電子の総数、積算状態数 N(E)です。これか ら、-15 ~ 0 eV までを占めている電子数 N(0) は 8 個であることがわかります。図 6 のバンド構造図だけからでは判断できませんが、}(A) の-12eV から 0eV のバンドは、 Si の 3sp3 混成軌道のバンドです。もともと Si 原子 1 つには、(3s)2(3p)2 の 4 つの電子 があったわけですが、N(0) = 8 というのは、基本格子に含まれる 2 つの Si 原子にある 71 電子とフォノンが結合して重い電子の様にふるまうもの(準粒子)をポーラロンといいます。多 かれ少なかれ、極性結晶中の電子はポーラロンとしてふるまいます。 72 結晶格子が格子振動などで変形することを考慮しないので、「rigid band model」と言います。 73 この言葉も、分野によって違う言葉で使われます。光電子分光などでは「一本のエネルギーバ ンド」を区別できないため、エネルギーの一番高い光電子シグナルのピーク全体の幅を「バン ド幅」ととることがありますので、混同しないようにしましょう。 74 「状態」とは state の和訳で、電子準位と同一と考えてください。電子数と言い換えた方がわか りやすいのですが、非占有準位には電子がいないのでこのような言い方になります。 75 E ~ E+δE の範囲にある状態の総数 N(E+δE) – N(E) が D(E) δE と定義されます。 3s, 3p 電子の総数になっているわけです。このようにして、図 9 右の N(E) が積算電 子数であることは納得してもらえたでしょう。図 9 左の D(E) というのは、N(E) のエ ネルギー微分、つまり D(E) = d N(E) / dE になります。実際に電気伝導度や磁化率な どを決めるのは状態密度 D(E)ですから、論文などでは D(E)だけを示します。 状態密度を各原子ごと、あるいは原子軌道ごとに分解すると、化学結合に関するイ メージを作ることができます。図 9 左では、結晶中の全状態密度(Total)、Si 原子の全 状態密度(Si total)、Si p 軌道だけの状態密度(Si p)と s 軌道だけの状態密度(Si s)をそれ ぞれ描いてあります。このような図を、投影状態密度 (projected density of state: PDOS) あるいは部分状態密度 (partial density of state: やはり PDOS)と呼びます76。 これから、0 eV より下の領域、つまり価電子帯は、p 軌道が主で、それより上の伝 導帯は、p 軌道と s 軌道が同程度の割合で混成していることが読み取れます。これは、 Si が四面体状の結合を作っており、sp3 混成軌道で化学結合が形成されていることを 示しています。これは、化学結合論で CH4 分子などを例に説明されるのと同じです。 ただし、実際に Si について DFT で厳密に計算してみると、様子が違うことにも気 がつくでしょう。まず、全状態密度は Si 原子の全状態密度の 2 倍近い値があります。 WIEN2k で計算される PDOS は、MT 球内の波動関数の積分値であり、MT 球外の波 動関数は、”Interstitial”として扱われます。つまり、”Total”と”Si total”の差は Interstitial 状態です。Si 結晶がもつダイヤモンド構造というのは、原子の充填率がもっとも低い 構造なのですが、これが、Interstitial の波動関数の割合が多い原因の一つです。他の 重要な理由として、Si が良い半導体であり、価電子帯側にある電子（つまり正孔）も 伝導帯の電子も、波動関数が原子の外まで拡がっていることを反映しています。特に この傾向は伝導帯で強く、3sp3 混成軌道だけでなく 4s 成分も入って s 軌道の混成が大 きくなります。この意味で Si の伝導帯は”s 性が強い”と言われます。 14 Total 1.5 12 DOS 1 N(E) 10 Si total 0.5 6 4 Si p Si p 0 -15 図9 8 Si s Si s -10 -5 0 Energy / eV 5 2 0 -15 -10 -5 Energy / eV 0 5 Si の(左) 状態密度 D(E)と、(右) 積算状態数 N(E)。 分子軌道法と比較すると、バンド計算で得られる PDOS の解釈にはより慎重になる 必要があります。分子軌道法のように LCAO を使っている場合、各原子軌道の分子軌 道への寄与は LCAO の係数で決まります。原子間の波動関数の寄与(非対角項)をどの 76 「投影」というのは数学用語で、例えば「ベクトル (0.7, 0.5, 0.2)の x 軸への投影ベクトルは (0.7, 0,0)である」というように、ある量を、ある特徴を持つ別の量の成分として表すときに使います。 ここでは、L/APW 法や PP-PW 法などで使われる平面波基底の波動関数を、特定の位置の原子 や、特定の対称性を持つ関数など (つまり s 軌道、p 軌道など) へ「投影」していることから、 このように呼ばれます。 ように各原子に分配するかは若干の任意性が残りますが、それほど大きな問題ではあ りません。ところが、平面波は空間全体に一様に広がっていますから、平面波基底の バンド計算では、原子ごとの PDOS を一意的に計算する方法がありません。そのため、 L/APW 法である WIEN2k は MT 球内の波動関数の積分値を原子の PDOS としていま す。VASP では MT 球がありませんので、各原子に半径を与え(Wigner-Seitz(WS)半径 と呼んでいます)、その中の波動関数の積分値を原子の PDOS として使っています。 Total DOS は変わりませんが、PDOS は MT 半径、WS 半径の選び方によって変わりま す。意味のある PDOS を得るためにはコツをつかむ必要がありますし、逆に、化学者 としての直観に合うように PDOS を故意に変えることもできますので、注意が必要で す77。 Absorption coefficient / 105 cm E1 E2 ε1 ε 80 70 60 50 40 30 20 10 0 -10 -20 -30 -40 -1 10) 光学スペクトル ε2 2 2.5 3 3.5 4 4.5 Photon energy / eV 5 7 6 5 4 3 measured@10K 125meV 2 1 calculated 0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 Photon energy / eV 図 10 (左) Si の複素誘電関数スペクトル。シンボルが実測値、実線が計算値 (CASPEP, PBE96)。 (右) LaCuOSe の光吸収スペクトル。実測値(上)と計算値(下) (WIEN2k (PBE96) + OPTICS)。 量子計算によって波動関数がわかると、摂動理論により誘電関数の実部スペクトル が計算できます。これを Kramers-Kronig 変換して虚部スペクトルも計算できます。た だし、遷移エネルギーの差が分母に入ってくるため、過小評価した DFT のバンドギ ャップをそのまま使うと絶対値に大きな影響を与えます。そこで、Scissors operator を用いてバンドギャップを実測値にあわせてから計算することで光学スペクトルを 計算します(図 10)。Si では、低エネルギーのε1 78の値は実験値とよく一致しています。 エネルギーが高く特徴的な構造が出る 3 eV 以上の構造の再現性はよくありませんが、 この原稿のために、えいやっ、と計算しただけの結果なので、きちんと計算すればも っとちゃんとした結果がでるかもしれません。LaCuOSe の場合も、矢印で対応させて いる特徴的な構造をよく再現できていることがわかります。ただし、DFT のエネルギ ー固有値はバンドギャップの補正をしても遷移エネルギーに対応しませんから、それ ぞれの構造のエネルギーはずれがあります。実測スペクトルで’125 meV’と書いて縦 線で示しているのは励起子79による鋭い吸収線です。DFT では励起子のような励起状 77 78 79 L/APW 法の場合、MT 球半径がとんでもない値だと、計算結果もとんでもなくなります。その ため、PDOS にもある程度の制限がかかりますが、VASP の場合には、電子構造の計算結果とは 全く関係なく WS 半径を決めることができますので、後者の方が”自由度”は大きくなります。 誘電関数 ε は、屈折率 n、消衰系数 κ と ε = (n - iκ)2 の関係があり、吸収係数 α は光の波長λ を用いて α = 4πκ/λ と求められます。 物質に光を入れると、価電子帯の電子を伝導帯に励起し、電子と正孔を作ります。バンド計算 では、このような電子も正孔も周期結晶中に無限に拡がった構造としてしか扱いません (もと 態は考慮していませんので、これに対応する構造は、計算した吸収スペクトルの肩と して見えているだけです。 11) 波動関数（電子密度） 量子計算をすると、波動関数も計算されます。ただし結晶軌道は、 Φ k (r ) = e ik ⋅r φ (r ) のように、Bloch の定理によって波数ベクトル k による位相 exp(ik･r)がかかります。 そのため、まともに Φ k (r ) の実数部と虚数部をプロットしても、k が 0 でない場合に は余計な位相項の変化が重畳して何がなんだかわかりません。そのため、その軌道の 電子密度になる Φ k (r ) をプロットすることがあります。 2 CBM Si GaN GaAs Ga Ga As N VBM N 図 11 異なる化学結合性状を持つ半導体の CBM と VBM の波動関数。図中の球は各原子、曲面 は波動関数の等高面。 図 11 には、Si, GaAs, GaN の CBM と VBM の波動関数 Φ k (r ) 2 (正しくは「波動関 数の自乗」です) を描いてあります。半導体の場合には、これらはそれぞれ、n 型半 導体での電子の通り道と、p 型半導体での正孔の通り道を可視化したものと捉えるこ とができます。さて、上の PDOS の項で見てきたように、Si の VBM は sp3 混成でで きていますが、p 軌道の寄与が大きくなっていましたね。それは図 11 でも、VBM 波 動関数が最近接の Si-Si 方向の結合上に大きい密度を持つことから見て取れます。一 方 CBM の方は s 性が強く、かなり丸みを帯びた形になっていますが、それでも、Si-Si 結合上に高い密度があることがわかります。Si-Si 結合の真ん中で波動関数の値がゼロ になる、いわゆる「ノード(節点, node)」がありますが、これは、CB が sp3 の反結合軌 道でできていることを反映しています。Si は一種類の原子だけで構成されているので 完全共有結合性ですが、GaAs や GaN では、Ga は電子親和力が小さくて陽イオンに なりやすいのに対し、As と N は Ga よりも電子親和力が大きく、陰イオンになりやす いですね。As と N を比べると、N の方が電子親和力が大きくなります。そのため、 Si => GaAs => GaN の順番で、化学結合のイオン結合性が強くなります。このことは、 もと周期構造を仮定しているため)、実際には、励起直後の電子と正孔は近傍の位置で局在して いる (相対論の「相互作用は光速以上で伝播しない」という要請を考えれば当たり前です) ため、 それぞれの電荷で引きつけあいます。このエネルギー(励起子の束縛エネルギー)が熱エネルギー よりも十分大きければ、LaCuOSe のように鋭い吸収線や発光線として観測されます。Si や GaAs の励起子束縛エネルギーは室温よりも小さいため、通常は問題にされません。 GaAs の VBM がほとんど As の軌道からできており、逆に CBM は As の軌道もあるも のの、Ga の s 軌道が強いことに現れています。GaN はこの傾向をもっと強くしてお り、CBM には N の軌道があまり現れていません。 全電子密度、価電子密度、up スピン密度、電位ポテンシャル分布なども同様に可視 化することで、さまざまな情報を得ることができます。 12) 3d 軌道の配位場分裂 (CASTEP, PBE96) ここでもう一つ、教科書の内容を第一原理計算で確認してみましょう。孤立してい る球対称原子では 5 つの d 軌道は縮退しています。結晶の中では、原子がおかれた環 境が球対称でないため、d 軌道がいくつかの組に分裂しますが、これを「配位子場分 裂 (ligand field splitting)」と呼びます。どのような組に分裂するのかは、原子がおかれ ている位置の点群とその点群の指標表からすぐにわかりますが、どのような順番で、 どれだけのエネルギーだけ分裂するかはわかりません。 Oh 点群の指標表 (一部抜粋) E 8C3 3C2 6C2 6C4 ’ A1g 1 1 1 1 1 Eg 2 -1 2 0 0 T1g 3 0 -1 1 -1 T2g 3 0 -1 -1 1 T1g 3 0 -1 1 -1 表2 自由 イオン Ti3d i 3 σ 6S4 6σ h d 1 2 3 3 -3 Oh配位子場 eg (4.3eV) 1 -1 0 0 0 1 2 -1 -1 1 1 0 1 -1 -1 1 0 -1 1 1 x2+y2+z2 80 (x2-y2,2z2-x2-y2) (Rx,Ry,Rz) (xy,xz,yz) (x,y,z) Ba Ti t2g (2.1eV) Ti 8S6 O バンドギャップ O Ba O2p (0eV) 図 12 立方結晶 BaTiO3 のバンドギャップ近傍の電子状態と波動関数。VBM は O2- 2p 軌道、伝導 帯は Ti4+ 3d 軌道でできており、3d 軌道が配位子場分裂で 2 組に分裂しています。 立方晶 BaTiO3 は図 12 に示すようなペロブスカイト型の結晶構造をもちます。ここ で大事なのは中心の Ti 原子で、上下左右に O2-イオンが配位しています。このような 対称性は Oh 点群に属しており、表 2 のような指標表で表されます。この右端の記述 から、d 軌道がどの組に分かれるかがわかります。ここで()で括っている組は、量子 計算では縮退した状態になる組です。ですから、(dxy, dxz, dyz)の組は T2g 表現に属して 三重縮退した軌道 (t2g 軌道などと書かれます) を作ること、 (dx2-y2, dz281)は Eg 表現に 属して二重縮退した eg 軌道を作ることが読み取れます。図 12 には CASTEP で計算し た波動関数を示していますが、一番エネルギーが高いのが eg 軌道で O2-方向に伸びた (dx2-y2 82, dz2)でできていること、その次が t2g 軌道で O2-のいない 8 つの方向へ(dxy, dxz, 80 81 82 x2+y2+z2 は r2 に等しく方位依存性がありません。つまり、球対称の意味です。 指標表では 2z2-x2-y2 になっていますが、x2+y2+z2 = r2 は球対称な成分なので、対称性の議論で は無視できます。2z2-x2-y2 は z2 と等価ということです。 dx2-y2 軌道は、図 12 では xy 平面の O2-の無い方向に向いていそうですが、これは原子軌道を表 dyz) 軌道が延びていることが判ります。VBM は O2- 2p 軌道でできていますが、ドー ナツのような形をしていて、よく知っている鉄アレイ型の p 軌道とは似ていません。 これは、2p 軌道が面内で縮退している (たとえば xy 面内の O2-イオンの 2px と 2py 軌 道が縮退している) ためです。縮退した波動関数を、例えば φ1 = x + y, φ2 = x – y とし てみると83、|φ1|2 + |φ2|2 = 2(x2 + y2) = 2 Rxy2 と、xy 面内の原点からの距離 Rxy だけの関 数になり、円対称になっていることが確認できます。このように、縮退した px, py 軌 道は面内で円対称になり、ドーナツ状の波動関数になります。 13) 欠陥とドーピング 量子計算の利点の一つは、実験では決定的なデータを得にくい問題に対して、補強 情報を与えてくれることです。そのような例に欠陥やドーピングの問題があります。 図 13 は、 結晶 Si を簡略化した Si5 クラスターについて、半経験分子軌道法ソフトウ ェアである MOPAC7 で計算したものです。ここで分子軌道法を使ったのは、バンド 計算で欠陥や不純物を扱うと非常に大きな super cell を組まなければならず、膨大な 計算時間がかかること、さらに、そうやって得られた結果は厳密であったとしても直 観的に理解しにくいことが多いためです。 4 LUMO(CBM) 2 Si5H12 Si5H11 Si5 HOMO (DB) 0 Si 3p -2 CB Eg Si 3s VB HOMO (VBM) Energy / eV Si CB -4 Si 3p VSi -6 -8 HOMO (DB) Si 3s VB -10 -12 -14 -16 -18 図 13 Si クラスターの電子構造 (MOPAC7)。Si (大きい球) を中心として 4 つの Si が配位した Si5 クラスターを基本に、外側の Si 原子の非結合軌道をすべて水素原子 (小さい球) で終端した Si5H12 クラスターと 1 つだけ非結合手を持つ Si5H11 クラスターについて比較している。両端には、Si の 電子構造と、Si 欠陥と非結合手があるときの電子構造を模式的に描いている。その内側は、HOMO, LUMO, ダングリングボンド(DB)に対応する波動関数。 図 13 中図のエネルギー準位図をみると、非結合手をすべて H で終端した Si5H12 ク ラスターは HOMO と LUMO の間に大きなエネルギーギャップがあるのに対し、非結 合手を一つ持つ Si5H11 や全く終端していない Si5 クラスターでは、Si5H12 クラスターの 現するときの座標軸と結晶の座標軸が違うためです。座標軸というのは人間が都合のいいもの を選んだだけですので、無限の任意性があります。量子計算を行うと、座標軸のとり方に関係 なく同一の計算結果が出てきます。この場合は、結晶の座標軸と Ti 原子の座標軸が xy 面内で 45o ずれているわけです。このような原子固有の座標軸を「局所座標」と呼びます。実際の計算 の際も、このように軸を取ることで計算を簡略化できることがあり、WIEN2k でも入力ファイ ル作成時に局所座標を計算しています。 83 縮退しているので、φ とφ が直交するなら、どうとっても構いません。一番簡単なのはφ = x, 1 2 1 φ2 = y ですが、簡単すぎるのでここの説明ではちょっと複雑な方で説明してみました。 エネルギーギャップの中央付近に余計な占有準位ができていることがわかります。こ れが「Si のダングリングボンド」で、HOMO(DB)と書いてある波動関数の絵から、こ れらが Si 上の H で終端されていない p 軌道であることがわかります。Si 中にこのよ うな状態がたくさんあると、光生成された電子－正孔対を非発光再結合させて LED や太陽電池の性能を落したり、トランジスタにゲート電圧を印加しても Fermi 準位が 動かない「フェルミ準位ピニング (Fermi level pinning)」を起こして Si デバイスの性 能を極端に低下させます。このような状況を改善するため、アモルファスシリコン中 のダングリングボンドは水素で、太陽電池の Si 表面・端面は酸化膜などを作って終 端し、ダングリングボンドを不活性化させています。 ただし、Si5H12 クラスターのエネルギーギャップの計算値は 10 eV 近くもあり、結 晶 Si の文献値 1.12eV とはぜんぜん合いません。DFT と同様の「バンドギャップ問題」 の寄与もありますが、それ以上に大きいのが、Si 原子 5 個だけの非常に小さなクラス ターで計算したためです。固体中の電子構造は、隣同士の波動関数が無限につながる ことでバンド分散をつくることはこれまで見てきたとおりです。特に良い半導体（有 効質量が小さい）ほどバンド分散は大きく、Si では 10 eV 近くあり、これが原子軌道 の上下に広がります。価電子帯、伝導帯とも、このバンド幅の半分ずつが上下に広が っているので、結局バンドギャップも、小さい Si クラスターよりも 10 eV 近く小さく なって、1.12 eV ほどになるわけです。実際に近い計算をするためには、クラスター のサイズを大きくすればいいのですが、Si の場合では、100 原子程度のクラスターで はエネルギーギャップが収束しないようです84。 E=ECBM EF 9 E = ECBM + 0.18 eV VB CB 8 P O 7 Fe DOS 6 P 5 4 B 3 VSi 2 1 P 0 -6 Si -4 -2 0 2 Energy / eV 4 図 14 VASP (PAW,PBE)で計算した、不純物原子１つを含む Si の電子状態。(左, 中) P ドープシリ コンの CBM 付近の波動関数。(右) Si 欠損とさまざまな不純物を含む Si の全状態密度。フェルミ エネルギーを基準に揃えてある。 もう一度バンド計算に戻りましょう。図 14 は、P ドープ Si の計算をした結果です が、2×2×2 の super cell (ブラベー格子が 8 原子をもちますから、この super cell は 64 原子を含んでいます) をつくり、角にある Si 原子を P 原子に置き換えてあります。構 造緩和させたあとの波動関数について、伝導帯側の最低エネルギーの準位(ECBM)と、 その 0.18 eV 上の波動関数を描いてあります。ECBM では波動関数が P 位置に集まって 84 Si は、BZ の半端な位置にエネルギーバンドの極値をもっていますから、中距離の波動関数の 位相が電子構造に大きな影響を与えていることは想像できます。一方で、酸化物結晶ではエネ ルギーバンドの極値はΓ点と BZ 境界にあるものも多く(BaTiO3 もその例)、この場合は最近接程 度のクラスターでも比較的良い結果が出ることもあります(BZ 境界の波動関数は、隣の単位格 子との反対称軌道であることを思い出しましょう。それよりも遠い波動関数との位相は+1､-1 が交互に続くだけなので、結局、最近接単位格子程度を考えても大差が無いことになります)。 いるものの、それよりも外側にも無視できない密度があります。つまり弱く局在化し た状態があり、これが、P ドープシリコンの 45 meV (文献値) の深さにできるドナー 準位に対応していると考えられます。その上の準位は純粋 Si の伝導帯に入っており、 super cell 全体に波動関数が拡がり、自由電子の状態になっていると考えられます。右 図には全状態密度をプロットしていますが、P ドープ Si のフェルミエネルギーは伝導 帯に近く85、また、伝導帯とつながった準位を作っていることがわかります86。P を取 り除いてダングリングボンドが 4 つある場合の計算(VSi)からは、図 13 と同様、VB と CB のほぼ中央に EF が来ることが確認できます。B の場合にアクセプターになること、 Fe や O 不純物が深い準位を作ることも読み取れます。 (A) (B) (C) (D) 200 4 2 150 0 酸素が抜けた跡に 捕まった電子(色中心) VO DOS Mg2+ O2- -2 100 Relaxed -4 50 Non-Relaxed -6 W L Γ X Z K 0 -6 -4 -2 0 2 Energy / eV 4 図 15 MgO 中の酸素欠損の状態 (VASP, PAW, PBE)。(A) 完全結晶から O2-イオンを除いた (白丸 で囲った部分) 構造。矢印はそれぞれの Mg2+イオンに加わっている力。(B) HOMO の波動関数と その模式図。(C) 緩和構造のバンド構造図。(D) 全状態密度。 もう一つ、典型的な酸化物である MgO の例を見てみましょう。酸化物の中には、 ZnO, SnO2, In2O3 などのように、透明でも電気伝導率が高い「透明導電性酸化物 (Transparent Coducting Oxide: TCO)」があり、平面 TV や太陽電池などの透明電極に使 われています。これらは Si と同様、それぞれ Ga3+, F-, Sn2+などの異価数不純物のドー ピングによってもキャリア濃度を制御できますが、不純物を添加しない純粋材料でも、 酸素欠損などによって電子キャリアが生成します。これは (4) ZnO (結晶) => ZnZn + VO·· + 2e·· (55) ZnO (結晶) => Zni + 2e 87 という反応式で表されます 。この考え方や、電気伝導度の酸素分圧依存性などから、 85 86 87 DFT のバンドギャップ問題のため、この議論には注意が必要です。Si での実験値とのずれは約 0.5eV ありますが、もしも P の準位はそのままで、本当の CBM が 0.5eV 高い位置にあるのであ れば、P の準位は浅いドナーではなく、ミッドギャップ準位と解釈する必要があります。この 問題を解決する方法は、より正確なバンドギャップがでる GW 近似などを使う方法です。しか しまだ一般的ではないため、バンドギャップなどを再現できるように DFT と HF の交換－相関 相互作用項を混ぜ合わせた B3LYP などの混合汎関数を使う方法もあります (バンド計算ソフト ウェ アでは CRYSTAL06 が対応しています)。 64 原子を含む super cell が、孤立した欠陥や不純物の議論するには十分な大きさとは限りませ ん。super cell が十分大きくないために P 軌道がバンドを作ってしまい、伝導帯とドナー準位が つながっているという計算結果が出ている可能性も高く、慎重に検討する必要があります。実 際、この super cell の P 濃度は 1/64 = 1.6% = 8×1020 cm-3 であり、実際の Si でも非常に高ドープ 状態をシミュレートしていることになります。 MA という書き方は Kröger-Vink の表記法と呼ばれ、理想結晶では A イオンが占めていた位置に M イオンが置換しているという意味です。つまり、ZnZn はもとのまま、VO は、酸素イオン位置 が空孔(vacancy)V に置き換わっている、Zni は格子間位置(interstitial, i)を Zn イオンが占めている これら TCO では、酸素欠損や格子間金属イオンの生成によって電子キャリアを発生 すると考えられています。 ところが、MgO などでは、同じように酸素欠損ができても電子キャリアは生成し ません。これは図 15 でわかるように、MgO 中の酸素欠損は、VBM から 3 eV ほど上、 ほぼバンドギャップの中央に電子を 2 つ捕獲した酸素欠損準位をつくり、伝導帯に電 子を残さないためです。このような準位を形成する原因の一つに、構造緩和が挙げら れます。電子が高いエネルギーの準位に存在する場合、構造を緩和させたり欠陥を作 ったりすることでより低いエネルギーに準位ができると、電子がそこへ移ることによ ってエネルギーを得します。構造緩和や欠陥生成には余計なエネルギーが必要ですが、 電子がバンドギャップの上にいれば、価電子帯近くまで準位を下げる余地があります から、エネルギー利得がバンドギャップエネルギー程度まで得られる可能性がありま す。TCO や MgO のバンドギャップ (TCO で 3 eV 以上、MgO で 7.8 eV) は十分大き く、構造緩和や欠陥形成を起こす可能性があります (自己補償効果88)。 図 15 の MgO の場合をみてみましょう。完全結晶から O2-イオンを除くと、もとも と Coulomb 力でひきつけあっていた周囲の M2+イオンが、O2-との引力がなくなった ために外向きの力を受けます。構造緩和計算をすると、この力により、酸素欠損周り の O2- – O2-間距離が 0.4216 nm と、平均値 0.4201 nm よりも大きくなります。HOMO 準位の波動関数(B)から、電子密度がほとんど酸素欠損位置に局在化していることが わかります。バンド構造(C)からは、この酸素欠損バンドに分散がほとんどなく、酸 素欠損は強く局在化していることがわかります89。上記のような自己補償効果が起こ っているかを調べるために、右図には、構造緩和前（完全結晶の構造から O2-イオン を取り除いただけの構造）と緩和後の状態密度を比較していますが、これから、緩和 させていない MgO 中の酸素欠損も、伝導帯に電子は残さず、もともとから深い局在 準位を作ることが確認できます。構造緩和の効果は、酸素欠損位置の空孔体積が増え たため、酸素欠損位置により多くの電子密度が捕獲されるようになっていることで確 認できます。 同様に、(4)の反応式で ZnO 中の酸素欠損はドナーになるように思えますが、これ は、VO や Zni がイオン化する場合には正しい反応式です。実際にこれらがイオン化す るかどうかは材料、ドーピング濃度、雰囲気によって異なるので、別途検証する必要 があります。実際、ZnO 中の酸素欠損を DFT で計算すると、VBM に近い位置に深い 準位を作ってしまい、ドナーにならないという報告があります90。もちろん、計算の 結果をそのまま信じるわけではありませんので、現在も実験結果を踏まえた議論が進 められています。 88 89 90 こと、を意味しています。上付き文字の ·· は、理想結晶の状態よりも+2|e|だけ正電荷が増えて いることを示しています。これは、O2-イオンが抜けたために相対的に正電荷が発生している、 Zn が格子間位置に来てイオン化して Zn2+になっている、と考えるために発生した正電荷です。 両辺での電気中性条件から自由電子 2e-の項が出てきます。もっと普通に、O2-イオンが抜ける ときには中性分子 O2 で抜けるため、電子 2 つを結晶に残していくと考えても構いません。 Self-comensation。ワイドギャップ半導体にドーピングが効きにくい主要因で、ドープされた電 子や正孔が欠陥を作ることで不活性化することです。 (B)のような波動関数の絵では、等高面の値をうまく調整すると、大抵の欠陥や不純物は局在化 しているように描けます。実際に局在化しているかどうかは、バンド分散があるかどうかを含 めて判断する必要があります。また、super cell のサイズを大きくしても結果が変わらないこと も確認することも必要です。 F. Oba et al., J. Appl. Phys. 90 (2001) 824.