修 士 論 文(題 目)

第一原理計算による異種金属界面の界面強度評価

平成 31 年度

岐阜大学大学院
自然科学技術研究科 博士前期課程
物質・ものづくり工学専攻

氏名 二村晟平目次

第1章

緒言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

第2章

第一原理計算の概要 . . . . . . . . . . . . . . . . . . . . . . . . .

32.1

原子単位 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32.2断熱近似と平均場近似 . . . . . . . . . . . . . . . . . . . . . . . . . .

42.3

密度汎関数理論 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52.4

局所密度近似 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

2.5

逆格子空間

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

2.6

ハミルトニアン . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.7

系のエネルギー . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.8

応力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.9

擬ポテンシャル法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.9.1

TM 型擬ポテンシャル . . . . . . . . . . . . . . . . . . . . . . 16

2.9.2

ウルトラソフト型擬ポテンシャル . . . . . . . . . . . . . . . . 20

2.10 電子占有数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.11 FFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.12 電子系の最適化手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.13 ノルム非保存形擬ポテンシャル法 . . . . . . . . . . . . . . . . . . . . 28
第3章

第 3 元素の表面エネルギーへの影響 . . . . . . . . . . . . . . . .

3.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

i

30

3.2

解析結果と考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

3.2.1

単元素の表面エネルギー . . . . . . . . . . . . . . . . . . . . . 35

3.2.2

第 3 元素による表面エネルギーへの影響 . . . . . . . . . . . . 36

第4章

異種金属界面に存在する第 3 元素の影響 . . . . . . . . . . . . .

4.1

解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4.2

第 3 元素による界面エネルギーへの影響 . . . . . . . . . . . . . . . . 45

第5章

結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

56

参考文献 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58ii

第 1 章 緒言
めっきは，単一材料では実現が難しい特性を材料に付加することを目的とした異種材料接合技術の 1 つである．自動車産業では部品の耐摩耗性の向上や防腐を目的に使用される．航空機産業では，比強度を高めた材料の疲労破壊や脆性破壊，腐食の対策として使用される．電子部品産業では導電性付加を目的に使用される．めっきがはく離した際には付加した材料の特性は失われるため，その界面強度を把握することが重要となる．現在，界面強度の評価方法としてスクラッチ試験 (1) や引張試
験 (2) ，テープ試験 (3) など様々な方法で界面強度評価が行われている．一方，アルミニウム合金に施したニッケルめっきへのボール摺動試験では，めっき後の熱処理によってはく離強度が上がることがわかっているものの，未処理材と比べてめっき界面の断面 TEM 像や XRD による結晶化度，XPS による成分分析などを行っても違いがわからず (4) ，界面のよりミクロな情報が必要とされている．
原子レベルでの界面強度の評価法として，原子スケールの微小試験片に対して精密な荷重を加え，破壊挙動を観察する研究 (5)−(9) などがされつつあるが，一般的には極めて困難である．そこで分子動力学法や第一原理計算などの計算材料科学的なアプローチが有効になる．分子動力学法 (Molecular Dynamics，MD) は，原子間の相互作用を簡単な関数を用いて表現し，ニュートンの運動方程式に基づいて原子の動きを追う方法である．一般的に実験のデータを基に原子間相互作用を計算しているため経験的手法とも呼ばれる．一方で第一原理計算は，原子の位置と種類を必要なデータとし，量子力学に基づいて電子レべルで原子間相互作用を計算する手法である．実験によるデータを必要としないことから非経験的手法と呼ばれる．異種金属
界面に関する研究として，Henz らは Ni，Al ナノ粒子の動的焼結プロセスを MD シ

1

ミュレーションにより検討している (10) ．Murch の研究グループでは，Al でコーティングをした Ni ナノ粒子の合金化反応の MD シミュレーションを行っている (11)(12) ．

A.Liu らは，銀の電気めっきに関して第一原理計算と MD による研究を行い，銀と銅表面への 5, 5-ジメチルヒダントイン (DMH) とニコチン酸の吸着挙動を報告している (14) ．W. Liu らは NiAl/Cr 界面の最適形状，熱力学的性質，電子構造を第一原
理計算により調べている (15) ．第一原理計算による界面の研究は他に，窒化アルミニ
ウムとアルミニウム界面の構造 (16) ，金属窒化物の表面エネルギーによる金型コー
ティング膜の離型性予測 (17) ，クロム層/シランカップリング剤界面の分子構造と接着強度の関係 (18) ，などがなされている．
我々の研究グループでは，Al 合金への Ni めっきを想定した Ni-Al 界面の MD シミュ
レーションを行っており，その結果 Ni-Al 界面が本質的に強固であること，表面エネルギーと弾性係数の差から界面ではなく Al 相側で必ず破断することなどを明らかにしている (19) ．一方で，実際のめっきでは P，Zn など種々の元素が含まれるが，MD シミュレーションでは適切なポテンシャルがなくその効果を議論することは難しい．また界面に存在する種々の元素の界面強度への影響をまとめた論文は少ない．
そこで本研究では，めっきの密着性予測に繋がることを目的として，第一原理計算
により様々な金属について表面エネルギーおよび界面エネルギーの評価や第 3 元素
による影響を検討した．

2

第 2 章 第一原理計算の概要
第一原理計算 (First principles calculation, Ab-initio calculation) とは，なんら実験データを参照せずに，対象とする物質の電子状態を原子番号と原子核の空間的配置を指定することのみで求めようとする解析手法である．実験で決めた原子間ポテ
ンシャルを用いないという意味で非経験的方法とも呼ばれる．そしてこの第一原理計算によって得られる電子状態から，エネルギー，原子に働く力，セルに働く応力などの諸物理量を高精度かつ定量的に求めることが可能となる．第一原理計算は大きく分けて，計算するモデルのサイズによってバンド計算とクラスター計算に分類される．バンド計算は結晶の周期性を利用して波数ベクトル空間で電子状態を解く方法である．それに対し，クラスター計算は有限サイズの原子集団の電子状態を実空間で解く方法であり，例えば分子軌道法などが挙げられる．固体材料の特性評価には主として前者のバンド計算が用いられる．本章では，第一原理バンド計算手法
として，局所密度汎関数法に基づく平面波基底疑ポテンシャル法による第一原理計
算手法について概説する．まず基礎として，一般的に広く用いられているノルム保存型擬ポテンシャルを用いた場合の系のエネルギー等の定式化について述べる．その後，ノルム保存型を用いた場合の定式化についても述べる．最後に，電子状態計
算の高速化手法についても述べる．

2.1

原子単位

バンド計算では，一般的に原子単位が用いられる．原子単位には Hartree 単位とRydberg 単位がある．本章では Hatree 単位を用いる．Table. 2.1.1 に Hatree 原子単位と SI 単位の関係を示す．

3

Table 2.1.1

Hatree atomic units and SI unit
HatreeSI
MASS
1.0
9.1095 × 10−31 kg
Length1.0
5.2918 × 10−11 m
Time
1.0
2.4189 × 10−17 s
Velocity
1.0
2.1877 × 106 m/s
Energy1.0
4.3598 × 10−18 J
Force
1.0
8.2388 × 10−8 N
Stress1.0
2.9417 × 1013 Pa

2.2断熱近似と平均場近似

通常，我々が扱う系は多数の原子核と電子からなる集合体である．そして電子間，原子核間，および電子と原子核との間の相互作用は多体問題であり，一般的に解くことができない．このような複雑な問題を実際に解くことが可能な問題へと帰着するために，通常，以下の 2 つの基本的な近似が導入される．

(a) 断熱近似
原子核は電子と比較すると非常に重く，電子よりもずっとゆっくりと運動する．このため，ある瞬間での原子配置に対して電子が速やかに基底状態をとると仮定する
ことができる．これを断熱近似 (Born-Oppenheimer 近似) という．この近似により，原子核は電子から見ると単なる外部のポテンシャル場とみなされ，原子系と電子系を独立に扱うことができる．

(b) 平均場近似
電子間相互の運動には Pauli の禁制による制約があり，またクーロン相互作用によって互いに避けあいながら運動するため，多電子系の運動を厳密に取り扱うことはき
わめて困難である．そこで，電子間の多体相互作用を一電子が感じる平均的な有効ポテンシャルで置き換える．この近似を平均場近似といい，バンド計算では通常，密度汎関数法が用いられる．

4密度汎関数理論

2.3

Hohenberg と Kohn は，外場ポテンシャル v(r)(原子核からの電場) 中における多電子系 (N 電子系) の基底状態の全エネルギー Etot が電子密度 ρ(r) の汎関数として∫Etot [ρ] =1 ∫∫ ρ(r′ )ρ(r) ′v(r)ρ(r) dr + T [ρ] +dr dr + E[ρ]2|r′ − r|(2.3.1)と表せることを明らかにした (20) ．右辺の各項はそれぞれ，原子核による電子ポテン
シャルエネルギー，相互作用する多電子系での電子の運動エネルギー，電子間クーロン相互作用エネルギー，他の全ての電子間多体相互作用を表す交換相関エネルギー
である．この Etot を最小にする ρ(r) が基底状態での電子密度分布となる．相互作用のない系での電子の状態を表す波動関数 (電子波動関数) を ψi とし，その運動エネ
ルギー Ts を

Ts [ρ(r)] =

occ
∑
i

