応用物理学会結晶工学分科会結晶工学スクールテキスト第 17 版（2024 年）

バンド構造を用いた材料開発（実践編）
東京工業大学 元素戦略 MDX 研究センター

神谷利夫

参 考 ： 関 連 す る python プ ロ グ ラ ム 等 を http://conf.msl.titech.ac.jp/Lecture/jsapcrystal/index.html にて公開しています。これらは教育目的で作成したもので、 間違いやバグ、計算精度の不足、プログラムの最適化などの保証は全くありませんので、自己責任でお使いくだ
さい。また、以下のチュートリアル資料と録画も公開しています。
・実空間像から理解するバンド理論 ・第一原理計算は何の役に立つか」
http://conf.msl.titech.ac.jp/D2MatE/2024Tutorial/tutorial2024-BandTheory.html

1. はじめに前の講義でバンド理論の基礎を学んだわけですが、それでは、この知識を実際の材料に応用するにはどうしたらよいでしょうか。私たちが物質の電子構造を知るためには量子理論に基づいた膨大な計算が必要となり、コンピュータで数値計算することになります。このようなソフトウェアを総称して「量子計算プログラム」と呼ぶことにしますが、現在では第一原理計算が使われる例が多くなっています。最近の第一原理
計算の精度は非常に高く、「材料設計」や「物性予測」にも広く使われるようになってきています。付録 1a, 2a で紹介しているように多くの量子計算プログラムが入手できるようになっていますが、それぞれに特徴があり、対象とする材料と物性によって使い分ける必要があります。本稿では、第一原理密度汎関数理論 (Density Functional
Theory: DFT) に基づくバンド計算プログラムである VASP 1 や WIEN2k 2 を使った結果を中心に、実際の材料をどのように理解するか、あるいは研究する上でこれらをどのように利用していけるかを説明していきます。
なお、ここで掲載するデータや図は、精査して論文に掲載されたものも、本原稿を書くために突貫で計算してみた、ほとんど検証していないものまで混在しています。
また、私の不勉強ゆえに不正確な記述もあります。書いてあることを鵜呑みにするのでなく、あくまでも、量子計算の可能性を知るための参考と考えていただければ幸いです。また、各計算結果の後に (VASP, PAW, PBE) などと書いてありますが、これは用いたプログラムと計算モデルを示しています。この場合は、VASP を用い、PAW(Projector Augmented Wave) 法 (付録 2a) で PBE 汎関数 (汎関数については付録 2b を参照) を用いたという意味です。波動関数などの視覚化には VESTA 3 や XCrySDen 4を使わせていただきました。
本講義では、テキストの本章を主に説明します。付録・脚注には、参考になりそうだけれども少し込み入った話や各論もまとめていますが、最初に本文を読む時には無視して読み進めていただければと思います。
また、本文を通して汎関数としては GGA PBE を主に使っていますが、最近では PBEの改良型の PBEsol や revPBE が使われる例が増えています。それぞれに特徴がありますので、私たち研究室で使っている代表的な汎関数を 12 章にまとめておきます。
量子計算の実際については、文献 5も参考になるかと思います。初学者向けの教科

1

https://www.vasp.at/
http://www.wien2k.at/3
http://jp-minerals.org/vesta/jp/
4
http://www.xcrysden.org/
5
佐々木泰造、末原茂共訳、密度汎関数理論入門 理論とその応用、吉岡書店 (2014).
2

書としては文献6 が挙げられます。文献7 が密度汎関数理論からより進んだ方法までを網羅しています。より広い内容を知りたい場合には 文献8 がまとまっています。
計算モデル・汎関数の解説は 文献9, 10 などに詳しく書かれています。特に文献 9 では、量子力学・密度汎関数法の歴史や各段階における密度汎関数理論の正当性 (運動量汎関数の微分への置き換え、スピン分極混成汎関数理論、混成汎関数) など、密度汎関数理論の基盤となる理論・参考文献も書かれています。文献11では、Quantum
Espresso 12というパッケージについて、実際に計算環境を作ってコマンドラインで計算をしていく手順が説明されています。
2. 量子計算で何ができるか？
第一原理計算の定義は研究者によってばらつきがありますが、「原子配列のみを入力することで精度の高い全エネルギーと信頼できる電子構造が出力される」という点でほぼ一致していると思われます。
本稿では「バンド計算」を扱いますが、これが与えてくれる情報はバンド構造だけではありません。あくまでも、結晶について効率的で精度が高い計算ができるがゆえに「バンド理論」が使われているだけであり、バンド構造は出力される結果の一部に過ぎません。ここで、第一原理バンド計算で得られる情報をまとめておきます。
1) 原子・電子の状態 (電子準位、波動関数、電位分布など) を可視化できる (10 章)
・教科書や講義で学んだ理論を直観的に理解できる
・自分が組み立てた直観的な理論・モデルが正しいかどうかを検証するツールになる
2) 精度の高い全エネルギー計算を行える
・絶対零度の基底状態 13 で安定な構造がわかる (15, 16 章)
最近では、フォノン分散計算を取り入れて自由エネルギーを計算し、有限温度での計算もある程度信頼できるようになっています (付録 2k, 2L)。
・生成エネルギーが計算できる (14 章)
・合成できるはずのない (準安定にもなっていない) 構造をスクリーニング
できる => 実験の構造解析結果が間違っていることを証明できる・標準状態で安定でない構造を安定化させる条件についてヒントをくれる
・合成条件 (化学ポテンシャル) に対して安定な相がわかる (付録 4b)
・与えられたフェルミ準位で安定な欠陥やドーパントの帯電状態がわかる(付録 4a)=>想定した欠陥種の範囲内での平衡フェルミ準位がわかる (付録 4a)

6

田中功、松永克志、大場史康、世古敦人 共著、材料電子論入門 第一原理計算の材料科学への応用、内田老鶴圃 (2017).
7
藤原毅夫著、固体電子構造論 密度汎関数理論から電子相関まで、内田老鶴圃 (2015).
8
R.M. マーチン著、寺倉清之、寺倉郁子、善甫康成訳、「物質の電子状態」、Springer Japan(2010).
9
M. Marsman et al, J. Phys.: Condens. Matter, 20, 064201 (2008).
10
常田貴夫著、密度汎関数法の基礎、講談社 (2012).
11
前園涼、市場友宏共著、動かして理解する第一原理電子状態計算、森北出版 (2020).
12
https://www.quantum-espresso.org/
13
一番全エネルギーが低くなる電子配置の状態

3) 電子構造を定量的に計算できる (バンド構造) (7 章)・直接遷移型半導体か間接遷移型半導体か (7 章)・光学スペクトル (電子の複素誘電関数 = 屈折率、吸収係数) (13 章)
・キャリア輸送特性 (有効質量、フェルミ準位、状態密度) (7, 8 章)
電子伝導度、電子比熱、電子熱伝導度、熱電係数 (別途散乱モデルを考慮する必要がある)・化学結合性状に関する定性的な知見 (9 章)
共有結合性、投影 (部分) 状態密度・フェルミ面の形状 (付録 7b)
・スピン分極した系 (磁性): スピン配置、スピン配列、自発分極など (付録 8)4) 格子系、電子－格子結合系の物性 (付録 2f, 付録 2g、付録 2k、付録 2L)・調和フォノン計算: フォノン分散、フォノンエントロピー
・非調和フォノン計算: フォノン－フォノン散乱時間、格子熱伝導度
・誘電率テンソル (Berry 位相)、有効電荷、圧電テンソル、弾性率テンソル
・赤外分光、Raman スペクトルなど
・電子－フォノン相互作用: 電子－フォノン散乱時間、移動度、電気伝導度など
3. 密度汎関数理論
相対論効果を無視すれば、私たちの世界で電子や原子の運動を支配している厳密な方程式は Schrödinger 方程式であり、イオンの位置を固定して計算する近似「BornOppenheimer の断熱近似」のもとで、N 電子系では以下の式で与えられます。
!
!− ∑# 𝛻# " + 𝑉$%& (𝒓𝟏 … 𝒓𝑵 ) + 𝑉$)$ (𝒓𝟏 … 𝒓𝑵 )+ Ψ(𝒓𝟏 … 𝒓𝑵 ) = 𝐸Ψ(𝒓𝟏 … 𝒓𝑵 )(3-1)"(3-2)式以降で出てくる、変数分離した一電子方程式と区別して、これを「全電子
Schrödinger 方程式」と呼びます。𝑉$%& (𝒓𝟏 … 𝒓𝑵 ) は原子核を含む外場との相互作用、
𝑉$)$ (𝒓𝟏 … 𝒓𝑵 )は電子－電子静電相互作用、E は系の全エネルギー (ハミルトニアンの期待値)、Ψ(𝒓𝟏 … 𝒓𝑵 ) は全電子波動関数です。
この方程式は各電子の 3N 個の座標に関する連立偏微分方程式であり、解くことが困難です。そのため、Pauli の排他律を考慮して変数分離した一電子 Hartree-Fock (HF)
方程式が導出され、分子の計算を主目的とする GAMMES、Gaussian などの非経験分子軌道計算プログラムなどで広く使われています。
!
!− 𝛻* " + 𝑉$%& (𝑟* ) + 𝑉$)$ (𝑟* )+ 𝜑* (𝑟* ) + 𝑣+* 2𝑟* , 𝜑# (𝑟* )4 = 𝜀* 𝜑* (𝑟* )
(3-2)
"

!

𝑣+* 2𝑟* , 𝜑# (𝑟* )4 = − ∑ 6∫ 𝜑 ∗ # (𝑟# ) - 𝜑. (𝑟# )𝑑𝑟# 9 𝜑# (𝑟* )
!"

(3-3)

ここで左辺第 1 項は相互作用のない粒子の運動エネルギー、Vext は外部ポテンシャル(原子核)と電子の相互作用、Ve-e は電子－電子静電相互作用 (Hartree ポテンシャルとも呼ばれる) です。vXl の X は交換相互作用 (eXchange interaction) の意味 (l は l 番目の電子を示す指標) であり、Pauli の排他律 (電子の交換により全波動関数が反対称になる) から導出される相互作用であるため、このように書かれます。HF 近似では他の電子との相互作用は全電子に共通な𝑉$)$ (𝑟* )で近似されており、平均場近似となっています。このことは、「他の電子の状態が変化すれば、注目している電子の状態も変化するはず」という電子相関相互作用が取り込まれていないことを意味しています。
一方、固体を扱う第一原理バンド計算プログラムでは、ほとんどが密度汎関数 14 理
論 (Density Functional Theory: DFT) を用いて計算を行っています。全電子 Schrödinger14

関数 g (x) が関数 f (x) の関数であるとき、g[f (x)] を f(x)の「汎関数」といいます。

方程式が、古典力学に各粒子の位置 xl と運動量 pl の量子力学的交換関係 15 を取り入れて導出されるのに対し、DFT は「あるポテンシャルの下での多電子系の基底状態は電子密度に対して一意的に決まる」という Hohenberg-Kohn の定理 16,17 に量子力学の要請を取り込んだ異なる形式の量子理論です。DFT の第 2 定理には「全エネルギーは正しい電子密度に対して最小の値を持つ汎関数である」という変分原理があり、これにより、電子密度を変数として全エネルギーを最小化するという実際の第一原理計算での計算手順が保証されています。このため、DFT における全エネルギーは𝜌(𝒓)の汎関数として𝐸 [𝜌(𝒓)]で与えられ、𝐸 [𝜌(𝒓)] = 𝑇/ [𝜌(𝒓)] + 𝑉$%& [𝜌(𝒓)] + 𝑉$)$ [𝜌(𝒓)] + 𝐸+0 [𝜌(𝒓)](3-4)とあらわすことができます。𝑇/ [𝜌(𝒓)] は相互作用が無いときの運動エネルギー、
𝑉$%& [𝜌(𝒓)]は外部ポテンシャルと電子の相互作用エネルギー、𝑉$)$ [𝜌(𝒓)]は電子－電子
静電相互作用エネルギーです。交換エネルギーを含むその他の寄与は全て𝐸+0 [𝜌(𝒓)]に取り込まれています。X は交換相互作用、C は電子相関 (electron Correlation) 18 相互作用を意味し、交換相関エネルギーと呼ばれます。
式(3-4)は具体的な汎関数が与えられれば、𝐸 [𝜌(𝒓)]を最小にする𝜌(𝒓)を求めることで解くことができます。しかしながら実際には𝐸 [𝜌(𝒓)]の形がわからないため、相互作用のない電子系で Schrödinger 方程式に一致するよう (独立電子近似) に一電子方程式
化した Kohn-Sham (KS) 方程式 19 が提案されました。KS 方程式は式(3-5)のように、
一電子 HF 方程式(3-2)とよく似た形をしています。
!
KS 方程式: !− 𝛻 " + 𝑉$%& [𝜌(𝒓)] + 𝑉$)$ [𝜌(𝒓)] + 𝑣+0 [𝜌(𝒓)]+ 𝜑* (𝒓) = 𝜀𝜑* (𝒓)(3-5)"

