講義スライドファイル更新版
http://conf.msl.titech.ac.jp/D2MatE/2023Tutorial/tutorial2023.html

チュートリアル

どのように第一原理バンド計算の条件を決めるか
神谷利夫
東京工業大学
国際先駆研究機構
元素戦略MDX研究センター

参考文献

参考Web: http://conf.msl.titech.ac.jp/Lecture/jsap-crystal/

参考文献

バンド計算法全般に関する参考文献
初心者、材料学者向け
材料電子論入門 第一原理計算の材料科学への応用
田中功、松永克志、大場史康、世古敦人 共著、内田老鶴圃 (2017).

量子計算の実際
密度汎関数理論入門 理論とその応用
佐々木泰造、末原茂共訳、 吉岡書店 (2014).

量子計算の物理的基礎
固体電子構造論 密度汎関数理論から電子相関まで
藤原毅夫著、内田老鶴圃 (2015)
物質の電子状態
R.M. マーチン著、寺倉清之、寺倉郁子、善甫康成訳、Springer Japan (2010).
密度汎関数法の基礎
常田貴夫著、 講談社 (2012)

特定の計算法・プログラムの参考文献
・ 第一原理シミュレータ入門 －PHASE & CIAO－
山本 武範、濱田 智之、山崎 隆浩、 岡本 政邦
アドバンスソフト発行、2004年初版
・ 固体の中の電子 WIEN2k 入門追加版
和光システム研究所、2006
・ (Xα – APW) スレーター分子軌道計算
菅野暁、足立裕彦、塚田捷、東京大学出版会 1982
・ (L/APW) Planewaves, pseudopotentials, and the LAPW Method
Ed. David J. Singh, Lars Nordstrom, Springer, 2006

・ (CRYSTAL) Hartree-Fock ab inito treatment of crystalline solids
C. Pisani, R. Dovesi, C. Roetti, Springer, 1988
・ The LMTO Method
H.L. Skriver, Springer, 1984
・ (Tight-Binding) 固体の電子構造と物性
W.A. ハリソン、現代工学社、1980

第一原理計算を始める方へ
無料プログラム• Quantum ESPRESSO (QE): https://www.quantum-espresso.org/
・ Linux版のみ配布。自力でコンパイル
・ GUI:PWgui http://www-k3.ijs.si/kokalj/pwgui/
最後のリリースが 2021/12。最近のLinuxで動くかどうか不明
・ その他関連ソフトウェア: https://www.quantumespresso.org/auxiliary-software/• ABINIT• Phase
有償プログラム• VASP: https://www.vasp.at/
アカデミック版 €4,000- ？
GUI: GoVASP, MedeAなど。高価

第一原理計算を始める方へ: Winmostar + QE
無償 (学生) ～ 60,000円 (教育機関向け) ～ 180,000円 (民間企業・官公庁)• Winmostar + Quantum ESPRESSO: Windows
https://winmostar.com/jp/
Quantum ESPRESSOのWindows版 (cygwin版) が入手可能
学生は在学中に限り無料の学生版を使用可能
https://winmostar.com/jp/editions/困ったときの検索
VASPの場合: VASP + キーワードで検索 => Vaspwiki が引っ掛かることが多い
理論的なこと: VASP + Theory で検索第一原理計算で何が計算できるか

第一原理計算とは
広義：
経験的パラメータを用いずに、物理の基礎方程式から望まれる計算結果を出力する
量子計算:
原子配列のみの入力から、量子方程式に基づいて精度の高い全エネルギーを出力する
 半経験法 (Tight-binding法)

全エネルギーから何が計算できるか
第一原理計算: 全エネルギー E を精度良く計算できる=> 原理的に、すべての物性が計算できる
・ 電子準位 (バンド構造): 𝑒𝑖 (𝒌) = 𝐸 𝑛𝒌,𝑖 − 𝐸(𝑛𝒌,𝑖 − 1)
(実際には一電子方程式の解)・ 安定構造:Etot が最小になる格子定数、原子座標を求める
13~14章・ 弾性率テンソル
13章, 付録2f~2g
1

𝑈 = 𝑈0 + 2 σ𝑖,𝑗,𝑘,𝑙 𝐶𝑖𝑗𝑘𝑙 𝑒𝑖𝑗 𝑒𝑘𝑙
𝜎𝑖𝑗 = σ𝑘,𝑙 𝐶𝑖𝑗𝑘𝑙 𝑒𝑘𝑙
・ 誘電率テンソル

歪 𝑒𝑖𝑗 を与えて U(𝑒𝑖𝑗 ) を計算
歪 𝑒𝑖𝑗 を与えて応力𝜎𝑖𝑗 を計算

付録2e~2g

𝑈 = 𝑈0 + 2 σ𝑖,𝑗 𝜀𝑖𝑗 𝐸𝑖 𝐸𝑗

1

電場 𝐸𝑖 を与えてU(𝐸𝑖 )を計算

𝐷𝑖 = 𝜀0 + 𝑃𝑖 = σ𝑗 𝜀𝑖𝑗 𝐸𝑗

分極 𝑃𝑖 はBerry位相から計算

エネルギー関数 U:
0 K, 定積: 内部エネルギー Etot (DFTの全エネルギー)
0 K, 定圧: エンタルピー H = U + PV
>0 K, 定積: Helmholzの自由エネルギー
F=U+
Felectron+
Fphonon

電子の状態密度から計算

フォノン状態密度
から計算

付録2k

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

エネルギー準位
E依存性
第1B.Z.で総和

k依存性

𝟏 𝝏𝟐 𝑬 𝐤
∗
𝒎 = 𝟐
ℏ 𝝏𝒌𝟐

有効質量

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

状態密度

バンド構造

全エネルギー

波動関数

光学スペクトル

Berry位相

|ψ|2

電気分極

電子分布

−𝟏
∞

𝒏𝒆 = න 𝒇𝒆 𝑬 𝑫𝒆 𝑬 𝒅𝑬

cicj<φｊ|H|φi>𝑬𝑪

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

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

磁気物性、半導体物性: 移動度、キャリア濃度、電子伝導度、Hall係数、熱電係数

第一原理計算による物性計算

第一原理計算全エネルギー U
自由エネルギーH = U – PV F = U + TS G = U – PV + TS

波動関数
Berry位相電気分極

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

誘電率
圧電定数

安定構造

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

弾性定数

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

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

反応エネルギー

欠陥生成
エネルギー

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

必要な計算精度は？

通常、計算精度はコストによって制限される
・ 計算時間・ 計算機費用・ 記憶容量

典型的な計算時間 (VASP、12core程度のPC 1台):
・ 100原子/単位格子、対称性が高い場合
電子構造: 1時間程度
構造緩和: 数時間程度
全部で1日程度の計算であれば、十分に高い精度で計算したっていいだろう通常、計算精度はコストによって制限される
欠陥、表面、接合などの計算
・ 300原子/単位格子、対称性が低い場合
構造緩和: 2~3日程度

予定しているすべての計算に数日以上かかるのであれば、精度パラメータを吟味して計算時間を節約する意味がある

今日の問: いかにして、必要とする精度を出せる
精度パラメータを決めるか

量子計算に必要な精度
1s軌道エネルギー (原子番号 Z)
H原子: 13.6 eV
Sn原子 (Z = 50) : 680 eV
全電子法の全エネルギー
擬ポテンシャル法の全エネルギー: 数十 keV: 数百 eV

物性にかかわる量:300Kの熱エネルギー : 26 meV
磁性相互作用: 数 meV
構造緩和の収束条件 : 数 meV

全エネルギーの有効数字計算過程での有効数字: 最低 6桁: 10桁？

64bit 浮動小数点型の有効数字: 丸め誤差

標準: IEEE 754 (binary64):
符号 : 1 bit
指数 : 11 bit (–1024 ~ +1023)
仮数 : 52 bit (4,503,599,627,370,495: 16桁)

64bit計算をする限り、変数型による誤差は無視できる計算過程で生じる誤差: 以後、計算誤差 と呼ぶ

・ 打切り誤差・ 無限個の平面波基底を使えば誤差を ０ にできる=> 有限個に制限する必要がある
・ 相互作用の和・ 積分誤差 (微分誤差は十分小さくできる)・ 収束誤差精度パラメータを適切に設定、必要な計算精度を確保
精度をいくらでもあげられる (変数型の丸め誤差内)

精度パラメータではどうしようもない誤差: モデル誤差
近似、汎関数・ 汎関数
・ Born-Oppenheimerの断熱近似 (原子核を固定)
フォノン－電子散乱
フォノン－フォノン散乱

・ 温度
熱膨張
自由エネルギー (相転移温度)

・ ゼロ点振動計算精度を上げても改善できない=> モデル誤差を勘案したうえで議論をする

可能な精度
計算誤差: < 10 meV/原子モデル誤差: 0 K、DFT
LDA: 300Kの格子定数の実測値に対して|誤差| < 2% (典型的には過小評価)
GGA PBE: 300Kの格子定数の実測値に対して|誤差| < 2% (典型的には過大評価)
バンドギャップ、エネルギー準位、有効質量などの絶対値:
実測値の絶対値を再現できると期待してはいけない

