平面波基底の第一原理計算法(第一回)
香

山

正

憲

で，平面波基底を用いた第一原理擬ポテンシャル法のハミル

.

は じ

め

に

トニアン，全エネルギー，電子密度分布，原子に働く力等の表式や演算の方法・手順を具体的に説明する(ここまで，第理論と計算手法の進歩，汎用コードの開発・普及により，二回，三回掲載分)．原子近傍のポテンシャルや波動関数の

物質・材料の研究開発で，第一原理計算が頻繁に用いられる取り扱いが鍵であり，また，フーリエ変換など逆空間表現をようになっている．結晶や表面・界面，合金や格子欠陥，ア

駆使して効率化を図る．第 7 章では，CarParrinello の第一モルファスなど，様々な系の安定構造やエネルギー，諸性質

原理分子動力学法に端を発する基底状態計算の「高速化技法」

を電子挙動まで掘り下げて探ることができる．実験観察と組について説明する(第四回掲載分)．

み合わせて物性発現の機構を解明したり，近年では機械学習やデータベースと組み合わせて，物質・材料の探索・予測・設計の可能性を広げるなど，第一原理計算は，材料科学・工学に広範なインパクトを与えている．材料研究者も，第一原
理計算の原理や汎用コードの計算の概要を理解し，精度や適用性，計算条件の知識を持つことが望ましい．

.


第一原理計算の基礎基本的近似と密度汎関数理論
第一原理計算とは

本章と次章では，「平面波基底の第一原理計算法」以外の手法も含めた，第一原理計算全体の基礎や方法論を概説す

筆者はこれまで，密度汎関数理論に基づく「平面波基底のる．物質・材料の安定構造や性質は，原子間結合の担い手で第一原理計算法」のプログラム開発や関連する計算手法開発

ある電子に支配される．第一原理計算とは，量子力学と統計を行い，それらを微視的観察と連携させ，エネルギー環境材

力学の理論のみにより，物質・材料中の電子や原子の挙動を料や構造材料の構造や性質の解明・設計に取り組んできた．

高精度に再現・予測する計算手法である．分子を扱う量子化第一原理計算では，電子の波動関数を何らかの「基底関数」学の理論や手法もあるが，結晶や固体には，固体物理で発展の線形結合で展開する．平面波基底の第一原理計算法は，基したバンド理論が用いられる(1)(4)．本講義では後者を扱う．

底関数に平面波を用いる手法であり，「第一原理擬ポテンシ

結晶のように周期性のある系(同じ構造単位(単位胞，ユニ

ャル法」と総称される手法である．原子間結合を担う価電子ットセル)が繰り返す系)の電子構造は，バンド理論に基づいのみを扱うので，全電子を扱う方法よりも負荷が少なく，巨て計算される(バンド計算という)．さらに電子構造に基づい大な固有値問題を新規アルゴリズムで効率的に解くので，比て，全エネルギー(凝集エネルギー)，安定原子配列，原子に較的大規模な構造が扱える． VASP, CASTEP, ABINIT な

働く力，応力，弾性定数，熱的・電子的・磁気的・光学的性

どの汎用ソフトは，全て，この手法を実行するものである．

質等，様々な物理量が計算できる．点欠陥・転位，表面・界

今回，平面波基底の第一原理計算法の原理，概要，計算技

面，アモルファスなどの諸構造や，拡散，偏析，変形・破術を講義する．大学院集中講義や勉強会で用いたノートを基

壊，吸着・反応といった諸現象も，大きなセル(スーパーセにしている．第 2 章，3 章で密度汎関数理論に基づく第一原ル)が繰り返すモデル構造を用いて扱える．結晶の単位胞

理計算の全体像を概説する(第一回掲載分)．第 4 章では，は，比較的単純な金属や化合物では通常 1～数原子，スーパ

第一原理擬ポテンシャル法の原理を説明し，第 5 章， 6 章

ーセルは数十～数百の原子を含む．計算時間は原子数(電子国立研究開発法人 産業技術総合研究所 エネルギー・環境領域 電池技術研究部門名誉リサーチャー(〒5638577 池田市緑が丘 18
31 産総研関西センター)
Lecture Notes on FirstPrinciples Methods Using a PlaneWave Basis Set (Part 1); Masanori Kohyama(Research Institute of Electrochemistry, Department of Energy and Environment, National Institute of Advanced Industrial Science andTechnology, Ikeda, Osaka)
Keywords: densityfunctional theory, KohnSham equation, Brillouin zone, Bloch's theorem, band structure, selfconsistent field
2021年12月 1 日受理[doi:10.2320/materia.61.588]




講義ノート

数)の 2 ～ 3 乗に比例して増えるので限界があるが，計算機

F ( x )が対応するが，汎関数は関数形 f ( x )に対し値 F [ f ]が能力と手法・アルゴリズムの進歩で扱える対象が拡大してい対応する．また，汎関数 F[ f ]は関数 f で微分できる(汎関数る．

微分dF[ f ]/df )．
全エネルギー Etot[ r]は次のように表される．