式(3-2)とは vX と vXC が違って見えますが、両者はみかけの計算式は全く同じになります。これはある意味当然の話で、どちらも同じ物理現象を扱っているので、得られる結果は本質的に同じにならないといけません。そのため、物理的な基盤が異なっても正しい方程式は数学的に似ることになります。KS 方程式を導出する際には、一電子
Schrödinger 方程式の解である電子密度 r(r) が、一電子系の DFT 方程式を満たすことを条件とし、r(r) を波動関数のような関数 f(r) で展開していますので、これらの帰結として式(3-2), (3-5) はよく似た形の一電子方程式になっています。
𝑣+0 [𝜌(𝒓)]は𝐸+0 [𝜌(𝒓)]の𝜌(𝒓)に関する変分として導出されます。これらは交換相関汎関数とも呼ばれますが、両者を区別する必要があるため、𝑣+0 [𝜌(𝒓)]は「ポテンシャル
汎関数」、𝐸+0 [𝜌(𝒓)]は「エネルギー汎関数」とも呼ばれます。11 章で出てくる modified
Beck-Johnson 汎関数のように、ポテンシャル汎関数しか与えられていないこともあります。
実際に量子計算を行って結果を解釈する際には、この両者を混同しても問題を起こすことはほとんどありません。しかし、違う量子理論ということを理解しておくに越
15

量子力学的交換関係を要請される物理量は解析力学に基づき決まっています。ラグランジアンが一般化座標𝑞! 、その時間微分𝑞̇ ! と時間 t の関数 𝐿(𝑞! , 𝑞̇ ! , 𝑡) として与えられる系では、𝑞!
に共役な一般化運動量は𝑝! = 𝜕𝐿/𝜕𝑞̇ ! で与えられ、𝑞! 𝑝! − 𝑝! 𝑞! = 𝑖ℏ が要請されます。
16 P. Hohenberg and W. Kohn: Phys. Rev. 136 (1964) B864.
17
Hohenberg-Kohn の定理では、基底状態が縮退していないこと、電子密度がポテンシャルで決められる関数に限定されるなどの芸減がありましたが、Levy, Lieb により拡張され、縮退している基底状態でも正しいことが証明されています。
18
本来の意味は「ある電子の状態がその他の電子の状態に与える影響」です。HF 近似は「他の電子の状態」を平均化しているため、一般に、HF 近似で取り入れられていない効果が電子相関と呼ばれます19 W. Kohn and L.J. Sham: Phys. Rev. 140 (1965) A1133.

したことはありません。一番重要な違いは、Schrödinger 方程式が各電子の座標 rl = (xl,yl, zl) (l = 1, 2, ･･･, N) を変数としているのに対し、DFT では電子密度 r(r) も波動関数に相当する f(r) も、空間座標 r = (x, y, z) だけを変数としています。そのため、N
電子系の Schrödinger 方程式が 3N 個 の変数の偏微分方程式であるのに対し、DFT では 3 個の変数しかありません。これが、DFT が多電子系に適した理論である理由の一つです。
また、前述したように、HF 方程式には電子相関が取り込まれていません。一電子
Schrödinger 方程式からスタートして電子相関を取り入れるためには配置間相互作用(Configuration Interaction: CI) や摂動法のような、HF 近似を越えた複雑なモデルを取り入れる必要があります。それに対して DFT の場合は、形式上、これらの効果もすべて vXC 項に含まれます。そのため、vXC 項の近似を上げることにより、理論や計算プログラムの枠組みを大きく変更せずに、多電子系のより高度な相互作用を含んだ計算が可能になります。
4. 第一原理計算で得られる「電子準位」は何を表しているか
HF 近似と DFT には他にも、似ているけれども違うことがいくつかあります。文献
11 の表現を借りると、「HF 方程式の𝜑* (𝒓* )は多電子 Schrödinger 方程式を一電子方程式化する際に全波動関数を変数分離するために導入した各電子の一電子波動関数」ですが、「KS 方程式の𝜑* (𝒓)は、一電子関数の自乗和が厳密な電子密度を表現するために導入された一電子関数」であり、本質的に異なります。ただし、HF 近似を「相関項ゼロの密度汎関数法ととらえれば、その極限で両者は接続」します (以上の「」中の表現は文献 11 の原文から改変しています)。

300

a-IGZO
Zn 3d

DOS (a.u.)400

状態密度 (PBE96)
Ga 3d
In 4d

500

O 2s

XPS

c-IGZO
200
100

c-IGZO(no In/Ga d)

ZnO

0-25 -20 -15 -10 -5 0

Binding energy / eV

5

Binding energy / eV

図 4-1 エネルギー準位に関する実験 (左図: X 線光電子分光、XPS) と DFT (PBE、全状態密度、DOS、と d 軌道の投影状態密度、PDOS) の比較。酸化物半導体である結晶
ZnO、結晶 InGaZnO4、アモルファス InGaZnO4 を例にした。エネルギー原点は XPS ではフェルミ準位 (EF)、DOS では価電子帯上端＋3.0 eV (ほぼ EF) にとって比較している。ZnO 3d 準位は実測で EF から -11 eV、DFT で -8 eV と、DFT の計算結果が小さくなっていることがわかる。c-IGZO (no In/Ga d)では、In 4d と Ga 3d 軌道を擬ポテンシ
ャルにとりこんでいるため、計算結果にはこれらの軌道が現れていない。

HF 近似で得られるエネルギー固有値 e i はイオン化ポテンシャル e i = E N - E N -1
(EN は電子が N 個あるときの全エネルギー) や電子親和力 e i = E N +1 - E N に対応します 20 が、DFT では e i = ¶E / ¶ni (E は全エネルギー、ni は i 番目の軌道を占有している電子数で、非整数として扱う) と、化学ポテンシャルのような量に対応します (10章、図 10-2 参照) 21。そのため、DFT で得られるエネルギー固有値は光電子分光法などで測定される束縛エネルギー (各軌道のイオン化ポテンシャル) よりも小さくなります 22 (図 4-1)。また、HF 近似では電子準位がイオン化ポテンシャルに対応すると述べましたが、実際の光電子分光法などで測定されるイオン化ポテンシャルは、イオン化による電子相関・構造緩和のため、HF 近似で得られるイオン化ポテンシャルよりもかなり浅くなります。
5. 全エネルギーから何が計算できるか

第一原理計算変分原理: 全エネルギーの期待値 <ψ|H|ψ> を最小化する固有状態をみつける

エネルギー準位
k依存性

電気双極子遷移確率<ψfin|x|ψini>

状態密度

バンド構造

全エネルギー

波動関数

E依存性
第1B.Z.で総和

光学スペクトル

Berry位相

|ψ|2

電気分極

電子分布

cicj<φｊ|H|φi>

有効質量

化学結合の知見
結合当たりエネルギー (iCOHP)

自由キャリア密度
イオン化不純物
状態密度有効質量

磁気物性、半導体物性: 移動度、キャリア濃度、電子伝導度、Hall係数、熱電係数第一原理計算全エネルギー U
自由エネルギーH = U – PV F = U + TS G = U – PV + TS

波動関数
Berry位相電気分極

フォノン分散微小変位モデルから微分(finite displacement method)

誘電率
圧電定数

安定構造

歪モデルから
全エネルギー微分

弾性定数

欠陥モデルとの全エネルギーの差

反応式中の物質の全エネルギーの差

反応エネルギー

欠陥生成
エネルギー

異なる帯電状態の欠陥モデルの比較
電荷遷移準位(ドナー／アクセプター準位)

図 5-1

20

第一原理計算と、各種物性計算の関係

Koopmans の定理: T. Koopmans, Physica 1 (1933) 104. 電子配置の変化によって軌道関数が変わらないことを前提としている。
21
Janak の定理: J.F. Janak, Phys. Rev. B 18 (1978) 7165. 関連したものに"Slater の遷移状態法"がある。ある軌道 i の電子数 ni が 1 から 0 に変化する際のイオン化ポテンシャルは、DFT では ni =0.5 のエネルギー固有値で近似できる。図 10-2 参照。
22
例えば W.R.L. Lambrecht et al., Phys. Rev. B 50 (1994) 14155.

2 章では、第一原理計算の重要な点は「精度の高い全エネルギーが出力される」と書きました。私たちが興味を持つ物性は全て、熱力学的（統計力学的）自由エネルギ
ーから導出されます。例えば定温・定圧条件では Gibbs エネルギー G = E + PV – TS,
定温・定積条件では Helmholtz エネルギー F = E – TS がわかれば、その微分等から、その温度・圧力における物性が計算できます。しかしながら、一般的に用いられる「断熱近似の第一原理計算」では、イオンの運動（フォノン）の効果を計算できませんので、多くの場合は DFT の計算で得られた全エネルギー E (熱力学の内部エネルギー)を使って、絶対零度における物性の計算を行います。最近ではフォノン計算も一般的に行われるようになっており、G を計算し、有限温度での安定構造・相転移などの計算もされるようになってきました。
ということで、「原理的には」第一原理計算により、2 章でまとめたような、私たちが興味を持つ多くの物性を計算できることになります。実際の計算では、DFT により第一原理計算を行った後、ポストプロセッサーと呼ばれる別のプログラムにより物性計算をしたり、結晶構造を変えて第一原理計算を行い、その結果を総合してフォノンなどの物性を計算します。つまり、第一原理計算を利用するということは、DFT 以外にも多くのプログラムを組み合わせることになります。図 5-1 に、第一原理計算と物性計算の関係図をまとめます。
6. バンド構造の読み方: フェルミエネルギー (EF)、バンドギャップ (Eg)
最初に、バンド構造図の読み方をまとめておきます。図 6-1(A) に WIEN2k で計算した Si のバンド構造図を示しますので、比べながら読んでみてください。

図 6-1 Si のバンド構造 (WIEN2k, PBE)。(A) DFT で計算したバンド構造。(B) 非占
有状態を平行移動させてバンドギャップ値を実測値にあわせたバンド構造。(C)第一 BZ と対称性の高い k 点の位置と記号。
1) 横軸は Bloch の波数ベクトル k を表しています。ただし k を数値で表記してもわ
かりにくいので、対称性が高い k 点にはそれぞれ記号が定義されています。記号の定義は脚注23の結晶構造データベースなどで調べられますが、知らなくても、重要な情報 (バンドギャップ、間接遷移型かどうかなど) はわかります。
23

例えば、bilbao crystallographic server (http://www.cryst.ehu.es/) を見てください。この Web ペー
ジから ”Space-group symmetry” => ”KVEC”を選び、空間群の番号 (たとえば Si の空間群はFd 3 m で、空間群番号は 227 です) を入力すると、図 5-1(C)のような第一 BZ の絵と、対応する点の記号、座標が表示されます。

1𝒌2) 加速定理により、Bloch の波数ベクトル k は𝑭 = ℏ 1& という古典力学に対応する関係を満たし24、k を電子が進む方向 (正しくは運動量の進む方向) とみなせます。
そのため、k はほとんどの場合、運動量のようにふるまいますが、脚注25のような例外もあります。
3) バンド構造は記号の書いてある k 点において縦線で区切られており、この縦線は、逆空間の原点 (G点) か、第一ブリルアンゾーン (Brillouin Zone: BZ) の境界のいずれかです。
4) 縦軸は電子のエネルギーを表し、下へ行くほど電子のエネルギーが深く、安定になります。
5) エネルギー固有値は k の連続関数 Ei(k) になり、これを「エネルギーバンド」と呼びます。そのため、バンド構造図は複数の複雑な曲線 E(k)で描かれています。
6) 波数ベクトル k をもつ電子は、飛び飛びのエネルギー Ei(k) しかとることができません 26。これは原子の中の電子と同じです。ただ、原子の中の電子は、たとえば 1s 軌道と 2s 軌道の中間のエネルギーを持つことはできませんが、結晶では k に
よってエネルギー固有値 E(k) が連続的に変わるため、「中間のエネルギー」も取ることができるようになります。
7) 電子のスピン分極を考慮しない場合は、一本のエネルギーバンド 27 には、単位格子あたり up スピンと down スピンの 2 つの電子が入ることができます。スピン分
極したバンドの場合は、これらのバンドを別々に考慮するため、それぞれ 1 つず
つが入ることができます 28。
24