第一原理法の計算結果の検証
第一原理計算であるから・・・・ どの計算方法、プログラムを使っても、物理的に意味のある量は必要な精度内で同じ結果である必要。
安定構造、全エネルギー差（凝集エネルギー）
全状態密度、占有状態のエネルギーバンド
例外: エネルギーの絶対値、部分状態密度

結果の確認方法
・ より高い精度のパラメータで計算を行い、結果が同一であることを確認する
・ 異なるプログラム、近似モデルで計算を行い、結果がモデルの精度差内で整合することを確認する
・ 実験結果を説明できるかどうか

自己無同着 (SCF) 計算
収束誤差

密度汎関数理論 (Density Function Theory: DFT)
数学的定理: Hohenberg-Kohnの定理
１．外部ポテンシャルVext(r)中で相互作用しあう電子ガスについて、
全電子密度(r) が与えられると、Vext(r)は一意的に決まる。
２．全エネルギーは(r)の汎関数 E[(r)] で与えられ、最小のE[(r)] を与える (r) が基底状態の電子密度
・ 形式的に、「交換相互作用」と「電子相関相互作用」の両方とも、
全電子密度の汎関数として厳密に扱える
簡単に電子相関を取り込める => 多体問題に適している・ ただし、汎関数の形を通して近似が入る
・ 実際の電子密度を再現する１電子系を解く問題に還元する「１自由電子」ではSchrödinger方程式の解に一致

Schrödinger方程式と密度汎関数理論
Hartree-Fock (HF) 方程式 (一電子Schrödinger方程式)
1 2− ∇𝑙 + 𝑉𝑒𝑥𝑡 𝐫𝑙 + 𝑉𝑒−𝑒 𝐫𝑙 𝜑𝑙 𝐫𝑙 + 𝑽𝑿𝒍 𝒓𝒍 , 𝝋𝒊 𝒓𝒍 = 𝜀𝑙 𝜑𝑙 𝐫𝑙
2
Kohn-Sham方程式 (密度汎関数法, DFT: Density Functional Theory)
1 2− ∇ + 𝑉𝑒𝑥𝑡 𝜌 𝐫 + 𝑉𝑒−𝑒 𝜌 𝐫 + 𝑽𝑿𝑪 𝝆 𝒓 𝜑 𝐫 = 𝜀𝜑 𝐫2

• 数式表現は似ている・Schrödinger方程式:
1. 古典的なハミルトニアンを量子化
2. 各電子の座標 rl に関する方程式
3. (HF近似) エネルギー固有値 はイオン化ポテンシャル
・密度汎関数法:
1. Hohenberg-Kohn定理（電子密度で系の基底状態が決まる)
2. 空間座標 r （電子密度 (r)の汎関数）だけに依存。
3. エネルギー固有値は軌道の化学ポテンシャル

密度汎関数法の一電子方程式
Kohn-Sham方程式
1 2− ∇ + 𝑉𝑒𝑥𝑡 𝜌 𝐫
2

・ 𝑽𝑿𝑪 𝝆 𝒓

+ 𝑉𝑒−𝑒 𝜌 𝐫+ 𝑽𝑿𝑪 𝝆 𝒓

𝜑 𝐫 = 𝜀𝜑 𝐫のモデル精度を上げれば、いくらでも厳密解に近づける・ プログラムの大きな変更を必要とせず、
新しい 𝑉𝑋𝐶 𝜌 𝐫 に対応できる

密度汎関数理論: 要するに

• DFTは任意の多体系で厳密に正確な量子理論である
全ハミルトニアンの固有値 (全エネルギー)、
電子密度は物理的意味が保証されている• ただし、基底状態のみしか保証されない• 励起状態が絡む量（伝導帯、バンドギャップ、
イオン化エネルギーなど） の妥当性は保証されない• 一電子方程式化した固有値、固有ベクトルの物理的意味は慎重に議論する必要がある原子単位 (a.u.): 方程式の規格化
 2 2e2 Z  −− (r ) = E (r )4 0 r 
 2me

r ' = arE ' = bE

 1 2 Z− 2  − r  (r ) = E (r )4 0 2−11a=2= 5.2918 10 m

me e 4me eb== 13.6 eV
2 22(4 0 ) 2e2 Z 2' −− (ar') = bE ' (ar')
2
4 0 a r '  2me a

共通単位: a.u.
単位: ボーア(bohr)
H 1s 軌道の半径

単位: リュードベリ (Rydberg)
H 1s 軌道のエネルギー準位

Z 2−  − 2 r  (r ) = E (r )b=me e 4(4 0 ) 2

2= 27.2 eV

単位: ハートリー (Hartree)

バンド計算: Blochの定理
一電子方程式: 𝐻 𝑟 𝜑 𝐫 = 𝜀𝜑 𝐫
Blochの定理
𝜑𝐤𝑙 (𝐫) = exp 𝑖𝐤 ⋅ 𝐫 𝑈𝐤 (𝐫) 結晶格子の周期関数= exp 𝑖𝐤 ⋅ 𝐫 σ𝑗 𝑐𝐤𝑗 𝑢𝐤𝑗 (𝐫 − 𝐫𝐣 ) 結晶格子の周期基底関数

𝑯𝒌 𝒓 𝒖𝒌 𝒓 = 𝜺𝒌 𝑼𝒌 𝒓 : k を変えながらこの方程式を解く
𝐻𝑘 𝑟

1
=
2𝑚

−𝑖ℏ𝛻 + ℏ𝑘 2 + 𝑉𝑒𝑥𝑡 𝜌 𝑟+ 𝑉𝑒−𝑒 𝜌 𝑟+ 𝑉𝑋𝐶 𝜌 𝑟

異なる k に関して完全に独立に解くことができる: k 並列化

一電子方程式の数値解法: 自己無撞着計算
Kohn-Sham方程式
𝐻(𝒓, 𝜌 𝐫 )𝜑 𝐫 = 𝜀𝜑 𝐫 : 固有値方程式を解くことに帰する
1 2
𝐻(𝒓, 𝝆 𝒓 ) = − ∇ + 𝑉𝑒𝑥𝑡2

𝜌 𝐫 + 𝑉𝑒−𝑒 𝜌 𝐫
ハミルトニアンに 𝝆 𝒓 が含まれている+ 𝑉𝑋𝐶 𝜌 𝐫・ 𝜌𝑖𝑛𝑖 𝐫 を仮定することで 𝐻(𝒓, 𝜌 𝐫 ) が計算できる

・ Kohn − Sham方程式を解いた結果、𝜌fin 𝐫 が得られる
𝝆𝒊𝒏𝒊 (𝒓) と 𝝆𝒇𝒊𝒏 (𝒓) は同じでなければいけない=> 繰り返し計算により、 𝝆𝒊𝒏𝒊 (𝒓) = 𝝆𝒇𝒊𝒏 (𝒓) を達成する

自己無撞着 (Self-Consistent) 計算

自己無撞着 (Self-Consistent) 計算の流れ
・一電子波動関数のハミルトニアンには電子密度として、波動関数が含まれている1 2− ∇ + 𝑉𝑒𝑥𝑡 𝜌 𝐫 + 𝑉𝑒−𝑒 𝜌 𝐫 + 𝑉𝑋𝐶 𝜌 𝐫 𝜑𝑙 𝒓 = 𝜀𝜑𝑙 𝒓
2

・最初の計算では、電子密度iniを何らかの方法で計算する必要がある： 原子の電子密度の和など・計算して得られた波動関数m(rm)から求めた「電子密度fin」は、波動関数を解く際に用いたものとは違うiniとfinは一致しないと、物理的意味がない

SCFサイクル・より整合性が高いだろうと考えられる電子密度newを推定し、
fin = iniになるまで繰り返す
同様の計算を行う
例： new = ini + kmix(fin + ini)
kmix： Mixing factor
SCFサイクルでの計算が発散しないようにかける係数。
1 に近い値を使うと発散しやすくなる。
収束性の高い方法が提案されている (Broydenの方法など)

SCF計算の収束判定条件
1. 波動関数 (場) [Self-Consistent Field: SCF]
|k(i) – k(i-1)| < 
2. 電荷密度 [Self-Consistent Charge: SCC]
|(r)(i) – (r)(i-1)| < 
3. 全エネルギー (VASP EDIFF)|Etot(i) – Etot(i-1)| < E
現在、この収束判定を使っているSC計算を
SCF計算と呼んでいる

全エネルギーによるSCF計算の収束判定

基本的な判定条件: |Etot(i) – Etot(i-1)| < E
量子力学の変分原理
ハミルトニアン H の正確な波動関数 Ψ について、全エネルギーの期待値 Ψ 𝐻 Ψ が最小値を取る
Ψ と異なる任意の関数 Ψ′ に関して Ψ′ 𝐻 Ψ′ > Ψ 𝐻 Ψ

(r) が自己無撞着に達していない場合、
𝛹′ 𝐻 𝛹′ が増加することがある

