付録
講義で触れない内容で、知っておくと役に立つことを記載しています。
(文献・脚注番号はリセットしてありますので、ご注意ください)
付録 1a: 分子軌道法とバンド計算
1) 分子軌道法
分子の電子状態を、原子軌道の一次結合 (Linear Combination of Atomic Orbitals:
LCAO) とし、量子力学に則って計算する方法です。構造モデルは有限の大きさです。電子構造は、縦軸にエネルギーをとって占有状態、非占有状態、場合によってはスピン状態を表示した「エネルギー準位図」で表わされます。横軸にエネルギ
ー、縦軸にそのエネルギーをもつ電子の数をプロットした「状態密度」も使われま
すが、分子軌道法の場合には、両者はほぼ同一です。代表的な計算プログラムに、
半経験的分子軌道法の MOPAC 1、第一原理法の GAMMESS 2、Gaussian 3 などがあり、これらは Windows の GUI で操作できる Winmostar 4 を使って計算できます。
2) バンド計算
周期的な構造をもつ結晶を扱います。Bloch の定理を利用して計算対象を単位格
子内の原子だけに減らしていますが、その代償として、バンド計算に固有のパラ
メータ「Bloch の波数ベクトル k」が現れます。エネルギー固有値は k の関数 E(k)になり、電子構造は、横軸を k、縦軸を E(k)とする「バンド構造」で表されます。
バンド構造を使うことにより電子構造を直観的に理解しやすくなりますが、多くの学生からはわけがわからないと嫌われ、「スパゲッティ」と呼ばれることもあります。バンド構造は対称性の高い特定の方向だけでしか表示できず、結晶の全体の電子構造を理解する場合には見落としが出る可能性があります。一方で状態密度はすべての k に関する積分 5 となるため、全電子状態の情報を判断することができます。
付録 1b: 第一原理法、半経験法、経験法、準経験法
1) 第一原理法 (first-principles method) 6 (非経験法)
量子論の基礎方程式をなるべく忠実にコンピュータを使って解く方法。基本的に入力は原子の種類と座標の「構造パラメータ」のみです。ただし、第一原理とは言っても、構造パラメータ以外に、計算のモデルの種類、精度を決めるパラメータ群を適切に入力する必要があります。これらのパラメータを多少変えたとしても計算結果に与えることは無く 7、本質的な結果に影響を与えるのは、構造パラメー
タだけです。

http://openmopac.net/index.html
http://www.msg.ameslab.gov/GAMESS/
3 http://gaussian.com/
4 http://winmostar.com/
5 正確には「第一ブリルアンゾーン内の積分」
6 ”最初から”という意味のラテン語から派生した”ab initio 法”とも呼ばれます。明確に区別されて使われてはいませんが、first-principles はバンド計算、ab inito は分子軌道法で使われることが多いようです。
7 というよりも、数値を変えても計算結果が変わらない程度に高い精度の計算をしていることが、正しい結果を得るための必要条件です。
122) 経験的分子軌道法、半経験法 (semi-empirical method)
物理方程式は Schrödinger 方程式を用いるものの、計算に時間がかかる項－エネ
ルギー積分や共鳴積分－を経験パラメータに置き換えることで計算時間を短くする方法。構造パラメータの他にこれらの実験などによって決めた経験パラメータを入力する必要があります。計算結果は経験パラメータによって大きく変わる（逆に言えば、都合の良い結果を出せるように変えることができる）ため、これらの選択が計算の信頼性を決定します。実験結果に合うようにパラメータを再設定しなおす必要がでることも多く、この時間の方が第一原理計算よりも長くかかることも珍しくありません。
分子軌道法で使われる (拡張) Hückel 法とバンド計算への拡張である Tight
Binding (TB) 法 8 が半経験法の例です。CNDO/MNDO 法 (MOPAC) や高度な TB
法では HF 近似などを取り込み、分子軌道法の係数を self-consistent に決定しますが、これらは経験的分子軌道法と呼ばれます。
3) 経験法 (empirical method)
量子論的な効果を明示的に考慮せず、原子間の相互作用ポテンシャルの形を実験から決め、古典力学の運動方程式 (Newton の運動方程式) を用いて材料中の原子の運動をシミュレーションする方法です。
格子力学 (Lattice dynamics: LD) 法、分子動力学 (Molecular dynamics: MD) 法、
モンテカルロ (Monte Carlo: MC) 法などの多くがこれを採用しています。これらには量子計算法 (Quantum MD, QMD, First-Principles MD, FP-MD などと呼ばれる)もあり、区別する場合には「古典 MD 法」などと呼ばれます。
4) 準経験法
半経験法で計算の軽減のために導入したパラメータを第一原理的に計算する方法が、準経験法と呼ばれることがあります。第一原理計算で構造を変えた多数のモデルの全エネルギーから原子間ポテンシャルを抽出して古典 MD を行う方法(例えば General Utility Lattice Program, GULP 9 にこの機能があります)、やはり第一原理計算から TB 法に必要な移送積分 (transfer integral) を計算する方法があります (DFTB 法 10、eHtuner 11 など)。最近では、量子 MD を行いながら経験ポテ
ンシャルを更新する On-the-fly potential や機械学習による machine-learning ポテン
シャル 12が研究されています。
付録 1c: 原子間ポテンシャル: MD 法との組み合わせ
準経験法の例として、第一原理計算で原子配置の違う複数の構造で全エネルギーを計算し、MD 法などに必要な原子間ポテンシャルを決める方法を図 A1c-1 に示します。
これは、MgO についてイオン間ポテンシャルを求めたものです (CRYSTAL88, HF)。

8
9

W.A. ハリソン、固体の電子構造と物性、現代工学社 (1984).

http://gulp.curtin.edu.au/gulp/

https://www.dftbplus.org/; http://www.dftb.org/
http://www.chem.wisc.edu/~danny/software/ehtuner/
12 https://compphysvienna.github.io/n2p2/
1011

-4.7835

Total energy / 10-15J/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.788

0

0.2

0.4

0.6

0.8

1.0

Lattice parameter / nm

-4.7840-4.7845-4.7850-4.7855
a=0.4nm-4.7860-4.7865

a=0.39nm-4.7870-4.7875

a=0.44nm
a=0.421nm

00.02 0.04 0.06 0.08 0.10 0.12 0.14

Ionic displacement / fractional coordinate

図 A1c-1 MgO について第一原理計算で求めたエネルギー曲線 (シンボル) と、それから求めたポテンシャル関数から求めたエネルギー曲線 (実線) 13
付録 2a: 計算法の分類、計算ソフト
第一原理法は量子方程式を数値的に解く「解析方法」と波動関数を近似する「基底関数」によって区別すると、特徴がわかりやすくなります。WIEN2k 14 は L/APW(Linearized/Augmented Plain Wave)+lo (local orbital) 法を使い、ある大きさの原子球(Muffin-Tin 球) の内側では数値的に解いた原子の波動関数を、外側では平面波 (Plain
Wave: PW) を基底関数として使います。VASP 15 , CASTEP 16 , Quantum ESPRESSO(PWScf) 17 , PHASE 18 , Advance/PHASE 19 , ABINIT 20 は、擬ポテンシャル (PseudoPotential: PP) 法と PW 基底を使うので、PP-PW 法とも略されます。以前は L/APW 法の計算精度が最も高いといわれましたが、最近では PP 法と同じプログラムコードで
実 装 が 可 能 で APW 法 の 長 所 を 擬 ポ テ ン シ ャ ル に 取 り 込 ん だ PAW (Projector
Augmented Wave) 法 21 が VASP などでも使われるようになり、L/APW 法と同等の計算精度が出ます。ATK/DFT 22 も PP 法を使いますが、基底関数は、数値的に解いた原子の波動関数を使います。CRYSTAL 23 はもともと Gauss 関数型基底 (GTO) を用いる全電子 HF 計算でしたが、現在では PP (Effective Core Potential, ECP, と呼んでいる)や DFT/HDFT 計算も可能です。
いずれも、安定構造や物性を計算するには十分高い精度を持ちますが、どのような物性が計算できるかはプログラム毎に違うので、目的に応じて使い分ける必要があります。
付録 2b: 汎関数
本文では DFT で使われているいくつかの密度汎関数 (交換相関相互作用汎関数)が出てきました。他にも、非常に多くの汎関数が提案されています。これらを区別するため、提案者の頭文字 (と発表年) を使った短縮名が使われます。たとえば PBE (あるいは PBE96) は PRL, 1996, p.3865 で Perdew, Bruke, Ernzerhof が提案したものです。
T. Kamiya, Jpn. J. Appl. Phys. 35 (1996) 3688
http://www.wien2k.at/
15 http://www.vasp.at/
16 Materials Studio, https://www.comtec.daikin.co.jp/SC/prd/castep/
17 http://www.pwscf.org/
18 https://azuma.nims.go.jp/software
13
14

19
20

http://www.advancesoft.jp/product/advance_phase/
http://www.abinit.org/

P. E. Blöchl, Phys. Rev. B 50, 17953 (1994); G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999).
https://www.quantumwise.com/products/atk
23 http://www.crystal.unito.it/
21
22

彼らは同じ年に混合汎関数を提案しており、こちらは PBE0 と呼ばれています。
以下に、よく使われる汎関数をまとめておきます。汎関数の効果については、例えば文献 24 などで比較されています。
LDA/LSDA:
CA (Ceperley-Alder)/PZ (Perdew-Zunger) [Perdew and Zunger, Phys. Rev. B 23 (1981)
5048]
PW92 (Perdew-Wang 92) [J.P. Perdew and Y. Wang, Phys. Rev. B 45 (1992) 13244]
GGA:
PBE あるいは PBE96 (Perdew-Burke-Ernzerhof) [J.P. Perdew, K. Burke, M. Ernzerhof,
Phys. Rev. Lett. 77 (1996) 3865; ibid., 78 (1997) 1396]
revPBE (Revised PBE) [Y. Zhang and W. Yang, Phys. Rev. Lett. 80 (1998) 890]
RPBE (Revised PBE) [B. Hammer, L. B. Hansen, and J. K. Nørskov, Phys. Rev. B 59 (1999)
7413]
PBEsol (PBE for solids) [J.P. Perdew, A. Ruzsinszky, G.I. Csonka, O.A. Vydrov, G.E.
Scuseria, L.A. Constantin, X. Zhou and K. Burke, Phys. Rev. Lett. 100 (2008) 136406]
WC (Wu-Cohen modification of PBE) [Z. Wu and R.E. Cohen, Phys. Rev. B 73 (2006)
235116]
付録 2c: どの汎関数が正確か: 最安定構造に関連して14 章で説明した構造緩和計算ですが、多形 25 や相転移が近傍にあるような微妙な構造では注意が必要です 26。たとえば鉄の場合には、LSDA 27 で計算すると面心立方構造、非磁性が最安定になるとの結論が出ますが、実際には体心立方構造で強磁性が安定です。これは GGA を使うことで実験と整合する結果が得られます 28。また、ルチル構造の TiO2 においては、PBE 汎関数を使うと安定格子体積が大きくなり、強誘電相が安定になってしまうという問題があります 29。鉄の場合には GGA、ルチルの場合は LDA で正しい結果が得られたわけです。必ずしも L(S)DA よりも GGA の方が正しいとは限りませんので、注意が必要です。
付録 2d: PDOS の曖昧さの問題
分子軌道法 (付録 1a) と比較すると、バンド計算で得られる PDOS の解釈にはより慎重になる必要があります。分子軌道法のように LCAO を使っている場合、分子軌道
F. Tran, J. Stelzl, and P. Blaha, J. Chem. Phys., 144 (2016) 204120.
polymorph。同じ化学組成で結晶構造が異なるもの (ルチル構造とアナターゼ構造の TiO2 など)。積層順序だけが違うものを特に「多型 (polytype)」という (閃亜鉛鉱構造とウルツ鉱構造の GaN など)。SiC や ZnSe のように、積層構造が異なる多型が多種存在する場合は、その積層
周期と対称性を使って、2H SiC, 3C SiC, 4H SiC などと書きます。FCC 構造では二次元六方格子の最密充填層が 3 層周期で配列することで 4 回軸が発生して立方対称となっているため、3C
SiC になります。偶数層周期では 4 回軸が現れないので、graphite のように六方最密構造となり、2H SiC や 4H SiC となります。
26 計算精度にも注意が必要です。平面波基底のカットオフエネルギーや k 点数を増やし、判断しようとしているエネルギー差よりもこれらの効果が小さくなることを確認する必要があります。
27 Local Spin Density Approximation。DFT のもっとも簡単な形式の一つである局所密度近似(LDA)をスピン分極計算に適用し、スピン密度の局所密度汎関数を VXC として計算する方法。
28 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.
29 B. Montanari and N.M. Harrison, Chem. Phys. Lett. 364 (2002) 528.
24
25

を各原子へ割り振る寄与率は LCAO の係数 |Ci|2 (密度行列の対角成分に対応) で決まります。原子間の結合に寄与している成分が非対角項 Ci*Cj+Cj*Ci に残っており、これをどのように各原子に分配するかに若干の任意性が残りますが、平面波基底の場合ほ
ど大きな問題ではありません 30。
ところが、バンド計算の基底関数としてよく使われている平面波は空間全体に一様に広がっていますから、原子ごとの PDOS を一意的に計算する方法は定義できません。
そのため、L/APW 法である WIEN2k は Muffin-Tin (MT) 球内の波動関数の積分値を原子の PDOS としています。VASP などで使う PW 法や PAW 法では MT 球を使いませんので、各原子に半径を与え (Wigner-Seitz (WS) 半径と呼んでいます)、その中の波動関数の積分値を原子の PDOS として使っています。Total DOS はこれらのパラメータによって変わることはありませんが、PDOS は MT/WS 半径の選び方によって大きく変わります。意味のある PDOS を得るためにはコツをつかむ必要がありますし、逆に、化学者としての直観に合うように PDOS を意図的に変えることもできますので、注意が必要です 31。結合状態により、原子半径やイオン半径に比例させるということも考えられますが、これらの半径は被占有軌道の半径を反映しているため、非占有軌道の
PDOS に不自然な結果が出ることがあります。別の考え方は、次の付録 2e で紹介します。
また、MT/WS 球を設定して PDOS を計算すると、MT/WS 球に属さない波動関数成分がでてきます。これは、現実の固体の波動関数が原子の波動関数よりも外部に拡がりをもつことを反映しており、特に、金属や半導体の伝導帯では無視できなくなります 32。このような、原子に帰属させにくい状態密度を “Interstitial (格子間)” 状態として表示することがあります。
付録 2e: イオン電荷の曖昧さの問題
上記のように、どの原子にどの空間領域の電子を帰属させるかを決めると、その中の電子数の積分値と原子核電荷の差から、そのイオンの価数を計算することができます。しかしながら、このようにして求めたイオン価数は、共有結合の寄与のため、化学的な直観によるイオンの形式電荷から比べると非常に小さくなります。表 A2e-1 では後述の 2 通りの方法でイオン電荷を計算していますが、例えば、イオン性が強いア
ルカリハライドでは形式電荷の 90%程度という結果が出ますが、酸化物ではイオン性が強いと考えられている MgO でも 80%程度に過ぎません。ZnO では 60%、共有結合性が大きくなる ZnSe だと、30%前後まで電荷が小さくなります。