1
< ψi | − ∇2 |ψi >2(2.3.2)と書くと，式 (2.3.1) は∫Etot [ρ] =1 ∫∫ ρ(r′ )ρ(r) ′dr dr + Exc [ρ]v(r)ρ(r) dr + Ts [ρ] +2|r′ − r|Exc [ρ] = T [ρ] − Ts [ρ] + E[ρ]∑occのように書ける．ここで，i(2.3.3)(2.3.4)は電子が占有している準位についての和をとることを表す．Exc は一電子近似のもとでの交換相関エネルギーであり，電子間相互作用を考慮した電子の運動エネルギー T [ρ] から，相互作用のない電子の運動エネルギー

Ts [ρ] を分離することによって，電子間の複雑な相互作用を全てこの項に押し込めている．
∫

電子密度に関する拘束条件 ρ(r)dr = N のもとで式 (2.3.3) に変分原理を適用することにより，以下の一電子シュレディンガー方程式 (Kohn–Sham 方程式) が得られる (21) ．
{}

1− ∇2 + veﬀ (r) ψi (r) = εi ψi (r)2(2.3.5)5ここで， veﬀ (r) は有効一電子ポテンシャルであり次式となる．
∫

veﬀ (r) = v(r) +ρ(r′ )
δExc [ρ]dr′ +′|r − r|δρ

(2.3.6)

第 2 項は電子間クーロン相互作用項，第 3 項は交換相関項である．
電子密度分布 ρ(r) は (2.3.5) 式の解から

ρ(r) =occ∑|ψi (r)|2(2.3.7)iとなる．以上のようにして，多電子問題は式 (2.3.5)∼(2.3.7) を Self–Consistent に解く問題に帰着される．

2.4

局所密度近似Kohn-Sham 方程式における，交換相関ポテンシャル ((2.3.6) 式第 3 項) には，多電子系を一電子近似したことによる複雑な相互作用が押し込められており，その汎関数の厳密な表現はわかっていない．そこで，電子密度の空間変化が十分緩やかであると仮定して，外場ポテンシャルが一定である一様電子ガスの交換相関エネルギー
密度 εxc を用い，
∫





εxc (r) ρ (r) drdδExc [ρ (r)]µxc (r) == εxc (r) + ρεxc (r) δρdρ

Exc [ρ] =(2.4.1)として計算する．つまり，電子密度 ρ(r) の点 r における交換相関エネルギーを同じ電子密度の一様電子ガス中のそれで代用する．これを局所密度近似 (Local Density

Approximation:LDA) という．
この εxc (r) の関数形についてはいくつか提案されている．以下に Perdew と Zungerの関数形 (22) を示す．

εxc (r) = εx + εcεx (r∗) = −(2.4.2)0.4582rs(2.4.3) −0.1423(rs ≥ 1)√εc (r) =  1 + 1.0529 rs + 0.3334rs
−0.0480 + 0.0311 ln rs − 0.0116rs + 0.0020rs ln rs (rs ≤ 1)
6

(2.4.4)

ここで，
(

rs =

3 1
4π ρ

)1
3

(2.4.5)である．交換相関ポテンシャル µxc は式 (2.4.1) よりµxc (r) = µx + µc

(2.4.6)

4µx (r) = εx3(2.4.7)

[

1



−0.1423√
1 + 1.0529 rs + 0.3334rs


(

)]

rs
1.0529(2.4.8)+1+(rs ≥ 1)√2
3(1 + 1.0529 rs + 0.3334rs )0.6668rs−0.0584 + 0.0311 ln rs − 0.0084rs + 0.00133rs ln rs (rs ≤ 1)

µc (r) = 

となる．

2.5

逆格子空間

第一原理バンド計算では，逆格子空間が用いられる．実空間における格子点の位置ベクトル R が，基本並進ベクトル a1 ,a2 ,a3 によって

R = n1 a1 + n2 a2 + n3 a3(n1 , n2 , n3 は整数)

(2.5.1)と表されるとすると，逆格子空間の基本並進ベクトル b1 ,b2 ,b3 は
a2 × a3 b1 = 2π
a1 · a2 × a3


a ×a

b = 2π31(2.5.2)2
a1 · a2 × a3

a1 × a2 b3 = 2π

a1 · a2 × a3と定義される．これらのベクトル b1 ,b2 ,b3 によって表される

G = m1 b1 + m2 b2 + m3 b3(m1 , m2 , m3 は整数)

(2.5.3)を位置ベクトルとする点の集合が逆格子であり，

G · R = 2π(m1 n1 + m2 n2 + m3 n3 )(2.5.4)7を満たす．結晶の並進対称性から，波動関数 ψ(r) と ψ(r + R) は同じ固有値をとる
関数となり，

ψ(r + R) = λψ(r) (|λ| = 1)

(2.5.5)の関係を満たす．式 (2.5.5) は Bloch の定理 (23) より

ψ(r + R) = exp(ik · R)ψ(r)(2.5.6)のように表される．ここで，k は波数ベクトル

k=

h1
h2
h3
b1 + b2 + b3
n1
n2
n3

(h1 , h2 , h3 は整数)

(2.5.7)である．式 (2.5.6) において，k → k + G としても (2.5.4) 式より同様に成立する．
したがって，G を全空間，つまり m1 , m2 , m3 を全ての整数についてとれば，k 点はG = 0 を中心とした Brillouin ゾーン (逆格子点を中心に近接する逆格子点へのベク
トルの垂直二等分線面で囲まれた空間) に限ってよいことになる．以上より，平面波
基底の第一原理計算では，無限の原子数の固有値問題を系の周期性により Brillouin
ゾーン内の各 k 点ごとの固有値問題に置き換えることができる．

2.6

ハミルトニアンk ベクトルについて n 番目の固有値をもつ波動関数 ψkn (r) を平面波で展開し，∑

Ψkn (r) =n
Ck+G| k+G>(2.6.1)Gと表す．ここで，|k + G >=1
exp[i(k + G) · r]Ω

(2.6.2)であり (Ω は全結晶体積)，規格直交条件

1∫exp [−i (k + G) · r] exp [i (k + G′ ) · r] dr<k+G | k+G > =Ω Ω
1∫=exp [i (G′ − G) · r] drΩ Ω′(2.6.3)= δGG′
8を満たす．式 (2.6.1) 中の ΣG は無限個の G についての和を表すが，実際の計算では
平面波の運動エネルギー |k + G|2 /2 がある一定の値 Ecut 以下のものについてのみ計算を行う．Ecut はカットオフエネルギーと呼ばれる．電子密度は

ρ (r) =

occ ∑
BZ
∑
n

=

fn fk |Ψkn (r)|2

k

occ ∑
BZ
∑∑∑
G

G′

n

1 n∗
′
n
fn fk Ck+G
′ Ck+G exp [i (G − G ) · r]Ωk

(2.6.4)で与えられる．ただし fn , fk はそれぞれエネルギー準位 n の占有数，k 点の重み付け∑BZ因子であり，k

は Brillouin ゾーン内の k 点についての和をとることを表す．以上のように平面波を基底関数として波動関数を展開すると，Kohn–Sham 方程式 (2.3.5)は次のように展開係数を固有ベクトルとする行列固有値問題となる．
∑

∑
1
n
n
= εkn
< k + G| − ∇2 + veﬀ |k + G′ > Ck+G
< k + G|k + G′ > Ck+G
′
′
2
G′
G′

=⇒

∑
G′

n
n
Hk+G,k+G′ Ck+G
= εkn Ck+G′

(2.6.5)以下にハミルトニアン行列要素 Hk+G,k+G′ =< k + G| − 12 ∇2 + veﬀ |k + G′ > の具体的な表現を示す. なお，各項の式変換において，|k + G′ > の具体的な表現を示す.
なお，各項の式変換において，

1∫
f (r) exp[−iG · r] exp[iG′ · r]drΩ Ω
1∫f (r) exp[−i(G − G′ ) · r]dr=Ω Ω

< k + G|f (r)|k + G′ > == f (G − G′ )(2.6.6)を用いる．

(a) 運動エネルギーの項
運動エネルギーの項は1
1< k + G | − ∇2 | k + G′ >= |k + G|2 δGG′
2
2

(2.6.7)となる．
一方，式 (2.3.6) に示したように veﬀ は原子核からのクーロン相互作用項 (v)，電子間クーロン相互作用項 (Vcoul )，交換相関項 (µxc ) からなる．平面波基底バンド計

9算では結晶結合に重要な役割を果たす価電子のバンド構造を効率的に計算するため，原子核からのクーロン項のかわりに内殻電子と原子核を正電荷をもったひとつのポ
テンシャルとして扱う擬ポテンシャル法が用いられることが多い．擬ポテンシャル法を用いることにより，膨大な平面波数を必要とする内殻電子の波動関数を直接扱うことなく価電子状態を正確に表すことができる (24)(25) ．擬ポテンシャルは 2.9 節
PPで後述するように，電子の角運動量に依存しない局所擬ポテンシャル Vloc,lと，依存PPする非局所擬ポテンシャル Vnloc
からなり，次式で表される．
PP
PP
VlPP (r − Ra ) P̂l = Vloc(r − Ra ) + Vnloc,l(r − Ra ) P̂l(2.6.8)

ここで，P̂l は角運動量 l への射影演算子，Ra は原子核の座標である．

(b) 局所項
局所擬ポテンシャルの行列要素は，PP< k + G |Vloc(r) | k + G′ >=1∫
V PP (r) exp [−i(k + G) · r] exp [i(k + G′ ) · r] drΩ Ω loc=

PP
Vloc(G − G′ )(2.6.9)である．結晶全体の局所擬ポテンシャルは格子周期関数であり，周期セル内の原子a からの距離 r に対する局所擬ポテンシャル VaPP,loc (r) を用いて
PP
Vloc
(r) =