上村洸、中尾憲司 共著、電子物性論、培風館 (1995); キッテル、固体物理学入門 第 8 版、丸善 (H17).
25
例えば BZ 境界では Bragg 条件 k’ = k + Ghkl によって電子が散乱されるため、散乱前後で k は保存されません。つまり、Bloch の波数 k は運動量保存則を満たさないので、厳密な意味での運動量ではなく、結晶運動量と呼ばれることもあります。このような運動量を保存しない散乱
は、フォノン散乱の場合はウムクラップ散乱と呼ばれています（運動量が保存される散乱は正常散乱と呼ばれる）。
26「結晶の固有状態としては」そのエネルギーをとりえない という意味。電子は結晶中でも原子による散乱を受け、進行波成分は減衰し、残りが散乱波となります。ところが、E(k) のエ
ネルギーをもつ場合に限り、3 次元周期構造による干渉が散乱波を打ち消し、進行波だけが残ります。そのおかげで、全く散乱を受けずに進行することができます (透過率が 100%の状態)。これが、固有エネルギーと固有状態、「エネルギーバンド」になります (薄膜による光の干渉と同じです。光は必ず表面で反射されますが、干渉条件を満たす場合は反射率が０にな
り、透過率が 100%になります)。(「チュートリアル： 実空間像から理解するバンド理論」(http://conf.msl.titech.ac.jp/D2MatE/2024Tutorial/tutorial2024-BandTheory.html) 参照)
固有エネルギー以外のエネルギーでも波動関数は存在しますが、波数ベクトルが虚数になります。そのような電子の進行波成分は、原子による散乱によって波数の虚部に反比例した減衰長で減衰します。厚さが減衰長程度の薄膜の場合にはこのような減衰があっても、対向電極にも電子が透過する確率が残るので、いわゆる「トンネル現象」がおこります。トンネル現象というのは、固有状態でない E(k) を持つことで起こる現象です。あるいは、「固有状態でのトンネル確率は 100%である」といっても構いません。
27
正確には、一つのバンドに「第一 BZ 全体」で入る電子の数が、単位格子当たりで電子 2 個になります。バンド構造図の一区画のバンド「一本」に入る数ではありません。
28
強磁性体の場合は、すべての磁性原子のスピンの空間的な向きが同じなので、up スピン、
down スピンと表現するのは自然に感じます。一方、反強磁性の場合は、同じバンドの電子であっても、隣同士の磁性原子の電子スピンの向きは逆に向くため、up/down という表現が適切でなくなります。このような場合、aスピン, bスピンなどと表記することもあります。

8) 特に明記されていない場合、縦軸のエネルギーの原点はフェルミエネルギー(Fermi energy: EF) とされています。つまり、0 K において、0 eV より下の軌道は電子が詰まっている被占有準位で、それより上の軌道は電子がいない非占有準位
です 29。
9) 金属であるか、絶縁体（半導体）であるかは、EF を交差するエネルギーバンド Ei(k)があるかどうかでわかります。ただし、バンド構造は特定の k 点だけしか表示しないため、BZ の他の k 点で EF と交差するバンドが無いかどうかを確認できません。そのため、8 章で説明する状態密度 D(E) が EF でゼロであるかないかで判断する方が信頼できます。
上でまとめたように、図 6-1(A) のようなバンド構造図の横軸は、第一 BZ 内の Blochの波数ベクトル k を表しています。ただし、3 次元に拡がる波数ベクトルのすべてを表示するのは無理がありますから、結晶格子 (正確には BZ) の中でも対称性の高い方向についてだけ抜き出して描きます。横軸の記号、W, L, G, X, K は逆空間における対称性の高い点、L, D、は対称性の高い逆空間の方位を表している (図 6-1(C)) のですが、脚注 23 のように、データベースなどで調べられますから、この記号を覚える必要はありません。ただ、G点は逆空間の原点 k = (0, 0, 0)、X 点は x 軸方向の BZ 境界(単純格子の場合は k = (1/2, 0, 0) 30) を表しているということくらいは覚えておきま
しょう (脚注31 を読むとその重要性がわかります)。
図 6-1(A)を縦軸のエネルギーの視点から説明していきましょう。}(a) と }(c) で示した領域はどのエネルギーでも Ei(k) 曲線に被りますが、}(b) の領域では、その範囲のエネルギーにかぶる Ei(k) 曲線がありません。つまり、どのような k をもつ電子でも、}(b) の範囲のエネルギーをとることはできないことを意味しています。このようなエネルギー範囲を禁制帯 (forbidden gap) やバンドギャップ (bandgap) とよびます。
7. バンド構造の読み方: 直接遷移、間接遷移
図 6-1 から、シリコンは G点 k = (0,0,0) に価電子帯上端 (Valence Band Maximum:
VBM) をもち、G–X 軸の途中に伝導帯下端 (Conduction Band Minimum: CBM) を持つことがわかります。バンド構造図が便利なのは、光学遷移における k 選択則 (13 章)
29

金属の場合は EF がバンドを横切るので EF は一意的に決まります。絶縁体や真性半導体ではバ
ンドギャップ中のどのエネルギーもこの条件を満たすので一意的に決まりません。半導体物理で使われる EF は電子の化学ポテンシャルの意味であり、Fermi-Dirac 分布関数、温度と電荷中性条件から決まり、真性半導体ではほぼバンドギャップの中央にあります。固体物理やバンド計算では、単純に「絶対零度での電子配置における最高の電子のエネルギー」をフェルミエネ
ルギーと呼び、価電子帯上端をエネルギー原点にとり、それを EF としたバンド構造図が良く使われます。
30
一般に k 点の座標は、逆格子ベクトル a*, b*, c*を単位とした内部座標で表されます。結晶学や回折学では逆格子ベクトルは a* = b ´ c / V のように定義されます (V は単位格子体積) が、物性理論やバンド理論では a* = 2p b ´ c / V と、係数が 2p だけ違うことに注意してください。この余計な係数のおかげで BZ 境界の k を 1/2 のような単純な分数で表現できるようになります。
31
LCAO での結晶軌道は Bloch の定理を使って
e ik×ri f (r - ri ) (f (r - ri)は原子の波動関数、ri はi 番目の単位格子の座標) で表されますが、G点では + f r - a + f r + f r + a + と簡単になります (a は格子ベクトル)。つまり、単位格子内の原子の波動関数を単純に足したもので、fが s 軌道であれば結合軌道になります。X 点は (1/2, 0, 0) ですから、a 軸方向の波動関数だけを抜き出すと、  - f r - a + f r - f r + a + になり、係数が +1, -1 と交互に変わることが判ります。これは分子軌道でいうところの反結合軌道です。つまり、G点から BZ 境界に kが変わるにつれ、波動関数は結合性から反結合性へと変わります (これは s 軌道の場合です。
s結合の p 軌道では反結合性から結合性へと変わります)。

å() () ()() () ()

を読み取ることができることです。光学遷移は一次近似では光の電場が誘起する電気双極子遷移で起こり、同じ k 点以外間の遷移確率はゼロになります。VBM と CBM が同じ k 点にある場合は、VBM => CBM の電子の遷移確率は高く、光学吸収係数も大きくなり得ますが、このような半導体を「直接遷移型半導体」と呼びます。Si の場合
は VBM を形成している電子準位の k 点と CBM を形成している k 点が異なりますので、一次近似では VBM => CBM 間の遷移では光学吸収係数がゼロとなる (禁制遷移)、「間接遷移型半導体」であることがわかります 32。間接遷移が起こるためには、運動量保存則を満たすために 33 波数ベクトルの差に等しいフォノンなどを吸収あるいは放出する必要があり、光学吸収係数が非常に小さくなります 34。また、発光ダイオー
ド (Light-emitting diode: LED) では、伝導帯に電子を、価電子帯に正孔を注入し、これらが再結合する際に出てくるエネルギーを光として放出させます。間接遷移型の場合
はバンド端間の電子の k が異なるために遷移確率が低く、他のフォノンや欠陥などを介した非発光遷移過程が先に起こってしまってほとんど発光せず 35、LED の発光層としては使えません。バンド構造図から、このようなことまで読み取れます。
図 7-1 に典型的なバンド構造をまとめますので、それぞれ、金属であるか半導体であるか、直接遷移型であるか間接遷移型であるか、考えてみてください。
6

(B)

(A)

(C)

(D)

エネルギー / eV

4
2
0
-2
-4
-6
G

M K

G AW

L

G

X WK WLG

X WK WLG

X WK

図 7-1 (A) Mg, (B) Si, (C) GaAs, (D) TiN のバンド構造 (WIEN2k, PBE)。
8. 有効質量と電子物性前の講義で「ほぼ自由な電子 (NFE) 近似」を学びました。NFE 近似はバンドギャップ以外にも「有効質量近似」というアイデアを提供してくれます。図 7-1(B,C) のそ
32

実験的に区別する場合には、吸収スペクトルa(E)が a (E ) = (E ± w ph - E g )2 / w 2 (間接遷移)であるか a (E ) µ (E - E g )1/ 2 / w 2(直接遷移) であるかで判断できることがあります。遷移確率のモデル化の違いにより、分母に𝜔" が入らない式を使う場合もあります。
脚注 26 で「k は運動量保存則を満たさないので厳密な意味の運動量ではない」と書きました
が、Bragg 条件を満たす場合などを除くと、通常の運動量と同じ保存則が出てきます。
34
一次摂動 (電気双極子遷移) での遷移確率はゼロですが、フォノンとの相互作用だけでなく高次の摂動によっても遷移が起こります (多光子吸収・放出過程)。ただし遷移確率は非常に小さくなります。
35
注入された電子や正孔が両方とも高いエネルギーをもって同じ k 点の状態にあれば、k 選択則を満たした許容遷移で発光再結合しそうです。しかし、高エネルギー電子 (ホットエレクトロン) は fs オーダーの短い時間でフォノンなどとの非弾性散乱によって一番エネルギーの低い状態に緩和します。一方、許容遷移でも再結合時間は ps 程度になりますから、ほとんどの場
合、ホットエレクトロンが CBM と VBM、あるいは局所的なエネルギー極小状態まで緩和した後で再結合が起こります。
33

れぞれのバンド構造について、CBM と VBM 付近の Ei(k) が二次関数的になっていることに気づいたでしょうか？ 36 NFE 近似ではなく LCAO (原子基底の一次結合,
Linear Combination of Atomic Orbitals。バンド計算の場合は Tight binding (TB) 法 37 が該当する) 理論に基づいてバンド理論を展開することも可能ですが 38、その場合、一次元結晶でのエネルギー分散 E(k)は(8-1)
E(k ) = e1 - 2 h12 cos(ka ) ~ e1 - 2 h12 + h12 a 2 k 2で与えられます (“TB バンド”と呼びましょう)。これを自由電子の(8-2)
E (k ) = E0 +  2 / 2me k 2と比較すると、あたかも電子が質量𝑚∗ = ℏ" (𝜕 " 𝐸/𝜕𝑘 " ))! = ℏ" /2|ℎ!" |𝑎"(8-3)
39を持っているとみなすことができます(図 8-1(A))。事実、半導体では、結晶中の電子を「誘電率 er の真空中を運動する質量 m*、電荷 –e の電子」、正孔を「電荷 +eの電子」として扱っても良い近似となることが知られています。そのため、m*を有効質量と呼び、このような近似を「有効質量近似」あるいは「有効媒質近似」と呼びま
す。

(

)

EnergyEnergy/ 2|h/12|h12|
|

2
1.5

自由電子
TBバンド

1
0.5

(A)
0
-0.5

-0.3

-0.1

2π/a

0.1

0.3

0.5