広く用いられている Mulliken 電荷解析では、非対角成分 CiCj を i 原子と j 原子に同量割り当てるという単純なやり方をしています。この場合でも、空間的広がりが大きく異なる基底関数を使うと、結果も変わってきます。
31 L/APW 法の場合、MT 球半径がとんでもない値だと、変な計算結果がでることがありますので、MT 球半径として選択できる値には制限があり、PDOS の任意性にもある程度の制限がかかります。PW-PP 法の場合には、電子構造の計算結果とは全く関係なく WS 半径を決めることができますので、後者の方が “自由度 (あるいは いいかげんさ)” は大きくなります。
32 このため、LCAO 法や TB 法では金属や半導体のバンド構造は正確に計算できないといわれることがあります。しかしながら、原子軌道を基底関数に使う場合も、主量子数の大きい軌道まで加えたり、拡がりの大きい diffuse 関数などを加えたりすることでこの問題は解決できます。現在の第一原理計算や TB 法ではこの指摘は当たらない場合が多いので注意が必要です。
30

表 A2e-1 MX 化合物のイオン電荷の計算値。WS 電荷は陽イオン M, 陰イオン X に異なる電荷値が計算されるが、Bader 電荷の場合は電荷総量が必ずゼロになるため、両者の電荷は同じとなる。
NaF
NaCl
LiF
LiCl

WS 電荷 M

WS 電荷 X

Bader 電荷

0.87
0.87
0.90
0.91

-0.55
-0.42
-0.53
-0.37

0.82
0.82
0.83
0.88

MgO
MgSe
ZnO
ZnSe

WS 電荷 M

WS 電荷 X

Bader 電荷

1.75
1.69
1.02
0.54

-1.12
-0.87
-1.11
-0.72

1.64
1.58
1.20
0.73

また、付録 2d で説明した PDOS の曖昧さのため、どの方法で PDOS を計算するかにより、イオン電荷も大きく変わります。原子間で電子密度の微分がゼロになる面
2(電流ゼロに対応する。zero flux surface, ∇ Ψ = 0 ) で電荷を分割 33 して曖昧さをなくした Bader charge analysis 34,35 が提案されていますが、この場合でも表 2e-1 のように形式イオン電荷よりかなり小さい値を与えます 36。
付録 2d で説明したように、PW 基底を使う VASP などの場合には、WS 半径を決めて有効電荷を計算することができます。ただし WS 半径の決め方には大きな任意性があります。この任意性をなくすための一つの方法として、私たちは、Bader 電荷解析
で求めた各原子の体積 (Bader volume) に相当するように WS 半径を決める方法を使うことがあります。図 2e-2 はこのようにして計算した WS 電荷と Bader 電荷を、周期表内の MX 型化合物 (仮想的な結晶構造を含む)で計算した結果です (VASP では WS
球の体積の和が単位格子の 90%程度になることがマニュアルで推奨されていますので、そのように WS 半径を決めています)。両者には±0.5 程度のばらつきがあります
が、良い相関が認められますので、任意性を排除した有効電荷の計算方法としては一定の利用価値があると考えられます。

WS charge / |e|WS charge / |e|Bader charge / |e|(B)

Bader charge / |e|(A)

図 2e-2 MX 型化合物の Bader 電荷と、Bader 体積の 90%から WS 球半径として計算した (A) M イオンと (B) X イオンの電荷 (WS charge) の相関関係。
その他には、イオン変位による誘電率から逆算される “ボルン有効電荷 (Born
Zero flux 面は、Lagrangian がゼロになる面 (Lagrangian の基準は無限遠点でゼロにとる)と等価であり、そのために zero flux 面で囲まれた体積は物理的に独立に扱えるという原理を基盤としています。
34 R. Bader, Atoms in Molecules: A Quantum Theory, Oxford University Press, New York, 1990.
35 http://theory.cm.utexas.edu/henkelman/code/bader/
36 Z. Xiao, H. Lei, X. Zhang, Y. Zhou, H. Hosono and T. Kamiya, Bull. Chem. Soc. Jpn. 88 (2015) 1250.
33

effective charge)” もよく使われます。誘電率に寄与する電荷をすべてイオンの点電荷qi に帰すると、電場 E 下においてイオンが δri だけ変位した場合に誘起される分極密度∆P は(A2e-1)𝚫𝚫𝑷𝑷 = 𝑉𝑉 −1 ∑ 𝑞𝑞𝑖𝑖 𝛿𝛿𝒓𝒓𝑖𝑖 = (𝜀𝜀 − 𝜀𝜀0 )Eとなります (V は和を取っている体積)。付録 2f で説明するように、δri はポテンシャ
ル曲率を用いて E に比例しますので、式(A2e-1) から誘電率 ε を計算できます。第一
原理計算ではイオンの変位 𝛿𝛿𝒓𝒓𝑖𝑖 に対応させて分極密度変化を計算しますが、このことは、式(A2e-1) から、イオン電荷 qi を逆算することができることを意味しています。
これが Born 有効電荷で、一般の場合には下記のようなイオンごとに与えられるテン
ソルであらわされます (𝛿𝛿𝛿𝛿𝑖𝑖,𝑘𝑘′は i 番目のイオンの変位の k’方位成分、𝑃𝑃𝑘𝑘 は電気分極のk 方位成分)。
𝑍𝑍 ∗ 𝑖𝑖,𝑘𝑘𝑘𝑘′ = 𝑉𝑉𝑉𝑉𝑃𝑃𝑘𝑘 /𝜕𝜕𝛿𝛿𝛿𝛿𝑖𝑖,𝑘𝑘′(A2e-2)電磁気学の教科書では、電気分極は∫ 𝒓𝒓𝜌𝜌(𝒓𝒓)𝑑𝑑𝒓𝒓で計算できると書いてありますが、これは電荷が有限の範囲に分布している場合でないと、一意的に計算できません。結晶のように無限周期構造があり、また、電子のように連続分布関数として電荷があらわされる場合は上記の式が使えません。この問題に対しては、原子を変位させたときの分極変化を Berry 位相で計算できるということがわかり、電気分極、誘電物性の計算にはこの方法が使われています (付録 2f を参照)。
ただし Born 有効電荷の場合は、誘電率にはイオン電荷の変位だけではなく、イオンの変位に伴う電子の移動により発生する電気分極変化の寄与が有効電荷に加わります。この効果は特に、誘電率の大きい物質で顕著です。そのため、誘電率が小さい
物質ではそれらしい値が得られます (表 A2e-2) が、TiO2 のような高誘電率材料、BaTiO3 などの強誘電体などでは、ボルン有効電荷はイオンの形式電荷よりかなり大きくなります (表 A2e-3 37)。電子移動の寄与が大きい極端な場合には、TTF-CA の強誘電分極で見られるように、電気分極から予期されるのとは逆方向にイオンが変位することもあります 38。
表 A2e-2
MgO
CaO
SrO
表 A2e-3
形式電荷BaTiO3
SrTiO3CaTiO3
PbTiO3
PbZrO3BaZrO337
38

MX 型化合物の Born 有効電荷 (VASP, PBE)。
Z*(M)
Z*(M)
1.98
NaCl
1.09
2.36
GaN
2.7
2.44
ZnO
2.2
立方晶ペロブスカイト型化合物の Born 有効電荷 36。
Z*(A)
Z*(B)
Z*1(O)
+2
+4
-2
2.75
7.16
-5.69
2.54
7.12
-5.66
2.58
7.08
-5.65
3.90
7.06
-5.83
3.92
5.85
-4.81
2.73
6.03
-4.74

Z*2(O)
-2
-2.11
-2.00
-2.00
-2.56
-2.48
-2.01

W. Zhong, R.D. King-Smith, and D. Vanderbilt, Phys. Rev. Lett. 72 (1994) 3618.
K. Kobayashi, S. Horiuchi, R. Kumai, F. Kagawa, Y. Murakami, and Y. Tokura, Phys. Rev. Lett. 108(2012) 237601.

イオンを隣の単位格子の等価位置に移動した際の電気分極を Berry phase を用いて計算してイオン電荷を求める方法も提案されています 39。断熱状態で (固有状態を変えずに、の意味。ここではバンドギャップを開いたまま、バンド間の電荷移動なく、と理解してください) イオンを移動させる経路が存在する場合には一意で整数のイ
オン価数が得られることが保証されるので、形式電荷と対応させやすい電荷値が得られます。一方、このことは、この方法が金属や非整数価数の系には使えないことを意味しています。
付録 2f: 弾性率、誘電率、圧電定数の計算
量子計算の教科書やプログラムのマニュアルでは、これら物性の計算を、波動関数を使って量子力学的に説明しているものが多いですが、見かけ上数式が複雑になり、理解しにくくなります。これらの物性は全エネルギー (正確には自由エネルギー。付
録 2k 参照) の微分から比較的簡単に説明できますので、ここではその流れで説明します。ここではバンド計算に対応できるよう、物質の全エネルギーを格子定数 ai と原子の内部座標 rl で 𝑈𝑈0 (𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 ) と表します。ここで応力 σij を加えた結果、格子定数が ai + δai に歪んだとすると (簡単のため、原子は内部的には変位しないとします) と、
この時の全エネルギーは1
𝜕𝜕𝜕𝜕
𝑈𝑈(𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 + 𝜹𝜹𝒂𝒂𝒊𝒊 ) = 𝑈𝑈0 (𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 ) + ∑𝑖𝑖,𝑗𝑗=𝑥𝑥,𝑦𝑦,𝑧𝑧𝛿𝛿𝑎𝑎𝑖𝑖 𝛿𝛿𝑎𝑎𝑗𝑗 + 𝑂𝑂(3) (𝛿𝛿𝑎𝑎𝑗𝑗 )(A2f-1)2

𝜕𝜕𝑎𝑎𝑖𝑖 𝜕𝜕𝑎𝑎𝑗𝑗

となります。ここで、応力のないときの平衡条件から、右辺では U の一次微分の項はゼロとしています。これと弾性率テンソルの定義
1𝛿𝛿𝑎𝑎 𝛿𝛿𝑎𝑎

1∆𝑈𝑈 = 2 ∑𝑖𝑖,𝑗𝑗=𝑥𝑥,𝑦𝑦,𝑧𝑧 𝐶𝐶𝑖𝑖𝑖𝑖 𝑠𝑠𝑖𝑖𝑖𝑖 = 2 ∑𝑖𝑖,𝑗𝑗=𝑥𝑥,𝑦𝑦,𝑧𝑧 𝐶𝐶𝑖𝑖𝑖𝑖 𝑎𝑎 𝑖𝑖 𝑎𝑎 𝑗𝑗(sij は歪テンソル) から、
𝜕𝜕𝜕𝜕

𝑖𝑖

(A2f-2)𝑗𝑗(A2f-3)𝐶𝐶𝑖𝑖𝑖𝑖 = 𝑎𝑎𝑖𝑖 𝑎𝑎𝑗𝑗 𝜕𝜕𝑎𝑎 𝜕𝜕𝑎𝑎𝑖𝑖𝑗𝑗が得られます。式(A2f-3) は、式(A2f-1) がなくても式(A2f-2) から直接導出できます
が、次の誘電率の説明のため、式(A2f-1) から説明しました。
次に、外部電界 E を引加した結果、原子が rl + δ rl に変位したとする (簡単のため、格子歪はないとします) と、原子の電荷を ql として、この時の全エネルギーは1
𝜕𝜕𝜕𝜕
𝑈𝑈(𝒓𝒓𝒍𝒍 + 𝜹𝜹𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 ) = 𝑈𝑈0 (𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 ) + ∑𝑙𝑙,𝑙𝑙′ ,𝑖𝑖,𝑖𝑖′
𝛿𝛿𝑥𝑥𝑙𝑙,𝑖𝑖 𝛿𝛿𝑥𝑥𝑙𝑙′ ,𝑖𝑖′ + ∑𝑙𝑙 𝑞𝑞𝑙𝑙 𝛿𝛿𝑥𝑥𝑙𝑙 𝑬𝑬 + 𝑂𝑂(3) (𝛿𝛿𝑥𝑥𝑙𝑙,𝑖𝑖 )2𝜕𝜕𝑥𝑥𝑙𝑙,𝑖𝑖 𝜕𝜕𝑥𝑥𝑙𝑙′ ,𝑖𝑖′

(A2f-4)となります。ここで平衡条件から、右辺では U の一次微分の項はゼロとしています。
式(A2f-1) を l 番目原子の内部座標の i 方位成分 xl,i で微分すると
𝑑𝑑𝑑𝑑
𝜕𝜕𝜕𝜕
= ∑𝑙𝑙′ ,𝑖𝑖′
𝛿𝛿𝑥𝑥𝑙𝑙′ + 𝑞𝑞𝑙𝑙 𝐸𝐸 = 0(A2f-5)
𝑑𝑑𝑥𝑥𝑙𝑙,𝑖𝑖