Eの選び方: 必要な全エネルギー精度より十分小さく
一般的: 10 meV/原子を狙う => 数100 meV ？
VASP default: EDIFF = 10-4 eV = 0.1 meVVASPのSCF計算過程の例
EDIFF = 1.0e-4
entering main loop
N
E
dE
d eps
ncg
rms
rms(c)
DAV: 1 0.323689510286E+03 0.32369E+03 -0.11911E+04 4536 0.710E+02
DAV: 2 0.184883082244E+02 -0.30520E+03 -0.25293E+03 3488 0.215E+02
DAV: 3 -0.202800298605E+02 -0.38768E+02 -0.37728E+02 3992 0.672E+01
DAV: 4 -0.210527298854E+02 -0.77270E+00 -0.77164E+00 4008 0.978E+00
DAV: 5 -0.210538423032E+02 -0.11124E-02 -0.11124E-02 3904 0.363E-01 0.153E+01
DAV: 6 -0.178804349254E+02 0.31734E+01 -0.82537E-01 4880 0.221E+01 0.981E+00
DAV: 7 -0.178340646623E+02 0.46370E-01 -0.53612E-01 4000 0.538E+00 0.530E+00
DAV: 8 -0.178062629120E+02 0.27802E-01 -0.82326E-02 3984 0.159E+00 0.977E-01
DAV: 9 -0.178065319620E+02 -0.26905E-03 -0.90429E-03 3952 0.390E-01 0.134E-01
DAV: 10 -0.178066921234E+02 -0.16016E-03 -0.11353E-04 3720 0.617E-02 0.590E-02
DAV: 11 -0.178069231601E+02 -0.23104E-03 -0.11731E-04 3864 0.465E-02 0.123E-02
DAV: 12 -0.178069229206E+02 0.23947E-06 -0.32970E-06 2312 0.107E-02
1 F= -.17806923E+02 E0= -.17806923E+02 d E =0.000000E+00

SCF計算が収束しない場合
1.

全エネルギーは単調に減少するが、減少の仕方が遅い
Mixing parameterを増やす

2.

全エネルギーは振動するが、収束していく
Mixing parameterを減らす

3.

全エネルギーが振動して収束しない
Mixing parameterを減らす
Smearing幅を増やす

SC計算の例: 直列抵抗を含むダイオードのI-V特性
SC-Diode.xlsx
i

I
等価回路
PV

1.

I の初期値 𝐼 (0) を仮定する

2.

𝐼 (𝑖) = 𝐼0 exp

I(I𝑰01, A 1)
𝐝𝐢𝐨𝐝𝐞𝑒𝑘𝑇

𝑉 − 𝑅𝐼 (𝑖−1)−1

により𝐼(𝑖) を計算

V3.

𝐼(𝑖) − 𝐼 (𝑖−1) が収束精度以下になったら計算終了4.

𝑰 𝒊 ∗ = 𝒌𝑰(𝒊) + 𝟏 − 𝒌 𝑰 𝒊−𝟏
で 𝐼 𝑖 ∗ を計算し、 𝐼 𝑖−1 に置き換えて
2 ~ 4 を繰り返す

直列抵抗 R

𝐼𝑑𝑖𝑜𝑑𝑒 = 𝐼0

𝑒𝑉
exp
−1
𝑘𝑇

𝐼 = 𝐼0 exp𝑒
𝑉 − 𝑹𝑰
𝑘𝑇−1

k: Mixing parameter

I
0
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28

Ical
|error|I0=
2-1E-122 n=
1.8-1E-12
1.8 T=
1.62-1E-12
1.62 R=
1.458-1E-12
1.458 V=1.3122-1E-12
1.3122
1.18098-1E-12 1.18098 k=
1.062882 -9.1E-13 1.06288
0.956594 4.31E-12 0.95659
0.860934 2.09E-10 0.86093
0.774841 5.77E-09 0.77484
0.697357 1.14E-07 0.69736
0.627621 1.66E-06 0.62762
0.564859 1.86E-05 0.56484
0.508375 0.000163 0.50821
0.457554 0.00115
0.45640.411914 0.006655 0.40526
0.371388 0.031631 0.33976
0.337412 0.116849 0.22056
0.315356 0.272927 0.04243
0.311113 0.321305 0.01019
0.312132 0.308953 0.00318
0.311814 0.312754 0.00094
0.311908 0.311626 0.00028
0.31188 0.311965 8.5E-05
0.311888 0.311863 2.5E-05
0.311886 0.311893 7.6E-06
0.311887 0.311884 2.3E-06
0.311886 0.311887 6.8E-07
0.311886 0.311886
2E-07

1.E-12 A
1
300 K
1 ohm
1
0.1

Density mixing methods

1. Straight mixing
𝝆𝒎𝒊𝒙 (𝑮) 𝒊 = 𝝆(𝑮)(𝒊−𝟏) + 𝑨 𝝆(𝑮)(𝒊) − 𝝆(𝑮) 𝒊−𝟏
A: mixing parameter (VASP:AMIX)
VASP default: 0.8 for ISPIN=1 & US-PP, 0.4 for ISPIN=2 or PAW

2. Kerker mixing (VASP:IMIX=1)

𝝆𝒎𝒊𝒙

(𝑮) 𝒊

= 𝝆(𝑮)(𝒊−𝟏) + 𝑨

𝑮𝟐
𝑮𝟐 +𝑩𝟐

𝝆(𝑮) 𝒊 − 𝝆(𝑮) 𝒊−𝟏

B = 0: Straight mixing (VASP:BMIX)
VASP default: 1.0
3. Broyden’s 2nd method and Pulay-mixing method (VASP:IMIX=4, default)
𝝆𝒎𝒊𝒙 (𝑮)𝒊= 𝝆(𝑮)(𝒊−𝟏)+ 𝑨 ∙ 𝐦𝐢𝐧

𝑮𝟐
, 𝑨𝐦𝐢𝐧𝑮𝟐 +𝑩𝟐

𝝆 𝑮

𝒊

−𝝆 𝑮

だいたいうまくいく: Amin = 0.4
金属の場合、Amin = 0.02のような小さい値が必要な場合もSCFサイクルが振動する => Aを小さく
SCFサイクルの収束が遅い => Aを大きく

𝒊−𝟏

バンド計算で使われるSmearing (ぼかし) の目的st
1. 1 BZの積分の高精度化 (補間)2. SCFの収束を安定化 (拡がった分布)
3. DOSの表示をみやすく (平滑化)SCF過程でρ(r)が不連続に変化する場合がある
AB結晶、HOMO－LUMOギャップが小さい (あるいは金属) の場合

ρini(r)
LUMO
EF
HOMO

ρfin(r)

i
A

B

i+1

Smearingがある場合
Smearing
関数
A

ρ(r)の
不連続な変化 全エネルギーの飛び
SCFの振動B

準位の反転

Iteration

準位の反転

Smearingがない場合

ρ(r)の変化が小さくなるバンド計算で使われるSmearingの目的

1. 1st BZの積分の高精度化 (補間)2. SCFの収束を安定化 (拡がった分布)
・ 金属の場合に特に重要
・ Smearing関数の幅を、反転する可能性がある準位をカバーする程度に大きくとる
・ SCFが収束したら関数幅を狭くできる場合もある
3. DOSの表示をみやすく (平滑化)構造緩和計算収束誤差

構造緩和計算と体積弾性率: Si-1160.137

Energy / Ry

-1160.138-1160.139

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

260

270

280

290

300

Volume / a.u.3

E = Emin + 1 / 2 B0 (V / V0 )2

B0 (GPa) = 87.57 GPa (exp: 97.88 GPa)

一般的な構造緩和： C12A7
VASP, PBE

VASPの構造緩和計算計算過程の例

ZnO
EDIFF = 1.0e-4 EDIFFG = -0.05 k = 6 x 6 x 4 ENCUT = 550 eV
N
E
dE
d eps
ncg rmsrms(c)
DAV: 1 -0.175746399246E+02 -0.17575E+02 -0.64567E+01 4208 0.259E+01 0.499E+00
DAV: 2 -0.186060451113E+02 -0.10314E+01 -0.13644E+01 3456 0.105E+01 0.314E+00
DAV: 3 -0.185976630745E+02 0.83820E-02 -0.18981E-01 3712 0.364E+00 0.134E+00
DAV: 4 -0.185980330337E+02 -0.36996E-03 -0.14890E-02 3416 0.420E-01 0.560E-01
DAV: 5 -0.185979882095E+02 0.44824E-04 -0.41890E-04 3672 0.155E-01
1 F= -.18597988E+02 E0= -.18597988E+02 d E =-.185980E+02
N
E
dE
d eps
ncg rmsrms(c)
DAV: 1 -0.273747994327E+02 -0.87768E+01 -0.40042E+01 3944 0.820E+01 0.223E+01
DAV: 8 -0.181697542307E+02 -0.34101E-04 -0.52887E-06 2992 0.602E-02
2 F= -.18169754E+02 E0= -.18169754E+02 d E =0.428234E+00
3 F= -.18598894E+02 E0= -.18598894E+02 d E =-.906219E-03
N
E
dE
d eps
ncg rmsrms(c)
DAV: 1 -0.186057936884E+02 -0.69047E-02 -0.23927E-02 4152 0.192E+00 0.595E-01
DAV: 2 -0.185999755658E+02 0.58181E-02 -0.66465E-03 3944 0.104E+00 0.204E-01
DAV: 5 -0.185993091877E+02 -0.67055E-05 -0.73909E-05 2072 0.270E-02
4 F= -.18599309E+02 E0= -.18599309E+02 d E =-.132098E-02
reached required accuracy - stopping structural energy minimization
収束したことを確認