∑∑
R

VaPP,loc (|r − ra − R|)(2.6.10)ra

PP(G) は以下より与えられる．
と表せることから，Vloc
PP
(G) =
Vloc

1 ∑
exp[−iG · ra ]VaPP,loc (G),Ωat a∫

VaPP,loc (G)=VaPP,loc (r) exp[−iG · r]dr
∫

= 2π
=

VaPP,loc (r) exp[−i|G|r cos ω]r2 sin ωdrdω

4π ∫ PP,loc(r)r sin(|G|r)dr
Va
|G|
10

(2.6.11)ここで，Ωat は周期セルの体積，ra はセル内の原子 a の位置ベクトル，R はセルの位置ベクトル，ω は G と r の間のなす角度である．

(c) 非局所項
非局所項の行列要素は，角運動量 l をもつ電子に対する原子 a からの非局所擬ポPP,nloc
テンシャル Va,l(r) により，PP< k + G|Vnloc(r)|k + G′ > =1 ∑exp[−i(G − G′ ) · ra ]VaPP,nloc (k + G, k + G′ )
Ωat a

PP
= Vnloc(k + G, k + G′ )VaPP,loc (k + G, k + G′ )
=

4π

∫

∑

(2l + 1)Pl (cos ω)(2.6.12)PP,nlocVa,l(r)jl (|k + G|r)jl (|k + G′ |r)r2 drl(2.6.13)となる (26) ．ここで，Pl は Legendre 多項式，jl は球 Bessel 関数であり，ω は k + Gと k + G′ との間の角度である．

(d) クーロンポテンシャルの項
電子密度分布 ρ(r) も格子周期関数であるのでフーリエ級数展開でき，

ρ(r) =∑

ρ(G) exp[iG · r]

(2.6.14)1∫
ρ(r) exp[−iG · r]Ω

(2.6.15)G

ρ(G) =となる．したがって，電子間クーロン項は Poisson 方程式 ∇2 Vcoul (r) = −4πρ(r)より，∇2 Vcoul (r) = −4π∑

ρ(G) exp[iG · r]

(2.6.16)

Gとなる．これを解いて，

Vcoul (r) = 4π

∑ ρ(G)
G|G|2

exp[iG · r]

(2.6.17)11

が得られる．これより，Vcoul (r) のフーリエ成分は

1∫
Vcoul (r) exp[−iG · r]drΩ Ω
∑ ρ(G′ )1∫′=4π′ 2 exp[iG · r] exp[−iG · r]dr
Ω Ω
|G
|
′
G

Vcoul (G) == 4π∑ ρ(G′ ) ∫
G′

= 4π

′ 2

|G |

1
exp[i(G − G′ ) · r]dr
Ω Ω

ρ(G)
|G|2

(2.6.18)であるから，電子間クーロン相互作用項のハミルトニアン行列要素は

1∫
Vcoul (r) exp [−iG · r] exp [iG′ · r] drΩ Ω1∫=
Vcoul (r) exp [−i (G − G′ ) · r] drΩ Ω< k + G |Vcoul (r) | k + G′ > == Vcoul (G − G′ )(2.6.19)となる．

(e) 交換相関ポテンシャルの項
交換相関項 µxc (r) も同様にフーリエ展開すると，µxc (r) =∑µxc (G) exp [iG · r]

(2.6.20)G1∫µxc (r) exp [−iG · r] drµxc (G) =Ω(2.6.21)となる．したがってハミルトニアン行列要素は (2.6.19) 式と同様に

1∫µxc (r) exp [−iG · r] exp [iG′ · r] drΩ Ω1∫µxc (r) exp [−i (G − G′ ) · r] dr=Ω Ω< k + G |µxc (r) | k + G′ > == µxc (G − G′ )となる．
以上により，ハミルトニアン行列要素は，

Hk+G,k+G′ =1|k + G|2 δGG′ + Vloc (G − G′ ) + Vnloc (k + G, k + G′ )2
+Vcoul (G − G′ ) + µxc (G − G′ )と逆空間での表式となる．

12(2.6.22)

2.7

系のエネルギー

全エネルギー Etot は，核（イオン）間相互作用エネルギー Ewald を加えて，Etot =∫
1∫εkn −Vcoul (r) ρ (r) dr + {εxc (r) − µxc (r)} ρ (r) dr + EEwald
2
n

BZ ∑
occ
∑
k

(2.7.1)と表される．εkn は式 (2.6.5) の固有値であり，Ewald は核間相互作用エネルギー（イオン間静電ポテンシャルエネルギー）を Ewald の方法 (27) によって表したもので，[

1 ∑ ∑ a a′ ∑4π
− |G|2
′
exp
[iG
·
(r
−
r
)]
exp
=
Zv Zv
a
a
2
2 a a′
4γ 2G̸=0 Ωat |G|

EEwald+−]1 ∑ ∑ a a′ ∑ erfc (|R + ra′ − ra | γ)dZv Zv
2 a a|R + ra′ − ra |R
∑ Zva2 γ

Z 2π
2πZ 2
√ −
+
lim
π
2Ωat γ 2 G→0 Ωat |G|2a(2.7.2)である．
ここで

ρ(r) =∑

ρ(−G) exp[−iG · r]

(2.7.3)

Gという関係を用いるとEtot =
+

BZ
occ
∑
∑
∑
1∑
n
PP
fk
fn
|k + G|2 |Ck+G|2 + Ωat
Vloc(G)ρ(−G)
2 k
n
G
G
BZ
∑
k

+

fk

occ
∑
n

fn

∑∑

G G′

′
n∗
n
PP
Ck+G
Ck+G
′ Vnloc (k + G, k + G )∑∑1
εxc (G)ρ(−G) + EEwald
Vcoul (G)ρ(−G) + Ωat
Ωat
2
G
G

(2.7.4)とフーリエ成分により表現できる．
PP(G) と Vcoul (G) は G = 0 で発散するが，これらの式 (2.6.11),(2.6.18) より，Vloc

発散成分は EEwald の発散項とうまく打ち消し合うため，次式のように表すことができる (28) ．

13

Etot =

BZ
occ
∑
∑
∑
2
1∑
n
PP
fk
fn
|k + G|2 Ck+G
+ Ωat
Vloc
(G) ρ (−G)
2 k
n
G=0
G̸=0

+

BZ
∑

fk

occ
∑

k

fn

n

∑ ∑
G=0 G′ =0

′
n∗
n
PP
Ck+G
Ck+G
′ Vnloc (k + G, k + G )∑∑ αa Z∑
1′
Vcoul (G) ρ (−G) + Ωat
εxc (G) ρ (−G) + EEwald
+
+ Ωat
2
a Ωat
G=0
G̸=0

(2.7.5)
ここで，EEwald は，式 (2.7.2) の第 5 項の発散成分を取り除いたものである．

2.8

応力

スーパーセルの平均応力 σαβ は，式 (2.7.5) に対称なひずみテンソル εαβ を用いて

r → (I + ε)r というスケーリングを適用し，それを対応するひずみテンソルの成分で微分することによって得られる (29)(30) ．Ωat ρ(G) や構造因子

Sa (G) = exp(−iG · ra )(2.8.1)はスケーリングの元のもとで不変であるから，平均応力は

∂Kγ
∂εαβ

= −δαγ Kβ( Kγ = (k + G)γ )∂Ωat= −δαβ Ωat∂εαβ(2.8.2)(2.8.3)という関係を用いることにより，

σαβ =1 ∂EtotΩat ∂εαβ

BZ ∑
occ ∑
1 ∑
n
|2 (k + G)α (k + G)β= −fk fn |Ck+GΩat k n G−

1 ∑∑∂V PP,loc (G)Sa (G){ a2Gα Gβ + VaPP,loc (G)δαβ }ρ(−G)
2Ωat G a∂(G )BZ ∑occ ∑ ∑ ∑ ∑
1 ∑
n
n∗
+
Ck+G
fk fn Sa (G − G′ )Ck+G′Ωat k n G G′ l a14∂
1 PP,nloc{V(k + G, k + G′ )}∂εαβ Ωat a,l2Gα Gβ
1∑
Vcoul (G)ρ(−G){− δαβ }
+
2 G
|G|2

×

+ δαβ

∑

(εxc (G) − µxc (G))ρ(−G)

G

+

Z ∑
1 ∂EEwald− δαβ 2αaΩat ∂εαβΩat a(2.8.4)と表すことができる．

2.9

擬ポテンシャル法

ブロッホの定理 (23) により，固体中の電子の波動関数は平面波基底により展開が可能である．しかし，平面波基底では原子核に強く引き付けられて局在している内殻電子の波動関数や，価電子密度の著しい変動を表現するには非常に多くの展開項数を要する．平面波数は解くべきハミルトニアンの次元数に比例し直接計算量に影響するので，これをできるかぎり少なくすることが望ましい．通常の固体材料では，内殻電子は原子核に強く引き付けられており，他の原子からの影響をほとんど受け
ず価電子がその特性を決定付けているといえるので，内殻電子と原子核をひとつのイオンと考え，原子間領域の価電子のみを取り扱うのが擬ポテンシャル法である．擬
ポテンシャル法は，その歴史の初期においては原子核付近で強い反発作用が現れた
り，原子核領域において真の波動関数と擬波動関数の 2 乗のノルムが一致していな
かったりしたため，self-consistent な計算には適用できなかった．そこで，Hamman
らは，これらの問題を解決した HSC 型 (BHS 型) と呼ばれるノルム保存型擬ポテン
シャルを開発した (24) ．しかし，第二周期元素や遷移金属では依然として非常に多くの平面波数が必要であったため，Troullier らはそれらの元素においても比較的少ない平面波数で扱える TM 型擬ポテンシャルを開発した (25) ．また，Vanderbilt らはノルム保存条件をはずすことにより，さらに少ない平面波数で計算を行えるウルトラ
ソフト型擬ポテンシャルを開発した (31) ．
本節では，まずノルム保存型擬ポテンシャルとして TM 型を説明する．その後，