𝜕𝜕𝑥𝑥𝑙𝑙′ ,𝑖𝑖′となり、この連立方程式から、電界 E 下での原子の変位 δ rl が決まります。電気分
極密度 𝑃𝑃𝑖𝑖 = 𝑉𝑉 −1 ∑𝑙𝑙 𝑞𝑞𝑙𝑙 𝛿𝛿𝑥𝑥𝑙𝑙,𝑖𝑖 (V は単位格子体積) と誘電率の定義
𝑫𝑫 = ε𝑬𝑬 = 𝜀𝜀0 𝑬𝑬 + 𝑷𝑷
(A2f-6)から、誘電率 ε が計算できます。圧電効果が無視できない場合は格子歪と原子の変位の両方を考慮し、𝑈𝑈(𝒓𝒓𝒍𝒍 + 𝜹𝜹𝒓𝒓𝒍𝒍 , +𝜹𝜹𝜹𝜹𝒊𝒊 ) − 𝑈𝑈0 (𝒓𝒓𝒍𝒍 , 𝒂𝒂𝒊𝒊 )
1
𝜕𝜕𝜕𝜕
1
𝜕𝜕𝜕𝜕
𝜕𝜕𝜕𝜕
= 2 ∑𝑙𝑙,𝑙𝑙′ ,𝑖𝑖,𝑖𝑖′ 𝜕𝜕𝑥𝑥 ,𝑖𝑖𝑖𝑖𝑥𝑥 𝛿𝛿𝑥𝑥𝑙𝑙,𝑖𝑖 𝛿𝛿𝑥𝑥𝑙𝑙′ ,𝑖𝑖′ + 2 ∑𝑖𝑖,𝑖𝑖′ 𝜕𝜕𝑎𝑎 𝜕𝜕𝑎𝑎 𝛿𝛿𝑎𝑎𝑖𝑖 𝛿𝛿𝑎𝑎𝑗𝑗 + ∑𝑖𝑖,𝑙𝑙 𝜕𝜕𝑎𝑎 𝜕𝜕𝑥𝑥 𝛿𝛿𝑎𝑎𝑖𝑖 𝛿𝛿𝑥𝑥𝑙𝑙 + ∑𝑙𝑙 𝑞𝑞𝑙𝑙 𝛿𝛿𝑥𝑥𝑙𝑙 𝑬𝑬
𝑙𝑙39

𝑙𝑙′ ,𝑖𝑖′

𝑖𝑖

𝑗𝑗

𝑖𝑖

𝑙𝑙

L. Jiang, S.V. Levchenko and A.M. Rappe, Phys. Rev. Lett. 108 (2012) 166403.

(A2f-7)

から、格子歪がない場合の誘電率 (clamped dielectric constant) と歪がある場合の誘電率が計算できます。逆圧電テンソル dtijk は定義より
𝑠𝑠𝑖𝑖𝑖𝑖 = 𝑑𝑑𝑡𝑡 𝑘𝑘𝑘𝑘𝑘𝑘 𝐸𝐸𝑘𝑘
(A2f-8)40ですので、電界印加下の格子歪を 式(A2f-7) から計算して求めることができます。
正圧電テンソル dijk は電束密度 D により𝐷𝐷𝑘𝑘 = ∑𝑖𝑖,𝑗𝑗 𝑑𝑑𝑘𝑘𝑘𝑘𝑘𝑘 𝜎𝜎𝑖𝑖𝑖𝑖 + ∑𝑘𝑘 𝜀𝜀𝑖𝑖𝑖𝑖 𝐸𝐸𝑖𝑖
(A2f-8)で定義されますが、この場合は Berry 位相で計算した電気分極密度 P を用いて計算をします。
以上の式をこれ以上具体的に書き下すことはしませんが、付録 2k のように U を自由エネルギーにとれば有限温度でも使えますし、バンド計算の場合は全エネルギーで置き換えることで絶対零度での弾性率、誘電率、圧電定数が計算できます。
𝜕𝜕𝜕𝜕式(A2f-4) で出てくる 𝜕𝜕𝑥𝑥 𝜕𝜕𝑥𝑥 は、Hessian (ヘッセ) 行列と呼ばれており、弾性率の𝑙𝑙𝑙𝑙′

計算だけでなく、構造緩和計算でもよく出てきます。Hessian 行列の計算をするには、摂動理論 (Density Functional Perturbation Theory: DFPT) を使って直接計算する方法と、原子を微少量変位させながら全エネルギーや力 (Hermann-Feynman Force) を計算して数値微分により求める方法 (Finite Displacement 法、Frozen Phonon 法などと呼ばれます) があります。
実測値の電気分極を対応させるには、試料形状による反分極場の考慮が必要ですが、他に局所場 (local field / Lorentz field) の補正が必要になります 41。VASP の場合は非補正値と補正値の両方が出力されます。
付録 2g: 電気分極の計算に関する注意

単位格子2

a+ r-

r-

単位格子1 r+

a+r+

s

r-–s

a+r-–sa+r+–s2a+r+–s

図 A2g-1 周期構造がある場合、分極の計算は一意的に行えない。単位格子 (図中四角形) 内で和をとると、単位格子の取り方によって、分極𝑷𝑷 = ∑𝑢𝑢𝑢𝑢𝑢𝑢𝑢𝑢 𝑐𝑐𝑐𝑐𝑐𝑐𝑐𝑐 𝑞𝑞𝑖𝑖 𝒓𝒓𝑖𝑖 の値が qa(a は任意の格子ベクトル) だけずれる可能性がある。
+q, r+

変位前-q, r-

変位後

r+

a

r+

+qは r+ から r+

に変位

r-

図 A2g-2 点電荷の変位∆𝒓𝒓+ による分極変化∆P。点電荷がは r+ から r++∆r+ に変位したと正しく追跡できれば正しい∆P が計算できるが、r++∆r+ − 𝒂𝒂 へ変位したとして計算すると、-qa だけ∆P を間違える。
40
41

池田拓郎著、圧電材料学の基礎、オーム社 (S59).
キッテル、固体物理学入門 第 8 版、丸善 (H17).

付録 2e の Born 有効電荷で説明したように、結晶の電気分極の絶対量は、一意的に定義できません。結晶のように周期構造がある場合は、δri に格子定数の整数倍の曖
昧さが発生します (図 A2g-1)。そのため、物理的に意味があるのは、原子核が変位した際の分極変化 ∆P になります。原子核の変位前後で各電荷がどれだけ変位したかを追跡できれば、∆P は一意的に決まります。逆に言えば、電荷を正確に追跡できないと、∆P を正しく計算できません (図 A2g-2)。しかしながら、電子のように連続的に電荷が分布している場合は、変位前後でどの点にあった電子がどこに動いたかを追跡することができません (追跡すること自体が物理的に無意味です)。この問題は、原子核が変位した際の分極変化 ∆P は Berry 位相を計算することで一意に決定できることがわかり、解決しました 42。現在の誘電物性関係の計算はすべて Berry 位相により計算されています。実際の計算では、付録 2f と同様、摂動理論や Frozen Phonon 法により Berry 位相の変化量を計算し、∆P を求めています。
付録 2h: 構造緩和計算の結果に関する注意表 13-1 の C12A7 については、文献にある結晶構造は立方晶ですが、緩和計算の結果は三斜晶になっています。これは、前者が X 線回折による構造解析で得られたもの
で、平均構造として立方晶となっているためです。実際の C12A7 では、12 個ある等価な原子位置のうち 2 つを O2-イオンが統計的に 43 占有しています。そのため、これらの O2- イオンの位置を決めて構造緩和計算をすると、対称性の低い正方晶や三斜晶になります。
以上のことは、原子が統計的にサイトを占めている固溶体 44では、緩和構造と実測の結晶構造の対称性を対応させるのが難しい場合があることを意味しています。一方で、固溶体や乱雑構造をもつ材料では、実験的には局所構造を知ることが困難ですが、それができるのは、構造緩和計算の大きな利点です。

ε1(0)付録 2i: 光学スペクトルの計算上の注意

図 A2i-1

15
14
13
12
11
10
9
8

HSE

PBE96
5

10

15
Nkx

20

25

Si の静的光学誘電率の計算値の k 点数依存性

11 章で光学スペクトルの計算について紹介しましたが、光学スペクトルの計算の際には、k 点をかなり多くとらないと計算結果が収束しないことに注意が必要です。図
R.D. King-Smith and D. Vanderbilt, Phys. Rev. B 47 (1993) 1651(R).
「ランダムに」を格好よく表現しているだけです。
44 半導体工学では「混晶」、金属工学では「不規則格子」と呼ばれる。
42
43

A2i-1 に、VASP で PBE と HSE で計算した静的光学誘電率 ε1(0) と一方位当たりの k
点数 Nkx (全 k 点数は Nkx3) をプロットしていますが、Nkx を 12 以上に取らないと十分に収束していないことがわかります (VASP では Scissors operator によるバンドギャッ
プ補正ができませんので、PBE では過小評価したバンドギャップで計算しており、
ε1(0)が実測値 11.9 よりもかなり小さくなっています)。VASP のマニュアルページでは、Si を例に精度を上げる手順が説明されています 45。
付録 2j: DFT で構造が正確に計算できない例
DFT ではバンドギャップが正確に求められないけれども、構造の精度は高いという話をしました。また、磁性体の構造は磁気構造 (スピン配置) にも大きく影響される
ため、正しい磁気構造で計算する必要があります (付録 8) し、電子相関が大きい系の計算結果の信頼性には注意する必要があることも説明してきました (10 章)。しかしながら、これらに該当しない "単純な” 系でも、DFT の計算結果が大きく外れる系があることが知られています。
一つの例が、希ガス原子や分子から構成されている結晶のように、構成原子間の相互作用が弱い場合です。分子間力の代表的なものが London 相互作用 (van der Walls 相互作用) ですが、この場合は、他の原子の影響により自身の電子分布が扁形して電気分極を形成する効果によるものであり、扁形しているという意味で励起状態の効果によって現れます。このため、基底状態を前提にしている DFT では分子間力を過小評価します。例えば、層状化合物で層間が分子間力で結合しているグラファイトの場合、格子定数の実験値は a = 0.2464, c = 0.6711 nm (ICSD#: 76767) ですが、PBE を使うと a= 0.2468, c = 0.7471 nm と、分子間力を弱く c 軸長を大きく計算してしまうことがわかります 46。この問題に対して、nonlocal correlation 項 Ecnl として分子間力項を取り入れる van der Waals density functional (vdW-DF) 法が提案されています 47,48。ただしこの場合、相関エネルギー、係数が異なる様々な vdW-DF が提案されており、対象とする系に対してどれが良い結果を与えるか、まだよく整理されていません。"第一原理計算的” に使うのはまだ難しいかもしれません。
また、BaTiO3 や PbTiO3 など強誘電体の場合も、構造精度が出ない場合があることが報告されています。これらは室温では正方晶となって強誘電性を示し、Curie 温度
Tc が格子軸率 c/a 比と良い相関があることが知られています。表 A2j-1 に PbTiO3 の計算結果を示していますが、PBE を使うと、c/a 比を非常に過大評価することがわかります。一方、HF 近似では、逆に c 軸長の方が a 軸長よりも短くなり、その結果、EE項の係数を適当に選ぶと、実測値に近い構造を再現できることがあります。BaTiO3 の場合には、実際の最安定構造は菱面体格子で、室温の強誘電体相とは違うので比較は難しいのですが、それでも c/a 比を過大評価することが確認できます。
付録 2b に挙げた Wu-Cohen 汎関数など、改善された汎関数も報告されていますが、割と身近で興味のある物質系に対しても、系統的に DFT の精度が悪い系が存在するということには注意が必要です。

https://cms.mpi.univie.ac.at/wiki/index.php/Dielectric_properties_of_Si_using_BSE
グラファイトの場合は、実は LDA で計算すると a = 0.2447, c = 6.70 nm と、実測値に近い結果が出ます。が、これは偶然の一致です。.
47 M. Dion, H. Rydberg, E. Schrö der, D.C. Langreth and B.I. Lundqvist, Phys. Rev. Lett. 92 (2004)
246401.
48 Review: D.C. Langreth et al., J. Phys. Condens. Matter 21 (2009) 084203
4546表 A2j-1 PbTiO3 の正方晶強誘電体層の安定構造の計算結果。表中上段から a 軸長、
c 軸長 (nm)、c/a 比。HSE 混合汎関数の場合は、EE 項の係数を変えた計算を行った。
HF
PBE0
HSE
HSE
HSE
実験値 PBE(CEE=0.25) (CEE=0.35) (CEE=0.50)
PbTiO3 0.3902 0.3875 0.3909 0.3855
0.38600.3856
0.3919a, c 0.4156 0.4510 0.3740 0.4353
0.43490.42890.3894c/a 比
1.065
1.164 0.9568
1.129
1.127
1.112
0.9936付録 2k: 自由エネルギー(c)

3R x 2

-1040

LT

S

HT

vib

LT

-1060

0

100 200
T (K)300(J mol-1 K-1)

100

F (kJ mol-1)

HT

-1020

0

400
200
T (K)

600

0

図 2k-1 YCu の高温相 (HT) と低温相 (LT) の(a)フォノン状態密度、(b)定積モル比
熱 CV、自由エネルギーF と格子エントロピー Svib 49。F の交差点が相転移温度。
Born-Oppenheimer の断熱近似を基盤とする第一原理計算では、全ハミルトニアンの固有値として内部エネルギーU が得られ、上記のように物性の計算も可能になります
が、原子核の運動を考慮しないため、基本的に絶対零度での計算になります。
一方、私たちが実験で得られる結果（相安定性、相転移、あるいは弾性率など自由エネルギーの微分で得られる物性など）は、通常は圧力一定、温度一定条件で測定されることから、Gibbs エネルギー G = U + PV – TS に基づいて決定されます。
このうち、系の体積 V は計算対象の単位格子体積として得られますので、大きな計算コストの増加なしに、自由エネルギーをエンタルピー H = U + PV として近似することができます。この場合は絶対零度、定圧条件に対応します。付録 5 のように、高圧下での相安定性・相転移などの計算が可能になります。
一方、体積 V 一定、有限温度での相安定性・相転移などを議論するには Helmholtz
エネルギーF の計算が必要になります。F は電子系 Felectrone と格子系 Sphonon の寄与分けられます。
𝐹𝐹(𝑉𝑉, 𝑇𝑇) = 𝐸𝐸0 (𝑉𝑉) + 𝐹𝐹𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 (𝑉𝑉, 𝑇𝑇) + 𝐹𝐹𝑝𝑝ℎ𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 (𝑉𝑉, 𝑇𝑇)
𝐹𝐹𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 は第一原理計算で得られた状態密度 D(e)と Fermi-Dirac 分布関数 f(e)から、𝐹𝐹𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 = 𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 − 𝑇𝑇𝑆𝑆𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒∞
𝐸𝐸𝐹𝐹 (0K)
𝐸𝐸𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 = ∫−∞ 𝐷𝐷(𝑒𝑒)𝑓𝑓(𝑒𝑒)𝑑𝑑𝑑𝑑 − ∫−∞𝑛𝑛(𝑒𝑒)𝑑𝑑𝑑𝑑∞

𝑆𝑆𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒𝑒 = −𝑘𝑘𝐵𝐵 ∫−∞ 𝐷𝐷(𝑒𝑒)[𝑓𝑓(𝑒𝑒)ln𝑓𝑓(𝑒𝑒) + (1 − 𝑓𝑓(𝑒𝑒)ln(1 − 𝑓𝑓(𝑒𝑒))]𝑑𝑑𝑑𝑑
で計算できますので、計算コストはほとんど増加しません。一方、𝐹𝐹𝑝𝑝ℎ𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 は第一 BZ
全体でフォノン分散ℏ𝜔𝜔𝑞𝑞 を計算し、以下の和を取る必要がありますので、非常に計算コストが高くなります。
49

H. Mizoguchi, J. Bang, T. Inoshita, T. Kamiya, and H. Hosono, Inorg. Chem. 58 (2019) 11819

1

𝐹𝐹𝑝𝑝ℎ𝑜𝑜𝑜𝑜𝑜𝑜𝑜𝑜 = 2 ∑ ℏ𝜔𝜔𝑞𝑞 + 𝑘𝑘𝐵𝐵 𝑇𝑇 ∑ 𝑙𝑙𝑙𝑙�1 − 𝑒𝑒 −ℏ𝜔𝜔𝑞𝑞/𝑘𝑘𝐵𝐵 𝑇𝑇 �

最近では図 2k-1 のように、有限温度相転移の計算の報告も増えています。
付録 3a: バンドギャップができるいろいろな成因
前のバンド理論の講義の「(ほぼ)自由な電子理論 (NFE)」では、BZ 境界では電子の進行波と反射波が干渉してバンドギャップを形成するということを学びました。しかしながら、私たちが普段問題にする「バンドギャップ」は、その下のエネルギー準位
が電子で完全に占有されていて、その上は電子が全く存在しないという意味の「バンドギャップ」のことです。後者の狭義のバンドギャップは NFE で出てくるバンドギ
ャップとは異なる成因で形成されています。
1) 共有結合（等極結合）Si など
結合、反結合軌道のエネルギー分裂によってバンドギャップが形成される。
2) イオン結合（異極結合）酸化物など 50
イオンのエネルギー準位の違いによってバンドギャップが形成される。
原子の 2p 軌道と 3s 軌道には大きなエネルギーギャップ (energy gap) があることと、本質的に同じ。
3) BZ 境界での干渉： Bragg 反射
NFE ででてくる広義のバンドギャップ。典型的な例は「パイエルス転移 51」。
4) 強電子相関系材料
上記の 1) ~ 3) によるバンドギャップは、DFT などの第一原理バンド計算でも再現できます 52。しかし、これらの計算では現実の材料で現れるバンドギャップを再現できないことがあります。典型的な原因に、電子相関を十分取り入れられていないことが挙げられ、特に d 電子、f 電子を含み、局在性の強い電子状態を作りやすい遷移金
属酸化物などで問題になります 53 。この場合、L(S)DA(GGA)+U などのいわゆる「beyond DFT」を使うことが必須になります (10 章)。
付録 3b: ブラベー格子 (Bravais lattice) と基本格子 (Primitive lattice)
GaAs の結晶構造は「閃亜鉛構造」と呼ばれます (図 A3b-1) が、Ga と As を同じ原子に置き換えると Si の結晶構造である「ダイヤモンド構造」になります。閃亜鉛構造
は 8 つの原子を含む立方体の単位格子で表されます（「立方晶に属する」といいます）
が、図 A3b-1(A) からわかるように、Ga だけでなる面心立方構造 (Face centered cubic:
FCC)と As だけでなる FCC 構造が (1/4, 1/4, 1/4) ずれて重なった構造です。しかし、
面心立方格子などは「ブラベー格子」といわれるもので、単位格子に含まれる原子の数を増やして見かけの対称性が高い形の単位格子をとったものです (図 A3b-1(B1))。
実際には、もっと小さい、原子数の少ない格子でも結晶の 3 次元周期性を表現するこ50