VASPの構造緩和計算過程の例-300
1.E-02
1.E-04
1.E-065

10 15 20 25

イオンステップ 4

1.E+00イオンステップ 3

25

エネルギー変化 / eV20

イオンステップ 4

15
イオンステップ 3

10
イオンステップ 2-25

5
イオンステップ 1

エネルギー / eV0

1.E+02

イオンステップ/SCFサイクル数

イオンステップ 2

イオンステップ/SCFサイクル数

-15
-20

EDIFFG = -0.05 k = 6 x 6 x 4 ENCUT = 550 eV

イオンステップ 1

ZnO
EDIFF = 1.0e-4

構造緩和計算の判定収束条件
内部座標のみの最適化• 座標変位 (Displacement)|xl,m (i) – xl,m(i-1)| < x
• Hellmann-Feynman theorem力 (VASPで EDIFFG < 0 のとき)
|Fk(i) – Fk(i-1)| < F
格子定数を含めた最適化 (Full relaxation)• 応力テンソル (Stress)|Plm(i) – Plm(i-1)| < P
• 格子定数変化|al(i) – al(i-1)| < a

いずれの場合も、全エネルギーが最小値を取る
全エネルギーで収束判定をすることも多い

精度 Ultra-fine で実行収束過程の確認 (CASTEP)

エネルギーは低くなる方向に変化しないといけない＝＞
そうでない場合は計算を中止し、初期値、アルゴリズムを再検討してやりなおす

変化量は必ずしも単調に減少しない

圧力差最大値
力最大値
原子位置変化

全エネルギー

横線は収束判定条件構造緩和の精度: エネルギー精度と座標精度
Mg – O ポテンシャルでの例

2
1.8
1.6
Kawamuraポテンシャル
𝑍𝑀𝑔 =2
1.4
𝑍𝑂 = –2
1.2
𝐴 = 2.53x10-14 J1
𝐵 = 0.165 Å
0.8
0.6
0.4
VASPの初期値
0.2
SCF (電子ステップ) の収束条件
0
EDIFF = 10-4 eV
構造緩和 (イオンステップ) の収束条件 1.555

𝑼(𝒓)𝑴𝒈−𝑶

𝟎. 𝟐𝑼(𝒓)𝑴𝒈−𝑶

エネルギー / eV

𝑒 2 𝑍𝑀𝑔 𝑍𝑂
𝑟
𝑈(𝑟)𝑀𝑔−𝑂 =+ 𝐴 exp −4𝜋𝜀0 𝑟𝐵

EDIFFG = 10-3 eVΔr/2 = 0.0024 Å
Δr/2 =0.0053 Å

EDIFFG = 1 meV
1.56
1.565
原子間距離 / Å

1.57

構造緩和の精度: 力精度と座標精度
1

Mg – O ポテンシャルでの例

0
1.555

F (eV/A)

1.56

1.565

1.57

1.56

1.565

1.57

-0.5

2
-1

エネルギー / eV

VASP: EDIFFG < 0 の時、Herman-Feymann力で収束判定
構造緩和 (イオンステップ) の収束条件
EDIFFG = –0.05=> |Fi| < 0.05 eV/Å

0.5

力 / eVÅ-1

𝜕𝑈 (𝑟)𝑀𝑔−𝑂
𝐹𝑖 (𝑟)𝑀𝑔−𝑂 = −
𝜕𝑥𝑖

F*0.2

1.5

1
0.5
0
1.555

原子間距離 / Å

構造緩和の精度: エネルギー精度、力精度と座標精度∆𝑈 𝑟

2 𝑈(𝑟 )
𝑑𝑈(𝑟0 )
1 𝑑 2 𝑈(𝑟0 )
1
𝑑
0
2
2
=
𝛿𝑟 +
𝛿𝑟
~
𝛿𝑟
𝑑𝑟
2 𝑑𝑟 22 𝑑𝑟 2

𝜹𝒓~ 𝟐𝐄𝐃𝐈𝐅𝐅𝐆/𝑼′′

𝐹 𝑟

𝑑𝑈 𝑟
𝑑𝑈 𝑟0
=−
~−
𝑑𝑟
𝑑𝑟

𝜹𝒓~|𝐄𝐃𝐈𝐅𝐅𝐆|/𝑼′′

𝑑 2 𝑈 𝑟0
𝑑 2 𝑈 𝑟0
−
𝛿𝑟~ −
𝛿𝑟
2
2
𝑑𝑟
𝑑𝑟

必要な精度低い精度でも大丈夫
1. 電子構造: バンド構造、状態密度、電荷ある程度高い精度が必要
1. 構造緩和: 格子定数、原子座標を実測値と比較する場合
格子定数の緩和計算には大きいカットオフエネルギーが必要
2. 光学スペクトル、キャリア輸送特性 (BoltzTraP2):
多くのk点数が必要 (BoltzTraP2には E(k) の内挿機能がある)
バンドギャップも合わせる必要がある

非常に高い精度が必要
1. フォノン計算、イオン誘電率など
特に 高誘電体、強誘電体、相転移が近い系など
格子定数の再現性も必要

イオン結晶の誘電率: 高精度の構造緩和計算が必要
・ 電子誘電率: 光学スペクトル（複素誘電率） 実部の 0 eVへの極限
VASPでは、INCARで LOPTICS=.TRUE. を設定
・ バンドギャップによって大きく影響を受けるので、実際の値に近いバンドギャップが得られる条件を使う
・ 空バンド数を十分入れる (VASP:NBANDS, 初期値の２～３倍)・ k点数を段階的に増やしながら、光学スペクトルが変化しない条件を探す・ イオン誘電率: LEPSILON = .FALSE. と設定
・ 結晶構造（原子座標、格子定数）に敏感な場合がある
構造緩和で、収束判定条件 (VASP:EDIFFG) を厳しめに設定する・ PBEsolのように、格子定数の再現性の高い汎関数を選ぶ

・ 荒い近似: 誘電率を電子誘電率とイオン誘電率の単純な和と近似
ε(total) ～ ε(optics) + ε(ion) - ε0複素誘電関数 (光学スペクトル)
ZnO, VASP, PBE, PREC=High, NBANDS=66 (デフォルト値22), CSHIFT=0.1, Tetrahdron法

14

k点数-NBANDS

12

ε1

10
8
6

8

7

ε1

4 4 3-66
6 6 4-66
8 8 6-66
10 10 8-66
12 12 10-66
14 14 10-66
16 16 12-66

8 8 6-66

6

10 10 8-66
12 12 10-66

5

4

14 14 10-66

2

16 16 12-66

4
0

0.2

0.4

0.6

E / eV

0
0

1

2

3

E / eV

4

5

0.8

1

Rutileの構造緩和計算
WIEN2k, PBE

実測の結晶構造
Experimental
a=
c=
x(O)=

緩和構造
Relaxed

4.594
2.959
0.306

Ti位置

a=
c=
x(O)=

4.64912.97230.30520-4018.060-4018.065-4018.070-4018.075

0.05

0.1

-0.1

-0.05 -4018.072 0

全エネルギー-0.05

全エネルギー-0.1-4018.074-4018.076-4018.078-4018.080-4018.082-4018.084-4018.080-4018.085

Ti位置-4018.070-4018.055-4018.086-4018.088

一般的な傾向： LDAは格子定数を過小評価
GGA(PBE) は過大評価

0.05

0.1

GGA 汎関数の効果
ZnO, VASP6, PREC=High, k mesh = 7 x 7 x 4
VASP
GGA tag

Functional
exp (300 K)
ICSD29272

Δa(%)

z

c

a
3.24275.19480.3826

Δc(%)

PE

PBE

3.28765.2821

0.38163

1.38

1.68

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

一次元のイオン結晶の誘電率-q+qEquilibrium

イオン変位の見積もり

lCase 1

外場が無い時の全エネルギー

l+x

l

1
2
U 0 = U ij (rij ) =  k (xi − xi −1 − l )2xi: i番目の原子の変位を部分座標で表す

lCase 2l’l’

外部電場Eによって電荷を持つイオンが変位したならば、系の全エネルギーは
1
2
U =  k (xi − xi −1 − l ) −  qi xi E
2

平衡条件

𝜕𝑈
=0
𝜕𝑥𝑗

𝑞𝑗

=> 𝛿𝑥𝑗 = 2𝑘 𝐸

イオン分極の誘電率𝒒𝟐𝜺 = 𝜺𝟎 +𝟒𝒍𝒌P= q jx j
2l

q2
=
E
4lk

D = E =  0 E + P

ポテンシャルの曲率 k で誘電率は決定される。

一次元のイオン結晶の誘電率2q = 0 +4lk

ポテンシャルの曲率 k で誘電率は決定される。

ルチル型TiO2のε33: ~180
一般的な酸化物の２０倍 => 曲率は１／２０
EDIFFG = –0.05 => δz(Ti) ~ 0.02
EDIFFG = –0.01 => δz(Ti) ~ 0.004

k点数 (k mesh / k spacing)