図 8-1 (A) TB 理論により求めたエネルギーバンドと自由電子近似したバンド。h12 は共鳴積分 (移送積分) で、バンド幅は 4|h12|になります。(B) SnO2 のバンド構造(VASP, PAW, PBE)。(C) SnO2 のG-Z の CBM エネルギーバンド (●)、自由電子近似バンド (実線)、電子の有効質量 (◇)。
図 8-1(B,C) には、SnO2 を例に有効質量の解析を示しています。SnO2 は TiO2 と同じルチル型の結晶構造を持ち、正方晶に属し、G–X とG–Z 方向で異なる有効質量 me^c,
me//c をもちます。図 8-1(C) はG–Z 方向の CBM バンドについて、E(k) の二階微分係数から計算した m*(E) を重ねてプロットしています。これから、Z 方向、つまり c 軸方向の有効質量 me//c は G点で 0.182me と得られます。文献値は 0.234me なので少し過36E(k)が二次関数的に振る舞わない材料も知られています。Graphene などでは直線的な正負の傾きを持つ 2 本のバンド E(k) の交差点が EF となり、Dirac 電子系などと呼ばれています。六回
対称性のために BZ 境界でバンドギャップがゼロになる k 点があり、そこでは直線バンドの交差点で有効質量がゼロとなります。バンドが直線のために電子の運動量が保存される、などの特徴があります。
37
W.A. ハリソン、固体の電子構造と物性、現代工学社 (1984).
38
化学をベースにしている人や化合物半導体を扱っている人の場合は、自由電子理論よりもむし
ろわかりやすいのではないでしょうか39
単純な場合は, a は該当する方向の格子定数になります。一般の場合には、kBD·a = pとなる長さ (kBD は BZ 境界の波数) になります。

小評価しています。G–X 方向、つまり a 軸方向では有効質量 me^c は 0.197me が得られ、やはり文献値 0.299me より過小評価しています。この点については後で説明します。また、正確な m*を求めるには、精度の高い微分公式を使うこと40、k 点のメッシ
ュを小さく変えながら収束した値を使う必要があることにも注意が必要です。
また、k がG点から Z 点方向へ変化するにつれ、有効質量が増加して発散し、その後で負の値になっています。これは E(k) の変曲点で曲率がゼロになるためで、これを境に m* が正の側を「電子バンド」、負の側を「正孔バンド」と呼ぶことがあります(負の有効質量の解釈については付録 7e を参照)。
有効質量を大雑把に見積もる方法もあります。一つのエネルギーバンドの最低エネ
ルギーと最大エネルギーの差を「バンド幅 41 (W)」と呼びますが、式(8-1) からは W= 4|h12|になることがわかります。単純なエネルギーバンドであればこの値はバンド構
造図から読み取れますので、 m* = 2 2 / Wa 2 から有効質量を計算できます。SnO2 の例
では W = 5 eV、格子定数 c = 0.3227 nm ですから、 m//* c = 2 2 / Wc 2 = 0.292me が得られます。E(k) から計算した値より結構過大評価しますので、大事なデータの場合はちゃんと E(k) 曲線から求めましょう。
7
PBE96

6
me* -1

5

PBE0

4

HSE
HSE*

3
2
1
0
0

0.5

1
Eg-1 (eV-1)

1.5

図 8-2 異なる汎関数で VASP により計算した ZnO のバンドギャップの逆数と電子の有効質量の逆数の関係。PBE0, HSE は、標準の混合パラメータ 1/4 と遮蔽係数を用いた結果。HSE*は標準の遮蔽係数の値を半分に減らした結果。実測値はキャリア濃度によって、0.25~0.3 me* 程度の範囲で変わります 42。
有効質量の過小評価の原因の一つは、イオン性結晶では、電子の電荷と光学フォノンが相互作用する (電子－格子相互作用) ことで実測される有効質量が重くなることが挙げられます。これを「ポーラロン」と呼びます。ポーラロンの有効質量は、Frölich
ポーラロンモデルで結合定数aを用いて、m* ~ m0*(1 + a/6) で表されます43。aは GaAs
では 0.068 と小さいですが、イオン性の強い SrTiO3 では 3.8 で、m*は m0*の 1.6 倍ほど
40#$

𝑘! が等間隔∆𝑘で並んでいるデータ (𝑘! , 𝐸! ) について 𝑘! における微分は、#"$(𝐸!'( − 𝐸!)( )/(2∆𝑘)、|= (𝐸!'( − 2𝐸! + 𝐸!)( )/∆𝑘#% " %&%!

"|=#% %&%!

を使うと、3 点だけのデータ点から短い計算時間で高い精度の微分の計算ができます。これより高次の数値微分公式もあります。
41
この言葉も、分野によって違う言葉で使われます。一般的な光電子分光測定などでは「エネルギーバンド」を一本ずつ区別することができないため、光電子シグナルのピーク全体の幅を「バンド幅」ととることがあり、一般に複数のバンドの全幅になります。
42
H. Fujiwara and M. Kondo, Phys. Rev. B 71 (2005) 075109.
43
H. Fröhlich, Adv. Phys. 3 (1954) 325.

になります。逆に言えば、イオン性の強い結晶で、バンド理論で計算した有効質量が実測値に一致したら、むしろどこかに大きい誤差が入っていることになります。
また、11 章で説明するように、DFT ではバンドギャップを過小評価し、その結果、伝導帯と価電子帯の波動関数の相互作用も強くなってバンド分散を過大評価し、有効質量を過小評価することがあります。例えば k･p 近似からは、有効質量 m*とバンドギャップ Eg の関係として
m / m* ~ 1+(2P2/m)Eg-1
(8-4)
44
-1
が得られます (P は運動量) 。図 7-2 は異なる汎関数で計算した ZnO の Eg と me* をプロットしています。PBE0 の計算結果に精度の問題が残っていますが、ほぼ 式(8-4)に従っていることが確認できます。
上記のように、有効質量ではいろいろな値が使われます。半導体で重要なのは、自由電子吸収解析などから得られた実験値で、単に「（キャリア）有効質量」と呼ばれます。これに対して、バンド構造 45 から求めた有効質量を「バンド有効質量」と呼んで区別することがあります。半導体・金属の熱電係数 (Seebeck 係数) は状態密度が大きな寄与を持ちますが、等方的な非縮重バンドを仮定して状態密度関数 D(E) から計算される有効質量は、バンドが縮重していたり異方性が大きい場合にはキャリア有効
質量とは大きく異なります。これを「状態密度有効質量」と呼ぶこともあります。第一原理計算で得られた状態密度関数から状態密度有効質量を求める方法については
付録 7d をご覧ください。
(B)

(C)

(A)

図 8-3 VASP で得られた𝐷(𝐸)から計算した半導体統計物理量。(A) EF = 0.4 eV, T =3000 K における状態密度関数 (黒線)、Fermi-Dirac 分布関数 (fe, fh)、DC(E)fe(E)とDV(E)fh(E)。(B) T = 300 K におけるキャリア密度の EF 依存性。バンドギャップ内の直線領域の傾きから有効状態密度が NC = 6.26´1017 cm-3, NV = 4.74´1018 cm-3 と得られます。(C) EF = 0.4 eV におけるキャリア密度、イオン化ドナー・アクセプター密度の温度依存性。ドナー／アクセプター状態は EC – ED = 0.28 eV, ND = 0.5´1017 cm-3, EA – EV
= 0.32 eV, NA = 0.8´1017 cm-3 としています。
キャリア密度、イオン化ドナー・アクセプター密度、有効状態密度などは、状態密度𝐷(𝐸)だけから計算できます。まず、状態密度を伝導帯状態密度𝐷0 (𝐸)と価電子帯状
態密度𝐷3 (𝐸)に分解します (図 8-3(A)で、縦鎖線より低エネルギー側が𝐷3 (𝐸)、高エネ
ルギー側が𝐷0 (𝐸))。
𝐷(𝐸) = 𝐷3 (𝐸) + 𝐷0 (𝐸)(8-5)
これらから電子密度 Ne、正孔密度 Nh は

44
45

P.Y. ユー、M. カルドナ、半導体の基礎、Springer (1999).
結晶格子が格子振動などで変形することを考慮しないので、「rigid band model」と言います。

5

𝑁$ = ∫6 𝐷0 (𝐸)𝑓$ (𝐸, 𝐸4 , 𝑇)𝑑𝐸

(8-6）

#

6$
𝑁7 = ∫)5
𝐷3 (𝐸)[1 − 𝑓$ (𝐸, 𝐸4 , 𝑇)]𝑑𝐸(8-7)
で計算できます(図 8-3(A)で、小さい青色ピーク領域と赤色ピーク領域の面積が𝑁$ 、
𝑁7 になる)。ドナー準位 (密度𝑁8 、エネルギー準位𝐸8 )、アクセプター準位 (𝑁9 、𝐸9 ) をパラメータとして与えれば、
𝑁8 : = 𝑁8 [1 − 𝑓$ (𝐸8 , 𝐸4 , 𝑇)](8-8）)
𝑁9 = 𝑁9 𝑓$ (𝐸9 , 𝐸4 , 𝑇)
(8-9）
でイオン化ドナー／アクセプター密度が得られます。図 8-3(B)は、𝐸4 を変えながら(85)~(8-9)を計算した結果です。ここで、EF がバンドギャップ中央付近では、非縮退領
域の式
𝑁$ ~𝑁0 exp[−(𝐸4 − 𝐸0 )/(𝑘; 𝑇)](8-10)[])/(𝑘
𝑁7 ~𝑁3 exp −(𝐸3 − 𝐸4
(8-11); 𝑇)に従いますので、この領域の傾きから伝導帯有効状態密度𝑁0 や価電子帯有効状態密度
𝑁3 を計算し、デバイスシミュレータに入力することもできます (ただし、計算精度を確保するためにはバンド端近傍の状態密度を非常に精度良く計算する必要があります。脚注46をご覧ください)。
実際の材料では、電荷中性条件
𝑁8 : + 𝑁7 = 𝑁9 ) + 𝑁$
(8-12)
を満たすように𝐸4 を決め、式(8-6)~(8-9)により𝑁$ 、𝑁7 、𝑁8 : 、𝑁9 ) を計算します。温度を変えながら計算した結果が図 8-3(C)になります。
金属の EF 近傍 (例: 図 7-3(A,D))や、半導体の価電子帯端近傍 (例: 図 8-3(B,C)) では、複数のバンドが交錯しています（マルチバンド）。この場合にはそれぞれのバンドが異なる有効質量を持ちますが、それを実験から得られる値と直接対応させることは困難です。このような場合は、バンド構造・状態密度と Boltzmann の輸送理論を組み合わせて、以下のように Hall 電圧や熱起電力などの計算が可能です。
𝜏(𝐸, 𝑇) = 𝜏/ 𝑇 < (𝐸 − 𝐸0 )-)!/"
r: 散乱因子(8-13)" 5?@% (6).
.
⟨𝜏 ⟩ = − ∫6 (𝐸 − 𝐸0 )𝜏(𝐸) 𝐷0 (𝐸 − 𝐸0 )𝑑𝐸T𝑛$
(8-14)
>

伝導度

#

$

𝜎% = 𝑒𝑛$ C∗ ⟨𝜏 ! ⟩
6

%

4

?6

!

Hall 係数: 𝑅D = E ;' = FG* = FG
( )
〈I+ 〉

*

(8-15)(q = –e for electron, +e for hole)(8-16)Hall 因子: 𝐹D = 〈I〉+%(8-17)!

Hall 移動度: 𝜇D = 𝐹D 𝜇1-#@&,-

Hall キャリア濃度: 𝑛D = 4 𝑛$

+

./.0

. ∫L),.M8(6)N IO !1 P16

熱電係数: 𝑆 = − $$𝜇1-#@& = C∗ ⟨𝜏 ! ⟩ドリフト移動度,∫L),.M8(6)N + I16

*

! ?6

+ $ ?Q0(8-18)(8-19)