固体中のイオンのエネルギー準位を正しく理解するには、「マーデルングポテンシャル」を考慮する必要があります。
51 フェルミ面が BZ 境界ではなく BZ 内部にある場合、結晶の周期を大きくしてフェルミ面に
BZ 境界が来るように周期構造が変形することがあります。そうすると、新しい BZ 境界でバンドギャップが開き、その分、電子系のエネルギーが低くなります。このエネルギー利得が結晶の周期性を変えることによる弾性エネルギーの損より大きいと、この周期の構造が安定化します。代表的な例がポリアセチレン。
52 バンドギャップ問題があるので、実験値を定量的に再現できるという意味ではありません。
53 バンドギャップを作る原子軌道が同一原子の場合に「ハバードギャップ」、異なる隣接原子の場合に「電荷移動ギャップ」などと区別されます。

とができ、そのような最小の格子を「基本格子」と呼びます（図 A3b-1(C1)）。文献な
どの格子定数は、通常はブラベー格子で表され、Si の場合は 296 K において aC = 0.5431
nm です (C は立方晶、cubic の意味)。一方、基本格子である菱面体格子 (Rhombohedral
lattice) では aR = 0.3840 nm, αR = 60o になります (R は菱面体の意味)。
実際の量子計算では、扱う原子数を減らして計算時間を短くするため、ブラベー格子ではなく「基本格子」を使います。計算結果も基本格子に対して出力されることが多いので注意が必要です。また、(B3) と (C3) に見られるように、単位格子をどうとるかによってバンド構造の見かけは変わってきます。ただしこれは第一 BZ 中の特定の高対称位置だけを抜き出したために見かけが変わっているだけで、状態密度、全エネルギーなどを比較すれば、単位格子の取り方によらずに厳密に同じになっていることが確認できます。
(A)(B1) Si,Ga

Si,As

Si,Ga

(C1)

Si,As

Si,Ga

cR

cC bC
aC
CASTEP Band Structure

(B2)

(B3)

(C3)

Energy (eV)
9

Si,As

bR
aR

CASTEP Band Structure

Energy (eV)
30

(C2)

8
7
6

5

5

20

4
3
2

10

1

0

10

0

-1
-2

0

-3

0

-4
-5
-6
-7

-10

-10

-8
-9

-10

-10
-11
-12

X

X

R

R

M

M

G

Γ

R

R

-20
W

W

L
L

Γ
G

X W K
X

W

K

図 A3b-1 Si, GaAs の単位格子。(A) 2 つの面心心立方構造が重なってできています。
(B1) ブラベー格子である面心立方格子、(B2) それに対応する第一 BZ と(B3) バンド構造。(C1) 基本格子である菱面体格子、(C2) それに対応する第一 BZ と(C3) バンド
構造。ブラベー格子では 8 つの Si 原子を含みますが、基本格子では 2 つに減っています 54。基本格子の体積はブラベー格子の体積の 1/4。
付録 3c: 対称性の表現
波動関数の対称性は、結晶・分子の対称性によって制限され、後者の対称操作によって変化しないか、符号が反転します。そのため、波動関数の対称性 (対称群の既約表現) を簡単な記号で表記しますが、その方法には主に、既約表現の次数に対応して(A,B), E, T を頭につける Mulliken の記法 55 と、Γを頭につけ連番を右下に振る Bethe
単位格子の角にある原子は 4 つの単位格子、面上にある原子は 2 つの単位格子で共有されていることを考慮して、それぞれ 1/4 と 1/2 の重みをかけて原子の数を数えてみましょう。
55 A などと書かれる。Z 軸回りの回転 C に対して符号が反転しない (偶である) 場合が A, す1gnる (奇) 場合が B であり、二重縮退、三重縮退の場合はそれぞれ E, T と表す。y 軸回りの回転に対して偶の場合は 1, 奇の場合は 2 を右下に。反転 i に対して偶の場合は g、奇の場合は u を
さらに右下につける。z 軸に垂直な面σh の反転に対して偶の場合は ’, 奇の場合は ’’ を右上に54

の記法などがあります（表 A3c-1 に Oh 点群の例を示す）。詳しくは文献 56, 57を参考にしてください。
表 A3c-1
A1g
Eg
T1g
T2g
T1u

Oh 点群の指標表 (一部抜粋)
Γ,R,H E 8C2 3C2 6C2
+
1 1
1
1
Γ1 Γ1
+
2 -1
2
0
Γ3 Γ12
+
’
3 0
-1
-1
Γ4 Γ15+’3
0-1
1
Γ5 Γ253 0-1-1
Γ4 Γ15

6C4
1
0
1
-1
1

i
1
2
3
3
-3

8S6
1
-1
0
0
0

3σh
1
2
-1
-1
1

6S4
1
0
1
-1
-1

6σd
1
0
-1
1
1