積分誤差バンド計算: Blochの定理
一電子方程式: 𝐻 𝑟 𝜑 𝐫 = 𝜀𝜑 𝐫
Blochの定理
𝜑𝐤𝑙 (𝐫) = exp 𝑖𝐤 ⋅ 𝐫 𝑢𝐤 (𝐫) 結晶格子の周期関数= exp 𝑖𝐤 ⋅ 𝐫 σ𝑗 𝑐𝐤𝑗 𝑢𝐤𝑗 (𝐫 − 𝐫𝐣 ) 結晶格子の周期基底関数

𝑯𝒌 𝒓 𝒖𝒌 𝒓 = 𝜺𝒌 𝒖𝒌 𝒓 : k を変えながらこの方程式を解く
1
2
𝐻𝒌 𝒓 =−𝑖ℏ∇ + ℏ𝒌 + 𝑉 𝒓
2𝑚異なる k に関して完全に独立に解くことができる: k 並列化

波数ベクトル k と第一原理計算

・ 異なる波数 k をもつ結晶軌道の形は異なる
・ 異なる波数 k をもつ電子密度は異なる
・ すべての k での電子密度を足し合わせる必要

・ 電子密度(波動関数)は k の連続関数なので、第一ブリルアンゾーン内の何点かを計算すればよい
第１B.Z.内で数値積分

バンド計算で使われるSmearingの目的

st1. 1 BZの積分の高精度化 (補間)2. SCFの収束を安定化 (分布)3. DOSの表示をみやすく (平滑化)数値積分: Rieman 積分

𝒈 𝒙
𝑥

𝑛

න 𝑔(𝑥′)𝑑𝑥′ = ෍ 𝒈 𝒙𝒊 ℎ𝑥0

x 0 x1 x2

xn

𝑖=0

x’
xi = x0 + ih

数値積分: 台形公式
渡部力ら監修、Fortran77による数値計算
ソフトウェア, 丸善株式会社 (平成1年)

h3 n
f ' ' (xi )誤差

12 i =0𝑛𝑥𝒈(𝒙𝒊 ) + 𝒈(𝒙𝒊+𝟏 )න 𝑔(𝑥′)𝑑𝑥′ = ෍
ℎ
𝟐
𝑥0

g(x’)

𝑖=0

x0 x1 x2

xnx’

Tetrahedron法: 線形補間 (台形公式に対応)
1. 第1ブリルアンゾーンを四面体に分割
2. 1つの四面体の頂点を(x0, y0, z0), (x1, y1, z1), (x2, y2, z2), (x3, y3, z3)とし, 頂点座標を [0, 1] に正規化する
3. 以下の一次式で補間する
𝐸 𝒌 = 𝐸000
+ 𝐸100 − 𝐸000 𝑘𝑥+ 𝐸010 − 𝐸000 𝑘𝑦+ 𝐸001 − 𝐸000 𝑘𝑧
𝐸𝑖𝑗𝑘 は 頂点 (i, j, k) におけるE(k)4．E(k) を積分する(0, 0, 1)(0, 1, 0)(0, 0, 0)(0, 1, 0)高次の積分: Methfessel-Paxton関数
M. Methfessel and A.T. Paxton, High-precision samplingfor Brillouin-zone integration in metals, Phys. Rev. B 40 (1989) 3616

デルタ関数をエルミート多項式で展開

n
(
− 1)
A =

 (x ) =  An H 2 n (x ) exp(− x 2 )


n

n =0

n!4 n DN(x)は (2N+1)次の多項式、2N次以下の多項式に直交(

N

DN ( x ) =  An H 2 n ( x ) exp − x 2

)

n =0

ステップ関数の近似
S N ( x ) = 1 −  DN (t )dtx−S 0 ( x ) = (1 / 2)(1 − erf ( x))

N=0
N = 10x

エルミート多項式 d2
d 2 − 2 x + 2n  H n ( x ) = 0 の解
dx
 dx


Hn(x)exp(-x2/2)は直交基底





−

(

)

H n ( x )H m ( x ) exp − x 2 dx =  mn 2 n  n!

m
(
− 1)
n−2m
(2 x )H n ( x ) = n! 
m = 0 m!(n − 2m )!
int( n / 2 )

H 0 ( x ) = 1, H1 ( x ) = 2 x H n ( x ) = 2 xH n −1 ( x ) − 2(n − 1)H n − 2 ( x )H n ' ( x ) = 2nH n −1 ( x ) = 2 xH n ( x ) − H n +1 ( x )調和振動子モデルの波動関数：(

n ( x ) = 2n) H (x )exp(− x / 2)

 n!

1/ 2

2

n

Methfessel Paxton関数の特徴
1. バンド構造、DOSが 2N 次以下の多項式で近似できる
場合、 N次のMP関数でSmearingしても積分誤差は発生しない
2. 単純なバンド構造、DOSの場合は、SMEAR幅をかなり大きく(Al で 1 eV以上) しても積分誤差は発生しない

3. バンド構造が複雑な場合 (d系など)、
SMEAR値の最適値はGaussianと同程度になるw1/2

4. 実際の Smearing幅 w1/2 は、
次数 N とSMEAR幅の値で変わる

SMEAR = 1

5. Occupancy に負の値や1.0 を超える値がでるNk点数と積分方法の比較
ZnO, PBE

Tetrahedron法
nx

ny

nz

Gaussian smearing

ISMEAR SIGMA

TOTEN

dE(meV)

nx

ny

nz

ISMEAR SIGMA

TOTEN

dE(meV)

10

10

8

-5

-17.808

0.000

1

1

1

0

0.1

-7.17010638.403

8

8

6

-5

-17.808

0.299

2

2

2

0

0.1

-17.027780.942

7

7

4

-5

-17.807

1.045

3

3

2

0

0.1

-17.684

124.285

5

5

3

-5

-17.800

7.600

3

3

3

0

0.1

-17.693114.540

4

4

3

-5

-17.784

24.175

4

4

3

0

0.1

-17.784

24.175

3

3

3

-5

-17.693114.540

5

5

3

0

0.1

-17.800

7.600

3

3

2

-5

-17.684

124.285

7

7

4

0

0.1

-17.807

1.045

2

2

2

-5

-17.027780.942

7

7

4

0

0.3

-17.807

0.946

7

7

4

0

0.5

-17.808-0.197

7

7

4

0

1

-17.821-12.766

単純な E(k) の半導体
• tetrahedron法でもGaussian smearingでも大きな違いはない• Gaussian smearingのガウス関数の幅 SIGMA も大きな違いはない
大きくとりすぎると伝導帯にも電子が入るので、誤差が大きくなるk点数と積分方法の比較
ZnO, PBE

Tetrahedron法
nx

ny

nz

Methfessel Paxton関数

ISMEAR SIGMA

TOTEN

dE(meV)

nx

ny

nz

ISMEAR SIGMA

TOTEN

dE(meV)

10

10

8

-5

-17.808

0.000

4

4

3

3

0.1

-17.784

24.279

8

8

6

-5

-17.808

0.299

4

4

3

2

0.1

-17.784

24.175

7

7

4

-5

-17.807

1.045

4

4

3

1

0.1

-17.784

24.432

5

5

3

-5

-17.800

7.600

4

4

3

5

0.1

-17.784

24.347

4

4

3

-5

-17.784

24.175

10

10

8

5

0.1

-17.808

0.007

3

3

3

-5

-17.693114.540

10

10

8

10

0.1

-17.808

0.000

3

3

2

-5

-17.684

124.285

7

7

4

3

0.1

-17.807

1.071

2

2

2

-5

-17.027780.942

7

7

4

5

0.1

-17.807

1.096

7

7

4

10

0.1

-17.807

1.062

6

6

3

5

0.1

-17.804

3.472

5

5

3

5

0.1

-17.800

7.710

単純な E(k) の半導体
• Methfessel Paxton関数のメリットははっきりしないk点数の効果: MgOの全エネルギー曲面
WIEN2k, k点数 200点

Exp.(RT)a = 0.422 nm (7.98 bohr)
V = 126.9bohr3 (primitive cell)Opt.
a = 0.4256 nm (8.046 bohr)
V = 130.2089

k点数の効果: MgOの全エネルギー曲面
WIEN2k, k点数 400点

-551.424

Energy / Ry-551.425

Exp.(RT)
ac = 0.422 nm
VR = 126.9 a.u.3 (primitive cell)-551.426-551.427

Opt.
ac = 0.4254 nm
VR = 130.0164 a.u.3-551.428-551.429-551.430-551.431

120

124

128

132
Volume / a.u.3

140

136

E = Emin + 1 / 2 B0 (V / V0 )2

B0 (GPa) = 156.9 GPa (exp: 162 GPa)

バンド計算における k点の選択
バンド理論において物理量

E(k), i(k)は k について連続に変化する＊実際の計算では、E(k) の特徴が十分反映できる程度のk点だけで計算すればよい only計算
・原子、分子・数百原子を含む結晶“良い固体” ではよくないバンド計算における k点の選択
 + BZ境界 (k mesh 222)・結合、反結合を考慮できる