自由電子で近似できる良い半導体のエネルギーは𝐸(𝒌) = (ℏ" /2𝑚∗ )𝑘 " のように k の平方に比例し、状態密度は𝐷+ (𝐸) = 𝐷+, 9𝐸 − 𝐸+ のように E の平方根に比例します。一方、状態密度の積分に
Tetrahedron 法を使うと、E は k に比例し、𝐷+ (𝐸)は𝐸 " に比例します。非縮退半導体のキャリア密度はバンド端の数 kBT 程度の状態密度によって決まるため、この領域で𝐷+ (𝐸) = 𝐷+, 9𝐸 − 𝐸+ の状態密度関数が再現できるように十分に細かい k 点メッシュでエネルギー固有値を計算し、積分を取る必要があります。これらの全てで DFT 計算をして固有値を求めると非常に時間がかかりますので、まず DFT で比較的荒い k メッシュで E(k)を計算し、BoltzTraP2 などの補間機能を使って計算する方が良いと思います。
46

これらのテンソル表記については、BoltzTraP の元論文47をご覧ください。図 8-4 に、
LaCuOSe という層状構造をもった半導体について、WIEN2k と BoltzTraP 48を使って計算した例を示します 49。

nHall,xx

1
-1

nDOS
-0.5

0

0.5

EF (eV)

1

1.5

2

30

5

(B)

0

µzz

Effective mass (m e)/Vs)(cm22/VS)
Mobility (cmMobility

21

100
10

10

(A)

nHall,zz-3cm-3))
Carrier density (x10
Carrier(1021 cm

1000

-5
-10
-15
-20
-1

µxx
0

1

21

2

-3 -3
Carrier))
Carrierdensity
density(10(1021cm
cm

3

(C)

20

mzz

10

mxx

0

-10
-1

0

1

2

3

21 cm-3
21-3)
Carrier
Carrierdensity
density (10(10 cm)

図 8-4 LaCuOSe について計算した (A) Hall 係数から計算したキャリア密度、(B)
動度、(C) 有効質量 (WIEN2k+BoltzTraP)。

移アモルファス物質のように周期性がない場合は、k が「良い量子数」とならないた
め、バンド理論に基づいて有効質量を定義することができません。しかしながら、実空間描像を使って局所的な有効質量が提案されており、酸化物半導体などのように式(8-2)が良い近似となるアモルファス物質もあります 50。
9. 状態密度

図 9-1 Si の(A) 状態密度 D(E) と (B) 積算状態数 N(E)。
図 9-1(A) は状態密度51 D(E) と呼ばれるもので、Si について計算したものです (図6-1(A) のバンド構造と対応させて見てください)。横軸に電子のエネルギー、縦軸にその付近のエネルギーを持つ電子が何個存在できるかをプロットしたものです。厳密には E ~ E+dE の範囲にある状態の総数 N(E+dE) – N(E) が D(E) dE と定義されますので、D(E) の単位は [状態数/eV] になります。ちょっとわかりにくいので、先に図9-1(B)を説明しましょう。これはエネルギー –15 eV から E までの準位を占有している電子の総数、積算状態数 N(E) です。これから、–15 ~ 0 eV (= EF) までを占めている電子数 N(0) は 8 個であることがわかります。もともと Si 原子 1 つには、(3s)2(3p)247
48

https://www.sciencedirect.com/science/article/pii/S0010465506001305

最新版は BoltzTraP2 になっています。
https://www.imc.tuwien.ac.at/forschungsbereich_theoretische_chemie/forschungsgruppen/prof_dr_gkh_madsen_theoretical_materials_chemistry/boltztrap/49
H. Hiramatsu, T. Kamiya, K. Ueda, M. Hirano and H. Hosono: Phys. Stat. Sol. A 207 (2010) 1636.
50
S. Kivelson and C.D. Gelatt, Jr., Phys. Rev. B 19 (1979) 5160.
51「状態」とは state の和訳で、電子準位と同一と考えてください。電子数と言い換えた方がわ
かりやすいのですが、非占有準位には電子がいないのでこのような呼び方になります。

の 4 つの電子があったわけですが、N(0) = 8 というのは、基本格子 (付録 3b 参照) に含まれる 2 つの Si 原子にある 3s, 3p 電子の総数になっています 52。このようにして、
図 9-1(B) の N(E) が積算電子数であることは納得してもらえると思います。図 9-1(A)の D(E) というのは、N(E) のエネルギー微分、D(E) = dN(E) / dE になります。実際に電気伝導度や磁化率などを決めるのは状態密度 D(E)ですから、論文などでは D(E) だ
けを示します。
負のCOHP: 結合性

図 9-2–iCOHP: COHPの積分のマイナス:
結合あたり結合エネルギー

AETiN2 (AE = Ca, Sr, Ba) の COHP53

状態密度を原子ごと、あるいは原子軌道ごとに分解すると、化学結合に関するイメージを作ることができます。図 9-1(A)では、結晶中の全状態密度 (Total)、Si 原子の全状態密度 (Si total)、Si p 軌道だけの状態密度 (Si p)と s 軌道だけの状態密度 (Si s)をそれぞれ描いています。このような図を、投影状態密度 (projected density of state: PDOS)
あるいは部分状態密度 (partial density of state: やはり PDOS)と呼びます 54。一方、位置の違う原子ごとに描いたものを局所状態密度 (local DOS: LDOS) と呼びます。Si のPDOS から、0 eV より下の領域、つまり価電子帯は p 軌道が主で、それより上の伝導帯は p 軌道と s 軌道が同程度の割合で混成していることが読み取れます。これは、Siが四面体状の結合を作っており、sp3 混成軌道で化学結合が形成されていることを反映しています。これは化学結合論で CH4 分子などを例に sp3 結合と正四面体の分子の形が説明されるのと同じです。一方で、Si の価電子帯は p 軌道性が強く、伝導帯端はs 軌道性が強く出ており、理想的な sp3 混成軌道でないことも読み取れます。後者は、伝導帯の波動関数が大きく空間的に拡がっているため、4s 成分も混成している結果で52

付録 3b 参照。Si は Bravais 格子として面心立方格子をとるので単位格子当たり 8 原子ありま
すが、物性計算ではで計算時間を短くするため、原子数が少なくなる基本格子で計算するのが一般的です。そのため、ここでの Si 原子数も基本格子内の 2 つになります。
53
X. He, T. Katase, K. Ide, H. Hosono, and T. Kamiya, Inorg. Chem. 2021, 60, 14, 10227.
54「投影」というのは数学用語で、例えば「ベクトル (0.7, 0.5, 0.2) の x 軸への投影ベクトルは(0.7, 0, 0) である」というように、ある量を、ある特徴を持つ別の量の成分として表すときに使います。ここでは、平面波基底から作られた波動関数を、特定の位置の原子や、特定の対称性を持つ関数など (つまり s 軌道、p 軌道など) へ「投影」していることから、このように呼ばれます。

す。この意味で Si の伝導帯は “s 性が強い” と言われます。
バンド計算、特に平面波 (Plane Wave: PW) 基底を使う場合には PDOS やイオン電荷の計算には大きな任意性が発生しますので、扱いや解釈には注意が必要です (付録2d, 2e 参照)。
また、化学結合の強さをより定量的に示す方法として、Crystal Orbital Hamilton
Populations (COHP) 55があります。これは、一電子状態𝜑# と𝜑R 間の化学結合のエネル
ギーを、Fock 行列の非対角項 𝐶𝑂𝐻𝑃#R = `𝜑# aℎa𝜑R b (ℎ: 一電子ハミルトニアン) を計算したものです。固体内の ℎ は負なので、𝐶𝑂𝐻𝑃#R が負であれば結合性、正であれば反結合性であると判断できます。また、関連する起動のエネルギー範囲で𝐶𝑂𝐻𝑃#R を積分
6CS%

して負の値を取った−𝑖𝐶𝑂𝐻𝑃#R = − ∫6C#G 𝐶𝑂𝐻𝑃#R 𝑑𝐸は結合エネルギーの目安になります (図 9-2)。
10. 波動関数（電子密度）

図 10-1 異なる化学結合性状を持つ半導体の CBM と VBM の波動関数。図中の球は各原子、曲面は波動関数の等高面。
量子計算をすると、波動関数も計算されます。ただし結晶軌道は、𝛷𝒌 (𝒓) = 𝑒 #𝒌⋅𝒓 𝑢(𝒓)のように、Bloch の定理によって波数ベクトル k による位相 𝑒 #𝒌⋅𝒓 がかかります。まともに𝛷𝒌 (𝒓) の実数部と虚数部を描いてみても、k が 0 でない場合には余計な位相項の変化が重畳して何がなんだかわかりません。そのため、多くの場合に、その軌道の電子密度に対応する|𝛷𝒌 (𝒓)|" をプロットします。
図 10-1 には、Si, GaAs, GaN の CBM と VBM の波動関数 |𝛷𝒌 (𝒓)|" (正しくは「一電子波動関数のノルムの自乗」です) を描いてあります。半導体の場合には、これらはそれぞれ、n 型半導体での電子の通り道と、p 型半導体での正孔の通り道を可視化したものと捉えることができます。さて、9 章の PDOS で見てきたように、Si の VBM
は sp3 混成軌道でできていますが、p 軌道の寄与が大きくなっています。このことは
図 10-1 でも、VBM の波動関数が最近接の Si–Si 方向の結合上に大きい密度を持つことから見て取れます。一方 CBM の方は s 性が強く、かなり丸みを帯びた形になっていますが、それでも、Si–Si 結合上に高い密度があることがわかります。Si–Si 結合の中央で波動関数の値がゼロになる、いわゆる「ノード (節点, node)」がありますが、
これは、CB が sp3 の反結合軌道でできていることを反映しています。Si は一種類の55

http://www.cohp.de/原子だけで構成されているので完全共有結合性ですが、GaAs や GaN では、Ga は電子親和力が小さくて陽イオンになりやすいのに対し、As と N は Ga よりも電子親和力が大きく、陰イオンになりやすいですね。As と N を比べると、N の電子親和力の方が大きくなります。そのため、Si => GaAs => GaN の順番で、化学結合のイオン性が強くなります。このことは、GaAs の VBM がほとんど As の軌道からできており、逆に CBM は As の軌道もあるものの、Ga の s 軌道が強いことに現れています。GaN
はこの傾向をもっと強くしており、CBM には N の軌道があまり現れていません。
全電子密度、価電子密度、up/down スピン密度、電位ポテンシャル分布なども同様に可視化することができ、様々な情報を得ることができます。
11. バンドギャップ問題とその起源: 密度汎関数の分類と特徴
さて、図 6-1(A) から読み取れる Si の最小のバンドギャップ（間接バンドギャッ
プ）は 0.6 eV でしたが、これは、300 K での実測値 1.1242 eV よりずいぶん小さい値です。絶対零度への外挿値ではさらに大きく 1.17 eV になりますから、測定温度の問
題ではありません。これが DFT の最大の問題と言われる、いわゆる「DFT のバンドギャップ問題 (Bandgap problem)」です。DFT ではバンドギャップをかなり過小評価しますので、バンドギャップの定量的な議論には使えません。逆に HF 近似には、バンドギャップを過大評価する「バンドギャップ問題」があります。

2.5
2.0

EVBM

ECBM
HSE0

80

PBE0
HF

60

1.5

GGA+U U=5eV

40

1.0

GGA+U U=3eV

20
US LDA
PBE96
HF
PBE0
HSE06
U=1 eV
U=2 eV
U=3 eV
U=5 eV
U=7 eV

0.5

(B)

100

DOS

Band gap (eV)

3.0

120

(A)

11.2 eV =>3.5

0 1 2 3 4 5 6 7 8 9 10 110-10

GGA(PBE96)
LDA

-5

0

5

Photon energy / eV1015

図 11-1 異なる汎関数を用いて計算した ZnO の(A)バンドギャップと(B)状態密度(VASP)。バンドギャップの実測値は 3.37 eV。
図 11-1 に、異なる汎関数 (vXC) を用いて計算した ZnO のバンドギャップと状態密度を示します。ZnO のバンドギャップの実測値は 3.37 eV ですが、HF 近似の結果は11.2 eV と、4 倍近く過大評価しています。一方で、左図左端の局所密度近似 (Local
Density Approximation: LDA) では 0.7 eV と、逆に過小評価しています。LDA のバンドギャップ問題が指摘された当初、その問題は空間の一点 r だけの関数 r(r) の汎関数
により計算している (これが ”局所密度” 近似と呼ばれる理由です) ためと考えられ、そのため、r の周囲の情報も含めるように r(r) の空間微分 Ñr (r ) を密度汎関数に取り 入 れ た 非 局 所 密 度 近 似 が 提 案 さ れ 、 現 在 の 一 般 化 密 度 勾 配 近 似 (Generalized
Gradient Approximation: GGA) へと至っています。GGA でもっとも広く使われている
PBE 汎関数 (汎関数については付録 2b 参照) を用いた結果を見ると、バンドギャッ
プは LDA よりは若干改善されているものの、依然として大きく過小評価しています。
このように、第一原理計算でバンドギャップを正確に計算するのは現在でも難しい問
題です。

DFT におけるバンドギャップ問題の起源の一つは、自己相互作用 (Self-interaction:SI) を補正できていないことにあります。(3-2)式の HF 方程式を水素原子のような一電子系で具体的に書き表すと、
V(𝒓 )𝑉$)$ (𝒓* ) = ∫ |𝒓 2 | 𝑑𝒓C(11-1)
2 )𝒓3となりますが、これは、自分自身の電荷が自分自身と相互作用する「自己相互作用」
であり、取り込んではいけないエネルギーです。一方で交換相互作用がV(𝒓 )𝑣+* (𝒓* ) = − ∫ |𝒓 2 | 𝑑𝒓C(11-2)2 )𝒓3となり、𝑉$)$ (𝒓* ) + 𝑣+* (𝒓* ) = 0と、自己相互作用をキャンセルし、結果として!
!− 𝛻* " + 𝑉$%& (𝑟)+ 𝜑* (𝒓* ) = 𝜀* 𝜑* (𝒓* )(11-3)"と、正確な Schrödinger 方程式に一致します。つまり、HF 方程式には自己相互作用の問題がありません。同様のことを (3-4) 式の KS 方程式で行うと、例えば基本的な LDAの交換相互作用を使って>!/>