15ノルム非保存型擬ポテンシャルとしてウルトラソフト型とそれを用いた場合の系のエネルギー等について説明する．

2.9.1

TM 型擬ポテンシャル

TM 型擬ポテンシャルは，まず擬波動関数の解析関数形を仮定し，これにノルム保存条件と少ない平面波数で収束させるための条件を課すことによりポテンシャルを構築する．以下にその手順を述べる．

1. まず，密度汎関数理論に基づき，孤立した原子に対して全電子計算を行う．具
体的には次式で表される動径方向の Kohn-Sham 方程式
[

]

1 d2
l (l + 1)−+ V (r) (rψnl (r)) = εnl (rψnl (r))
2
2 dr
2r2

(2.9.1)を解くことにより，各角運動量成分 l の動径方向の電子の感じる真のポテン
シャル VlAE (r) と真の波動関数 ψlAE (r)，および，その固有値 εAEnl を求める．

2. 内殻領域で節を持たない擬波動関数 ψlPP (r) を次式のような解析関数形で表す．
{
PP
ψscl,l(r) =ψlAE (r)(r ≥ rcl )lr exp [p(r)] (r ≤ rcl )(2.9.2)p (r) = c0 + c2 r2 + c4 r4 + c6 r6 + c8 r8 + c10 r10 + c12 r12(2.9.3)

ここで，rcl は角運動量 l に対する内殻領域の半径である．このようにおくと式
PP(r) が(2.9.1) より，価電子によって遮蔽 (screening) された擬ポテンシャル Vscr,l次式で表される．
AE Vl (r)PP(r) = Vscr,l εl +

′

′′

l+1 ′
p (r) + [p (r)]p (r) +
r2

(r ≥ rcl )2(r ≤ rcl )(2.9.4)

3. ここで，ノルム保存型擬ポテンシャルが満たすべき各種の条件を課す．
(a) ノルム保存条件16∫ rcl0|ψlPP (r)|2 r2 dr =

∫ rcl
0

|ψlAE (r)|2 r2 dr(2.9.5)

より，
[∫ r

]

cl

2c0 + lnr2(l+1)

0exp [2p(r) − 2c0 ] dr = ln

[∫ r

cl

0

|ψlAE (r)|2 r2](2.9.6)(b) 式 (2.9.4) の 2 次微分までが rcl で連続である条件[]

P(rcl )p(rcl ) = ln
rcll+1P′ (rcl ) l + 1
p′ (rcl ) =−P(rcl )
rcl

(2.9.7)(2.9.8)2(l + 1) ′2
p (rcl ) − [p′ (rcl )]rcl
2(l + 1) ′p (rcl )p′′′ (rcl ) = 2VlAE′ (rcl ) +rcl2
2(l + 1) ′′−p (rcl ) − 2p′ (rcl )p′′ (rcl )rcl
4(l + 1) ′p (rcl )p′′′′ (rcl ) = 2VlAE′′ (rcl ) −rcl2
4(l + 1) ′′2(l + 1) ′′′+p (rcl ) −p (rcl )2
rcl
rcl(2.9.10)−2 [p′′ (rcl )] − 2p′ (rcl )p′′′ (rcl )(2.9.11)p′′ (rcl ) = 2VlAE (r) − 2εl −

2(2.9.9)

ここで，′ は r による微分を表し，P (r) ＝ rψlAE (r) である．
′′

PP
PP
(c) Vscr,l(r) の r = 0 における曲率が 0 である条件 (Vscr,l(r) = 0)c22 + c4 (2l + 5) = 0(2.9.12)

4. これらの非線形連立方程式を解く．まず c2 を仮定し，式 (2.9.12) から c4 を決める．残りの 5 個の係数は式 (2.9.7)〜式 (2.9.11) の連立一次方程式であり，ガウス消去法により求める．最後に求まった係数を用いて c2 が妥当であるか式(2.9.5) により判断する．c2 の決定には bisection 法を用いる．
5. 以上により求まった擬ポテンシャルから，価電子による遮蔽効果を取り除くことにより内殻電子を含めたイオンの裸のポテンシャルを得る．

17

PP
PP
PP
Vion,l(r) = Vscr,l(r) − Vcoul(r) − µPP
xc (r)(2.9.13)

PP
ここで，Vcoul(r) はクーロンポテンシャル，µPP
xc (r) は交換相関ポテンシャルである．

6. 擬ポテンシャルを局所成分と非局所成分に分解する．
PP
PP
Vion,l(r) = Vion,loc
(r) +

∑

Vnloc,l (r)P̂l(2.9.14)l

ここで，P̂l は角運動量 l への射影演算子である．

擬ポテンシャルの KB 分離型表現
平面波展開による第一原理分子動力学法では，大きなハミルトニアン行列を繰り返し解く必要があるため，その繰り返しの中で変化しない量はメモリー上に記憶しておくことが高速化の基本となる．特に式 (2.6.13) の非局所項は，平面波の 2 乗のループを含んでおり計算時間がかかるとともに，記憶する量も平面波数の増加に対してその 2 乗で増える．そのため，大規模な計算ではすぐにメモリー容量に破綻をきたす．そこで，非局所項に次式で表される KB 分離型表現 (32) を用いれば，平面波の 2 乗のループは 1 乗のループとなる．
KB
Vnloc,lk(r) =

PP|Vnloc,l(r)ψlPP (r) >< ψlPP (r)Vnloc,l (r)|P̂l< ψlP P (r)|Vnloc,l (r)|ψPPscl,l (r) >(2.9.15)

これを用いると，行列要素の非局所項は，
KBk(r)|k + G′ > =< k + G|Vnloc,l∑ (4π)2
l

×
×

ΩCl

{∫ ∞

{∫0∞
0
l
∑

×

}

PP
ψlPP (r)Vnloc,l(r)jl (|k + G|r)r2 dr}PP(r)jl (|k + Gd|r)r2 dr
ψlPP (r)Vnloc,l∗Ylm (k + G)Ylm(k + G′ )m=−l(2.9.16)となる．ここで，18

Cl =< ψlPP (r)|Vnloc,l (r)|ψlPP (r) >Vnloc,l (r) =∑∑R(2.9.17)

PP,nloc(|r − ta − R|)Vl,a(2.9.18)aである．したがって，PP,nloc
PPPP(r) >(r)|Vl,a(r)|ψlaCla =< ψlaAla (k + G) =∫ ∞0(2.9.19)

P P,nlocPPψla(r)Vl,a(r)jl (|k + G|r)r2 dr(2.9.20)とおくと，
KB< k + G|Vnloc,l(r)|k + G′ > =×(4π)2 ∑ ∑ 1Ωat l a Clal∑{exp[−iG · ra ]Ala (k + G)Ylm (k + G)}

m=−l∗
× {exp[iG′ · ra ]Ala (k + Gd)Ylm(k + G′ )}(2.9.21)と書ける．平面波展開係数との積は，
∑
G′

KB
n
< k + G|Vnloc,l(r)|k + G′ > Ck+G′=(4π)2 ∑ ∑ 1Ωat l a Cla×l∑{exp[−iG · ra ]Ala (k + G)Ylm (k + G)}

m=−l

{

×

∑
G′

}
′
′
′
n
∗
Ck+G
′ exp[iG · ra ]Ala (k + G )Ylm (k + G )(2.9.22)となり，

AYlam (k + G) = Ala (k + G)Ylm (k + G)∗(2.9.23)

をあらかじめ記憶しておけば計算が速くなる．また，この行列要素を計算した際に，

CAYnalkm (k + G) =

∑

′
′
′
n
∗
Ck+G
′ exp[iG · ra ]Ala (k + G )Ylm (k + G )Gを記憶しておけば後のエネルギーや原子に働く力の計算が高速化できる

19(2.9.24)

2.9.2

ウルトラソフト型擬ポテンシャル

Vanderbilt らは，擬ポテンシャルの作成時にノルム保存条件をはずすことによってさらなるソフト化を達成したウルトラソフト型擬ポテンシャルを開発している．しかしながら，それをはずしたことによって生じるノルムのずれを補う計算が系の全エネルギーや電子密度等に必要となる．
全電子計算により求められた真のポテンシャルを VAE とすると，真のシュレーディ
ンガー方程式は，真の波動関数 Φi を用いて(T + VAE − εi )|Φi >= 0

(2.9.25)と書ける．ここで，r > rloc で VAE と一致するように局所ポテンシャル Vloc を r < rlocの領域で適当に決める．また，r > rcl で Φi と一致し，r < rcl で節を持たない擬波
動関数を Ψi とすると，擬波動関数の満たすべきシュレーディンガー方程式は以下のようになる．
′(T + Vloc + VNL− εi )|Ψi >= 0 ,

′
VNL
=

|χi >< χi |< χi |Ψi >(2.9.26)′
ここで，VNL
は非局所ポテンシャルであり，関数 χi は

|χi >= (εi − T − Vloc )|Ψi >(2.9.27)と定義する．χi は r > R = Max(rcl , rloc ) では 0 となる局在した関数である．非局所′
ポテンシャル VNLは次のように変形できる．
′
VNL= Σi,j Bij |βi >< βj |(2.9.28)

ただし

Bij =< Ψi |χj > ,|βi >= Σj (B −1 )ji |χj >(2.9.29)また，< Ψi |βj >= δij(2.9.30)である．擬ポテンシャルにノルムの保存条件を課さなかったことにより，内殻領域において電子密度が20

Qij (r) = Φ∗i (r)Φj (r) − Ψ∗i (r)Ψj (r)(2.9.31)だけ不足している．また求められた波動関数も，ノルムが
∫

Qij =

r<rcl

Qij (r)dr(2.9.32)だけ不足している．これを考慮して重なり積分演算子 S を

S =1+

∑

Qij |βi >< βj |(2.9.33)ijと定義すれば，規格直交条件が以下のように満足される．

< Ψi |S|Ψj >= δij(2.9.34)′
これを (2.13.12) 式に含めるためには，非局所ポテンシャル VNL
も変形を加える必要がある．よって，(T + Vloc + VNL )|Ψi >= εi S|Ψi >VNL =(2.9.35)∑(Bij + εj Qij )|βi >< βi |(2.9.36)

ij

となる．

2.10

電子占有数

金属では Fermi エネルギー εF の近傍に多くのエネルギー準位が存在するため，整数の占有値では問題が生じる (33) ．たとえば時間とともに Fermi エネルギー近傍の２つの準位が交差してしまうと，電子密度が不連続に変化してしまう．このような問題を避けるために，Gaussian Broadening(34) という方法を用い，fn のかわりに非整数の占有数 fi[(1
εi − εFfi =1 − erf2
σ)](2.10.1)を導入し．フェルミレベルに対して σ の幅で占有状態をぼかしてある程度の非占有状態も計算する．実際の数値計算では

21

2

∑

fi = Z(2.10.2)iとなるように εF を決定する．Z はセル内の総価電子数である．このとき，fi に関する自由度が増えるので，全エネルギー Etot のかわりに自由エネルギー Ef

Ef = Etot − T S
S = −kB(2.10.3)∑{fi ln fi + (1 − fi ) ln(1 − fi )}(2.10.4)iを考えなければならない．

2.11

FFT

n
固有方程式を解いて求めた固有値 Ck+G
を∑

ukn (G) =n
Ck+G(2.11.1)

Gとおけば，フーリエ逆変換より

ukn (r) =

∑

n
Ck+G
exp[iG · r]

(2.11.2)′
n∗Ck+G′ exp[−iG · r]

(2.11.3)

Gとなる．同様にu∗kn (r) =∑Gであるから，

ukn (r)u∗kn (r) =

∑∑
G

G′

′
n
n∗
Ck+G
Ck+G
′ exp[i(G − G ) · r]

(2.11.4)

したがって，式 (2.6.4) より

ρ(r) =

occ ∑
BZ
∑
n

1
fn fk (ukn (r)u∗kn (r))Ω
k

(2.11.5)22

となり電子密度分布が得られる．すなわちハミルトニアンから求められる固有ベクnトル Ck+Gをフーリエ変換することにより，実空間の電子密度分布 ρ(r) を式 (2.6.4)

に従って直接評価するより高速に計算できる．ρ(r) が求められれば交換相関エネル
ギー，交換相関ポテンシャルの実空間における値が得られ，フーリエ変換によって逆空間での値も求められる．このように実際の計算ではフーリエ変換を多用するため，一般に高速フーリエ変換 (Fast Fourier Transformation:FFT) のプログラムが用
いられる．

2.12

電子系の最適化手法

平面波基底による電子状態計算では，前節で定式化された Kohn-Sham 方程式をセルフコンシストに解くことによって固定した原子配置に対する電子の基底状態を求める．オーソドックスな収束計算手法は，ハミルトニアン行列 (式 (2.6.22)) の対角化を繰り返す方法であるが，この方法では対象とする系によっては多大な計算労力を必要とする．そこで，近年電子状態計算を効率的に行う方法が開発された (35)−(38) ．
本節では共役勾配法についてその概要を示す．

共役勾配法の原理
共役勾配法は，一般には正定な係数行列をもつ連立１次方程式を最適化の考えに立って解くために，あるいは，多次元空間の２次関数 F (X) の最小化問題を解くために用いられる計算手法である．共役勾配法では，前者の問題は結局後者の問題に帰着され，適当な初期値 X0 から出発して順次修正を加えながら · · · ，Xm−1 ，Xm ，Xm+1，· · · と変化させて F (X) を最小にする X を探索する．
密度汎関数法に基づく電子状態計算では系の全エネルギー Etot は，電子密度すな
わち波動関数の汎関数で表され正しい波動関数によって最小化される．したがって，
平面波基底の波動関数を用いた場合には，系の全エネルギーを最小にする係数ベク
トル Cnk を規格直交条件のもとで求める計算を行えばよい．
（ここで，Cnk は平面波

23n展開係数 Ck+G
を成分に持つベクトルである．
）すなわち，′Etot= Etot −∑λmn (< Ψmk |Ψnk > −δmn )mn= Etot −

∑

(

λmn

mn

∑

)
m∗
n
Ck+G
Ck+G
− δmn

(2.12.1)

Gの最小化を考える．ここで，λmn =< Ψmk |Ĥ|Ψnk >=

∑∑
G

G′

n
m∗
Ck+G
Ck+G
′ Hk+G,k+G′(2.12.2)

である．共役勾配法では，次式を残差ベクトル（G の数だけの成分を持つ）として各バンド n の各 k 点ごとに最適化を行う．
[

]

[

]

′
∑
∂Etot
n
Rnk = −=−
(Hk+G,k+G′ − λnn ) Ck+G
(G=G1 ,G2 ,···,Gmax )(2.12.3)
′ ,
m∗
∂Ck+G
′
G

以下に金属の電子状態計算において代表的な BKL 法 (39) について解説する．

BKL 法
BKL 法と並んで共役勾配法のもう 1 つの代表的な手法であり，全エネルギーの最小化を行う TPA 法 (40) は，絶縁体と半導体には有効であるが，金属には適さない．
これは，金属ではフェルミ面がぼやけるために非占有状態も考慮しなければならな
いことによる．このため，占有状態にしか依存しない全エネルギーを最小化する方法では適切な電子状態計算を行うことができない．そこで，BKL 法では占有状態と非占有状態の両方について計算できるエネルギー期待値 εkn =< Ψkn |H|Ψkn > の最小化を行う．したがって，BKL 法は，金属はもちろん絶縁体と半導体についても有効な方法である．
具体的な手法としては，まず波動関数の展開係数を成分とする係数ベクトルの残差ベクトルを求める．次に preconditioning という処理を施し，共役方向ベクトル (探索方向) を求める．それをもとにして εkn を最小にするような新たな係数ベクトルを求める．以上の手順を εkn が収束するまで繰り返した後に，電子密度とハミルトニ
アンの更新を行い全エネルギーを計算する．

< 残差ベクトル >24′Etot を εkn に置き換えることによって，式 (2.12.1) の Etot
は ε′kn に置き換わるとす

ると，残差ベクトルは式 (2.12.3) より次式で表される．
[
i
Rnk
=−

∂ε′kn
m∗
∂Ck+G

]i

(

)

i
= − H − λin I · Cnk

ただし，式中の i は， i 回目のステップにおける[(2.12.4)という意味を表し，

]

n,i
i
i
i
Ĥ Ψnk>= Ck+GH = [Hk+G,k+G′ ] , λin =< Ψnk
Cnk′ ,(2.12.5)である．これは，i のステップにおいて εkn を最小にする方向 (最急降下方向) を示す
ベクトルを表している．
iRnkには，最終的に得られる次のステップの波動関数 Ψi+1
nk が同じ k 点における n以外の全バンドの波動関数 Ψmk (m ̸= n) と直交するように，直交化処理が施される．
∑ (

i′
i
Rnk
= Rnk
−

)

i∗
i
Cmk
· R ink Cmk(2.12.6)m̸=n<preconditioning>′i
残差ベクトル Rnk
に対して preconditioning という処理を施す．大きな逆格子ベク

トルについては平面波の運動エネルギーが大きくなるが，このことが残差ベクトルに影響して収束性を悪化させる．preconditioning は，この問題を回避して収束を速i
めるために行われる．preconditioning された残差ベクトルを Gnk
とすると
i
i′
Gnk
= K i · Rnk(2.12.7)と表される．ここで(27 + 18x + 12x2 + 8x3 )(27 + 18x + 12x2 + 8x3 + 16x4 )Ekin (G)
x =
i
Ekin
1
Ekin (G) =|k + G|221iii| − ∇2 |Ψnk
= < Ψnk>
Ekin
2
KGG′ = δGG′

(2.12.8)(2.12.9)(2.12.10)(2.12.11)である．式 (2.12.8) は，経験的にそれがよいとされている式である．最後に直交化処理が施される．

25

(

)

i′
i
i∗
i
i
Gnk
= Gnk
− Cnk
· Gnk
Cnk
−

∑ (

)

i∗
i
i
Cmk
· Gnk
Cmk

(2.12.12)m̸=nii
ここで，Gnk
は Cnkと直交しなければならないことに注意が必要である．< 探索方向>
探索方向は，次のようにして定められる．
′

i−1
i
i
+ γ i Fnk
= Gnk
Fnk

γi =






(2.12.13)i′∗
Gnk· R ink′i−1 ′ ∗′
· R i−1Gnknk(i > 1) 0(2.12.14)(i = 1)

さらに，直交化処理と規格化処理を施す．
(

)

i′
i
i∗
i
i
Fnk
= Fnk
− Cnk
· Fnk
Cnk
i
Dnk
= (

i′
Fnk
i′∗
i′
Fnk
· Fnk

(2.12.15)(2.12.16))12< 新たな係数ベクトルの組み立て >新たな係数ベクトルの組立ては次のように行われる．
i
i
C i+1
nk = αC nk + βDnk

(2.12.17)

結合係数 α と β は，エネルギー期待値 εkn を最小化するように決定される．すなわ
ち，C ink ，Dink を基底とする 2 × 2 ハミルトニアン行列，[

C ink∗ HC ink C ink∗ HD ink
D ink∗ HC ink D ink∗ HD ink

]

[

=

ε11 ε12
∗
ε12
ε22

]

(2.12.18)

を組立て，この行列の小さい方の固有値 γ ，{ε11 + ε22
(ε11 − ε22 )2γ=−+ ε12 ε∗12
2
4

} 12

(2.12.19)に対応した固有ベクトルによって次式で与えられる．

α = {

ε12
ε12 ε∗12 + (ε11 − γ)2

β = −{

ε11 − γ(2.12.20)}12(2.12.21)}1

ε12 ε∗21 + (ε11 − γ)2 226

以上の手順を εkn が収束するまで繰り返せばよい．計算の全体的な手順を以下に示す．

1. 係数ベクトル Cnk の適当な初期値を，行列計算などによって，全 k 点の全状態について作成する．
i
i+1
2. 各 k 点の各状態について，Cnk
から Cnk
を組立てる一連の計算を反復し，適当な条件で打ち切る．打ち切り条件は，例えば，1 回のステップでの εkn の減少値が，最初のステップでの減少値の 30%以下や一定値以下になることである．反復計算が打ち切られれば，同じ k 点における次の状態についての計算へと移る．

3. 全 k 点の全状態について，1,2 の計算が終了したら，この時点で初めて電子密度とそれに伴うハミルトニアンの更新を行い，全エネルギーを求める．

4. 全エネルギーが収束すれば，計算を終了し，そうでなければ再び 1 ∼ 3 を行う．

27

2.13

ノルム非保存形擬ポテンシャル法

Vanderbilt らは，擬ポテンシャルの作成時にノルム保存条件をはずすことによってさらなるソフト化を達成したウルトラソフト型擬ポテンシャルを開発している．しかしながら，それをはずしたことによってノルムのずれを補う計算が系の全エネル
ギーや電子密度等に必要となる．
全電子計算により求められた真のポテンシャルを VAE とすると，真のシュレーディ
ンガー方程式は，真の波動関数 Φi を用いて(T + VAE − εi )|Φi >= 0

(2.13.1)と書ける．ここで，r > rloc で VAE と一致するように局所ポテンシャル Vloc を r < rlocの領域で適当に決める．また，r > rcl で Φi と一致し，r < rcl で節を持たない擬波
動関数を Ψi とすると，擬波動関数の満たすべきシュレーディンガー方程式は，以下のようになる．
′′(T + Vloc + VNL− εi )|Φi >= 0, VNL=|χi >< χi |< χi |Ψi >(2.13.2)関数 χi は

|χi >= (εi − T − Vloc )|Ψ >(2.13.3)と定義する．χi は r > R = Max(rcl , rloc ) では 0 となる局在した関数である．非局所′
ポテンシャル VNLは次のように変形できる．
′
VNL= Σi,j Bij |βi >< βj |(2.13.4)

ただし

Bij =< Ψi |χj >, |βi >= Σj (B −1 )ji |χj >(2.13.5)

また，< Ψi |βj >= δij(2.13.6)である．擬ポテンシャルにノルムの保存条件を課さなかったことにより，内殻領域において電子密度が28

Qij (r) = Φ∗i (r)Φj (r) − Ψ∗i (r)Ψj (r)(2.13.7)だけ不足している．また求められた波動関数も，ノルムが
∫

Qij =

r<rcl

Qij (r)dr(2.13.8)だけ不足している．これを考慮して重なり積分演算子 S を

S =1+

∑

Qij |βi >< βj |(2.13.9)ijと定義すれば，規格直交条件が以下のように満足される．

< Ψi |S|Ψj >= δij(2.13.10)′
これを (2.13.12) 式に含めるためには，非局所ポテンシャル VNL
も変形を加える必要がある．よって，′(T + Vloc + VNL)|Φi >= εi S|Φi >′VNL=(2.13.11)∑(Bij + εj Qij )|βi >< βi |(2.13.12)

ij

となる．

29

第 3 章 第 3 元素の表面エネルギーへの影響
本章では Fe，Al，Mg，Ni 表面に H，O，Si，P，Cr，Co，Cu，Zn，Mo，W が存在した際の表面エネルギーの変化について計算し，表面エネルギーの大小から密着性の予測を立てる．

3.1

解析条件

本研究での解析はすべて Kresse らにより開発された平面波基底ウルトラソフト擬
ポテンシャル法に基づく第一原理バンド計算コード VASP（Vienna Ab-initio Sim-ulation Package）を用いて行った．交換相関項には局所密度近似（Local Density
Approximation, LDA）に勾配を考慮した一般化密度勾配近似（Generalized Gradient Approximation, GGA）を用いた．また波動関数の収束計算には共役勾配法を採用した．カットオフエネルギー，バンド数，波動関数と電子密度の FFT メッシュ
は，原子種と電子数から VASP が定める値を用いた．Table 3.1.1 に対象とした元素の結晶構造，格子長さ，スーパーセル内の原子数を示す．また Fig. 3.1.1 に周期境界条件を設けた各結晶構造のバルクのスーパーセルを示す．これらのバルクモデルに対して構造緩和計算を行い，バルクモデルでの全エネルギー Ebulk を求めた．次にFig. 3.1.2 に示すようにバルクモデルの c 軸方向に 1nm 以上の真空層を設け，表面を 2 つ有するスラブモデルを作成し，構造緩和計算は行わずスラブモデルの全エネ
ルギー Eslab を求めた．バルクモデル，スラブモデルいずれの計算においても，逆格子空間の k 点メッシュは Monkhorst-Pack 法に従い 12 × 12 × 2 とした．表面エネル

30

ギー Es は以下の式で求める．

Es =Eslab − Ebulk
2S(3.1.1)

ここで S はスーパーセルの断面積である．
次に第 3 元素による表面エネルギーへの影響を考えるため，c 方向の原子面上で原子を第 3 元素に置換したバルクモデルおよびスラブモデルを作製した（Fig. 3.1.3，

Fig. 3.1.4）．bcc 構造の最小周期では (001) 面に 1 つしか原子が存在しないため，a ，
b 方向 2 格子分に拡張し，2 原子置換して 50%置換とした．他の構造では 2 原子ずつ（図では境界上の原子がイメージ粒子として表示されているため 5 原子）あるう
ちの 1 原子を置換して 50%とした．また図では上下のイメージ粒子も表示されている．それぞれのモデルの全エネルギーを前述と同様に求めた．Fig. 3.1.4 に示すようにスラブモデルでは第 3 元素のある面とない面が存在するので置換した表面のエネ
ルギー Es′ を次式で求めた．
′
′
Eslab
− Ebulk
′
− Es
Es =

(3.1.2)S

31

Material
Mg
Al
Fe
Ni
Si
P
Cr
Co
Cu
Zn
Mo
W

Table 3.1.1 Calculation conditions
Crystal structure Lattice constants(nm) Number of atoms
a=0.321hcp
20c=5.210
fcc
0.405
20
bcc
0.287
10
fcc
0.352
20
diamond0.543
40
a=0.331
black phosphorus20b=0.438
bcc
0.288
10
a=0.251hcp
20
c=0.407
fcc
0.361
20
a=0.266hcp
20
c=0.495
bcc
0.314
10
bcc
0.316
10

32

Fig. 3.1.1

Supercells for bulk model (a)bcc, (b)fcc, (c)hcp, (d)diamond, (e)black
phosphorus

Fig. 3.1.2

Supercells for slab model (a)bcc, (b)fcc, (c)hcp, (d)diamond, (e)black
phosphorus

33

Fig. 3.1.3

Supercells for bulk model including third elements (a)Fe, (b)Ni, (c)Mg，(d)AlFig. 3.1.4

Supercells for slab model including third elements (a)Fe，(b)Ni, (c)Mg,
(d)Al

34

3.2

解析結果と考察3.2.1

単元素の表面エネルギー

バルクモデルの構造緩和計算により得られた格子定数を Table 3.2.1 に示す．構造緩和計算により得られた格子定数はいずれも既知の格子定数と 1%程度の誤差である．次に表面エネルギーおよび文献値 (41)−(48) ，誤差を Table. 3.2.2 に示す．表面エ
ネルギーは文献値にもばらつきがあるが同程度の値を示している．Cr が文献値より低く，Co は多数の文献値に比べると高い．また P は c 軸方向は金属結合より弱い
ファンデルワールス力によって結ばれているため，金属原子より値が著しく小さい
表面エネルギーを示したものと考える．

Table 3.2.1 Lattice constants of various materials
Material Lattice constants(nm) Literature values(nm) Error(%)
a=0.317a=0.321a:-1.25c=5.20c=5.210c:-0.19

Al

0.402

0.405

-0.74

Fe

0.284

0.287

-1.05

Ni

0.351

0.352

-0.28

Si

0.545

0.543

0.37

a=0.328a=0.331

a:0.90b=0.449b=0.438

b:2.51

0.285

0.288

-1.04

a=0.249a=0.251

a:-0.80c=0.404c=0.407c:-0.73

0.362

0.361

0.28

a=0.267

a=0.266

a:0.38c=0.495c=0.495

c:0

Mo

0.315

0.314

0.32

W

0.317

0.316

0.32

Mg

P
Cr
Co
Cu
Zn

35

Table 3.2.2
Material

3.2.2(001) Surface energy of various materials(001)Surface energy

Literature values(J/m2 )

0.79,(41) 0.29,(45) 0.64(46)

Mg

0.53

Al

0.92

1.35,(41) 1.04(42) 1.04,(43) 0.67(44)

Fe

2.40

2.43(41)Ni

2.21

2.43,(41) 2.59,(42) 2.60,(43) 1.44(44)

Si

2.16

2〜2.5(48)

P

0.02

-

Cr

2.93

3.98(41)

Co

3.93

Cu

1.44

2.17,(41) 1.65,(42) 1.92,(43) 1.27,(44) 2.06(46)

Zn

0.34

0.99,(41) 1.23,(45) 0.24,(47) 1.06,(47) 0.62(47)

Mo

3.58

3.84,(41) 3.53(42)

W

4.36

4.64(41)

0.67,(47) 0.70,(47) 0.69(47)2.78,(41) 2.80,(43) 2.86,(45) 3.18(46)2.15,(47) 2.30,(47) 2.22(47)

第 3 元素による表面エネルギーへの影響

Fe，Al，Mg，Ni 単元系モデルの表面元素を 50%第 3 元素に置換した系の表面エネルギーを Fig. 3.2.1〜Fig. 3.2.4 にそれぞれグラフとして示した．いずれのグラフも左端の棒グラフが第 3 元素に置換する前の表面エネルギーであり，これより上ならば表面エネルギーが高く，下になれば低くなる効果がある．Mg 表面はどの元素に置換しても表面エネルギーが上昇し，特に Si，P，Cr，Co，Zn に置換すると表面エネ
ルギーの変化は著しく上昇する．Fe 表面も Si に置換すると著しく表面エネルギーが増加し，O，P，Cu，Zn で減少した．Al 表面は Cr，Co，Cu，Mo，W に置換すると表面エネルギーが増加し，O，Si でわずかに低下する．Zn ではほとんど変わらず，

36

P 添加で最も表面エネルギーが小さくなる．Ni 表面も添加元素によって表面エネル
ギーの増減が変わるが全体的に低下する傾向にある．H，O，P，Zn で表面エネル
ギーが著しく減少し，特に O 原子では 0 に近い値となった．Al と Ni は，Si 以外は単元素の表面エネルギーが Al と Ni より大きい元素（Cr，Co，Mo，W，Al については Cu も）に置換すると表面エネルギーが増加している．Si の表面エネルギーはAl より小さいが，Si 単元素の表面エネルギーはダイヤモンド構造に対するものであること，Si や P といった非金属元素は金属結合ではない結合形態となったためなどが考えられる．

Fig. 3.2.5〜Fig. 3.2.8 は，各元素の完全結晶のエネルギーから，元素置換を行った
バルクのエネルギーの変化（青色の棒グラフ）と，表面をつけたスラブモデル同士の，第 3 元素置換によるエネルギーの変化（橙色の棒グラフ）を示したものである．
例えば Fig. 3.2.5 の左端の H では，バルクモデルよりもスラブモデルの方が置換によるエネルギー増加が上回っているため表面エネルギーは無置換のそれに比べてわずかに上がる．また右端 W ではバルク，スラブともに置換によってエネルギーが下がるが，バルクモデルよりもスラブモデルの方が置換によるエネルギー減少量が小さいため表面エネルギーは無置換のそれに比べてやはり上がる．2 つの棒グラフの差が著しい元素は，第 3 元素が表面にあるか内部にあるかで効果が著しく異なることを意味する．これで第 3 元素がバルク内，表面のどちらで大きく影響を及ぼすか議論する．

Fe では，Mo，W 以外は Ebulk および Eslab が増加する．Si は表面に存在すること
で著しくエネルギーが上昇していることがわかる．表面エネルギーが大きければその表面を形成しにくくなることから (19) ，Si の添加で密着性が上がることが予測される．逆に P は橙色の棒グラフが青色の棒グラフに比べて半分以下であることから
表面エネルギーが減少し，P を含む面で破断しやすくなる．

Al では H と Zn 以外は Ebulk および Eslab を減少させる．Ebulk と Eslab の差が比較的大きいのが，先の Fig. 3.2.2 でエネルギーが増加した Cr，Co，Mo，W である．Cu
は表面でほとんどエネルギー変化がなく，バルクでエネルギー減少したため表面エ37

ネルギーが上昇した．

Mg では Si と H がバルク内でエネルギーを減少させるが表面ではエネルギーを増加させる，Fe と Al にない傾向を示した．これは Si，H が Mg 内部に拡散しやすい可能性を示している．また Zn は表面エネルギーを著しく増加させ，他の元素はバルク，表面ともにエネルギーを減少させる効果があるが減少量が後者で少ないために，
表面エネルギーは増加する．

Ni の場合，O はバルクのエネルギーの増加がスラブモデルの増加量より大きいため，内部への侵入がないものと予測される．P もスラブモデルでのエネルギー減少が大きいので Ni 内部に拡散しにくい．

Fig. 3.2.1

Change in surface enrgy of Fe by third element

38

Fig. 3.2.2

Change in surface enrgy of Al by third element

Fig. 3.2.3

Change in surface energy of Mg by third element

39

Fig. 3.2.4

Change in surface energy of Ni by third element

Fig. 3.2.5

Change in Ebulk and Eslab of Fe by third element

40

Fig. 3.2.6

Change in Ebulk and Eslab of Al by third element

Fig. 3.2.7

Change in Ebulk and Eslab of Mg by third element

41

Fig. 3.2.8

Change in Ebulk and Eslab of Ni by third element42

第 4 章 異種金属界面に存在する第 3 元素の影響
本章では，Fe，Al，Mg 上への Ni めっきを想定した計算を行う．Fe/Ni，Al/Ni，

Mg/Ni それぞれの界面に存在する第 3 元素による界面エネルギーの変化から密着性について議論する．

4.1

解析条件

3 章と同じく VASP を用いて解析を行い，諸条件も 3 章と同様にした．Fig. 4.1.1に Fe/Ni，Al/Ni，Mg/Ni のスーパーセルを示す．Fe/Ni 界面は (100) 面で接合させた正方形断面の角柱セルである．正方形断面一辺の長さは Fe の格子長さ 0.284nm お
よび Ni の格子長さ 0.351nm の平均値である 0.3175nm×2 としている．c 方向はセル
中央と上下境界の界面が同じになるように Fe は 5 格子分，Ni は 5.5 格子分とした．

Al/Ni 界面は断面 1 格子分の角柱セルで計算した．断面 1 辺の長さは Al の格子長さ0.402nm と Ni の格子長さ 0.351nm の平均値 0.3765nm としている．c 方向は Al，Ni
それぞれ 5 格子である．Mg/Ni 界面は (0001) 面と (111) 面を整合させた長方形断面の角柱セルを用いた．寸法は同様に平均値から設定した．c 方向は Mg10 層分，Ni 側
は 2 つの界面が同じになるように 9 層とした．計算量の制約からいずれも平均格子長さを用いざるを得ないため，各単結晶の全エネルギーは 3 章と異なる．そこでそれぞれの界面モデルに対応したひずみを与えたバルクモデルでの計算も行った．以上の界面モデル，バルクモデルの全エネルギー Eab ，Ebulk から以下の式より界面エネルギー Ei を求めた．

Ei =Eab − Ea − Eb
2S

(4.1.1)43

(4.1.2)
ここで Ea および Eb は接合させた元素 a，b のひずみを考慮したバルクエネルギー

Ebulk である．次に第 3 元素による界面エネルギーへの影響を調べるため，3 章と同様に，界面の 50%の原子を第 3 元素に置換したモデルを作成した（Fig.4.1.2）．第 3元素のある界面とない界面が存在するので以下の式より置換した界面のエネルギーEi′ を求めた．
Ei′ =

′
Eab
− Ea′ − Eb′− Ei
S(4.1.3)(4.1.4)

Fig. 4.1.1

Supercells for interface model (a)Fe/Ni，(b)Al/Ni，(c)Mg/Ni

44

Fig. 4.1.2

4.2

Supercells for interface model (a)Fe/Ni，(b)Al/Ni，(c)Mg/Ni

第 3 元素による界面エネルギーへの影響

Fig. 4.2.1〜Fig. 4.2.3 に Fe/Ni，Al/Ni，Mg/Ni の第 3 元素を含まない完全な界面の界面エネルギーおよび第 3 元素に置換した界面の界面エネルギーを示す．Fe/Ni 界
面では Mo，P を添加した場合以外は，正の値を示したことから凝着性を示さず不安定である．Mo，P は Fe/Ni 界面において界面の密着性をあげる効果があると予測できる．

Al/Ni 界面は第 3 元素がない界面で最も高い凝着性を示しており，MD シミュレー
ションで報告されている「本質的に強固」という事実に一致する．Si，P を除いて，第 3 元素が存在すると凝着性は著しく減少し，O，Mo では正の値を示した．Si，Pは Al バルク，Ni バルク中でのバルクエネルギーに与える影響が小さかったので界面の凝着性が維持されたと考えられる．一方 O，Cu は Al バルクに与える影響は少な45

かったが，Ni バルクに与える影響が大きいため凝着性が著しく低下した．

Mg/Ni 界面では，いずれの値も正値を示した．この理由として今回は様々な第 3元素による影響に着目し，構造緩和を行わなかったため，Mg と Ni 間の距離が適切ではなく真空層ができたことによる．したがって以降では Fe/Ni，Al/Ni 界面についてのみ議論する．

Fig. 4.2.4〜Fig. 4.2.11 に Fe/Ni，Al/Ni の電子密度分布を，原子を置換する (100)
面と (110) 面 (セル対角線方向の断面) について示す．スケールは界面近傍の電子密度に合わせ 0〜1/Å としたためそれより密な部分はすべて赤色に表示されている．また図中には電子密度が 0.05 変化する毎に結んだ等電子密度曲線を示している．

Fe/Ni では界面の原子を H，O に置換すると著しく電子密度が変化し，Fe，Ni 間より電子が疎になった．Fe/Ni の界面エネルギーを減少させた Mo，W，Cr はいず
れも Fe-X（X: Mo，W，Cr）間では Fe 原子間の電子密度と同程度の電子密度となり，Ni-X 間にも電子密度が中程度で幅広く共有された．W は正の界面エネルギーを示したが，界面距離を変えれば凝着性を示す可能性がある．

Al/Ni も H，O は著しく電子密度を変化させた．Al は Fe と違い原子核まわりに価電子が局在しておらず Al 内で広がっている．第 3 元素は程度の差があるが価電子が局在しているため反力を生じた可能性がある．Si，P は局在の度合が低くまた Al側への影響が少ないため界面の凝着性があまり低下しなかった．

46

Fig. 4.2.1

Change in interface energy of Fe/Ni by third element

Fig. 4.2.2

Change in interface energy of Al/Ni by third element

47

Fig. 4.2.3

Fig. 4.2.4

Change in interface energy of Mg/Ni by third element

Electron density distribution of Fe/Ni interface on the (100) plane

48

H

O

Si

P

Cr

Co

Cu

Zn

Mo

W
Fig. 4.2.5

Electron density distribution of Fe/Ni interface with third element on
the (100) plane49

Fig. 4.2.6

Electron density distribution of Fe/Ni interface on the (110) plane

50

H

O

Si

P

Cr

Co

Cu

Zn

Mo

W
Fig. 4.2.7

Electron density distribution of Fe/Ni interface with third element on
the (110) plane51

Fig. 4.2.8

Electron density distribution of Al/Ni interface on the (100) plane

52

H

O

Si

P

Cr

Co

Cu

Zn

Mo

W

Fig. 4.2.9

Electron density distribution of Al/Ni interface with third element on
the (100) plane

53

Fig. 4.2.10

Electron density distribution of Al/Ni interface on the (110) plane

54

H

O

Si

P

Cr

Co

Cu

Zn

Mo

W

Fig. 4.2.11

Electron density distribution of Al/Ni interface with third element on
the (110) plane55

第 5 章 結言
めっきの密着性予測に繋がることを目的とし，第一原理計算を用いて表面エネ
ルギー，界面エネルギーの評価や第 3 元素による影響を議論した．

2 章では，平面波基底での第一原理計算の基礎理論および電子状態計算の高速化について説明した．

3 章では，まず Fe，Al，Mg，Ni，Si，P，Cr，Co，Cu，Zn，Mo，W の 12 種類の元素について単結晶のバルクと，真空層を設けたスラブモデルの全エネルギーを求め，各単元系の表面エネルギーを求めた．次に Fe，Al，Mg，Ni の 4 金属について，バルクモデルとスラブモデルそれぞれで第 3 元素に置換して全エネルギーを求め，第 3 元素の表面，バルク内での影響および表面エネルギーの変化を調べた．第3 元素は H，O，Si，P，Cr，Co，Cu，Zn，Mo，W の 10 種類を対象とした．Fe では O，P，Cu，Zn の添加で表面エネルギーが下がり，Si は著しく上昇させた．Al は

H，O，Si，P が表面エネルギーを下げたがその効果は小さい．Mg ではいずれの元素も表面エネルギーを著しく上昇させた．Ni は多くの元素で表面エネルギーが下がり，特に O では 0 に近い値となった．

4 章では Fe/Ni，Al/Ni，Mg/Ni 界面を有するスーパーセルを用いて，界面の Fe，Al，Mg を第 3 元素に置換した時の界面エネルギーに与える影響および電子密度分布の変化を調べた．Mg/Ni 界面は Mg/Ni 間の格子長さが不適切だったためいずれも正の界面エネルギーとなった．Fe/Ni 界面も正の界面エネルギーだったが，Mo，Pを置換すると負となり凝着性のある界面となった．Al/Ni 界面では第 3 元素のない
界面エネルギーが大きな負値を示し，本質的に強固であるとの報告を裏付けた．第3 元素に置換するといずれも負の界面エネルギーが小さくなり密着性を低下させる．

56

ただし Si と P は比較的影響が小さかった．

57

参考文献

(1) Y. Mutoh, 日本機械学会論文集 A 編 68-670 (2002), 909-915(2) K. Nakasa, 材料 44-498 (1995), 321-326(3) Y. Nakamura, 鉄と鋼 92-11 (2006), 676-682(4) Non published private data(2018)(5) Y. Yan, Eng Fract Mech 193 (2018), 137-150(6) Y. Okura, 日本機械学会論文集 A 編 77-775 (2011), 495-504(7) T. Sumigawa, Mater. Sci. Eng. A 527 (2010), 6518-6523(8) T. Kitamura, 日本機械学会論文集 A 編 69-684 (2003), 1216-1221(9) H. Hirakata, 材料 55-12 (2006), 1073-1080(10) B.J. Henz, Mol. Simulation 35 (2009), 804-811(11) E. V. Levchenko, Comput. Mater. Sci. 47 (2010), 712-720(12) A. V. Evteev, Intermetallics 19 (2011), 934-941(13) H. Z. Wu J. Alloys Compd. 581 (2013), 519-525(14) A. Liu, Sci. Rep. 4 (2014), 3837
(15) W. Liu, Phys. Rev. B 73 (2006)58

(16) S. Ogata, 日本金属学会誌

60-11(1996), 1079-1083(17) R. Nakanishi, 表面技術

64-10 (2013), 547-550(18) M. Minamiyama, 材料

67-10 (2018), 930-936(19) K.Yashiro, Mater. Trans. 59-11 (2018), 1753-1760(20) P. Hohenberg, Phys. Rev. B 136 (1964), 864-871(21) W. Kohn, Phys. Rev. A 140 (1965), 1133-1138(22) J. P. Perdew, Phys. Rev. B 23 (1981), 5048-5079(23) C. キッテル, 固体物理学入門 上 第 7 版, 丸善 (1998)(24) G. B. Bachelet, Phys. Rev. B 26, (1982), 4199
(25) N. Troullier, Phys. Rev. B 43 (1991), 1993(26) L. D. ランダウ, E. M. リフシッツ, 量子力学 1, 東京図書 (1992)(27) 上田顕, コンピュータシミュレーション, 朝倉書店 (1990)(28) 香山正憲, 固体材料の電子状態の基礎, 材料学会勉強会資料 (1993)(29) O. H. Nielsen, Phys. Rev. B 32 (1985), 3780
(30) P. G. Dacosta, J. Phys. C: Solid State Physics 19 (1986), 3163
(31) D. Vanderbilt, Phys. Rev. B 4 (1990), 7892
(32) L. Kleinman, Phys. Rev. Lett. 48 (1982), 1425
(33) 尾口多美夫，佐々木泰造，固体物理, 25, (1990), 857(34) C.L. Fu, Phys. Rev. B 28 (1983), 5480

59

(35) M.P. Teter, Phys. Rev. B 40 (1989), 12255.
(36) R.D. King-Smith, Phys. Rev. B 49 (1994), 5828.
(37) M.J. Gillan, J. Phys. Condens. Matter. 1 (1989), 689(38) I. Stich, Phys. Rev. B 39 (1989), 4997
(39) Teter, Phys. Rev.B 40 (1989),12255(40) D.M. Bylander, Phys. Rev. B 42, (1990), 1394
(41) L. Vitos, Appl. Surf. Sci. 411 (1998), 186-202(42) S. G. Wang, J. Phys. Chem. Solids 61 (2000), 1295-1300(43) Baoqin Fu, Appl. Surf. Sci. 256 (2010),6899-6907(44) 浦野靖久, J. Chem. Software 3-2 (1996), 91-96(45) Baoqin Fu, Appl. Surf. Sci. 255 (2009) , 9348-9357(46) H. L. Skriver, Phys. Rev. B 46 (1992), 7157-7168(47) E. Aghemenloh, Com. Mat. Sci. 46 (2009), 524-530(48) 水口隆, 日本機械学会誌 68-2 (2004), 86-89

60

謝辞

本研究をするにあたり浅学非才な著者に対し，丁寧な指導をしていただいた．屋
代如月教授に心よりお礼を申し上げます．本論文を完成させるにあたり，多くの助言をしていただいた内藤圭史助教にも心より感謝いたします．ともに切磋琢磨し合った研究室メンバーや先輩方，同じ専攻の友人達にも御礼申し上げます．

61