・単純な酸化物では比較的特徴を表す“波動関数が拡がっている”(BZ内部で曲率が反転する)と場合はよくないバンド計算における k点の選択
実際には・・・
1. 最初に荒い k点メッシュ で計算
2. E(k)の特徴が反映できる
k間隔 k を求め、BZ内のk点数を計算する
3. kx, ky, kz をk 間隔ごとに変えながら KS方程式を解く
4. k を減らしながら(全k点数を増やしながら)計算し、計算結果 (全エネル
ギー、バンド構造、状態密度)が変わらなくなる最小のk点が望ましいカットオフエネルギー(平面波基底を使う場合)

打切り誤差

電子は波である: 平面波が良い基底関数である理由
Schrödinger方程式
 2 2 + V (r ) (r ) = E(r )
−
 2m


(r ) = C exp (ik (r)  r )e-, Ek (r ) =2m(E − V (r ))
2

E=5eV

5

V (r )

Psi(real part)2m(E − V (r ))k (r ) =
2


6

転送行列法

4eV

4

3eV

3

2eV

2
V(x)

1eV

1
0
-3

-2

-1

0

Position / nm12

平面波基底Fourier変換 (数学の定理)
2𝜋
周期関数は波数 𝑘𝑙 = 𝑙 の平面波の無限和で正確に表現できる
𝑎1次元：𝑓 𝑥= σ∞𝑙=−∞ 𝐴𝑙 exp

2𝜋
𝑖 𝑙𝑥
𝑎

3次元：𝑓 𝐫 = σ∞ℎ,𝑘,𝑙=−∞ 𝐴ℎ,𝑘,𝑙 exp 𝑖𝐆ℎ𝑘𝑙 ⋅ 𝐫

無限個の平面波基底を使えれば、どのような関数も正確に表現できる=>計算時間、メモリーの制限から、有限個の平面波だけを使う: 打切り誤差

平面波法
一次結合の基底関数として平面波を使う

k (r ) = exp(ik  r ) Chkl uhkl (r )uhkl (r ) = expiG hkl  r 

波数 Ghkl の平面波は格子周期の関数の完全基底系:
すべての hkl について和を取れば、完全に正しい解になる=> 実際の計算では |Ghkl| < Gmax の範囲で近似する
H 11 − ES11

H 12 − ES12

H 21 − ES 21

H 22 − ES ss

H 1n − ES1n
H 2 n − ES 2 n

H n1 − ES n1

H n 2 − ES n 2  H nn − ES nn

uh 'k 'l ' H uhkl =  e−i (k + G h'k'l' )r=0 2 2
 i (k +G hkl )r + V (r )e
dr
−
 2m


2 2
=  hkl ,h 'k ',l 'k + V * (G hkl − G h 'k 'l ' )2m
実際の計算のほとんどがポテンシャルのフーリエ変換=> GPUで高速化が容易平面波基底の打切り• 𝒌ℎ𝑘𝑙 = 𝑮ℎ𝑘𝑙 の最大値で指定
カットオフ波数2ℏ2ℏ
2• 𝐸𝒌 =
𝑘𝑐𝑢𝑡 =
2𝑚𝑒
2𝑚𝑒

𝐆ℎ𝑘𝑙,𝑐𝑢𝑡2の最大値で指定

カットオフエネルギー
WIEN2kの場合: Rkmax=Min(RMT)*kcut

波動関数の変化スケールと Ecut: Na原子の3s動径関数0.8

DV-X法で計算

Psi / arb. unit

0.6

2kcut ~
L

0.4

L

0.2
a ~ 0.5nm L ~ 0.01nm
必要な平面波の数
2𝑎/𝐿 3 ~ (2*50)3 = 1,000,000
VASP: Ecut = 400 eV L ~ 0.064 nm

0

-0.2
-0.4

0

0.5

1
1.5
Radius / A

2

2.5

3

カットオフエネルギーの効果
ZnO, VASP, PBE, PREC=High, NBANDS=22, k = 7x7x4

20

ΔE / meV

0
-20

-40
-60

300

400

500
Energy / eV

600

700

EDIFF ENCUT TOTENdE(meV)
1.00E-04
350 -17.880609 -61.4945
1.00E-04
400 -17.839466 -20.3519
1.00E-04
440 -17.816869 2.24513
1.00E-04
480 -17.808055 11.05938
1.00E-04
520 -17.805917 13.19678
1.00E-04
550 -17.810213 8.90098
1.00E-04
600 -17.81168 7.43439
1.00E-04
650 -17.814328 4.78603
1.00E-04
700 -17.817961 1.15269
1.00E-04
800 -17.819114
01.00E-05
1.00E-05

800 -17.819114 0.0003
480 -17.808053 11.06101-80
10 meV/原子 (40 meV/cell) の精度を出すには 400 eV で十分
構造緩和などで 数meV/cell の精度を出すためには 650 eV 以上が必要

計算精度の決め方: 基本方針

基本方針
1. 必要な全エネルギー精度 εE を決める
2. 精度パラメータを大きくしながら、全エネルギーの変化が εE 以下になる
最小パラメータを見つける
どのように精度パラメータを増やすか
・ 計算時間を見積もる
対角化は基底関数の数 N に対して O(N3) の計算時間がかかるN ∝ kmax3 ∝ Ecut3/2: Ecut を 400 eV から 450 eVに増やすと、N は 1.2倍、対角化の計算時間は 1.7倍になる
k点数は nk * ny * nz ~ nk3 に比例する
nk を 4 から 5に増やすと、N は 1.95倍、対角化の計算時間は7.4倍になる

平面波・擬ポテンシャル法におけるEcut
擬ポテンシャルの種類・作り方に依存するので、マニュアルを読むこと

一般論：
1. ウルトラソフト擬ポテンシャル : 250 eV～？
2. ノルム保存擬ポテンシャル : 500 eV～？
3. VASP (POTCARにデフォルト値が書かれている):
ほとんどのポテンシャル : 400 eV～
ハードポテンシャル: 500 eV～
4. 圧力計算の場合は、Ecut を1.3倍程度にする(Pulay stressの誤差をなくす)数値積分の精度: ak積
𝒈 𝒙x0x1x2x’ 𝑥𝑖 = 𝑥0 + 𝑖𝒉

xn

変数のメッシュ幅 h で積分精度が決まる

B.Z. の積分では、
𝟐𝝅 𝟏
∆𝒌 =
(KSPACING)が同じなら同程度の積分精度になる
𝒂 𝒏𝒌

神谷の場合:
𝟐𝝅
ak積 𝒂𝒏𝒌 = , 一定 で k点数の目安をつける∆𝒌k点数の選び方: 神谷のやり方
ak積: a は格子定数、k は a 方向の k 点数

荒い計算
高精度計算
非常に高精度a [nm]
0.4
0.4
0.4

ak [nm]
1.5
2.0
2.5

nk
4
5
7 (6)

Δk [nm-1]
0.625
0.500
0.357 (0.417)

LaCo2B2: EF でバンドが交錯する例
LUMO:
0.116 eV

La
B
Co

HOMO:
-0.341 eV

多くのk点数を取る必要があるnk = 7 x 7 x 3
a = b = 3.62 Å
c = 10.23 Å
ak = 2.53, 2.53, 3.10 nm
MA

Z 

X

M



バンド計算における k点 の選択
1. k点間隔はバンド構造の変化を反映できるだけ細かく
・ イオン結晶、直接遷移半導体:
荒いメッシュでも比較的精度が出る
・ 金属、BZ内バンドの折れ曲がり、
Fermi準位で複数のバンドが交錯している場合など
細かくとらないと精度が出ない
2. 点を含めないほうが (若干)精度が上がる(Monkhorst-pack)Gamma-centered vs Monkhorst pack
H.J. Monkhorst, J.D.Pack, Phys. Rev. B 13, 5188 (1976)
山内淳、表面科学 Vol. 28, pp. 135-143 (2007)

k点は対称性により Nk 個の等価位置を持つ: その k点の計算結果は、 Nk 個の k点の重み=> 対称要素上に無い k 点を使った方が得
Monkhorst-PackΓ-centered mesh𝑛𝑖 ∗𝑛𝑖 +𝟏/𝟐 ∗σ𝒌=σ𝒌{𝑛𝑖} = 𝑖=𝑥,𝑦,𝑧𝒂 𝒊 {𝑛𝑖 }
𝑖=𝑥,𝑦,𝑧 𝑁 𝒂 𝒊𝑁𝑖𝑖𝑛𝑖 = 0, 1,2, ⋯ , 𝑁𝑖

注意:
Tetrahedron法の計算では
Γ-centeredの方が簡単
x, y鏡映面があれば、A点だけの計算で良い

G, E, F, H点 4点の計算が必要

VASPでの精度確認: 神谷のやり方
Ecut (VASP:ENCUT)
1.
2.

400 (350) eV から開始
圧力計算をする場合は、1.3倍、450 eVまで視野に入れる
Ecut を 10%程度ずつ増やす(350, 380, ) 400, 440 (, 480, 520)

k点数 (VASP:KPOINTS)
1. ak = 2.0 (1.5) の k点数から開始。a = 0.4 nm なら nk = 5 (4)2. nk が小さい場合は1ずつ増やす
nk が大きい （6以上） 場合は 2ずつ増やす

全エネルギーだけでなく、安定構造、バンド構造、状態密度も収束していることを確認

汎関数、擬ポテンシャル: モデル誤差

汎関数

Kohn-Sham方程式
 1 2