𝑣+0 (𝒓) = −3𝛼 !XY 𝜌(𝒓)+(11-4)となり、𝑉$)$ (𝒓) + 𝑣+0 (𝒓) ≠ 0となり、自己相互作用が完全にキャンセルされず、誤差として残っていることがわかります。この効果は図 11-2 の Ne = 1.0 にみられるように、軌道準位を浅く過小評価します(Slater の遷移状態法 (脚注 21) では、軌道の占有電子数が Ne = 0.5 で、Koopmans の定理のイオン化エネルギーが得られると近似されます)。このため、自己相互作用を補正することでより正しい電子構造が得られる「自己相互作用補正 (Self-interaction correction: SIC) 法」も使われています。

図 11-2 電子数 Ne を変えながら、水素原子の 1s 軌道エネルギーE1s を LDA (式(11-1),(11-4)で計算した結果。Ne ~ 0.5 で正しい E1s = 13.6 eV に一致しますが、Ne ~ 1 では𝑉$)$ (𝒓) > 𝑣+0 (𝒓)となり、軌道エネルギーを浅く過少評価します。”(optimized)”
は、H 1s 動径関数の指数係数を変分パラメータとして最適化した結果ですが、その場合でも全体の傾向は変わりません。
バンドギャップ問題の他の原因としては、非占有準位の問題があります。被占有準位というのは、現実に電子がいる状態そのものを計算した結果です。ところが非占有状態というのは、量子理論の方程式を解いた結果として現れた、数学的な解にすぎません。計算しているときには存在していない電子の状態ですから、その結果に物理的意味はありません。さらに、密度汎関数法の基礎理論である Hohenberg-Kohn 理論では、基底状態の正確さしか保証していません。そのため、被占有軌道を real state、非占有軌道を「仮想軌道 (virtual state)」と呼ぶことがあります。ですから、物理的意味

の無い仮想軌道とのエネルギー差であるバンドギャップの計算値が実測値に合わないからといって、DFT や HF 近似の精度が悪いということではありません。非占有軌
道やバンドギャップは DFT の前提外であるということです。
ただし、「仮想軌道」とはいっても、経験的に、非占有準位の構造が励起状態の電子構造をかなり良く反映していることが確認されています。そのため、これを「伝導帯」の電子構造として解釈するということが一般的に認知されています。「バンドギャッ
プ問題」に対しては、実験値のバンドギャップを再現するように、「非占有準位全部を適当なエネルギー Dsc だけ平行移動してしまう」という対応がされ、Dsc を「Scissors
operator」と呼びます 56 (図 6-1(B))。
また、実験的にバンドギャップを測定する場合と計算モデルでの条件の違いも考慮する必要があります。半導体のバンドギャップは多くの場合、光学吸収スペクトルから決められる「光学ギャップ Eopt」です。この場合は、価電子帯端から電子を一つ取り除き、伝導帯端に電子を一つ加えた場合の全エネルギー変化をバンドギャップとして測定していることになります。つまり、
dE
dE
Eopt = {E (M + 1) - E (M )}- {E (M ) - E (M - 1)} =dn N = M +0 dn N = M -0
N =M +0(11-5)dE{}= e M +1 ( M + 0) - e M ( M - 0) = e M +1 ( M ) - e M ( M ) DFT +dn N = M -0
57となります 。ここで、E(N)は電子が N 個ある場合の全エネルギー、N = M は半導体
で価電子帯が完全に占有されている電子数、N = M + 0 は CBM に微少量の電子がある
状態、N = M – 0 は VBM に微少量の正孔がある状態、eN(M) は M 個の電子があるときの N 番目のエネルギー準位で、Janak の定理より最後の式の第一項が DFT で計算したバンド構造から得られるバンドギャップになります。これに対して、最高被占有準位が VBM から CBM へ飛ぶ際に dE / dn が不連続に変わると、光学ギャップは DFT で計算されるバンドギャップよりもこの値だけずれることになります。このような電子相関効果が「バンドギャップ問題」の原因の 1 つです。その電子相関の効果は、被占有準位を深く、非占有準位を浅くしますが、実際の DFT の計算結果では、被占有準位を浅く計算していることがバンドギャップを過小評価する主要因であるといわれています 58,59。
Beyond DFT 法 と し て 、 電 子 密 度 の 二 階 微 分 Ñ 2 r (r ) や 運 動 エ ネ ル ギ ー
t = (- 2 / 2m)Ñ 2f (r ) を取り入れた meta-GGA 法も使われています。Modified BeckJohnson (mBJ) という汎関数は、比較的実測値に近いバンドギャップを出してくれま
す60。
また、HF 近似と DFT がバンドギャップの値を逆方向に間違えることを利用し、これらを適当な割合で混合した混成汎関数が使われることが多くなっています。式(3-2)の HF 項 𝑣+* は波動関数の汎関数であり、電子密度の汎関数ではありません。この点
56

物理学者というのは、普段は厳密に数式を扱っているように思えますが、必要なときには極端な単純化や、非常に荒っぽい手段も平気で使います。
57
文献 5 の 9 章、文献 8 の 7.3 章
58
N. Hamada and S. Ohnishi, Phsy. Rev. B 34 (1986) 9042.
59
M.S. Hybertsen and S.G. Louie, Phys. Rev. B 34 (1986) 5390.
60
一電子方程式に含まれるポテンシャル汎関数𝑣+0 しか与えられておらず、エネルギー汎関数
𝐸+0 がわかっていません。そのため、全エネルギーに関して自己無撞着計算ができないという欠点があります。

は、KS 法を電子密度よりも一般的な汎関数に拡張した一般化 KS 法61 により条件付
きで正当化されています (文献 8)。HF 項は Exact Exchange (EE) 項とも呼ばれます。
平面波基底を用いたバンド計算では EE 項のため、DFT に比較して計算時間が非常に長くなるという問題がありますが、それでも後述の GW 近似よりは軽いことから、広く使われるようになっています。最近使われている実験値に近いバンドギャップを与える汎関数は、だいたい次のグループに分けられます。
(i) EE 項に遮蔽パラメータを入れる方法 (Screened Exchange: sX)(ii) 実測のバンドギャップに合うように EE 項と密度汎関数を適当な比で混合した混成汎関数。B3PW91, B3LYP など (Becke の 3 パラメータ法で PW91, LYP 汎関数を用いたもの)。混合係数が経験的に決められている。
(iii) (ii) に似ているが、係数が理論的に決定されている混成汎関数。PBE0 62 や、
その EE 項に sX を用いた HSE 63 など。
(iv) その他。meta-GGA 汎関数である modified Becke Johnson (mBJ) は多くの半導体において実測値に近いバンドギャップを与える 58,64,65。
このうち、特に (ii) のグループは、混合係数に任意性がでてくるため、第一原理法の概念からは一歩後退したように思えます。しかしながらそもそも、DFT では汎関数の選択自体に任意性があるため、任意の汎関数の和は常に DFT 理論の枠内で正しい汎関数になります (ここで “正しい” というのは、理論の枠組みのことであり、正確な計算結果を与えるという意味ではありません)。そのため、汎関数の係数に経験パラ
メータがある場合も、現在の一般的な第一原理計算とは “第一原理” の程度は変わりませんので、これらすべてが “第一原理法” として分類されることがあります。ただし、図 11-1 を見てわかるように、(iii) の場合も含め、混成汎関数を使ってもバンドギ
ャップの計算値やイオン化エネルギーが合わないことも珍しくありません。(iii) のPBE0 や HSE の場合でも、混合係数や遮蔽定数などを変えて実測のバンドギャップに合わせることもされています 66。そのため、これらを “混成汎関数法” と呼び、第一原理法と呼ばない研究者もいます。
任意パラメータを使う簡便な方法に DFT+U 法があります。これは、d 軌道あるいは f 軌道のような局在性の強い軌道のみに、同じ軌道内の電子―電子反発ポテンシャ
ル U (Coulomb パラメータ) と交換ポテンシャル J (Exchange パラメータ) を補正として考慮する方法で、これらの軌道エネルギー値や、遷移金属化合物のように d 軌道がバンドギャップを作っている場合にはバンドギャップの値を実測値に合わせることができます。この場合は任意パラメータが入ってくるので、厳密な意味で第一原理法と呼ぶには抵抗があります。一方、U, J パラメータを第一原理計算から見積もる方
法もあります (Constraint LDA 法 67)。
バンドギャップや励起状態も含めた計算として最も精度が高い計算法として GW
61

A. Seidl, A. Görling, P. Vogl, J.A. Majewski, and M. Levy, Phys. Rev. B, 53 (1986) 3764.
J.P. Perdew, M. Ernzerhof and K. Burke, J. Chem. Phys. 105 (1996) 9982.
63
J. Heyd and G. E. Scuseria, J. Chem. Phys. 118 (2003) 8207 の ERRATUM が 124 (2006) 219906 に出ているため、HSE03/HSE06/HSE などの呼称が使われています。
64
A.D. Becke and E.R. Johnson, J. Chem. Phys. 124 (2006) 221101.
65
F. Tran and P. Blaha, Phys. Rev. Lett. 102 (2009) 226401.
66
例えば F. Oba, A. Togo, I. Tanaka, . Paier and G. Kresse, Phys. Rev. B 77 (2008) 245202 では HSE のパラメータを標準の a = 0.25 から 0.375 に増やしています。
67
O. Gunnarsson, O.K. Andersen, O. Jepsen and J. Zaanen, Phys. Rev. B 39 (1989) 1789; I.V. Solovyev
and P.H. Dederichs, Phys. Rev. B 49 (1994) 6736.
62近似が使われるようになっています。しかしながら、励起状態を正確に計算するために、空のバンド数、平面波のカットオフエネルギー、誘電関数の周波数などを大きくとる必要があり、混成汎関数法よりもはるかに長い計算時間がかかります。また、基底状態の計算モデル (DFT、DFT+U、混成汎関数などのどれを使うか) でも結果が変わりますが、なるべく実測に近いバンドギャップを与えるモデルで計算した結果を基底状態に使う方が良いようです。
分子などの孤立系の励起エネルギーの計算では時間依存 DFT (TDDFT) 法も使われています (例えば Gaussian 68 などの分子軌道法プログラム) 。比較的短い時間で計算ができますが、結晶のバンド計算プログラムで実装しているものは多くありません。
交換相関カーネルの構築の工夫によりバンドギャップの計算値は改善できるということです (例えば Advance/PHASE 69 には post-SCF 計算として Linear-response TDDFT(LR-TDDFT) が実装されており、誘電関数や光吸収断面積の計算が可能です)。
+ 非占有状態+ 占有状態

一般化乱雑位相近似 (RPA)
HF Exact Exchange
混合汎関数法+ eXhyper-GGA2
2
2+ t = (- / 2m)Ñ f (r ), Ñ r (r ) meta-GGA+ Ñr (r )GGA
LDA/LSDAr (r )
図 11-3 Jacob’s ladder による汎関数の階層的構造。
5 ∆ PBE96
◊ LDA

-1
0

NnO
NiO

LaCuOSZnO

LaZnOP

0

LaZnOAs

1FeO
CoO
SrGeO3
LaCuOSe

2

LaMnOPLaMnOAs

3● HSE▼ PBE0♦ B3LYPLaMnOSb