こうした第一原理計算は，基礎である密度汎関数理論の確
立，第一原理擬ポテンシャル法など高精度且つ高効率の手法

f

Etot[ r]＝T[ r]＋ V(r… )r(r… )dr… ＋の開発，基底状態計算の高速化技法の開発など， 1980 年代

f
f…… … …

)r(r )r(r ′dr… dr… ′||r－r ′e22＋Exc[ r]＋EI－I

から 90 年代に飛躍的に発展し，基幹部分が開発された．原(1)

子の動きのたび毎に大規模電子構造計算を行い，原子の挙動

第一項は電子系の運動エネルギー，第二項は核(イオン)電を高精度に扱う第一原理分子動力学法も実現している．材料

子間のポテンシャルエネルギー，第三項は電子間静電相互作の機械的性質に関する第一原理計算も急増している(5)(7)．

用エネルギー，第四項が複雑な多電子間相互作用を全て含む



交換相関エネルギー Exc，第五項は原子核間の静電相互作用

断熱近似と平均場近似

エネルギーである．第三項に 1 / 2 が付くのは，この積分形

原子核に比べて桁違いに小さな質量の電子は，原子核の動では二重に数えるからである．最初の四項が r の汎関数できに速やかに追随し，その時の原子配置に対し瞬時に最安定ある．平面波基底法では，価電子と価電子を取り去った正イ状態をとると見なせる．従って，通常の第一原理計算では，オンの集団を扱うので， r(r… )は価電子密度分布， V (r… )は正固定した原子配置の電子構造を計算する(断熱近似， Born 

イオンからのポテンシャル， EI－I は正イオン間静電相互作

Oppenheimer 近似)．多くの現象はこの立場で扱えるが，電用になる．

子系と原子系の間に直接のエネルギー移動が生じる現象(格KohnSham 方程式



子振動による電子の散乱など，電子格子相互作用)は扱えな密度汎関数理論から，全エネルギー汎関数 Etot [ r ]を最小

い．
一方，物質中の電子集団は，静電相互作用とパウリの原理にする r ( r… )が求まれば，電子構造を解明したことになる．

の元で複雑な多体相互作用により運動しており，単純な解析

それを求める方程式が KohnSham 方程式(9)である．r(r… )が的取り扱いは不可能である．通常，各電子が他の電子からの占有された一電子波動関数ψi(r… )の電子密度分布の重ね合わ

平均的ポテンシャル場の中で独立に運動しているという描像

せを用いる(平均場近似)．

occr(r… )＝2∑|ψi(r… )|2

パウリの原理からフェルミ粒子である電子の集団の全体の(2)i多電子波動関数は，スピンを含めた電子の交換に対しマイナで表されるとする(2 はスピンの和，i の和は占有状態のスになる(反対称)．この性質を満たす多電子波動関数の形の和)．運動エネルギーを一つがスレーター行列式で，ここから平均場近似として，多

f (

occ

T＝2∑ ψi(r… ) － 静電相互作用と◯ 交換相互作用を持って平均電子が互いに◯i)

2 2: ψi(r… )dr…
2m

(3)

場中を運動する描像が構築される(HartreeFock 近似)．こ

とすると

 ，◯
 以外の電子間のダイナミックな相互作用の段階では，◯
 相関相互作用(電子相関))は取り入れられない．次節の密(◯(－2m : は運動エネルギー演算子，:＝&&x ＋&&y ＋&&z )，

 ，◯
 ，◯
 の多電
度汎関数理論は，実質は平均場近似だが，◯

Etot を最小化するための KohnSham 方程式が導出される．

子間相互作用を全て取り入れる試みである．



2

2

2

2

2

2

2

2

[ 2m :＋V (r… )]ψ(r… )＝E ψ(r… )Hψi(r… )＝ －密度汎関数理論

2

2

eff

i

i

i

Veff(r… )＝V(r… )＋VH(r… )＋mxc(r… )多電子間相互作用の取り扱いが計算の精度を決定する．バンド計算では密度汎関数理論( density2functional

theory;DFT)(8)(9)を用いる．分子を扱う量子化学では多電子波動関

f… … …

VH(r… )＝e 2)r( r ′
dr… ′||r－r ′

mxc(r… )＝dExc[ r]/dr

(4)
(5)
(6)
(7)

数の組み立て(スレーター行列式の線形結合など)を工夫する

これは一電子波動関数のセット{ψi( r… )}を決める方程式であ理論が発達したが，密度汎関数理論は，多電子波動関数でなる．Veff は，原子核(イオン)からのポテンシャル V，電子かく電子密度分布 r ( r… )に着目する．これは，場所 r… に電子がらの静電ポテンシャル VH，複雑な多体相互作用を表す交換存在する確率，電子雲の密度分布である．

相関ポテンシャル mxc の和である． mxc は汎関数 Exc の関数

密度汎関数理論の結論は，(1)原子核やイオンからのポテ
ンシャル場 V ( r… )の下での電子集団の基底状態(絶対零度の微分で，それ自体も r の汎関数である．
この方程式は Lagrange の未定係数法から導出される．式

最安定状態)は，ただ一つの電子密度分布 r(r… )を持つ，( 2)( 1 )の汎関数 Etot は，式( 2 )から一電子波動関数のセット

基 底 状 態 の r ( r… ) は ， r ( r… ) の 汎 関 数 で あ る 全 エ ネ ル ギ ー{ψi}の汎関数 Etot [{ψi}]と見なせる． r についての最小化を

Etot [ r ]を最小にする，ということである．ここで汎関数と{ψi}の規格直交条件付きの最小化で実行する．そのためには，関数の関数という意味．通常の関数は，値 x に対し値

Qtot[{ψi}]＝Etot[{ψi}]－∑i jli j(＜ψi|ψj＞－di j)

ま て り あ 第61巻 第 9 号(2022)
Materia Japan

(8)




の 変 分 dQtot[{ψi}]/dψ
i ＝0 を 扱 う (注＜ψi|ψj＞＝∫ψi (r… )ているのは，式( 1 )の全エネルギー汎関数 Etot を最小化す

ψj(r… )d r… )．直交化は別過程で付加するとすれば

る電子密度分布 r が基底状態の正しい r と Etot を与えると

Qtot[{ψi}]＝Etot[{ψi}]－∑ili(＜ψi|ψi＞－1)(9)

いうことだけである．一電子波動関数のセット{ψi}は，(極で，変分は dEtot [{ψi}]/ dψ
i－ lψ
i i＝ 0 である．式( 1 )の汎関論すれば)式( 2 )を通じて式( 1 )の最小化を実行する手段

数微分 dEtot [{ψi}]/ dψi は，第二項～第四項について r での(パラメータ)である．KohnSham 方程式の導出時に用いた微分を先に行い，式( 2 ), ( 3 )から KohnSham 方程式が導運動エネルギーの一電子波動関数による表現式( 3 )も，T＋出される．

Exc に封じ込められた多体相互作用を扱う際の一つの選択でV(r… )r(r… )dr…
{f)er(r )r(r ′＋f|r…… －r… ′… | dr… dr… ′＋E [r]}/dr･dr/dψ2f＝ －[ 2m :＋V(r… )＋e f|rr… －(r… r′… )′| dr… ′＋dE [r]/dr]ψ

dEtot[{ψi}]/dψ
i ＝dT[{ψi}]/dψ
i ＋d
2

xc

i

2

2

2

xc

＝Hψi

ある．つまり，KohnSham 方程式の一電子波動関数ψi とその準位 Ei の物理的な意味が明確でないという弱点がある．
一電子波動関数とその準位という概念自体がモデルであるためである．固有値 Ei は最高占有準位が，系からその電子を取り去るエネルギーに対応することが知られているが，他のi準位の意味は明確でない．固有値 Ei は，観測されるバンド(10)構造(多電子系のスペクトル)を少なくとも定性的には再現す式( 2 ), ( 3 )はスピンを区別して扱わないので 2 が係数についている．式( 10 )にも係数 2 がつくかと思われるが，一つのスピンの状態についての汎関数微分なので係数 2 は無くる．


局所密度近似と密度勾配近似てよい．ここで dEtot [{ψi}]/ dψ
i＝ Hψi の関係は重要である

密度汎関数理論の最大の問題点は， r ( r… )の汎関数である(第 7 章で用いる)．式( 1 )の第五項 EI－I は r を含まないの

交換相関エネルギー Exc [ r ]とその汎関数微分である交換相で Kohn Sham 方程式には関わらない．なお，ψi でなくψi

関ポテンシャル mxc[r]の厳密な形がわからないことである．

についての汎関数微分を行っても，ψi についての同様の方

一方，一様電子ガス(一様な正電荷バックグラウンド下の自程式になり同じである(第 6 章(5)節で扱う)．

由電子集団)については，理論や数値解析から関数形がわか式( 5 )の Veff は系の全電子に共通である．静電ポテンシっている．そこで，r(r… )の空間変動が小さいとして，r… 毎に

ャル VH は，全電子からの寄与で，各状態ψi の自分自身か

その地点の電子密度を一様電子ガスの表式に代入して用いるらの寄与(自己相互作用)も含むが， mxc 内の自己相互作用で簡便法が局所密度近似(local density approximation; LDA)で打ち消される．電子相関の強い系では，この点がうまくいかある(10)．exc(r… )を交換相関エネルギー密度としてず精度を悪くする(文献(3)の第 9 章に詳しい)．一方，電子
スピンについて区別せずに扱ってきたが，スピン分極する系
等では，up と down の各スピンの波動関数ψsi (s＝↑or↓)でExc[r]＝ exc(r… )r(r… )dr…

f

(12)

mxc(r… )＝dExc[r]/dr＝exc(r… )＋dexc/dr･r(r… )(13)

別々に KohnSham 方程式が組み立てられ，交換相関ポテンとなり，実空間メッシュ点 r… m 毎に r(r… m)の値から exc(r… m )，シャル mxc が msxc でスピン毎に異なる．

mxc(r… m)の値を計算して与える．

式( 4 )はシュレディンガー方程式の形 Hψi＝Eψ
i i をしてお実際には r ( r… )の空間変動は小さくない．そこで，電子密り，固有値 Ei ，固有状態ψi を求める問題となる．異なる固度に加えて密度勾配 :r(r… )(厳密には|:r(r… )|)も取り入れて有状態の波動関数は互いに直交化している(∫ψi(r… )ψj(r… )dr… ＝Exc, mxc を表す一般化密度勾配近似(generalized gradient ap-di j)．Veff は r の寄与({ψi(r… )}の寄与)を含むので，方程式はproximation; GGA)(11)(12)も提案されている．

入力と出力の r ( r… )が一致するように解く．これを自己無撞局所密度近似を用いた計算は，物質の凝集エネルギーを高着(selfconsistent field; SCF)計算という．最終解の{ψi(r… )}めに，ボンド長を短めに再現する傾向がある(実験値との誤からの r ( r… )が Etot を最小にする基底状態の電子密度分布で差 2 ～ 5 内外)．密度勾配近似ではかなり改善される． Fe

ある．

の最安定構造が強磁性の bcc 構造であることも密度勾配近ここで，式( 4 )の両辺に左からψi を作用させて積分し，

似で初めて再現される(文献(1)の第 5 章)．しかし，密度汎占有状態の和をとると，右辺は 2∑occ
Ei となり( 2 はスピンi

関数理論は基底状態を扱う理論であり，励起状態に関わる物の和)，左辺の各項は r やψ
i をかけた積分となり，式( 1 )と性値(光学的性質等)の精度は良くない(LDA, GGA 共にバンドギャップ値は実験値の 1/2～2/3)．また，フォンデルワーの比較から次式が導出される．
occ

Etot[r]＝2∑Ei－
i

1
2

ルス相互作用は動的な電子相関に起因し，従来の LDA ，

f…

VH(r )r(r… )dr… ＋Exc[r]f－ mxc(r… )r(r… )dr… ＋EI－I

GGA では扱えなかったが，最近，密度汎関数理論の枠内で(11)扱う手法が開発されている(13)．一方，強相関系といわれる物質群(例えば NiO など)について，密度汎関数理論では，

右辺第二項は，式( 6 )からの∫VH(r… )r(r… )dr… が式( 1 )の第三

金属か非金属かなど基本的な電子状態の再現性も難しい場合

項の 2 倍になっているためである．

がある．密度汎関数理論の限界を超える理論や手法(量子モ

ところで，密度汎関数理論で厳密に正しいことが証明され

ンテカルロ法等(14))も研究されている．

講義ノート

.




第一原理計算の基礎周期的ポテンシャル場における固有値・固有関数
格子と逆格子結晶やスーパーセルの周期系では格子点毎に同じ単位胞が繰り返す(図)．格子点(単位胞)の位置を示す格子ベクトルR• は，基本並進ベクトル a… 1, a… 2, a… 3 を用いて R• ＝ l1a… 1 ＋ l2a… 2 ＋
l3a… 3(l1, l2, l3 は整数)である．単位胞は，a… 1, a… 2, a… 3 を稜とする平行六面体を取ることができ，原子位置は単位胞内相対位置ベクトル t• a を用いて R• ＋t• a である．
結晶やスーパーセルでは，式( 5 )のポテンシャル Veff( r… )が格子周期性(並進対称性)を持ち，任意の R• に対しVeff(r… ＋R• )＝Veff(r… )

(14)

図1

結晶の単位胞の並びとボルンフォンカルマンの周期境界条件．単位胞(体積 Qc )が基本並進ベクトル …a1, …a2, …a3の周期で繰り返す．一方，三方向に N1 …a1, N2 …a2, N3 …a3 のサイズの体積 Q の結晶部分が，さらに外側に繰り返す
ボルンフォンカルマンの境界条件を考える．

図2

ブリルアンゾーンと既約領域(正方晶の例)．黒丸が G•
点．図に垂直方向が kz 軸，ブリルアンゾーンは kx － ky
面が正方形，kz 軸に沿って正と負の正四角柱． •k1～ •k7 点
は，既約領域の •k 点を系の対称操作で回転した S •k 点．
既約領域は kz が正の領域の三角柱( 1 / 16 の部分)．反転
対称のある正方晶( D4h など)では対称要素は 16 個， S •k
により既約領域がブリルアンゾーン全体を埋める．反転対称のない場合(C4v など)，対称要素数は 8 個，時間反転対称で－ •k や－ S •k も同じ固有値，既約領域は同じ(1/16)．

となる．ポテンシャルは単位胞毎に同じものが繰り返す．
一方，基本逆格子ベクトルが，実空間の格子に対応して，b• 1＝2pa… 2×a… 3/a… 1･a… 2×a… 3, b• 2＝2pa… 3×a… 1/a… 1･a… 2×a… 3,b• 3＝2pa… 1×a… 2/a… 1･a… 2×a… 3(15)で定義され，a… i･ •bj＝2pdi j を満たす．逆格子点(逆格子ベクトル)が G• ＝ m1b• 1＋ m2b• 2＋ m3b• 3( m1, m2, m3 は整数)のように周期的に繰り返す逆格子空間(k• 空間)が設定できる．格子ベク
トルとの内積がG• ･R• ＝2pM(M は整数)

(16)である．実格子の単位胞に対応して，逆格子空間の単位胞がブリルアンゾーン(Brillouin zone)である(図)．b• 1, b• 2, b• 3 を稜とする平行六面体でもよいが，逆格子点を中心に近接する
逆格子点へのベクトルの垂直二等分面で区切られた同じ体積の領域が取られる．実空間の単位胞の体積は a… 1, a… 2, a… 3 による平行六面体の体積 Qc ＝| a… 1 ･ a… 2 × a… 3 |で，逆格子空間のブリル ア ン ゾ ー ン の 体 積 は b• 1, b• 2, b• 3 に よ る 平 行 六 面 体 の 体 積1である．単位胞の体積 Qc が大きいほ|b• 1･b• 2×b• 3|＝(2p)3Q－c

どブリルアンゾーンの体積は小さくなる．


ブロッホの定理，k• 点，ブリルアンゾーン内積分

格子ベクトル R• ，逆格子ベクトル G• で特徴付けられる周期系(並進対称性のある系)において KohnSham 方程式を解く計算が，バンド計算である．KohnSham 方程式は 2 次の偏微分方程式(波動方程式)で，解は波動の性質を持つ．ブロッホの定理( Bloch's theorem )から，周期系の固有関数は，以下の性質を持つ(厳密には群論を用いて証明される(15)(16))．固有関数は，必ず波数ベクトル k• を持つ波(進行方向 k• ，波長 2p /|k• |)であり， k• はブリルアンゾーン内に限

2|Ukn• (r… )| であり，式(18)から，その分布も格子周期関数でられる．さらに固有関数は以下の形を持つ(図)．

ある．

ψkn
• (r… )＝Ukn• (r… )exp[ik• ･r… ]

(17)

Ukn
• (r… )は格子周期関数で，任意の R• の並進についてUkn• (r… ＋R• )＝Ukn• (r… )本稿でのブリルアンゾーンは，原点を中心とした第一ブリ
ルアンゾーンのことである．図 2 に示すように，第一ブリ(18)を持つ固有状態，例えば，理ルアンゾーンから外に出た k• ′を満たす．n はバンド指標である．単位胞当たりの電子数に想的な自由電子として，波動関数が exp [ ik• ′･ r… ]の電子状態

応じて k• 点毎に数個の固有状態を求めるが，エネルギーの低＝k• ＋G• として，を考えよう．何らかの逆格子ベクトル G• で k• ′い順に番号 n を付ける．式( 17 )から，固有状態ψkn• ( r… )の電k• が第一ブリルアンゾーンに入るようにできる．波動関数は2子密度分布の式( 2 )への寄与|ψkn• ( r… )| は， exp 項が消えて･ r… ]＝ exp [ iG• ･ r… ] exp [ ik• ･ r… ]＝ UG• ( r… ) exp [ ik• ･ r… ]で，式exp [ ik• ′ま て り あ 第61巻 第 9 号(2022)
Materia Japanれらの和として実行する．例えば，式( 2 )の電子密度分布は，
occr(r… )＝∑2
n

f

Qc
(2p)3

BZ

occ

2
fkn
• |ψkn• (r… )| dk• ＝∑∑wk• i n|ψ•ki n(r… )|2
n

•ki

occ

＝∑∑wk• in|Uk• in(r… )|2
n

(19)

•ki

で与えられる．一行目の表現は占有されたバンド毎のブリル
アンゾーン( BZ )内の占有状態部分の積分で， 2 はスピン，
図3

結晶中の固有関数の概念図．(a)格子の周期ポテンシャ
ル V ( r… )のもとでは，ブロッホの定理から，( b )格子周
期関数 u •k(r… )と(c)平面波 ei •k･r… の積の形を(d)電子の固有関数ψ•k(r… ) が有する．u •k(r… )や ei •k･r… , ψ•k(r… )は複素数であるが，簡単のため実数成分のみを模式的に示す．

Qc/(2p )3 はブリルアンゾーン体積で割る規格化因子， fkn• は

バンド毎の k• 点について占有部分で 1，非占有で 0 に取る．
金属の場合，同じ n でもブリルアンゾーン内で k• 点によって，占有状態と非占有状態があり得る．式( 19 )の一行目の右側と二行目は，ブリルアンゾーン内積分をメッシュの k• i
点の値の重み付きの和で実行する．スピンや占有，非占有含( 17 )のように第一ブリルアンゾーン内の k• の固有状態と見なせる．exp[ik• ･r… ]の状態より上のバンドの状態(n＞1)になる．

めて，重み wk• in で調整する．
積分の精度と効率の観点から，ブリルアンゾーン内の k• 点のサンプル法，メッシュ法が提案されている(17) ．金属の場

格子ベクトルから逆格子ベクトルやブリルアンゾーンを組

合，フェルミ面近傍の占有部分と非占有部分の区別のため，

み立てる理由は，周期系の電子の固有状態は k• で識別され，細かな k• 点メッシュ点が必要になる．限られた k• 点数で積分

ブリルアンゾーン内の全 k• 点で固有値，固有関数を求めれ精度を上げる方法として，フェルミ面付近で Gaussian 等でば，電子構造が完全に解けたことになるからである．式

ボカした占有率を導入する方法等がある(18)(19)．

( 14 ) , ( 17 ) , ( 18 ) を 式 ( 4 ) と 組 み 合 わ せ る と ， 周 期 系 のKohn Sham 方程式は k• 点に依存した形となり， k• 点毎に固
系の対称性とブリルアンゾーンの既約領域

有値 Ekn• ，固有関数ψkn• を求める．求まった Ekn• ，ψkn• は，バメッシュの k• 点を抽出する領域は，ブリルアンゾーン全域

ンド指標 n を固定し k• を変化させると少しずつ変化し，k• にではなく，時間反転対称性と結晶系の持つ対称要素{S|t• S}かついての連続関数とみなせる．ブリルアンゾーン内の連続なら，既約領域( irreducible part )に絞ることができる(15) (図固有状態の一集団を「バンド」と呼ぶ(複数のバンドが重な

2)．既約領域の k• 点について固有値，固有関数を求めれば，る場合， n の値の入れ変わりがある)．ブリルアンゾーン内

Sk• や－ k• の固有値や固有関数との関係から，ブリルアンゾの k• に応じた Ekn• の変化の様子をバンド構造(バンド分散)とーン全体から k• 点を抽出したことになるのである．

いう．

既約領域とは，その領域の k• 点を，系の全ての対称要素のブリルアンゾーン内には，稠密に k• 点が存在する．周期系

回転行列 S で回転した Sk• で，ブリルアンゾーンの全体(反全体の単位胞(格子)の数 N とブリルアンゾーン内の k• 点の転対称性を持つ系)または半分(反転対称性を持たない系)が数は理念上，同じである．図 1 のマクロの周期境界条件(ボ埋め尽くせる領域である．例えば，対称要素が48個ある dia-

ルン・フォンカルマン(Bornvon Karman)の周期条件)で，mond 構造では，既約領域はブリルアンゾーン全体の 1 / 48

単位胞が a… 1, a… 2, a… 3 方向に各々 N1a… 1, N2a… 2, N3a… 3 のサイズでの部分である．反転対称のない zincblende 構造では，対称繰り返す体積 Q＝NQc(N＝N1N2N3，N は巨視的な数)の結晶

要素は24個だが，既約領域は同じく 1/48である．

部分を考え，そこでの並進自由度の総数 N に対応して k• 点

まず，時間反転対称から，－ k• の固有値 E－kn
• は k
• の固有

が定義されるからである(16) ．従って，ブリルアンゾーン

値 Ekn• と同じ値で，固有関数(波動関数)も互いに複素共役のは ， バ ン ド 毎 に ス ピ ン 含 め て 2N の 状 態 ( 電 子 ) を 収 容 す関係がある(15) ．従って，必要な領域は，ブリルアンゾーン

る(16)．単位胞当たりでは 2 電子である．bccNa 金属の価電の半分に絞られる．上記の diamond 構造と zinc blende 構子バンドを考えると，bcc 単位胞に 1 原子，価電子 1 個なの造の比較で，反転対称のない系でも反転対称のある系と同じで， bcc のブリルアンゾーンは半分だけ埋まる． Diamond

既約領域になる理由は，この時間反転対称の効果である(た構造の Si は，単位胞 2 原子，価電子が 8 個， diamond 構造

だし，スピンが関わると成り立たない場合もある)．

のブリルアンゾーンでは，重なった 4 バンド分が全て埋ま

次に対称要素の効果を考える．{ S | t• S }は，行列 S で回転る．単位胞当たり 2 電子が常に一つのバンド(ブリルアンゾ(鏡映含む)した後， t• S の並進を施す操作(t• S は基本並進ベクーン全体)を埋める(スピン分極した場合は 1 電子が 1 バン

トルより小さい並進で「らせん」「映進」の場合にのみ存在)

ド)．

である(15)．結晶が{S|t• S}を対称要素として持つ(その操作でさて，式( 2 )，( 3 )，( 11 )などの占有状態についての和結晶が不変)なら，Sk• での固有値が k• での固有値と同じ値で，は，ブリルアンゾーン内の k• 点の占有部分の積分になる．実Sk• での固有関数は k• の固有関数を{ S | t• S }で移したものにな

際には，メッシュ点など離散的な k• 点で Ekn• , ψkn• を求め，そ

る．




講義ノート

こ の 証 明 は 以 下 で あ る ． 実 空 間 で の 対 称 操 作 は r… ′＝が関わる原子位置 R• a ＝ R• ＋ t• a について， R• a′＝{ S | t• S} R• a＝ SR•{ S | t• S} r… ＝ Sr… ＋ t• S である．関数 f( r… )に作用させると，括弧内＋( St• a ＋ t• S )＝ R• ″＋ t• a ′のように，原子座標の調整にのみ t• S がの r… には逆変換で

関わる．
(20)

以上のように， k• での固有値 Ekn• ，固有関数ψkn• ( r… )が求まと表される．k• での固有値，固有関数の KohnSham 方程式

•

…

•

…

…

•

れ ば ， Sk• で の 固 有 値 は 同 じ 値 で ， Sk• で の 固 有 関 数 は ，ψkn• ( r… )を実空間で対称操作したものになる． k• と Sk• の固有{S|tS}f(r )＝f({S|tS}－1r )＝f(S－1(r－tS))の両辺に対称操作を作用させると{S|t• S}Hψkn• ＝{S|t• S}Ekn
• ψkn
• ,
H{S|t• S}ψkn• ＝Ekn• {S|t• S}ψkn•

関 数 の ノ ル ム の 二 乗 ( 電 子 密 度 ) の 分 布 は ， 式 ( 22 ) か ら(21)－1 ( － ))|2 で，後者は前者を実空間で
2|Ukn
r… t• S• ( r… )| と| Ukn• (Sとなる．ハミルトニアンと対称操作が可換であるためであ{S|t• S}の操作で移したもの．ψSkn• (r… )を求めなくても，ψkn• ( r… )る．これから，波動関数{S|t• S}ψkn• は， k• の固有状態ψkn• と同の電子分布の対称操作でその寄与を得ることができる．なじ固有エネルギー Ekn• を持つ固有状態であると言える．

お，対称性(空間群)の議論は文献(15)を参照のこと．

一方，{S|t• S}ψkn• は，式(20)と式(17)を用いて



{S|t• S}ψkn• (r… )＝{S|t• S}Ukn• (r… )exp[ik• ･r… ]
－1(r－t ))exp[ik･S－1(r－t )]
＝Ukn
• (S
… •S
•
… •S
－1(r－t ))exp[iSk• ･r… ]
＝exp[－iSk• ･t• S]Ukn
• (S
… •S

第一原理計算の手順

第一原理計算は，以下のような手順となる(図)．
(22)となる．内積の関係式 k• ･S－1r… ＝Sk• ･r… を用いている．式(22)－1 ( r － t ) ) の 部 分 に つ い て ，の exp [ － iSk• ･ t• S ] Ukn
• (S
… •S
exp[－iSk• ･t• S]は定数なので無視し，残りの部分は

 実格子から逆格子，ブリルアンゾーンを組み立て，規約領◯域と k• 点メッシュ，各種計算条件を決定する．
 単位胞内初期原子座標の設定，初期の input の電子密度分◯

布 rin(r… )の設定(自由原子の電子密度分布の重ね合わせ等)．

－1(r＋R－t ))＝U (S－1(r－t )＋S－1R )
Ukn
• (S
•
… • •S… •S•kn

 原子配列と rin(r… )を式( 4 )～( 7 )に入れて，ポテンシャル◯－1(r－t )＋R′
＝Ukn
(23)
• (S
… • S • )＝Ukn• (S－1(r… －t• S))－1 ( r － t ))が格子周期関数と言となり， exp[－ iSk• ･tS]Ukn
• (S
… •S

Veff(r… )＝V(r… )＋VH(r… )＋mxc(r… )，ハミルトニアンを組み立て••る． k• 点毎に Kohn Sham 方程式を解き，固有値 Ekn• ，固有える． R ′
＝ S－1R が同じ格子の格子ベクトルであること(回関数ψkn• を求める(この部分の詳細は次章以降で論じる)．

転 S や S－1 で同じ格子系に移る)，Ukn• ( r… )が格子周期関数で 求めた全 k• 点の全固有状態をエネルギーの低い順に並べ，◯あることを使っている．こうして，{S|t• S}ψkn• ( r… )はブロッホの定理を満たす Sk• の固有状態ψSkn• ( r… )と言える(複素数の位

フェルミ準位を決める．各状態の占有率を決める．

相の不定性はある)．

ンゾーン内積分から，output の電子密度分布 rout(r… )を計算．

2
 占有された固有関数の電子密度分布|ψkn◯• ( r… )| のブリルア

なお，系の対称要素{ S | t• S }について， k• と Sk• の関係には

 rin(r… )と rout(r… )を比較，前回の rin(r… )に rout(r… )との差の成◯並進成分 t• S は関わらない．また，上記で出てきた格子ベク に戻分を加え(chargemixing )，次回の rin ( r… )を構築し，◯トルの対称操作にも t• S は直接には関わらない．内部座標 t• aる．

図4

第一原理計算の手順．SCF ループと構造緩和ループ．右下の二重四角部分が狭義のバンド計算で，KohnSham 方程式の固有値・固有関数を •k 点毎に求める．最近のアルゴリズムでは，入力した初期波動関数(破線の四角部分)の繰り返し法による最適化計算として解く．

ま て り あ 第61巻 第 9 号(2022)
Materia Japan 以上の過程を rin ( r… )と rout ( r… )が一致するまで繰り返す◯( SCF ループ)．収束条件は，例えば，単位胞内全メッシュ
点 r… m で|rout(r… m)－rin(r… m)|D(24)のように入出力差が閾値 D 以下になること．実空間メッシ
 節)．
ュについては後述(第 6 章

簡単な chargemixing 法としては，例えば＋1N
N(r… m)＝rN
rN
in
in(r… m)＋a(rout(r… m)－rin(r… m))

(25)

＋1
を前回の rNのように，次回の rNin
in に入出力差の a 倍( a ＜1)を r… m 毎に加え作成する．金属では振動が起きやすいので

a を小さくする．後述(第 7 章)のように r のフーリエ成分を扱ったり，振動を防ぐ工夫が提案されている(20)．
SCF 計算で{ψkn• ( r… )}と r ( r… )が求まれば，式( 1 )から全エ 節)や応力(21) がネルギーが決まる．原子に働く力(第 6 章計算できる．力に従って原子変位を与える過程を繰り返せば安定原子配列が求まる．応力に従って格子ベクトルを変えて計算を繰り返せば，格子定数や単位胞(格子)の形が最適化できる．


バンド構造図と状態密度バンド構造図や状態密度の計算は， SCF 計算が終了した
後，確定した r(r… )による Veff(r… )のもと，様々な k• 点についての Kohn Sham 方程式の固有値計算で行う． k• 点は SCF計算で用いた k• 点とは異なる．バンド構造図では，特別の対称性を持つブリルアンゾーン内のエッジや頂点の k• 点で計算

( 4 ) R.M. マ ー チ ン  物 質 の 電 子 状 態 ( 上 ， 下 ) ， 丸 善 出 版 ，(2012).
( 5 ) S. Ogata, Y. Umeno and M. Kohyama: Model. Simul. Mater.
Sci. Eng., 17(2009), 013001(133).
( 6 ) J. Pokluda, M. Cerny, M. ½Sob and Y. Umeno: Prog. Mater.
Sci., 73(2015), 127158.
( 7 ) R. F. Zhang, S. H. Zhang, Y. Q. Guo, Z. H. Fu, D. Legut, T. C.
Germann and S. Veprek: Phys. Rep., 826(2019), 149.
( 8 ) P. Hohenberg and W. Kohn: Phys. Rev., 136(1964), B864
B871.
( 9 ) W. Kohn and L. S. Sham: Phys. Rev., 140(1965), A1133
A1138.
(10) J. P. Perdew and A. Zunger: Phys. Rev. B, 23(1981), 5048
5079.
(11) J. P. Perdew, K. Burke and M. Ernzerhof: Phys. Rev. Lett., 77(1996), 38653868.
(12) J. A. While and D. M. Bird: Phys. Rev. B, 50(1994), 4954
4957.
(13) I. Hamada: Phys. Rev. B, 89(2014), 121103(14).
(14) W. M. C. Foulkes, L. Mitas, R. J. Needs and G. Rajagopal:
Rev. Mod. Phys., 73(2001), 3383.
(15) G. バーンズ物性物理学のための群論入門，培風館，(1983).
(16) 香山正憲まてりあ，60(2021), 717721.
(17) H. J. Monkhorst and J. D. Pack: Phys. Rev. B, 13(1976),
51885192.
(18) C.
L. Fu and K.M. Ho: Phys. Rev. B, 28(1983), 54805486.
(19) M. Mathfessel and T. Paxton: Phys. Rev. B, 40 (1989), 3616
3621.
(20) G. P. Kerker: Phys. Rev. B, 23(1981), 30823084.
(21) O. H. Nielsen and R. M. Martin: Phys. Rev. B, 32(1985),
37803791.
(22) G. Lehmann and M. Taut: Phys. Stat. Sol. B, 54(1972), 469477.
(次号へつづく)され， Ekn• の変化の様子(バンドの分散)を探る．状態密度★★★★★★★★★★★★★★★★★★★★★

は，既約領域内の緻密な均一メッシュの k• 点で固有値計算を

1985年 東京大学大学院工学系研究科(金属材料学専攻)博士課程中退
同年 工業技術院 大阪工業技術試験所(現産業技術総合研究所関西センター)入所
2004年 産業技術総合研究所 ユビキタスエネルギー研究部門 グループ長
2015年 産業技術総合研究所 エネルギー・環境領域
電池技術研究部門 首席研究員
2021年～現職
専門計算材料科学，粒界・界面・ナノ構造の材料科
学，新規手法・コードの開発
Modelling and Simulation in Materials Science andEngineering 誌 Editorial Board memberE
mail: mkohyama＠aist.go.jp

行い，エネルギー軸でのヒストグラム(各エネルギー刻みに存在する固有状態の数)から求める．内挿補間から連続曲線の状態密度を求める方法もある(四面体(tetrahedron)法(22))．

文献( 1 ) 金森順次郎，米沢富美子，川村 清，寺倉清之固体構造と物性，岩波書店，(1994).
( 2 ) 小口多美夫バンド理論，内田老鶴圃，(1999).
( 3 ) 藤原毅夫固体電子構造論，内田老鶴圃，(2015).




香山正憲

★★★★★★★★★★★★★★★★★★★★★

講義ノート