−  + Vext ( (r )) + Ve−e ( (r )) + VXC ( (r )) (r ) = (r ) 2VXCの形がわからない=> いろいろな近似でいろいろな VXC が提案されている

汎関数
L(S)DA: Local (Spin) Density Approximation
局所密度近似1
VXc = −3 ((3 / 8 )  (r ))3
局所: 一点 r だけで決まる

非局所: r 以外の情報も考慮する ＝＞ 微分で取り入れる

GGA: Generalized Gradient Approximation
一般化密度勾配近似x =  / 

4
3

E XC = E XC

LSDA

− b   


4
3

2

x
NL
dv
+
E
X
1 + 6bx sinh −1 xExact exchange:
Hartree-Fock交換相互作用

単純金属の安定構造 (PBE vs LDA)
実験 LDA
Li b
HCP
Be h
HCP
B
Diamond
C d
Diamond
Na b
FCC
Mg h
FCC
Al f
FCC
Si d
Diamond
K b
HCP
Ca f
FCC
Sc h
HCP
Ti h
HCP
Ni f
FCC
Cu f
HCP
Zn h
HCP
Ga
HCP
Ge d
Diamond
Rb b
FCC

PBE96
BCC
HCP
Diamond
Diamond
FCC
FCC
FCC
Diamond
HCP
FCC
HCP
FCC
FCC
FCC
HCP
Diamond
HCP

Sr
Y
Zr
Nb
Mo
Rh
Pd
Ag
Cd
In
Sn
Sb
Te
Cs
Ba
La
Hf
Ta

実験
f
h
h
b
b
h
f
f
h
d

b
b
h
h
b

LDA
BCC
HCP
HCP
BCC
BCC
FCC
FCC
HCP
HCP
HCP
Diamond
HCP
SC
FCC
BCC
FCC
HCP
BCC

PBE96
BCC
HCP
HCP
BCC
BCC
FCC
FCC
FCC
HCP
HCP
HCP
SC
HCP
FCC
BCC
FCC
HCP

実験 LDA
Au f
FCC
Hg
SC
Tl h
HCP
Bi
HCP

PBE96
BCC
FCC
BCC
HCP

鉄の基底状態: LSDAとGGAFCC FeE (mRy)BCC FeVolume (a.u.3)
C. S. Wang, B. M. Klein, H. Krakauer, “Theory of Magnetic and Structural Ordering in Iron”, PRL 54, 1852 (1985)
T. C. Leung, C. T. Chan, and B. N. Harmon, “Ground-state properties of Fe, Co, Ni, and their monoxides: Results of thegeneralized gradient approximation”, PRB 44, 2923 (1991)

Lattice dynamics of TiO2 rutile: influence of gradient corrections in
density functional calculations
B. Montanari, N.M. Harrison, Chem. Phys. Lett. 364, 528 (2002)

GGA/LSDA ground state of iron （VASP,FM)
3.5

-6

GGA

Energy per Fe / eV

(US-LDA)

Spin moment

-7

FCC

-7.5

2
1.5
1

BCC
2

2.5
Fe - Fe distance / nm

FCC

0
3

-0.5

-6

3.5

-6.5

3

2

2.5

-7
FCC

-7.5
-8
-8.5

BCC

-9
-9.5

BCC

0.5

-8
-8.5

LSDA

2.5

Spin moment

Energy per Fe / eV

(PAW-PBE)

3
-6.5

2

2.5
Fe - Fe distance / nm

2.5
Fe - Fe distance / nm

3

BCC

2
1.5
1
0.5

FCC

0
3

-0.5

2

2.5
Fe - Fe distance / nm3

Rutileの構造緩和計算
WIEN2k, PBE96

実測の結晶構造
Experimental
a=
c=
x(O)=

緩和構造
Relaxed

4.594
2.959
0.306

Ti位置

a=
c=
x(O)=

4.64912.97230.30520-4018.060-4018.065-4018.070-4018.075

0.05

0.1

-0.1

-0.05 -4018.072 0

全エネルギー-0.05

全エネルギー-0.1-4018.074-4018.076-4018.078-4018.080-4018.082-4018.084-4018.080-4018.085

Ti位置-4018.070-4018.055-4018.086-4018.088

一般的な傾向： LDAは格子定数を過小評価
GGA(PBE96) は過大評価

0.05

0.1

GGA 汎関数の効果
ZnO, VASP6, PREC=High, k mesh = 7 x 7 x 4
VASP
GGA tag

Functional
exp (300 K)
ICSD29272

Δa(%)

z

c

a
3.24275.19480.3826

Δc(%)

PE

PBE

3.28765.2821

0.38163

1.38

1.68

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

近似法の精度
一般論: 全電子・全ポテンシャル法 > 擬ポテンシャル (PP) 法(OPW, LAPW/APW) > PP
擬ポテンシャル:ノルム保存PP > ウルトラソフトPP(大きなEcutが必要)(小さいEcutで精度を出せる)

PAW (Projector-Augmented Wave) 法
擬ポテンシャルの形に LAPW法を取り込む
L/APW法と同程度のエネルギー精度

化合物の生成エネルギーの比較:
LDA (US-PP) vs. PBE (PAW)

ZnO
Ga2O3
In2O3

Euc / 化学式量Eform
(kJ/mol)

Ef0 (exp)(kJ/mol)

5.921 (5.714)-21.27 / 2-422.9348.3

PBE

5.466 (5.714)-18.20 / 2-367.1

348.3

LDA

5.921 (5.943)-140.67 / 4-719.1

544.6

PBE

5.602 (5.943)-120.19 / 4-574.7544.6

LDA

7.298(7.123)-522.45 / 16-627.6462.9

PBE

6.647 (7.123)-449.06 / 16-518.2462.9

VXC

D (DX)

LDA

重要なこと1. (64bit 浮動小数点の精度以内で) 計算精度はいくらでもあげられる
計算精度: 再現性を担保するのに必要な精度
議論する内容に必要な精度を確保すること

2. 計算精度をいくら上げても、モデル精度以上に正確な計算はできない「バンドギャップ、エネルギー準位、有効質量が実測値に一致した」というのは、多くの場合、ただの偶然である

3. 構造モデル依存性と結果の精度、一般性を判断
わずかな格子定数の違いにより結果が本質的に変わる場合がある(特に高誘電率、ソフト材料、相転移近傍)
本当に、計算結果が唯一・普遍的な結果といえるのかバンドギャップ問題

電子のエネルギー (eV)

Si (WIEN2k, PBE) 実測値: Eg = 1.12 eV (300K)
8

8

6

6
(C)

4
2

4

2
(B)

0

0

-2

-2

-4

-4
(A)

-6

-6

-8

-8

-10

-10

-12

-12

-14
W

L    X WK
波数ベクトル↑0.53 eV-14
W

L    X WK
波数ベクトルバンドギャップ問題: HF近似とDFT
HF近似 (CRYSTAL06, 3(6)-21G)

DFT (WIEN2k)
PBE

10

10

5

5

0

0

-5

-5

-10

-10
-14 W

-14
L


X
波数ベクトルK

L    X WK W
波数ベクトル

LDA

L    X W K
波数ベクトルよりもっともらしい
バンドギャップの計算 － beyond DFT －(i) 配置間相互作用 (Configuration Interaction: CI)
分子軌道法 (Gaussianなど) で広く使われている
バンド計算で実行するのは難しい
バンド計算の場合(i) 自己相互作用補正 (Self-Interaction Correction: SIC)
(ii) EE項に遮蔽パラメータを入れる方法
Screened Exchange: sX近似(iii) HF近似とDFTが Eg の値を逆方向に間違える=> これらを適当な割合で混合した混合汎関数 (Hybrid DFT)(a) 実測のバンドギャップに合うようにEE項と密度汎関数を
適当な比で混合したもの。
B3PW91, B3LYPなど。混合係数が経験的に決められている。
(b) (a)に似ているが、係数が理論的に決定されている。
PBE0 や、その遮蔽版 (HSE)(iv) GW近似 (準粒子近似)

混合汎関数: PBE0, HSE03/06
PBE0混合汎関数
E xcSR , HF PBE 0 = aE xHF , SR + (1 − a )E xPBE + EcPBE
混合パラメータ a = 1/4
HSE混合汎関数
E xcPBEh = aE xHF , SR ( ) + (1 − a )E xPBE , SR ( ) + aE xPBE , LR ( ) + EcPBE

混合パラメータ a = 1/4
遮蔽パラメータ  → 0: PBE0

 → : PBE (GGA)(バンドギャップが合うように調節されたりすることもある)
HSE03:  = 0.15
HSE06 (HSE03のERRATA):  = 0.15/21/2 = 0.106 (HF part) = 0.1521/3 = 0.189 (PBE part)John P. Perdew, Matthias Ernzerhof and Kieron Burke
J. Chem. Phys. 105 (1996) 9982
Jochen Heyd, Gustavo E. Scuseria, Matthias Ernzerhof
J. Chem. Phys 118 (2003) 8207; 124 (2006) 219906