Calculated Eg (eV)

4 ○ PBE96+U(1.5 eV)12
34Experimental Eg (eV)

図 11-4 異なる計算モデルによるバンドギャップの計算値と実測値の関係。左図は文献 70 より。右図は文献 70 の遷移金属酸化物 (FeO, CoO, MnO, NiO) のデータと表 11-1 のデータより作成。
上述の SIC 法も良いバンドギャップやイオン化エネルギーを与えると言われていますが、一電子軌道の取り方で計算結果が変わり、電子軌道の物理的意味が曖昧になります。
以上のように、現在量子計算で使われている汎関数・近似モデルには様々なものが68
69

http://gaussian.com/
http://www.advancesoft.jp/product/advance_phase/

あり、”Jacob’s ladder 70”として図 11-3 のようにまとめられています。
図 11-4 に、異なる計算モデルによるバンドギャップの計算値と実測値を比較した
結果をまとめています。左図 71 からは、バンドギャップが小さい場合には HSE より
も PBE0 の方が実測に近い値を与え、大きい場合には PBE0 がより良い値を与える傾向を示しています。
また、一般的に、DFT で計算したバンドギャップの大小の傾向は、類似した物質系
では直接比較できるとして議論されることもあります。しかしながら、文献 70 と表
11-1 のデータから作成した右図を見ると、DFT (LDA, PBE) で得られるバンドギャッ
プが実測値と単調な相関があるわけではないことが読み取れます。特に立方晶 SrGeO3の例にみられるように、実測値が 2.7 eV と大きいものでも、DFT では負のバンドギャップを与える場合もあることもわかります。表 11-1 からは、遷移金属を含む場合には混成汎関数ではバンドギャップを過大評価し、むしろ EE 項の比率を下げないといけない場合があることもわかります。このように、構造や構成元素が大きく異なる物質の DFT のバンドギャップの比較はできないと考えるべきです。
表 11-1 密度汎関数法 (GGA PBE) と混成汎関数法 (HSE) で計算したバンドギャッ
プと実測値 (Eg(obs)) の比較 (単位 eV)。*で示した Mn 化合物のバンドギャップ Eg(PBE) として書いてある値は、Mn 3d 軌道に U = 1.5 eV を用いた PBE+U 計算の結果である。**で示した ZnO、SrGeO3 の Eg(HSE) は PBE0 を使った結果。CEE は、HSE 計
算における EE 項の比率。
CEE EgEg
Eg
CEE Eg
Eg
Eg(HSE) (PBE)(obs)(HSE) (PBE) (obs)
LaCuOS
0.25 3.00
1.72
3.1
LaMnOP
0.25 1.88
0.871* 1.3
0.30 3.26
0.15 1.37
LaCuOSe 0.25 2.78
1.55
2.8
LaMnOAs 0.25 1.74
0.878* 1.4
0.30 2.97
0.15 1.53
LaCuOTe 0.25 2.27
1.21
LaMnOSb 0.25 1.17
0.586* 1.0
LaZnOP
0.25 1.50
0.63
1.7
0.15 0.89
72
LaZnOAs 0.25 1.38
0.55
1.5
SrGeO30.25 2.7**
-0.4
2.7
**
LaZnOSb 0.25 0.57
0.03
ZnO
0.25 3.17
0.7
3.37
12. 密度汎関数の選択
11 章でみたように、安定構造、電子構造、物性の計算結果は、密度汎関数の選択に左右され、特に半導体で重要なバンドギャップでは密度汎関数の選択は重要になります。一般的には金属には LDA が適しており、半導体の場合は GGA の方がより適切な場合が多いです。安定構造について汎関数を比較したのが表 12-1 です。
また、一般に LDA では 300 K の実測格子定数を< 2%過小評価します。逆に、PBEは 300 K の格子定数を同程度過大評価しますが、構造緩和計算が絶対零度における安定構造であることを考えると、フォノン計算など構造敏感な計算に使うには、PBE の安定構造では精度が悪いことがあります。PBEsol 汎関数は、300 K での実測構造の再現性が高いと言えますが、バンドギャップは PBE よりも小さくなる傾向にあります。

70

J.P. Perdew and K. Schmidt, AIP Conf. Proc. 577 (2001) 1; J.P. Perdew, A. Ruzsinszky, J. Tao, V.N.
Staroverov, G.E. Scuseria, G.I. Csonka, J. Chem. Phys. 123 (2005) 062201.
71
M. Marsman et al, J. Phys.: Condens. Matter, 20 (2008) 064201.
72
H. Mizoguchi, T. Kamiya, S. Matsuishi, H. Hosono, Nat. Commun. 2 (2011) 470.

表 12-1 VASP で使える汎関数で計算した ZnO の構造緩和計算結果と、300K の実測構造との誤差。
VASP
GGA tag

Functional

a

c

z

Δa(%)

Δc(%)

PE

exp (300 K)
ICSD29272
PBE

PS

PBEsol3.24195.19330.38163

-0.03

-0.03

RE

revPBE

3.28805.28120.38155

1.40

1.66

RP

RPBE

3.32685.37000.38158

2.59

3.37

AM

AM05

3.25515.2236

0.38202

0.38

0.55

91

PW91

3.28285.27570.38162

1.24

1.56

B3

B3LYP

3.27235.26240.38179

0.91

1.30

3.24275.19480.3826

3.28765.2821

0.38163

1.38

1.68

図 12-1 VASP で使える汎関数で計算した ZnO のバンドギャップの値。MBJ が mBJ、
LMBJ は派生版の local mBJ の結果であり、PBE0 で得られるのと同程度の値が得られる。
上述のように、PBEsol は構造の再現性が高い一方でバンドギャップは小さく計算
し、mBJ は実測のバンドギャップの再現性が高いですが、エネルギー計算ができません。そのため、私たち研究室では、以下のような汎関数の選択をベースに考えています。
(i) まず、あたりをつける。計算の速い DFT で構造緩和、電子構造計算を行う。
金属の場合は LDA、バンドギャップが開いている場合は PBE や PBEsol など。
(ii) 以下はバンドギャップが開いている場合。
(ii) エネルギー計算（構造緩和、反応エネルギーなど）: PBEsol(iii) 電子構造 (バンド構造、電子構造): mBJ、あるいは電子配置等によって構造に影響が出る場合などは混成汎関数 (バンドギャップが 1.0 eV 程度であれば HSE、大きければ PBE0)。場合によっては実測値を再現する汎関数を選択する。
(iv) キャリア輸送特性: (iii)で計算した電子構造を使い、BoltzTraP2 を使って電子構造

を補間して計算(vi) フォノン計算: (ii)で計算した安定構造を使い、VASP、phyonopy (調和フォノン近似)73、ALAMODE (非調和フォノン計算)74 などを使って計算。汎関数は(ii)に合わせる
13. 光学スペクトル
量子計算によって波動関数がわかると、外部交流電界下の摂動理論により誘電関数の虚部スペクトル e2(hn) が計算できます。電子遷移については、バンド間 (interband)と バンド内 (intraband) 遷移の効果があります。バンド内遷移とは、同じバンド内で電子の k が変わる遷移で、要するに、電界などで電子が加速される効果です。半導体の評価では自由キャリア吸収と呼ばれ、Drude の自由電子モデルで良く記述されます。
半導体で特に興味がもたれるのがバンド間遷移で、この場合は e1(hn), e2(hn) は以下の式で表されます 75。
@ 4 4

XY$ +

𝜀!,## 4 (ℎ𝜈) = 1 + 3C ℏ" ∑*,*[ 6 +33)6,"" +(13-1)