x2+y2+z2 58(x2-y2,2z2-x2-y2)(Rx,Ry,Rz)(xy,xz,yz)(x,y,z)波動関数の既約表現がわかると、いくつか重要な情報が即座にわかります。移送積
分 (共鳴積分) Hij = <φi|H|φj> では、ハミルトニアン H は全対称である (すべての対称操作に対して符号が変化しない、つまり指標表の値がすべて 1 であるため、φi とφj が異なる既約表現に属する場合は Hij = 0 になります (既約表現が異なるため、いずれかの対称操作によりφiHφj の符号が反転し、全空間積分は必ずゼロになる)。つまり、既約表現の異なる波動関数間では相互作用はなく、それらの状態を別々に考慮して構わないことがわかります。
また、光励起を考える場合は、双極子遷移確率 χij = <φi|x|φj> がゼロでないことが必要になりますが、そのためには、φi かφj のいずれかが、x と同じ既約表現に属する必要があります。ここでは x 偏光の場合について説明しましたが、y 偏光、z 偏光の場合も、x を y, z に置き換えて読むことができます。
付録 3d: 3d 軌道の配位子場分裂 (CASTEP, PBE96)

図 A3d-1 立方結晶 BaTiO3 のバンドギャップ近傍の電子状態と波動関数。VBM は O22p 軌道、伝導帯は Ti4+ 3d 軌道でできており、3d 軌道が配位子場分裂で 2 組に分裂しています。
9 章で説明した波動関数の可視化を、教科書によく出てくる内容で確認してみましょう。孤立している球対称原子では 5 つの d 軌道は縮退しています。結晶の中では原子がおかれた環境が球対称でないため、d 軌道がいくつかの組に分裂します。これをつける。
犬井哲郎, 田辺行人, 小野寺嘉孝 著, 応用群論 －群表現と物理学－, 裳華房 (2007 年 増補版)
57 ユー、カルドナ、半導体の基礎 (Springer) p. 38 に、Koster によるものと、BSW (Bouckaert,
Smoluchowski, Wigner) によるものがあります。Koster 表記ではΓ3 が、BSW 表記ではΓ12 になります。
58 x2+y2+z2 の対称性は方位依存性がなく、r2 と同等です。つまり、球対称の意味です。
56

「配位子場分裂 (ligand field splitting)」と呼びます。どのような組に分裂するのかは、原子がおかれている位置の点群とその点群の指標表からすぐにわかりますが、どのような順番で、どれだけのエネルギーだけ分裂するかはわかりません。
立方晶 BaTiO3 は図 A3d-1 に示すようなペロブスカイト型の結晶構造をもちます。
ここで大事なのは中心の Ti 原子で、上下左右前後に 6 つの O2-イオンが配位しています。このような対称性は Oh 点群に属しており、表 A3c-1 のような指標表で表されます。この右端列の記述から、d 軌道がどの組に分かれるかがわかります。ここで () で括っている組は、量子計算では縮退した状態になる組です。ですから、(dxy, dxz, dyz) の組は T2g 表現に属して三重縮退した軌道 (t2g 軌道などと書かれます) を作ること、(dx2-y2, dz2 59) は Eg 表現に属して二重縮退した eg 軌道を作ることが読み取れます。
図 A3d-1 には CASTEP で計算した波動関数を示していますが、一番エネルギーが高いのが eg 軌道で O2-方向に伸びた (dx2-y2 60, dz2) でできていること、その次が t2g 軌
道で O2- のいない 8 つの方向へ (dxy, dxz, dyz) 軌道が延びていることが判ります。VBM
は O2- 2p 軌道でできていますが、ドーナツのような形をしていて、よく知っている鉄アレイ型の p 軌道とは似ていません。これは、2p 軌道が面内で縮退している (たとえば xy 面内の O2-イオンの 2px と 2py 軌道が縮退している) ためです。縮退している波動関数を、例えば φ1 = x + y, φ2 = x – y とすると 61、これらの電子密度の和は |φ1|2 + |φ2|2= 2(x2 + y2) = 2 Rxy2 と、xy 面内の原点からの距離 Rxy だけの関数になり、軸対称になっていることが確認できます。このように、縮退した px, py 軌道は面内で軸対称になり、ドーナツ状の波動関数になります。
付録 4a: 欠陥の帯電状態とフェルミ準位、ドーピング
反応エネルギー、相安定性や、半導体中の欠陥の安定性とドーピングなどの計算をする際には、熱力学で出てくる Gibbs-Duhem の式が基盤になります。
(A4a-0)
𝐺𝐺 = ∑𝑖𝑖 𝑛𝑛𝑖𝑖 𝜇𝜇𝑖𝑖
ここで、𝜇𝜇𝑖𝑖 は化学種 i の化学ポテンシャル、𝜇𝜇𝑖𝑖 は i の数、G は自由エネルギーであり、i には原子の他、電子の化学ポテンシャル＝フェルミ準位𝐸𝐸𝐹𝐹 を考慮する必要があります。自由エネルギーとしては、付録 3k で説明したように、対応する圧力や温度条件、あるいは必要とする近似に依り、内部エネルギーU、エンタルピー H、Helmholtz エネ
ルギーF、Gibbs エネルギーG などを選び、以下の議論を行います。
以下は、U にもとづき、文献 62 に報告されている ZnO の結果に基づいて説明します。図 4a-1 は、VASP を用いて HSE 混合汎関数 (ただし、バンドギャップを実測に合指標表では 2z2-x2-y2 になっていますが、x2+y2+z2 = r2 は球対称な成分なので、対称性の議論では無視できます。つまり、対称性の観点からは 2z2-x2-y2 は z2 と等価です。
60 d2x2-y2 軌道は、その座標からは図 A3d-1 では xy 平面の O の無い方向に向きそうですが、これは結晶の座標軸と原子軌道を表現するときの座標軸が違うためです。座標軸というのは人間が都合のいいものを選んだだけですので、無限の任意性があります。量子計算を行うと、座標軸のとり方に関係なく同一の計算結果が出てきますが、特定の座標を取ることにより、計算が簡
単になります（弾性テンソルなどが固有軸をとることで対角化されるのと同じ）。この場合は、結晶の座標軸と Ti 原子の座標軸が xy 面内で 45o ずれているわけです。このような原子固有の座標軸を「局所座標」と呼びます。実際の計算の際も、このように軸を取ることで計算を簡略化できることがあり、WIEN2k では入力ファイル作成時に局所座標を指定します。
61 縮退しているので、φ とφ が直交するのであれば、どうとっても構いません。一番簡単なの
1
2
はφ1 = x, φ2 = y ですが、簡単すぎるのでここの説明ではちょっと複雑な方で説明してみました。
62 F. Oba, A. Togo and I. Tanaka, Phys. Rev. B 77 (2008) 245202.
59

わせるため、EE 項の混合割合を 0.375 にしている) で計算した結果を示しています。
ここで VO は酸素空孔、Zni は格子間 Zn など、欠陥種を表しています 63。半導体中では EF が変化することにより、それぞれの欠陥 D の帯電状態が変化します。理想結晶
における電荷状態からの相対的な電荷変化を q とすると、その欠陥 Dq の生成エネル
ギー EfD,q(EF,µ) は
𝐸𝐸 𝑓𝑓 𝐷𝐷,𝑞𝑞 (𝐸𝐸𝐹𝐹 , 𝜇𝜇) = 𝐸𝐸D,q − 𝐸𝐸0 − 𝑛𝑛Zn 𝜇𝜇Zn − 𝑛𝑛𝑂𝑂 𝜇𝜇𝑂𝑂 + 𝑞𝑞�𝐸𝐸𝐹𝐹 − 𝐸𝐸𝑉𝑉𝑉𝑉𝑉𝑉 0 �(A4a-1)で計算されます。ここで、E0 は理想結晶、ED,q は欠陥 Dq を含むモデルの全エネルギーです。ここで、欠陥の全エネルギーは十分大きな supercell モデルで計算します。nZn,
nO は、欠陥モデル生成時追加あるいは削除された原子数。例えば(A4a-2)ZnO (結晶) => ZnZn + VO2+ + O + 2eという VO2+ 形成の欠陥反応を考えると、nO = -1、nZn = 0、q = +2 となります。µM は原子 M の化学ポテンシャルで、反応雰囲気によって変化します。図 A4a-1 で “Oxygenpoor” と “Oxygen-rich” と 2 種類のグラフが描かれているのは、ZnO が安定に存在す
る酸化極限と還元極限では µO が異なり、 EfD,q(EF,µ) も変わるためです。例えば、ZnOを還元しすぎると金属 Zn が析出しますので、Zn が析出し始める条件、ZnO と Zn が化学平衡にあり、ZnO 中のµZn が金属 Zn 中のµZn 値に一致する条件が還元極限です。
これと、式(A4a-0)にもとづく式(A4b-1)から、µO が決まります。酸化極限では O2 が過剰であるため、その境界では ZnO 中のµO は O2 のµO に一致します (いずれの場合も、
析出境界では相が共存すること、そのため、異なる相内の同種原子の化学ポテンシャ
ルは同じであること、に注意。µM の条件については付録 4b で説明します)。式(A4a-1)
より、“Oxygen-poor” と “Oxygen-rich” のグラフは、それぞれの欠陥種について nM･µM の違いだけ平行移動しているということがわかります。式(A4a-1) 右辺の最後は、帯電欠陥中の電子のエネルギーが EF によって変化する項です。この種の計算では、単純化のためにエネルギーの基準を理想結晶の EVBM (EVBM0) にとることが多いため、
図 A4a-1 のエネルギー原点も EVBM0 になっています。
以上の計算は、式(A4a-1) を見ると単純に見えますが、実際には以下のようにいく
つもの補正計算が必要になります。詳しくは文献 64, 65, 66 などをご覧ください。
(1) DFT のバンドギャップ問題。
混合汎関数法など実測のバンドギャップを再現できるモデルを使う。
(2) エネルギー基準 (EVBM) の補正
一般に、欠陥モデルと理想結晶モデルで計算した EVBM は一致しませんので、内殻軌道準位などを用いて EVBM を合わせる(3) DFT の欠陥モデルは一般には小さく、現実の欠陥濃度よりも数桁高い濃度でMA という書き方は Kröger-Vink の表記法と呼ばれ、理想結晶では A イオンが占めていた位置に M イオンが置換しているという意味です。つまり、ZnZn はもとのまま、VO は酸素イオン位
置が空孔 (vacancy) V に置き換わっている、Zni は格子間位置 (interstitial, i) を Zn イオンが占めていること、を意味しています。上付き文字の ·· は、理想結晶の状態よりも+2|e|だけ正電荷が増えていることを示していますが、最近では 2+ と書くことも増えています。これは、O2-イ
オンが抜けたために相対的に正電荷が発生している、Zn が格子間位置に来てイオン化してZn2+になっている、と考えるために発生した正電荷です。両辺での電気中性条件から自由電子
2e-の項が出てきます。もっと普通に、O2-イオンが抜けるときには中性分子 O2 で抜けるため、
電子 2 つを結晶に残していくと考えても構いません。
64 F. Oba, M. Choi, A. Togo, and I. Tanaka, Sci. Tech. Adv. Mater. 12 (2011) 034302.
65 S. Lany and A. Zunger, Phys. Rev. B 78 (2008). 235104.
66 S. Lany and A. Zunger, Model. Simul. Mater. Sci. Eng. 17 (2009) 084002.
63

の計算しかできない。そのため、希薄極限への補正計算を行う必要がある。
・高密度ドナーにより EF が高くなっている場合、それらの電子のエネルギー総和を ED,q から差し引く・欠陥が帯電している場合、欠陥間の静電・分極相互作用が発生するので、
このエネルギーを ED,q から差し引く 67

VZn0 => VZn(ε0/-)

ε-/2-

真性ドーピングリミット

図 A4a-1

ZnO 中の欠陥の生成エネルギーの帯電状態と EF 依存性。
(左図) 酸化極限 (酸素過多)、(右図) 還元極限 (酸素不足)式(A4a-1) からわかるように、q に帯電した欠陥の生成エネルギーは、EF に対して
q の傾きを持った直線になります。逆に言えば、図 A4a-1 の傾きから、欠陥の帯電状態がわかります。また、各欠陥で異なる q 毎に図の EF 全域にわたって直線が引けますが、図を見やすくするため、それぞれの EF において EfD,q(EF,µ) が最低になるもの
のみをプロットしています。例えば、VZn は p 型極限 (EF ~ EVBM0 = 0) では q = 0 が最も安定であり、EF ~ 0.7 eV で正孔を捕獲して q = –1 に帯電状態が変わります。つまり、電荷遷移準位 ε0/- が ~0.7 eV であり、これが第一アクセプター準位に対応します。
VZn は、さら EF が上がるとε-/2- ~ 2.4 eV でもう一つ正孔を捕獲して q = –2 の状態になります。VO は EF が深い場合では q = +2 価にイオン化して電子 2 つを放出している ”double donor” であり、~2.2 eV で電子を 2 つ捕獲して 中性欠陥 (q = 0) になることがわかります。ここで、q = +1 の VO+ は、VO0、VO2+ のいずれよりも EfD,q(EF,µ) が大きいために安定な帯電状態としては現れず、電荷遷移準位は ε0/2+ だけが現れ、これがドナー準位に対応します。
ここで注意が必要なのは、VO のドナー準位は ECBM ~ 3.3 eV よりも 1.1 eV も深く、浅いドナーとはならないことです。一般に、ZnO, SnO2, In2O3 などの酸化物半導体では、酸素欠損が浅いドナーとなって非ドープ試料でも電子伝導を示すと考えられていました。式(A4a-2) はこのような状況を想定し、VO が double donor となっていると仮定した反応式です。ところが、図 A4a-1 の結果も、他の計算モデルによる第一原理計
算の結果 68 もすべて、ZnO 中の VO は深いドナー準位を形成する点で一致しています。そのため現在では、ZnO をはじめとする非ドープ酸化物半導体のドナーの起源が見直されています。
図 A4a-1 では、Zni が ε+2/0 ~ 3.3 eV で浅いドナーになるように見えますが、この場
67
68

Y. Kumagai and F. Oba, Phys. Rev. B 89 (2014) 195205.
S. J. Clark, J. Robertson, S. Lany and A. Zunger: Phys. Rev. B 81 (2010) 115311.

合は生成エネルギーが EfD,q(EF,µ) ~4.2 eV と非常に大きいことが問題になります。温
度 T における欠陥 Dq の平衡濃度 ND,q は𝑁𝑁𝐷𝐷,𝑞𝑞 = 𝑁𝑁sites exp�−𝐺𝐺 𝑓𝑓 𝐷𝐷,𝑞𝑞 (𝐸𝐸𝐹𝐹 , 𝜇𝜇)/𝑘𝑘𝐵𝐵 𝑇𝑇�
(A4a-3)ffで計算できます (G D,q は欠陥生成の自由エネルギー。ここでは E D,q で近似する) が、
ZnO の作製温度 (欠陥の凍結温度) を 1,300 K としても、濃度はわずか ND,q / Nsites =6×10-17 であり、Nsites = 6×1022 cm-3 とすると ND,q = 2×106 cm-3 と、Si の真性キャリア密度 ~1010 cm-3 のわずか 1/10000 しかありません。つまり、仮に Zni が生成すれば浅いドナーになるものの、ZnO を作製する通常の条件では、無視できる濃度しかできないということです。そのため、図中に Hi や HO で示しているように不純物水素が浅いドナーの起源ではないかと考えられています。
図 A4a-1 からはまた、ドーピング限界も見積もることができます。VO+2 の EfD,q(EF,µ)は EF が~1.8 eV 以下で負の値になります。つまり、p 型 ZnO を作ると自発的に VO を生成するという意味です。しかしながら、VO2+ を生成すると電子を放出して EF を上昇させることになりますので、結果として、このように EfD,q(EF,µ) が負になる欠陥が存在すると、EF はこの付近で固定され、これが ZnO の p 型ドーピング限界となります。ただしこれは、「何をしてもこれ以上の範囲に EF が動くことは（まず）無い」という意味での限界です。実際の限界値はこのドーピング限界よりもはるかに低い可能性があります。例えば、Oxygen-poor 極限での HO, Hi は CBM 下 1 eV 付近に EF を固定するように見えます (このような計算結果は複数の酸化物の計算で報告されています) が、負の生成エネルギーが現れる場合には、水酸化物などの異相が生成することが多く、それによる主相の安定条件も考慮する必要があります。また、酸化物中ではHO, Hi は ZnO と同様に浅いドナーになりそうですが、例えば PbO では、水素ドープ
により Pb が還元されるために n 型ドーピングが実現できません 69。
以上では一つの欠陥種に着目してドーピング限界を議論しましたが、厳密には、すべての欠陥を考え、ある EF における平衡濃度を 式(A4a-3) で計算すると、放出される全電荷 は Qe = ∑ qN D ,q で計算できます。第一原理計算によって伝導帯と価電子帯の状態密度 DC(E)、DV(E)が計算できるので、EF を与えると、伝導帯中の自由電子濃度
ne と価電子帯中の自由正孔濃度 nh も計算できます。電荷中性条件 Qe + nh = ne を満足する EF が、「平衡フェルミ準位 EF,e」となります。このような計算をすると、純粋なSnS では VS2+, VSn2-, SnS- が主要な欠陥ですが、それぞれの平衡濃度は 5.0×1015,
3.9×1015, 5.2×1015 cm-3 となります。この際の正孔と電子の濃度がそれぞれ 1.3×1016,
1.0×1016 cm-3 であり、EF,e はほぼギャップ中央の EF,e – EV = 0.41 eV (Eg ~ 1.1 eV) とな
ります。このことから、真性 SnS は強く電荷補償された結果として、少し p 型にドープされた半導体であることを説明できます 70,71。
付録 4b: 構成原子の化学ポテンシャルと相の安定条件
式(A4a-1) には構成原子 M の化学ポテンシャル µM がでてきますが、この上限、下限は次のように母相である ZnO が安定になる条件で決められます。実際の合成実験
では化学ポテンシャルを意識することはありませんが、雰囲気ガスの分圧などを制御
パラメータにとるため、物理化学的には化学ポテンシャルを制御して合成を行ってい
M. Liao, S. Takemoto, Z. Xiao, Y. Toda, T. Tada, S. Ueda, T. Kamiya, and H. Hosono, J. Appl. Phys. 119(2016) 165701.
70 Vidal, J. et al., Appl. Phys. Lett. 100 (2012) 032104.
71 F.-Y. Ran, Z. Xiao, Y. Toda, H. Hiramatsu, H. Hosono and T. Kamiya, Sci. Rep. 5 (2015) 10428.
69

ることになるということに注意してください。そのため、構成成分の化学ポテンシャ
ルが与えられたときに、どのような相が安定となるかを判断することが、この章の目的になります。

図 4b-1Cs2SnI6 の安定条件。

ZnO のような二元系化合物では、平衡状態では、∆μZn + ∆μO = ∆GfZnO(A4b-1)f
を満足します。ここで∆G ZnO は ZnO の生成自由エネルギー (DFT で計算した全エネ
ルギーにもとづく生成エンタルピー∆H(ZnO)で近似)、∆μZn, ∆μO は、Zn, O の単体の化
学ポテンシャルからの差であり、実験条件で決まるパラメータです。ただし、式(A4b1) の制約により、∆μZn と∆μO のいずれかを決めれば、もう一つの化学ポテンシャルは自動的に決まるので、自由度は１しかありません。ZnO が単体よりも安定になるためには、構成成分の化学ポテンシャルが担体中の値より低くないといけませんので、∆μZn ≤ 0, ∆μO ≤ 0 を満足する必要があります。"O-poor/Zn-rich" 極限は単体の Zn が析出しはじめる境界、つまり単体の Zn と ZnO が平衡状態にあり、両者で∆μZn が等しい条件となります。よって、∆μZn = 0, ∆μO = ∆GfZnO と、全ての化学ポテンシャルが決まります。同様に "O-rich/Zn-poor" 極限では単体の O2 と ZnO が平衡状態にあるという条件から、∆μZn = ∆GfZnO, ∆μO = 0 となります。
以上のことを、三元系化合物であるダブルペロブスカイト型化合物 Cs2SnI6 で見てみましょう 72。式(A4b-1)に対応する母相 Cs2SnI6 の平衡条件は2∆µCs + ∆µSn + 6∆µI = ∆H(Cs2SnI6) = −10.49 eV(A4b-2)となります。Cs-Sn-I 系には化学組成の異なる多くの化合物が存在し、それらが異相として析出する可能性があります。それらすべての化合物の生成エンタルピーを計算した結果、CsI, SnI2, SnI4, CsSnI3 が Cs2SnI6 の安定条件境界で共存することが確認され、これらの異相よりも Cs2SnI6 が安定となる条件として、以下の 4 式が導かれます。
∆µCs + ∆µI< ∆H(CsI)= −3.80 eV(A4b-3)∆µSn + 2∆µI
< ∆H(SnI2) = −1.69 eV(A4b-4)∆µSn + 4∆µI
< ∆H(SnI4) = −2.40 eV(A4b-5)∆µCs + ∆µSn + 3∆µI < ∆H(CsSnI3) = −5.77 eV(A4b-6)式(A4b-2) より、実験条件に関する自由パラメータは 2 つに絞られますので、図 A4b1 では、∆µSn と∆µI をパラメータとし、式(A4b-3)～(A4b-6)の 4 本の直線を描いています。例えば式(A4b-3)の条件∆µCs + ∆µI = −5.25 eV – 0.5 ∆µSn – 2∆µI < −3.80 eV(A4b-3’)0.5∆µSn + 2∆µI > 1.45 eV(A4b-3’’)から、A-D の直線よりも上でないと Cs2SnI6 が安定でないことがわかります。式(A4b3)～(A4b-6) についても同様の検討を行うことで、Cs2SnI6 は図 A4b-1 の四角形 ABCD
72

Z. Xiao, Y. Zhou, H. Hosono and T. Kamiya, Phys. Chem. Chem. Phys. 17 (2015) 18900.

の範囲内だけで安定であることがわかります。∆µI が小さい、つまりヨウ素分圧 (活量) の高い条件では CsI が安定になります。∆µCs が高い、Cs-poor な条件では CsSnI3が析出すると判断できます。
付録 5: 圧力誘起相転移への応用例
体積－全エネルギー曲線の応用例として、圧力相転移の理論計算の例を挙げておきます。図 A5-1(A) には、異なる結晶構造について分子単位あたりの体積（つまり単位格子体積）を変えながら全エネルギーを計算した結果を示します。常温常圧では BaS
は NaCl 型構造 (B1 構造) をとりますが、高圧下で CsCl 型 (B2 構造) へと構造転移を起こします。熱力学から、圧力は P = ∂E / ∂V で計算できるので、図 A5-1(A) のよ
うな E–V 図から、図 A5-1(B) の E–P 図が得られます。P, T = 一定の条件でどの相が安定になるかはギブスエネルギーG = E + PV – ST で決まります。高圧では ST 項は無視できますが、PV 項が無視できません。そこで、同じ圧力で構造転移に伴う体積変化∆V を考慮すると、図 A5-1(B) の ”B1, E” と ”B2, ∆H” の大小が交差する点が構造転移を起こす転移圧力であり、約 9 GPa であることがわかります。これは、文献値 ~6
GPa と良い一致をしています。同様の計算を NaCl についてしてみると、計算値は 28
GPa と、報告値 30 GPa をよく再現することが確認できます。

Energy (eV)

-2

B1

-4

B2

(B) -2-3
Energy (eV/molecule)

(A) 0
-6
-8

-10
-12

0

50

100

150

Volume (A3 )

200

250

-4
-5
-6
-7
-8
-9
-10
-11

P transition

B1, E
B2,E
B2, ∆H0

10 20 30 40 50 60 70 80 90 100

Pressure (GPa)

図 A5-1 BaS の圧力誘起相転移の計算。(A) E–V 図、(B) E–P 図。高圧化では PV 項が無視できないため、転移圧力は ∆H = ∆E + P∆V がゼロになる点できまります。これは、図(B)では “B1,E” と “B2, ∆H” が交差する点となります。
付録 6a: 実験で決めがたい構造の決定: アモルファス物質の局所構造
第一原理計算の利点の一つは、実験的には決定の難しい局所構造に関する情報を与えてくれることです。付録 6a~6d ではいくつかの例について紹介します。
結晶では有限個のパラメータで原子の位置を決められ、これらを X 線/中性子回折法などの構造解析で決定できます。しかしながら、アモルファス物質の場合は、電子密度の動径分布関数のように、平均的な構造情報しか実験では得られません。
そこで「局所構造を決められる」という第一原理計算の利点が活きるわけですが、構造緩和計算では、初期構造モデルの近くで全エネルギーが極小値を持つ構造 (local
minimum) に収束してしまいます。そのため、なるべく実際に近い初期構造モデルを作らないといけません。有限温度の動力学シミュレーションができる MD 法で、モデ
ル構造を一旦融解させ、再固化させてアモルファス構造を作るのが直接的で単純な方法です。しかし、第一原理 MD 法は計算が遅いため、時系列を追って熔融－再固化をシミュレーションするには非常に長い時間がかかります 73。そこでまず、同じ構成原73

現在のコンピュータと第一原理計算法では、10,000 ステップ以上の第一原理 MD 計算も現実的に実行可能で、実際にアモルファス構造を第一原理計算だけで作ることも一般的になってい子をもち構造が近い結晶の構造や、アモルファス材料の実験データに合うように、大雑把に原子間ポテンシャルを決めます。これを用いて古典 MD 法で溶融－再固化シミュレーションを行い、アモルファスの初期構造モデルを作ります (melt-quench 法)。
この構造から第一原理構造緩和計算をすることで、アモルファス材料の安定構造や電子構造の計算ができます (図 A6a-1) 74。

図 A6a-1 MD 法と VASP (PAW, LDA)で求めたアモルファス InGaZnO4 の原子配列(左)。84 原子を持つ supercell 75について計算した。右図は、最低非占有準位の波動関
数 |ψ|2 を等高面で描いている。
付録 6b: 実験で決めがたい構造の決定: 水素の位置
表 A6b-1 水酸化物結晶の生成エネルギーの計算値と文献値
計算値(eV/molecule) 文献値(eV/molecule) 文献値(kJ/mol)
Mg(OH)2 8.79
8.74
924.66Ca(OH)2
9.52
9.32
986.09Sr(OH)2
9.31
9.16
968.89Ba(OH)2
8.85
8.95
946.3
Fe(OH)2
5.25
5.43
574.04水素の位置は X 線構造解析では決められないため、水素を含む物質の構造データ
ベースに水素の位置が載っていないものが多くあります。このような場合にも、基本結晶構造から水素位置の見当がつく場合があり、それを基にして構造緩和計算させる
ことで水素の位置を決定できます。そうでない場合にも、適当に水素の初期位置を入力し、高温での第一原理 MD 計算などを行うことで最安定構造を見つけられることもあります。ただし、このようにして得られた構造が妥当であるかどうかは、慎重に判断する必要があります。簡単に入手できる比較データとして、上記の方法で水素位置を決めた水酸化物の生成エネルギーの計算例を表 A6b-1 にまとめました。いずれも高い精度で生成エネルギーを再現していることがわかり、得られた構造モデルの妥当性を支持しています。ここには挙げていませんが、ICSD 76 #28275 の水素のデータが含ます。
K. Nomura et al., Phys. Rev. B 75 (2007) 035212.
75 材料シミュレーションでは、基本となる単位格子を複数個つなげて大きい格子の基本構造をつくることが多い。これを supercell と呼びますが、結晶学ででてくる本来の超格子 (superlattice)
や、薄膜などででてくる人工超格子 (artificial superlattice)とは意味が異なります。
76 Inorganic Crystal Structure Database, http://www2.fiz-karlsruhe.de/icsd_home.html
74まれていない Mg(OH)2 から構造緩和をすると、表 A6b-2 のように、実測値よりも大きく c 軸長が縮みます (水素が抜けているので当たり前ですが)。報告されている水素のない構造中の大きな格子間位置に化学組成に合うように水素を入れたモデルから構造緩和をすると、表 A6b-2 のように、実測値と格子定数がよく一致しますので、この構造が確からしいということが確認できます。実際に、中性子回折による水素位置を報告している ICSD #34401 と比較しても、水素を入れた構造緩和結果がよく実測を再現できていることがわかります。
表 A6b-2

Mg(OH)2 で水素を入れない場合と入れた場合の構造緩和結果。
a (Å)
c (Å)
z(O)
z(H)
ICSD (#28275)
3.147
4.768
0.217
3.235
3.477
0.252
構造緩和 (H 無し)
3.162
4.721
0.223
0.429
構造緩和 (H 有り)
ICSD (#34401)
3.142
4.766
0.22160.4303

付録 6c: 実験で決めがたい情報: 簡便な欠陥計算
付録 4a では ZnO を例に、欠陥のできやすさとドーピングの EF 依存性の計算について説明しましたが、この計算はかなり手間がかかります。ここでは、欠陥を含めたモデルの単純な計算からどのようなことがわかるかについて説明します。

図 A6c-1 Si クラスターの電子構造 (MOPAC7)。Si (大きい球) を中心として 4 つのSi が配位した Si5 クラスターを基本に、外側の Si 原子の非結合軌道をすべて水素原子(小さい球) で終端した Si5H12 クラスターと 1 つだけ非結合手を持つ Si5H11 クラスターについて比較している。両端には、Si の電子構造と、Si 欠陥と非結合手があるときの電子構造を模式的に描いている。その内側は、HOMO, LUMO, ダングリングボンド(DB)に対応する波動関数。
図 A6c-1 は、 結晶 Si を簡略化した Si5 クラスターの電子構造を半経験分子軌道法
ソフトウェアである MOPAC で計算したものです。ここで分子軌道法を使ったのは、
バンド計算で欠陥や不純物を扱うと非常に大きな supercell を組まなければならず膨
大な計算時間がかかること、さらに、そうやって得られた結果は厳密であったとしても直観的に理解しにくいことが多いためです。図 A6c-1 中図のエネルギー準位図をみると、非結合手をすべて H で終端した Si5H12 クラスターは HOMO と LUMO の間に
大きなエネルギーギャップがあるのに対し、非結合手を一つ持つ Si5H11 や全く終端していない Si5 クラスターでは、Si5H12 クラスターのエネルギーギャップの中央付近に余計な占有準位ができていることがわかります。これが「Si のダングリングボンド」
で、HOMO (DB)と書いてある波動関数の絵から、これらが H で終端されていない Si
上の p 軌道であることがわかります。Si 中にこのような状態がたくさんあると、LEDや太陽電池の効率を落したり、トランジスタの移動度を下げたりします。このような状況を改善するため、アモルファスシリコン中のダングリングボンドは水素で、太陽電池の Si 表面・端面は酸化膜などを作って終端し、ダングリングボンドを不活性化させています。
ただし、Si5H12 クラスターのエネルギーギャップの計算値は 10 eV 近くもあり、結晶 Si の文献値 1.12eV とは全く合いません。DFT と同様の「バンドギャップ問題」の寄与もありますが、それ以上に大きい原因が、Si 原子 5 個だけの非常に小さなクラス
ターで計算したことにあります。固体中の電子構造は、隣同士の波動関数が無限につながることでバンド分散をつくっています。特に良い半導体（有効質量が小さい）ほどバンド分散は大きく、Si ではバンド幅が 10 eV 近くあり、これが原子軌道の上下に広がります。価電子帯、伝導帯とも、このバンド幅の半分ずつが上下に広がっている
ので、結局バンドギャップも、小さい Si クラスターよりも 10 eV 近く小さくなって、
1.12 eV ほどになるわけです。実際に近い計算をするためにはクラスターのサイズを大きくすればいいのですが、Si の場合では、100 原子程度のクラスターではエネルギーギャップが収束しないようです 77。
付録 6d: 実験で決めがたい情報: 簡便なドーピング計算

図 A6d-1 不純物原子１つを含む Si の電子状態 (VASP, PAW, PBE96)。 (左, 中) P ドープシリコンの CBM 付近の波動関数。(右) Si 欠損とさまざまな不純物を含む Si の全状態密度。フェルミエネルギーを基準に揃えてある。
付録 4a で説明したように、欠陥やドーパントの帯電状態は EF によって変わります。
77

Si は、第一 BZ の半端な位置にエネルギーバンドの極小値をもっていますから、中距離の波動関数の位相が電子構造に大きな影響を与えていることは想像できます。一方で、酸化物結晶ではエネルギーバンドの極値はΓ点と BZ 境界にあるものも多く(BaTiO3 もその例)、この場合は最近接程度のクラスターでも比較的良い結果が出ることもあります (BZ 境界の波動関数は、隣の単位格子との反対称軌道であることを思い出しましょう。それよりも遠い波動関数との位相は+1 と-1 が交互に続くだけなので、結局、最近接単位格子程度を考えてもΓ点と BZ 境界の電子構造をよく表現できることになります)。

そのため、ドーピングできるかどうか、また、ドナー準位やアクセプター準位がどこに来るかは、ドーピングモデルの単純な計算では正しいことはわかりません。
しかしながらそれでも、ある不純物や欠陥がドナー／アクセプターになるかならな
いか、について定性的な知見を与えてくれます。図 A6d-1 は、P ドープ Si の計算をした結果ですが、2×2×2 の supercell (ブラベー格子が 8 原子をもちますから、この supercell
は 64 原子を含んでいます) をつくり、角にある Si 原子を P 原子に置き換えてあります。構造緩和させたあとの波動関数について、伝導帯側の最低エネルギーの準位(ECBM) と、その 0.18 eV 上の波動関数を描いてあります。ECBM では波動関数が P 位置に集まっているものの、それよりも外側にも無視できない密度があります。つまり弱く局在化した状態があり、これが、P ドープシリコンの 45 meV (文献値) の深さにできるドナー準位に対応していると考えられます。その上の準位は純粋 Si の伝導帯に入っており、supercell 全体に波動関数が拡がり、自由電子の状態になっていると考えられます。右図には全状態密度を示していますが、P ドープ Si のフェルミエネルギーは伝導帯に近く 78、また、伝導帯とつながった準位を作っていることがわかります 79。
P を取り除いてダングリングボンドが 4 つある場合の計算 (VSi) からは、図 A9d-1 と同様、バンドギャップのほぼ中央に EF が来ることが確認できます。B ドープの場合にアクセプターになること、Fe や O 不純物が深い準位を作ることも読み取れます。
付録 7a: スピン軌道相互作用
以下に示すように、半導体の価電子帯上端はスピン軌道 (SO) 相互作用により分裂し、スピン分裂バンドを形成します。SO 相互作用は、スピンの磁気モーメントを仮定すれば、係数を 2 倍間違えることを除いて、古典的な電磁気学でも説明できます。
しかしながら、スピンがちゃんと量子方程式に現れるのは、Schrödinger 方程式ではな
く相対論的量子力学方程式である Dirac 方程式であるため、第一原理計算でも、相対論補正を考慮した場合のみ、このバンド分裂が出てきます。
第一原理計算で相対論効果をどのように取り入れるかはプログラムによって異なります。多く使われているのは、原子の波動関数や擬ポテンシャルは Dirac-Coulomb
方程式を解き (Full relativistic)、結晶軌道の計算では相対論補正を摂動項として計算する方法 (Semi-relativistic) です。この場合に一番大きな補正項が SO 相互作用になり、次式で表されます。
H SO =

(∇V (r ) × p ) ⋅ σ
4m 2 c 2(A11-1)
ここでσは Pauli のスピン行列、V(r)は他の電荷が作る静電ポテンシャル、p は運動量です。V(r)が中心場の場合 (角度によらず r = |r| のみで決まる)、

DFT のバンドギャップ問題のため、この議論には注意が必要です。Si での実験値とのずれは
約 0.5 eV ありますが、もしも P の準位はそのままで、本当の CBM が 0.5 eV 高い位置にあるのであれば、P の準位は浅いドナーではなく、ミッドギャップ準位と解釈する必要があります。
より正確なバンドギャップが得られる GW 近似や混合汎関数法などを使うと、この問題を軽減できます。
79 64 原子を含む supercell が、孤立した欠陥や不純物の議論をするには十分な大きさとは限りません。supercell が十分大きくないために P 軌道がバンドを作ってしまい、伝導帯とドナー準位がつながっているという計算結果が出ている可能性も高く、慎重に検討する必要があります。
実際、この supercell の P 濃度は 1/64 = 1.6% = 8×1020 cm-3 にもなっており、実際の Si でも非常に高ドープな状態をシミュレートしていることになります。
78

H SO =Ze 2
s ⋅ l = λs ⋅ l
2m 2 c 2 r 3(A11-2)となり、このλを「スピン軌道結合係数」と呼び、全角量子数 j = l + s で jz = 1/2 と 3/2の場合には、バンドの分裂∆0 は 3λ/2 になることが知られています。
図 A7a-1 には SO 相互作用を入れて計算した Si と GaAs のバンド構造を示します。
GaAs の M–Γ上のバンドを見るとはっきりわかりますが、VBM にはΓ点で重なった、
バンド分散が違う 2 つのバンドがあり、軽い正孔 (light hole: lh)と重い正孔 (heavy
hole: hh) バンドと呼ばれます。その下に SO 相互作用で∆0 離れた位置にもう一つバンドがあり、これは「スピン軌道分裂 (soh, split off hole) バンド」と呼ばれます。Si の場合にも、VBM はΓ点で重なった 2 つのバンドからなっているように読み取れますが、実験的には 44 meV ほど分裂しています。表 A7a-1 に見られるように、高い精度
で∆0 が計算できることがわかります。

図 A7a-1 スピン軌道相互作用を考慮した Si, GaAs のバンド構造 (VASP, PAW, PBE96)。
X 点が<100>、M 点が<110>方位。
表 A7a-1

代表的な半導体の SO 分裂値∆0 。
括弧内は VASP, PAW, PBE96 での計算値。
∆0 [meV]∆0 [meV]Si
44 (47)
Ge
295
(341)
GaAs 341 (328) InAs 380
GaSb 750
InSb 810

付録 7b: フェルミ面バンド構造図では BZ の一部のエネルギー準位しか描画しないため、正しいバンド
ギャップの値を読み取れない可能性があります。この問題については、状態密度を使
えば BZ 全体の和が読み取れます。一方、電子伝導率、磁化率などの物性に直接関与するのはフェルミ準位付近の電子ですが、バンド構造を見ただけでは、全体像がつかみにくいという問題があります。このような場合、フェルミ面 (Fermi Surface: FS) を見ることが有効です。
一つのバンドに着目してフェルミ準位を横切る k 点 (kF,フェルミ波数) を集めると、
図 A7b-1 のように、第一 BZ 内で曲面を描きます。これが FS で、たとえば金属 Na の場合、かなりきれいな球面になります。自由電子モデルから k F = (E F 2me )1 / 2 /  ですから、
このことは、Na のバンドがいずれの方位にも同じ有効質量をもち、自由電子近似がよく成立することを示しています。図 A7b-1 の白線の枠が第一 BZ の境界を示していますが、Na は体心立方構造をもつため、FCC の Si とは第一 BZ の形が異なります。
Mg も六方稠密構造をもつため、六方柱状の第一 BZ をもちます。Mg の FS には、中心に平べったい円盤、上下の辺に半円状、それと、六角形周囲に広がった 3 種類の FSが見えます（ギザギザが見えるのは計算しているメッシュが荒いため）が、これは、複数のバンドが Mg のフェルミ準位を横切っているためです。

図 A7b-1 金属のフェルミ面。
バンド理論では、「一つのバンドに 2 つの電子が入れるので、Na は半分しかバンドを埋めずに金属になる」という説明がされますが、それだと、Mg では電子がバンドを埋めて絶縁体になるはずです。実際にそうならないのは、2 次元以上では、多面体の第一 BZ と同じ体積を球体で実現すると、球体の一部が多面体の外にはみ出すためです。つまり、Mg が持つ 2 つの価電子は第一 BZ の体積をぴったり埋め尽くす数ですが、FS は球に近い形のため、平行六面体の第一 BZ を全部埋め尽くす前に、一部が第二 BZ にはみだします。そうなると、第一 BZ にも第二 BZ にも非占有準位が残る
ため、Mg は絶縁体にはなりません 80。BZ は逆格子ベクトルの整数分の並進移動によって第一 BZ 内に還元できるので、第二 BZ にはみだした FS は、第一 BZ では負の曲率を持つ曲面になります。これは Al の FS で見ることができます。
半導体の EF はバンドギャップの中央付近にあるので、金属と同じ意味での FS を描くことはできません。しかし、CBM, VBM 付近の任意のエネルギーでバンドが横切る面を描くことができ、これによってバンドの対称性や異方性がわかりやすくなります。
これらは n 型、p 型半導体の FS と考えることもできますが、「等エネルギー面」と呼ばれることもあります。

図 A7b-2 Si の CB と VB のフェルミ面。VB のフェルミ面はわかりやすいように二次元でも表示した 81。
図 A7b-2 に Si の CB, VB の FS を描いていますが、これからいくつか重要なことがわかります。まず、Si の CBM は第一 BZ 内に 6 か所あります。Si は立方晶ですから、
x, y, z を入れ替えても全く同じ状態があります。CBM はバンド構造のΓ-X 線上の (kCBM, 0, 0) にあるので、同じ状態が (±kCBM, 0, 0), (0, ±kCBM, 0), (0, 0, ±kCBM) の 6 か所に80

81

第二 BZ に拡張されたエネルギーの高いバンドは第一 BZ の一つ量子数の高いバンドに対応するため、Mg の 3s バンドを埋め尽くす前に 3p バンドに電子が入ると言うのと同じである。
P.Y. ユー, M. カルドナ著, 半導体の基礎, Springer (1999).

あることになります。FS を描くと、このことがはっきりわかります。これらの CBMは同じエネルギーを持ちますが、波数ベクトルは異なるので、フォノンによる散乱を受けやすくなります。このような散乱が「谷間 (inter-valley) 散乱」と呼ばれます。
付録 7a の SO 分裂で説明したように図 A7b-2 の Si の VBM には hh バンドと lh バ
ンドが重なっています。VB の FS も 2 つの面がありますが、内側の FS が隠れて見えないため、図 A7b-2 右図に二次元断面図も示しています。これから、内側の FS (lh バンド) は少し歪んでいますが球形に近く、いずれの方位へも同じ有効質量をもつことがわかります。一方外側の FS (hh バンド) は金平糖のように大きく球形から歪んでいます。これは、図 A7a-1 で hh バンドの曲率が M–Γ上とΓ–X 上で大きく異なっている、
つまり、<110>方位の正孔の有効質量が<100>方位よりも重いことに対応しています。
付録 7c: Si の有効質量
Si は、広く使われている物質であるにもかかわらず、複雑なバンド構造を持ちます。
表 A7c-1 に、Si の有効質量の文献値をまとめました。有効質量は、一般的な関係として(me* −1 )ij = (d 2 E(k ) / dki dk j )/  2(A7c-1)で与えられます。図 A7a-1 にみられたように、Si の CBM はΓ–X 線上にあります。つまり、立方格子の a 軸方向に走っているある運動量を持った電子が最低エネルギー
Ekx(k)をもちます。この時 d2Ekx(k)/dkx2 で決まる mle*と、d2Ekx(k)/dky2, d2Ekx(k)/dkz2 で決まる mte*は異なり、それぞれ、「縦有効質量」「横有効質量」と呼ばれ、実験的に mle* =0.98me と mte* = 0.19me と得られています。第一原理計算の結果からは、mle* = 0.96me,
mte* = 0.09me が得られます (VASP, PAW, PBE)。mle*の一致は良いのですが、mte*は 1/2程度に過小評価しています。k·p 摂動法という近似から(

me* = 1 + 2 P 2 / me E g) m−1(A7c-2)と、Eg が小さくなるほど有効質量が小さく傾向がわかりますが、このことは、DFT がバンドギャップを過小評価することにより、me*もまた過小評価されていることを示
唆しています。
Si の正孔に関しては、VBM 付近では 3 つのバンドがあり、lh, hh, sh バンドに対して異なる有効質量があります。図 A7b-2 の FS で見たように軽い正孔の計算値はいず
れの方位にも mlh* ~ 0.19me であり、異方性は小さいことがわかります。重い正孔、分裂バンドについては<110>方位では mhh* = 0.83me, msoh* = 0.14me であるのに対し、<100>
方位では mhh* = 0.26me, msoh* = 0.22me と大きく異なります。
表 A7c-1

82
83e

Si の有効質量 82,83。単位は電子の静止質量 me。
電子
縦質量
横質量
状態密度
mle* = 0.98
mte* = 0.19
mde* = 0.33
正孔
重い正孔
軽い正孔
mhh* = 0.49
mlh* = 0.16
SO 分裂バンド
状態密度msoh* = 0.29
mdh* = 0.55

S.M. Sze, Physics of Semiconductor Devices, 2nd ed., Wiley (1981).
http://ecee.colorado.edu/~bart/book/effmass.htm

付録 7d: 状態密度有効質量と有効状態密度半導体の物性解析やデバイスシミュレーションでは、状態密度をパラメータで入力することがあります。NFE の状態密度は
𝐷𝐷𝐶𝐶 (𝐸𝐸) = 𝐷𝐷𝐶𝐶0 �𝐸𝐸 − 𝐸𝐸𝐶𝐶
(A7d-1)
3/2
3/2
3
3
∗
𝐷𝐷𝐶𝐶0 = �√2/𝜋𝜋 2 ��𝑚𝑚𝑒𝑒∗ /ℏ �𝑀𝑀𝐶𝐶 = �√2/𝜋𝜋 2 ��𝑚𝑚𝑑𝑑𝑑𝑑/ℏ �(A7d-2)*で与えられ、一般的な「キャリア有効質量」me と「CBM の多重度 (縮重度)」MC で決まりますが、多重度を含んだ有効質量 mde*を「状態密度有効質量」と呼びます。
ここで Si の場合について考えてみます。すでにキャリア有効質量で見てきてきた
ように、Si の CBM 付近では有効質量に異方性があり、エネルギーは
2
2
2
2
2
2
ECBM ,kx (k ) = ( 2 / 2 ) k x / mte* + k y / mle* + k z / mle*
(A7d-3)

(

となり、

(

mde* = mle* mte* mte*)

) =0.33me
1/ 2(A7d-4)が得られます。既に FS のところで説明したように、Si では CBM の状態が 6 か所にあるので、多重度 MC は 6 になります。式(A7d-2) の MC を mde*に繰り込むことができますが、そうすると
1/ 2
2/3(A7d-4)mde* ' = (mle* mte* mte* ) M C =1.08meが得られます。上述のように、デバイスシミュレーションなどでは伝導帯の状態密度の代わりに、有効状態密度 NC を入力することがあります。NC も
3/ 2*(A7d-5)N C = 2(2πmde
k BT / h 2 ) M Cと、mde*がわかっていれば計算できます。(A7d-4) 式の値を使うと、2.8×1019 cm-3 が得
られ、文献 76 に載っている値と一致することが確認できます。
正孔の場合は、(A7d-6)
DV (E ) = DV 0 EV − E* 3/ 2(A7d-7)/ 3 MV
DV 0 = ( 2 / π 2 ) mdh((*
NV = 2 2πmdh
k BT / h)) M

2 3/ 2V(A7d-8)となりますが、Si の VBM は lh, hh バンドが重なっています。このような場合 D(E)はそれぞれの和として、
3/ 2* 3/ 2* 3/ 2(A7d-9)mdh= mlh* + mhhで表され、mdh* = 0.57me が得られます。ここで、VBM があるΓ点は第一 BZ 内に一つしかないので、MV = 1 です。ただし、hh バンドは球対称ではないため、これを考慮す
ると 0.81me になります。この mdh*の値を 式(A7d-8) に入れると NV = 1.04×1019 cm-3 となり、これも文献 76 の表と一致します。
式(A7d-1), (7d-2), (7d-6), (7d-7) によれば、第一原理計算から得た DC/V(E)を使い、DC/V(E)2 をプロットして傾きから mde*, mdh*が求められることになります。しかしながら実際にはこのようなプロットから一意的に状態密度有効質量が決まるとは限りません (実際に Si の VBM ではこの関係に乗らない)。そのため、mde*や NV などを第一
原理計算から得る場合には、式(A7d-1), (7d-2), (7d-6), (7d-7) が計算で得られた DC/V(E)を再現するように mde*, mdh*を決めるほうが良いかもしれません。このような計算を
SnO TFT のシミュレーションに応用した例を図 A7d-1 に示します 84。
図 A7d-1 では、異方性のある VB バンドでも自由電子的な状態密度 (𝐷𝐷(𝐸𝐸) ∝84

Y. Ogo, H. Hiramatsu, K. Nomura, H. Yanagi, T. Kamiya, M. Kimura, H. Hirano and H. Hosono, Phys.
Stat. Sol. A 206 (2009) 2187.

�𝐸𝐸 − 𝐸𝐸0 ) であるためにデバイスシミュレーションに使っても精度の出せる mdh* を出せましたが、一般の𝐷𝐷(𝐸𝐸)はもっと複雑です。そのような場合は、7 章で説明したような方法で有効状態密度や、それに対応する状態密度有効質量を計算することができます。

図 A7d-1 VASP, PAW, PBE で計算した SnO の状態密度と、状態密度有効質量のフィッティング結果。これから mdh* = 2.05me を得た。
付録 7e: 電子と正孔、正の有効質量と負の有効質量9 章で、電子の有効質量は 𝑚𝑚𝑒𝑒 ∗ = ℏ2 (𝜕𝜕 2 𝐸𝐸/𝜕𝜕𝜕𝜕 2 )−1 と表せ、バンドの曲率 𝜕𝜕 2 𝐸𝐸/𝜕𝜕𝜕𝜕 2が負になると質量も負とみなせることを説明しました。半導体の場合はこのような状
況は価電子帯上端付近で発生します。外部電界 Fx がある場合には電荷 –e の電子は
d2x / d t2 = –eFx / me* と、me* が正で電荷 –e が負の電子は Fx と逆方向に加速されます。正孔を電子の非占有準位として定義すると、電子の移動と正孔の移動は逆方向になりますので、正孔の電荷を正の e とし、質量の符号も逆に mh* = -me* と定義すると、d2x / dt2 = eFx / mh* と、mh* が正で電荷 e が正の正孔は Fx と順方向に加速されることになり、私たちの直観に自然な定義となります。一方、運動量は、電子、正孔のそれぞれに me*ve, mh*vh となり、電界 Fx とはそれぞれ逆方向、順方向となり、やはり、直観に自然な定義であることが確認できます。
4. 負曲率の正孔
3. 負曲率の電子

2. 正曲率の正孔
1. 正曲率の電子

力 F = -|e|E が印加された時
1. 正曲率の電子
運動量、速度とも、F の向きに加速
2. 正曲率の正孔 (正の電荷)
運動量はFの向き、速度は逆向きに加速
3. 負曲率 (負の質量) の電子運動量はFの向き、速度は F と逆の向きに加速
4. 負曲率 (正の質量) の正孔 (正の電荷)
運動量、速度ともに F と逆の向きに加速

図 A7e-1 電子と正孔の有効質量と運動量

一方、例えば超伝導では、運動量が+P の電子と –P の電子が Cooper 対を作ると説明されることがあります。この説明をそのまま読むと、対を作る電子が実空間イメー
ジではどんどん離れ、対を作った状態を維持できないように思えます。ここで重要な

のは、一般的に Cooper 対も電子と正孔が対を作っていることで、上記の –P をもつ電子というは正孔に対応します。量子力学における粒子の速度は、群速度から 𝑣𝑣 =𝑑𝑑𝑑𝑑(𝑘𝑘)/𝑑𝑑𝑑𝑑 = ℏ−1 𝑑𝑑𝑑𝑑(𝑘𝑘)/𝑑𝑑𝑑𝑑 と導出され、これから vx = Px / m* という、古典力学に対応する自然な結果が出てきます。一方、この式からは、正孔を本来の「負の有効質量をもつ電子」と解釈すると、電子と正孔とは速度の方向が逆転します。実際、正孔は上述のように電子の空いた孔ですので、電子の移動方向と正孔の移動方向は逆になるので、これは直観にも一致する正しい結果です。つまり、+P の電子と –P の電子 (正孔) は同じ方向の速度を持つので、実空間イメージでも対が崩壊することはない、という自然な結論が出ます (半導体や絶縁体中の励起子が、運動量 P の電子と正孔が対を作っているのと同じ)。
以上のような疑問は、正孔という単純化したモデルによっておこる見かけ上の矛盾です。このようなことが生じたら、本来の物理モデルに立ち戻って考えることが必要です。以上の状況を図 7e-1 にまとめました。この 3, 4 は同じ電子状態なので、起こる物理現象も同じです。しかしながら、電子とみるか、正孔とみるかで、速度は逆向きになります。
付録 8: 磁気構造
磁性体中の磁気配列は、中性子回折により実験的に決定することができます。一方、第一原理計算からも予測が可能です。今までの結晶構造計算と同様にこの場合も、適当なスピン配列を仮定して構造緩和計算を行い、各スピン配列での最小エネルギーを求めます。その結果、一番エネルギーの低いスピン配置が基底状態であり、他のスピン配置とのエネルギー差 (∆E) がスピン－スピン相互作用の大きさに相当することになります。表 A8-1 に示すように、LaMnOP という層状結晶では、面内反強磁性配列が強磁性配列よりも 1.8 eV も安定ですが、面間のスピン配列によるエネルギー差はなく、面間のスピン間相互作用は無視できることがわかります 85。一方、BaMn2P2 のように面間隔が小さくなると、スピン相互作用も 80 meV ほどになり、AFM/FM 配列が安定でなります。
表 A8-1 Mn 化合物の磁性構造と結晶構造、Mn イオンあたりのスピンモーメントの計算結果。FM, AFM, NM はそれぞれ、強磁性配列、反強磁性配列、非磁性を示す。ス
ピン配置の a-b, c はそれぞれ、a-b 面内および c 軸方向での配列を表す。
a
b
c
化合物
スピン配列
∆E
µMn
a-b
/
c
式量 Z
LaMnOP FM/FM--0.4091
0.40910.82592.1
Z=2
AFM/AFM-1.7810.41100.41100.89403.9
AFM/FM-1.7820.41110.41110.8955
3.9
AFM/FM
0.40580.40580.8843
2.26
実測
BaMn2P2 FM/FM--0.4139
0.41391.254
3.5/3.1
Z=2
AFM/AFM-2.6490.40710.4071
1.302
3.9
AFM/FM-2.5670.38930.4072
1.309
3.9
AFM/AFM 0.4037
0.40371.30524.2
実測
付録 9a: 接合界面と輸送特性
多くの第一原理計算プログラムでは 3 次元周期をもつ結晶しか扱えないため、表面
85

T. Kamiya, H. Yanagia, T. Watanabe, M. Hirano and H. Hosono: Mater. Sci. Eng. B 173 (2010) 239.

や界面の計算をする場合は、一次元方向に大きな supercell を組み、界面や表面がこの方向に周期的に現れる構造について計算を行います (ATK や CRYSTAL 86 のようにスラブなどの計算ができるプログラムもあります)。
図 A9a-1(a,b) は、異なる構造の Au/ZnO/Au 界面についての緩和構造と静電ポテン
シャル分布を表しています 87。この結果から、図(a)では ZnO の両端が Au 電極と良い電気接触をしており、静電ポテンシャルが ZnO 中で均一になっていることがわかります。一方で図(b) の場合は静電ポテンシャルが大きく変化しており、左端電極界面の ZnO の Zn 面と電極の電気的接触が悪く浮遊電位に近い状態になっています。
ATK では非平衡グリーン関数法を用いることで、このような系の電子の透過率スペ
クトルと、電子散乱がない場合の電流－電圧特性を計算することができます。その結果を図 A9a-1(c)に示しますが、非オーミック特性が現れており、理想界面を持つ
Au/ZnO 接触は Schottky 接触であることが理論的に確認できます。
表面・界面の原子配列を実験的に決めることは難しく、また、決められたとしても、その配列が物性測定にかかっているマクロな範囲で理想的に保たれていることは稀です。そのため、理想的なモデルで計算した結果が測定結果にどれだけ対応するかは常に注意する必要があります。

図 A9a-1 Au/ZnO/Au 界面についての (a,b) 緩和構造と静電ポテンシャル分布(VASP,
PAW, PBE96)、および、(c) 非平衡グリーン関数法で計算した電流-電圧特性 (ATK/DFT)。
(a),(b)では界面における原子配置が異なる。それぞれについて下の図は、静電ポテン
シャルの等高面を表している。
付録 9b: 表面とイオン化ポテンシャル、仕事関数
表面と界面の計算モデルの違いは、表面の場合には原子を置かない空間を supercell
内に作ることだけです。図 A9b-1 は異なる表面を持つ SnO2 (Rutile 構造) の supercell
モデルで計算した静電ポテンシャル分布です。ここで界面の DFT 計算に問題がでます。付録 9a のような界面の計算では VBM のエネルギー差を直接局所状態密度などから求めることができますが、表面の場合は固体と真空のエネルギー準位 (真空準位,
Evac) との差が問題になります。しかしながら、バンドギャップ問題と同じ理由で、非占有状態である Evac は DFT では大きく過小評価されてしまいます。一方、静電ポテ
ンシャル Evac,pot は Poisson 方程式により電荷密度から正確に計算できます。しかも、
孤立自由電子の最低エネルギーは静電ポテンシャルに等しいため、Evac,pot を Evac とみ

8687

http://www.crystal.unito.it/
T. Kamiya, K. Tajima, K. Nomura, H. Yanagi and H. Hosono: Phys. Stat. Sol. A 205, 1929 (2008).

なすことができます 88,89,90。
SnO2 の (001)面は電荷中性であり、かつ、電気分極も持たないため、図 9b-1 にみられるように、真空中の静電ポテンシャル (Evac,pot) の計算結果は一定になります。仕事関数(φ) は定義から、φ = (Evac + EN-1) – EN で計算できます。ここで Evac は無限遠点での真空準位です。EN は系に N 個の電子がある場合の全エネルギーであり、電子数を変えることで EN-1 も計算できます。このようにして得られるφは、この場合はイオン化エネルギーに相当し、8.0 eV が得られます。この値は光電子分光による実測値 8.6
eV よりも 0.6 eV 小さいだけです。仕事関数を求める場合には、物質中のキャリア濃度で変わる分を考慮し、バルク中の計算結果と仕事関数の相対値を使って比較する必要があります。
図 A9b-1 の SnO2 の(100)面や(110)面の計算結果からわかるように、表面が帯電している場合、真空中にも電束が漏れ出し、電界が発生します。実際の測定でこのような状況が起こると、いわゆるチャージアップ効果と区別がつかず、また、測定装置の電位境界条件で結果が変わるため、well-defined な測定は困難です。また、実際の物質表面では環境中からの分子吸着、電子の取り込みや、バルク中の電子伝導によりこのような電束の漏れは抑えられているはずであり、このようなモデル計算に現実的な意味があるかは慎重に考える必要があります。プログラムによってはこのような"artefact" な電界の補正ができますが、その場合も補正の物理的意味と測定条件・結果にどう対応するのかを十分考慮しなければなりません
また、10 章 で、DFT では被占有準位を浅く過小評価することを説明しました。これはイオン化ポテンシャルや仕事関数の計算でも問題になり、特に半導体のイオン化
ポテンシャルを上記の方法で計算する場合も、混合汎関数法や GW 近似でないと実測値を再現できないことも多くあります 91。

図 A9b-1 SnO2/表面構造モデルの静電ポテンシャル分布。(001)表面については、真空領域の厚さを 5, 10, 20Å と変えて計算した結果を示している。振動しているのが SnO2
固体領域であり、原子がある位置かどうかで静電ポテンシャルが変わる。平坦あるいは直線的に変わっている領域が真空領域である。

C.J. Fall, N. Binggeli and A. Baldereschi: J. Phys.: Condens. Matter 11, 2689 (1999).
A. Baldereschi, S. Baroni and R. Resta: Phys. Rev. Lett. 61, 734 (1988).
90 H. Zhu, M. Aindow and R. Ramprasad: Phys. Rev. B 80, 201406 (2009).
91 Y. Hinuma, A. Grüneis, G. Kresse and F. Oba, Phys. Rev. B 90 (2014) 155405.
8889