汎関数の選択に関する神谷のアドバイス
１．まず、あたりをつける ／ 構造緩和計算
計算の速い DFT を使う
金属: LDA
半導体、絶縁体: GGA (PBEsol, PBE, revPBE)
２．より正確な計算 (バンドギャップ、欠陥電荷遷移準位)
狭バンドギャップ半導体: HSE 混合汎関数
広バンドギャップ半導体: PBE0 混合汎関数
バンドギャップの再現性の高い汎関数を選択(場合によっては 混合係数や遮蔽係数を調整)

積分・Smearing

バンド計算で使われるSmearingの目的

st1. 1 BZの積分の高精度化 (補間)2. SCFの収束を安定化 (分布)3. DOSの表示をみやすく (平滑化)Smearing関数 𝒇𝒔 (𝑬)1/ 2

w1/2 = 1 eV 0.5

0.4

 ln 2 G ( x ) =  w1/ 2 

2




x
 
exp − ln 2

w1/ 2  




0.3

1
1
L( x ) =w1/ 2 1 + ( x / w1/ 2 )20.2

𝑓′ 𝑥 = exp −𝑥 2 / 1 + exp 𝑥

0.1
0-3 -2 -1 0 1 2 3
畳み込み (convolution) 演算により D(E) を平滑化:∞
𝑫 𝑬 ∗ 𝒇𝒔 (𝑬) = ‫׬‬−∞ 𝑫 𝑬′ ∗ 𝒇𝒔 (𝑬′ − 𝑬) 𝒅𝑬

2コンボリューションで平滑化した場合の問題
200

Calc. DOS
計算結果
Convoluted
畳込後

180
160

140
120
100
80
60

40
20
0
-8

-7

-6

-5

-4

-3

-2

-1

0

1

Smearing関数の幅だけ、バンド端より外側に裾がはみ出す=> 正確にバンド端、バンドギャップを決めることができない
半導体・絶縁体であれば、Tetrahedron法を使うほうが良い

状態密度のSmearing D(E)*fS(E)
D(E ) = DV 0 (EV − E )1/ 2+ DC 0 (E − EC ) * (1 + rand[−0.5,0.5])
1/ 2

=0

D(E )* f s (E ) =  D( E ' ) f s ( E '− E )dE '-2  E  4eV
E  -2 , 4  E

w = 0.2 eV

108

10次のMethfessel Paxton関数
5次のMethfessel Paxton関数

6

2次のMethfessel Paxton関数
1次のMethfessel Paxton関数
1
1
L( x ) =w 1 + (x / w)2
1/ 2
2
G ( x ) = (1 / w) exp − ( x / w)

4

(

2

f ' ( x ) = exp( x / w) / (1 + exp( x / w)))

2

D (E ) = D (E ) + D (E )0-3-2-1 0 1 2 3 4 5VEnergy (eV)C

対角化に取り入れるバンドの数(空バンドの数)
VASP:NBANDS

平面波法: 原理的には平面波が基底関数
𝜑𝐤 (𝐫) = exp 𝑖𝐤 ⋅ 𝐫 σ 𝑐ℎ𝑘𝑙 𝑢ℎ𝑘𝑙 (𝐫), 𝑢ℎ𝑘𝑙 (𝐫) = exp 𝑖𝐆𝐡𝐤𝐥 ⋅ 𝐫
H 11 − ES11

H 12 − ES12

H 21 − ES 21

H 22 − ES ss

H 1n − ES1n
H 2 n − ES 2 n

H n1 − ES n1

H n 2 − ES n 2  H nn − ES nn=0uh 'k 'l ' H uhkl =  e−i (k + G h'k'l' )r 2 2
 i (k +G hkl )r
−

+
V
(
r
)
dr

e
 2m


2 2
=  hkl ,h 'k ',l 'k + V * (G hkl − G h 'k 'l ' )2m•基底関数として平面波を使うと、固有値方程式が巨大になる
ZnO, VASP, ENCUT=500 eV の場合で23,040個の平面波基底を使う•

基底関数として固有値方程式の解 σ 𝑐ℎ𝑘𝑙 exp 𝑖𝐆𝐡𝐤𝐥 ⋅ 𝐫 を用いる=> 対称化適合軌道になっているので、対称性別に対角化を実行できるエネルギーの高い非占有軌道を対角化に入れる必要がなくなる平面波法: 固有値方程式の解を基底関数に使う• 基底関数として固有値方程式の解 σ 𝑐ℎ𝑘𝑙 exp 𝑖𝐆𝐡𝐤𝐥 ⋅ 𝐫 を用いる=> 対称化適合軌道になっているので、対称性別に対角化を実行できる• SCFに到達していれば、占有状態だけの計算で良い• SCFに到達していない場合、「空のバンド」 を入れて対角化・SCF計算を実行しないといけない: VASP:NBANDSで設定

平面波法: 空バンドの効果• SCFに到達していれば、全エネルギーへの影響はない• 非占有状態の電子構造に大きく影響
エネルギーの高い伝導帯バンド・状態密度の形
NBANDSの決め方• 一度SCF計算などを実行して、出力からデフォルト値を確認(価電子軌道数を計算してもいい)• NBANDSを1.2倍程度しながら、バンド構造、状態密度が変わらなくなることを確認する (全エネルギーは変わらない場合が多い)• 光学スペクトルの計算では (VASP)デフォルト値の２～３倍にする
神谷の場合:
小さいモデルの場合、デフォルト値の2倍を使うことも多い伝導帯の形がおかしかったら
例： ••伝導帯端から 5 eV以上上のバンド数が少ない、
スパイクや不連続がある
伝導帯の状態密度が連続的に変化せず、ピーク形状の和になっている
良い半導体であれば、伝導帯の状態密度は 𝐸 − 𝐸𝐶 に比例するあり得る原因:• 非占有軌道の数 (NBANDS) が足りない• B.Z.の積分点数が足りない
解決方法:• NBANDSを増やす• Tetrahedron法を使う

PDOSの問題

鉄系超伝導体 LaFeOP の投影状態密度

PDOSの曖昧さの問題
波動関数を 原子ごとに分割する一意的な方法はない
PDOSの計算方法
1.

2.
3.

Mulliken Population解析: LCAOでよく使われる
LCAOの係数の対角項|Ci|2
非対角項 (化学結合) CiCj + CjCi を半分ずつ
WIEN2k (L/APW法)
Muffin-Tin 球内の電子数
VASP (PW法)
各原子に半径 RWS の球を与え、その球内の電子数
2,3: 原子に割り当てられていない電子: Interstitial

4. Wigner-Seitz胞 に分割し、その中の電子数: イオン/波動関数の大きさを無視(WS胞: 原子間結合の垂直二等分面が作る最小の多面体)

5. Bader charge解析
波動関数のフラックス (電流) がゼロになる多面体で分割する
全電荷密度の谷を境界にすることと同義
・ 一意的 ・ 各原子の電子数の総和が全電子数に一致する

おまけ Wigner-Seitz (WS) 格子:
面心直方格子空間を同じ多面体で埋め尽くすにはどうしたらよいか？
Wigner-Seitz (WS) 胞:
2つの格子点の垂直二等分面が作る最小の多面体

a2a1

おまけ 逆格子のWS胞: ブリルアンゾーン
第１B.Z.: 逆格子の原点 (Γ点) を中心とする最小のWigner-Seitz胞
実格子のWS胞

a2

逆格子のWZ胞: 第一BZ

a*2
a1

Γ
a*1

Bader電荷を利用した一意的なPDOSの計算例
1. Bader charge解析を実行=> 各原子 i に割り当てられた多面体の体積 Vi (Bader volume)
2. Vi に相当する球の半径を計算 ri = (3Vi / 4π)1/3

3. kri の半径の球の体積の和 Vsphere が単位格子体積 Vcell の 90% 程度になるように、
ri を規格化: ri’ = kri を WS半径に使う
Bader電荷解析の問題点:
結晶構造中に大きな空隙があると、空隙内の電荷も最近接の原子に割り当てられる
例: WO3 のAイオン位置、蛍石型構造の中心位置、12CaO・7Al2O3など、
イオン空孔に対応する空隙を持つ結晶は多いイオン価数の曖昧さの問題
イオンの価数 ＝ 原子核電荷 － イオンに属する電子数 (PDOSの積分)・ PDOSの曖昧さと全く同じ問題
・ 共有結合の寄与のため、イオンの形式電荷より小さくなる
イオン性が強いMgOでも 形式電荷 +2 の60% 程度 (LCAO)
70~90% 程度 (イオン半径に比例)
PDOSを使わない方法 (電気分極・Berry位相)・ ボルン有効電荷 (Born effective charge)
イオン変位による誘電率から逆算
イオンの変位に伴う電子の移動による効果が加わる
誘電率が小さい物質ではそれらしい値が得られる
高／強誘電率材料でははるかに大きくなる
・ イオンを隣の単位格子に移動した際の電気分極から計算
イオン価数が整数で、バンドギャップが常に開いている
場合にしか使えないまとめ• 時間とお金を大切に• でも、外部に発表しても恥ずかしくない精度を
1. 必要な計算精度をはっきりさせる
2. そのために必要なパラメータを決定
3. 計算結果を２重・３重に確認する
・ 実験結果との整合性
・ 科学常識との整合性
常識に合わない結果が得られたら、理由を考える

・ 異なるアルゴリズムのプログラムと比較
4. モデル精度以上の議論をしない