𝜀",## 4 (ℎ𝜈) = 3C ("Y\) ∑*,*[ 𝑓**4 ,## 4 2𝑓(𝐸* ) − 𝑓(𝐸*[ )4𝛿(𝐸*[ − 𝐸* − ℎ𝜈)(13-2)

𝑓**4 ,## 4 = ℏ+ (𝐸* − 𝐸*[ ) < 𝑙′|𝑥# |𝑙 >𝒌 < 𝑙|𝑥#[ |𝑙 [ >𝒌(13-3)"Y+ $ +

%

3

34

%

"C%

ここで𝑓**4 ,## 4 は振動子強度と呼ばれます。WIEN2k の誘電関数計算をしている OPTICSなどでは、まず e2(hn) を計算してから、Kramers-Kronig 変換することで e1(hn) を計算してます。式(13-1) のように、これらの式では分母に遷移エネルギーが入っており、DFT の過小評価したバンドギャップを使うと、吸収端エネルギーだけでなく、誘電関数の値にも大きな影響を与えます。そこで、WIEN2k や CASTEP などでは、Scissors
operator を用いてバンドギャップを実測値にあわせてから計算することで光学スペク
トルを計算します。

(A)

E1

E2

e1

e

80
70
60
50
40
30
20
10
0
-10
-20
-30
-40

e2
2

2.5

3
3.5
4
4.5
Photon energy / eV5

図 13-1 (A) Si の複素誘電関数スペクトル。シンボルが実測値、実線が計算値(CASPEP, PBE)。(B) LaCuOSe の光吸収スペクトル。実測値(上)と計算値(下) (WIEN2k(PBE) + OPTICS)。
図 13-1(A) には CASTEP で Scissors operator の補正をした Si の結果を示しますが、
e1(0) 76 の値は実験値とよく一致しています。図 13-1(B) の LaCuOSe の場合も、矢印
73

https://phonopy.github.io/phonopy/
https://alamode.readthedocs.io/en/latest/75
花村榮一著、固体物理学 第２版、裳華房 (1987).
76
誘電関数 e は、屈折率 n、消衰系数 k と e = (n - ik)2 の関係があり、吸収係数 a は光の波長74で対応させている特徴的な構造をよく再現できていることがわかります。ただし、
DFT のエネルギー固有値はバンドギャップの補正をしても遷移エネルギーに対応しませんから、それぞれのスペクトル構造のエネルギーには ずれがあります。実測ス
ペクトルで “125 meV” と書いて縦線で示しているのは励起子 77 による鋭い吸収線です。DFT では励起子のような励起状態は考慮していませんので、これに対応する構造は、計算した吸収スペクトルの肩として見えているだけです。
式(13-3) に出てくる積分 < 𝑙|𝑥#[ |𝑙 [ >𝒌 は遷移双極子モーメント Tll’ で、波数 k が軌道間の遷移の場合は𝑇**4 = ∫ Ψ*[.[ ∗ (𝒓)𝒓Ψ*,𝒌 (𝒓)𝑑𝒓 = ∫ 𝑢*[ ∗ (𝒓)𝒓𝑢(𝒓)exp[𝒊(𝒌 − 𝒌[ ) ∙ 𝒓]𝑑𝒓
(13-4)となります。ここで Bloch 関数の条件 Ψ*,𝒌 (𝒓) = 𝑢* (𝒓)exp (𝑖𝒌 ∙ 𝒓) を使っています。右辺で𝑢* (𝒓) は結晶格子と同じ周期性を持ちますので、前半の被積分関数も周期関数になります。𝒌 ≠ 𝒌[ の場合には exp の項は周期的に正負が逆転するため、この積分は必ずゼロになり、この電子遷移は双極子近似の範囲内では起こらないことになります。
これがいわゆる「k 選択側」です。
14. 全エネルギー: 生成熱など
上述のように、現在の第一原理計算では、入力した構造モデルの全エネルギーを高い精度で出力してくれます。WIEN2k で PBE を使うと、Si の全エネルギーは-1160.141
Ry 78 が得られます。これは、基本格子に含まれる Si 原子 2 つ分のエネルギーです。
1 原子あたり 8 keV になりますから、普段私たちが研究で接しているエネルギーから比べると数桁大きいことがわかります。これは、WIEN2k では全電子 L/APW+lo 法を用いており、1s 軌道からのすべての電子軌道までエネルギーを加算しているためです79。私たちが普段測定するような物性量、例えば生成エネルギーなどと対応させるには、化学反応式の両辺にでてくる物質の全エネルギーを計算し、そのエネルギー差を求める必要があります。
例 14-1: Na の昇華熱 (Na (結晶) => Na (原子)のエネルギー差) (VASP, PAW, PBE96)※ Na(結晶) の全エネルギー: E = -2.62 eV/cell※ Na(原子) の全エネルギー: E = -0.22 eV/atom※ Na(結晶) => Na(原子)のエネルギー差: DE = 1.09 eV = 105 kJ/mol※ Na(原子) 気体の PDV 項 RT = 2.49 kJ/mol (300 K)を足してエンタルピーにする:
l を用いて a = 4pk / l と求められます。
77

物質に光を入射すると、価電子帯の電子を伝導帯に励起し、電子と正孔を作ります。バンド計
算では、このような電子も正孔も周期結晶中に無限に拡がった構造としてしか扱いません (もともと周期構造を仮定しているため) が、実際には、励起直後の電子と正孔は近傍の位置で局在している (相対論の「いかなる相互作用も光速以上では伝播しない」という要請を考えれば当たり前です) ため、それぞれの電荷で引きつけあいます。このエネルギー (励起子の束縛エネルギー) が熱エネルギーよりも十分大きければ、束縛状態の励起子が安定になり、LaCuOSeのように鋭い吸収線や発光線として観測されます。Si や GaAs の励起子束縛エネルギーは室温よりも小さいため、通常は問題にされません。
78
Ry は Rydberg の原子単位で、エネルギーについては 1 Ry = 13.6 eV (水素原子の 1s 軌道エネル
ギー)。他の量子計算プログラムでも出力、入力が原子単位になっているものが多く、また、
Hartree の原子単位も使われているので注意が必要です (1 Hartree = 27.2 eV)。
79
エネルギーの浅い価電子状態しか考慮しない擬ポテンシャル法では、全エネルギーの計算値はかなり小さくなります。例えば VASP (PAW, PBE)で Si の計算をすると-10.844 eV が得られます。

DH = 107 kJ/mol

文献値: 108 kJ/mol 80

例 14-2: NaCl の生成エネルギー 81※ NaCl(結晶) の全エネルギー : E = -27.26 eV/cell (4NaCl)※ Na(結晶) の全エネルギー: E = -2.62 eV/cell (2Na)82※ Cl2(分子) の全エネルギー: E = -3.57 eV/cell (2Cl)※ Na(原子) の全エネルギー: E = -0.220 eV/atom※ Cl(原子) の全エネルギー: E = -0.296 eV/atom83※ 凝集エネルギー
NaCl(結晶) => Na(原子) + Cl(原子):-7.33 eV/NaCl = -608 kJ/mol 文献値 641 kJ/mol※ Cl-Cl の結合エネルギー: E = -0.290 eV/Cl-Cl bond = 288 kJ/mol
文献値 242.8 kJ/mol※ 生成エネルギー NaCl(結晶) => Na(結晶) + 1/2 Cl2(分子):-3.720 eV/Na = 358.7 kJ/mol
文献値 411 kJ/mol
例 14-3: Si の凝集エネルギー (結合解離エネルギー)
NaCl の例と同様、Si(結晶)と Si(原子)の SCF 計算を行い、全エネルギー差をとります。
※ Si(結晶) の全エネルギー: E = -43.37 eV = -5.42 eV/atomSi(原子) の全エネルギー: E = -0.86 eV
凝集エネルギー Si(結晶) => Si(原子)-4.56 eV/Si = 440 kJ/mol
文献値 446 kJ/mol 84※ 結合解離エネルギーは、これを結合数で割ればよい。Si 原子１つは 4 つの SiSi 結合をもちますが、結合は 2 つの原子で共有されているので、
Si-Si の結合解離エネルギー: E = 220 kJ/mol 文献値 224 kJ/mol 85

15. 単純な結晶の安定構造（構造緩和計算）と体積弾性率
さて、第一原理計算で全エネルギーが求まるということは、それが最小になるような格子定数や原子座標を探すことで最安定な構造を計算できることを意味しています。このような計算を構造緩和計算と呼びます 86。ただしここで注意が必要なのは、DFT で求められる電子状態は基底状態であることと、原則として電子分布は絶対零度の Fermi-Dirac 分布を仮定しますし、原子核の運動は考慮していないので、得られた80『バーロー物理化学、藤代亮一訳、東京化学同人』より。
生成エネルギー: 一般的には化学反応式の右辺と左辺の全エネルギーの差であるが、ここでは、化合物をその原料の安定な単体から合成するのに必要なエネルギーの意味（標準生成エンタル
ピーの定義）。
82
分子や原子の全エネルギーは、大きい単位格子内に一つの原子・分子を入れた構造をつくり、
k 点を G点だけに限定した計算 (Gamma-Only 計算) から求めます。
83
凝集エネルギー: 0 K の凝集体を、構成する原子に分けてばらばらにするのに必要なエネルギー。似た単語に格子エネルギーがある: 0 K の結晶をイオンに分けてばらばらにするのに必要なエネルギー。
84『理化学辞典 (岩波書店)』より。
85『化学便覧基礎編 II (第 3 版, 丸善) 』より: (CH3)6Si2 で測定した値。
86
収束の安定性確保や計算時間の関係で、格子定数を固定して原子座標のみの緩和計算がされることも多い (Fixed-cell relaxation などと呼ばれる)。これと区別して、すべての格子定数まで最適化する計算を Variable-cell relaxation、Full relaxation と呼ぶことがあります。
81結果は、「基底状態、絶対零度」での構造になります。

Energy / Ry

-1160.138-1160.139

Exp.(RT)
aC = 0.5431 nm
VR = 270.5 a.u.3 (primitive cell)Opt.
aC = 0.5472 nm
VR = 276.67 a.u.3-1160.142

図 15-1(B)
Exp.(RT)
ac = 0.422 nm
VR = 126.9 a.u.3 (primitive cell)-551.427-551.428

Opt.
ac = 0.4254 nm
VR = 130.0164 a.u.3-551.429-1160.143-1160.144-551.425-551.426-1160.140-1160.141-551.424(A)
Energy / Ry

-1160.137-551.430

260

270

280
290
Volume / a.u. 3300-551.431

120

124

128
132
Volume / a.u. 3

136

140

Si (A)と MgO (B) の単位格子体積に対する全エネルギー変化 (WIEN2k, PBE)。

Si の例を見てみましょう。図 15-1(A) のように、V0 = 276.69 a.u.3 87 に極小を持つエ
ネルギー曲線が得られます。横軸の体積は基本格子である菱面体格子のものですので、計算しやすいように、4 倍してブラベー格子である面心立方格子の体積に変換してから格子定数を求めてみましょう。a = (4V0)1/3 = 0.5472 nm が得られます。文献値は 296
K で 0.5431 nm ですから、計算誤差は 0.8% にすぎません。ただ、実験値は 296 K での値ですので、0 K での外挿値 (ac = 0.5426 nm) と比較するのが自然ですが、いずれにしても大きな違いはありません。また、全エネルギーE と体積 V は、
2
E = Emin + B0 (V /V0 ) の関係がありますから、極小点の曲率から体積弾性率が B0 = 87.6

GPa と求まります。実験値は 298 K において 97.88 GPa ということですから、11% 程度の誤差になります。
MgO の場合は、V0 = 130.02 a.u.3 に極小を持ちます。MgO も FCC 構造を持ち、このV0 は基本格子である菱面体格子の体積なので、Si の場合と同様の手順で a = 0.4254 nmが得られます。文献値は室温で 0.422 nm ですから、この場合も計算誤差は 0.8% にすぎません。体積弾性率は B0 = 157 GPa であり、室温の実験値 162 GPa と 3%程度の誤差で一致します。
ここでもう一つ、量子計算の特徴と注意を書いておきます。図 15-1 をみてわかるよ
うに、計算で得られる全エネルギーは 1000 Ry 程度 (WIEN2k は全電子 L/APW 法なので、内殻電子の大きなエネルギーの和となるために、全エネルギーは非常に大きな値になります。重原子を含む場合はさらに桁で大きくなります) ですが、力学物性に効くのは mRy 以下の変化です。量子計算で安定構造や物性を正確に計算するには 8~10桁以上の計算精度が必要といわれる理由です。
16. 一般的な構造緩和計算: 安定構造と局所構造
Si や MgO の結晶構造は対称性が高く、自由パラメータが格子定数の a だけしかありません。そのため、a を横軸にとった図 15-1 のようなグラフだけで最安定な構造を決められますが、一般の結晶では自由パラメータが多くなり、グラフからは最安定な構造を決められなくなります。一般の場合は、繰り返し逐次計算による非線形最適化
88によって最小エネルギーを実現する自由パラメータを探し、最安定な構造を決めま87

ここでも原子単位が使われています。長さの単位は Bohr 半径 a0 = 0.05292 nm で、単位として ’bohr’ も使われます。
88
共役勾配 (Conjugate Gradient: CG) 法、最速降下 (Steepest Descent: SD) 法などの最適化アルゴ
リズムや、焼きなまし (Simulated Annealing) 法、熱運動を抑えた分子動力学 (dumped MD) 法す。
表 16-1 に代表的な構造緩和計算の結果をまとめていますが、ほとんどの場合に1~3% 以内の精度で構造パラメータが決まっていることがわかります。現在の DFT では、基底状態の電子構造と結晶構造については、この程度の精度が得られます。逆に、これよりも実測値との一致が悪い場合は、(i) 構造モデルが間違っている、(ii) 計算条
件 (Ecut、k 点数 89、L/APW 法の RMT など) が不適切、(iii) 密度汎関数が適切でない、(iv) スピン分極を考慮しないといけない、(v) スピン配置構造が間違っている、などをまず疑います。これらが問題ない場合、例えば強相関電子系である、分子間力結合(van der Walls 相互作用) などの励起状態を含む摂動項が重要な系であるなど、DFT を適用することの妥当性について再検討する必要があります 90 (系統的に DFT の精度がよくない系については、付録 2j を参照)。また、一般的な傾向として、LDA では格子定数を過小評価し、GGA では過大評価します。構造パラメータの計算精度は、
DFT+U、混成汎関数や GW 近似を用いることでさらに改善されるのが一般的です。
表 16-1 種々の結晶の構造パラメータ。括弧内は緩和計算で得られた値 (VASP, PAW,
PBE)。括弧外が文献値。z(M)などは、M 原子の座標のうち可変パラメータであるもの。
軸長の単位は Å。
Al (FCC)a = 4.0498 (4.0462)
Ca (FCC)a = 5.5884 (5.5194 )Mg (HCP)a = 3.2094 (3.1869) c = 5.2103 (5.1978)
Na (BCC)a = 4.235 (4.2044)
Si
a = 5.41985 (5.4663)
GaAs
a = 5.65359 (5.7605)
GaN (wurzite) a = 3.186 (3.2454) c = 5.176 (5.2897) z(N) = 0.375 (0.3758)
NaCl
a = 5.62 (5.6506)
MgO
a = 4.2109 (4.2362)
CaO
a = 4.8112 (4.8378)
ZnO
a = 3.2427 (3.2545) c = 5.1948 (5.2141) z(O) = 0.3826 (0.3816)
In2O3
a = 10.117 (10.032)
SnO2
a = 4.738 (4.7154) c = 3.1865 (3.1836)
b-Ga2O3
a = 12.23 (12.026)
b = 3.04 (2.9927)
c = 5.8 (5.7185) b =103.7 (103.86)
InGaO3(ZnO)1 a = 3.299 (3.2949) b = 5.714 (5.7042) c = 26.101 (25.4037)
12CaO·7Al2O3 a = 11.989 (12.0284, 11.997, 11.9884)(C12A7)
a = 90 (a=89.9895, b=89.9334, g=89.9619)
17. 未知構造探索
第一原理計算では、実際には存在しない仮想構造の計算ができることも大きな長所です。そのため、安定な未知構造を探索することで新材料探索が行われています。しなどが使われます。
バンド計算では第一 BZ 内の異なる k ベクトルについて波動関数を計算し、全電子密度などを求めます。この際に使う k 点の数を「k 点数」と呼ぶことがあり、精度に重大な影響を与えます。
90
実験の解析やモデルの方が間違えていることを見つけることもあります。DFT の結果と実験に大きな違いが出た場合、その理由を突き詰めることで新しい知見が得られるかもしれません。
89

かしながら、第一原理計算による安定構造の計算は、逐次計算により安定構造を求めるため、初期構造を与える必要があること、初期構造近傍で連続変化で到達する構造しか得られないという制限があります。
そのため、化学組成だけを与えて、いくつかの初期構造から不連続に構造を変化させて最安定構造を調べる方法がとられています。遺伝的アルゴリズムを用いる
USPEX 91や、粒子群最適化を用いる CARYPSO 92。などがあります。

91
92

https://uspex-team.org/en/uspex/overview
http://calypso.cn/