要約
酸化イットリウム (Y2 O3 ) が酸化物分散型強化鋼 (ODS 鋼) の高温での機械的特性を向上させることは経験的に知られているが，そのメカニズムについては不明な点が多い．本研究では，ODS 鋼開発の一助とすべく，Y2 O3 ならびに Al,Ti を添加した Fe の第一原理計算を行い，全自由エネルギー，弾性係数，O 原子の溶解熱等から検討した．
まず，bcc の Fe 単位格子を 2 × 2 × 2 並べたスーパーセルを用いて，その中の bcc の第 3 近接位置の 2 つの Fe 原子を Y 原子に置換し，Y-Y 結合を挟む 3 つの O 原子の添加
サイト (bcc の八面体中心) の組み合わせを 5 通り考慮した解析を行った．弾性係数の行列式は今回考慮した全ての組み合わせで正値をとり，力学的には存在し得る系であることが明らかにされたが，いずれの系も弾性係数の値は Fe 単体より低下し，Y2 O3の添加が直接剛性アップするわけではないことが分かった．また，5 つの組み合わせの中で全自由エネルギーが低い構造では O 原子は大きな負の溶解熱を持ち，強く結合することが示唆された．
次に，第 3 元素の役割を検討するために，先の組み合わせにおいて Y-Y を Y-Al や

Al-Y のように Al, または Ti に置換した構造について全自由エネルギー，弾性係数，O原子の溶解熱を調べた．その結果，Al,Ti の添加によって Y2 O3 だけの場合よりも弾性係数が増加し剛性がアップするケースが多く見られた．O 原子の溶解熱は，Y2 O3 だけの解析で全自由エネルギーが低く溶解熱が大きいことが示されたサイトでは，添加によって溶解熱が小さくなり (結合が弱まる)，逆に溶解熱が小さかったサイトでは添加によって結合が強められていた．このことは，第 3 元素の添加によって酸素が Y 原子の近傍のおいて様々な形態で結合し得る可能性を示唆している．

2

Summary
It is known that Yttrium oxide improves the mechanical properties of oxide dispersion strengthened (ODS) ferritic steels at elevated temperature. In the present study,
we have evaluated the elastic coeﬃcients, CIJ , total free energy, and Oxygen solution
heat by using the ab-initio DFT calculation to obtain new insight for the developmentof ODS steels. First, we analyzed bcc-Fe containing Y2 O3 using a supercell consist of
2 × 2 × 2 bcc unit lattices. Two yttrium atoms are substituted in the third neighborsite and three oxygen atoms are set in the center of bcc octahedrons around Y-Y connection. Here, we have considered 5 cases for the combination of the three octahedral
site. It is shown that all the models can be exist since the minor determinants of CIJ is
positive(lattice stability analysis). However, the components of each elastic coeﬃcientare lower than those in pure Fe. That is, Y2 O3 never improve the elastic stiﬀness. On
the oxygen solution heat, it is suggested that the Y2 O3 -conﬁguration which has lowest
total free energy has large negative solution heat in its oxygen site. Second, we havesubstituted Al or Ti with Y in the previous Y2 O3 -Fe calculations. It is noteworthy
that the elastic coeﬃcients are increased from the Y2 O3 -Fe systems in many cases. On
oxygen solution heat, the heat solution is weaken in the doping site above mentioned
in the Y2 O3 -Fe system while it is strengthened in the other oxygen doping site. Thatis, the addition of Al or Ti can lead the possibility of the oxygen solution in various
conﬁguration around Y, Al and Ti atoms.

修士論文

第一原理格子不安定解析による酸化物分散強化機構の基礎的研究

指導教員：屋代 如月

山本 智2010 年 2 月
神戸大学大学院 工学研究科 博士過程前期課程 機械工学専攻

Master thesis

Fundamental Study for Oxide Dispersion
Stlengthen Mechanism:
Ab-initio Lattice Instability Analysis

Tomo YAMAMOTO
February 2010
Division of Mechanical Engineering,
Graduate School of Technology
Kobe University,Kobe,Japan

目次
1

緒論

1

2

第一原理分子動力学法の概要

3

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

3

2.2 密度汎関数法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

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

6

2.4 逆格子空間 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

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

8

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

12

2.7 応力 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

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

14

2.8.1

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

15

2.8.2

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

19

2.9 電子占有数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.10 FFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.11 電子系の最適化手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22格子不安定性解析の概要

27

3.1 不安定条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

3.2 応力と弾性係数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

29

3.3 弾性係数による無負荷平衡点での格子不安定性評価 . . . . . . . . . . .

30

3

i

目次

4

5

ii

酸化イットリウム含有 Fe の第一原理格子不安定解析

32

4.1 解析方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

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

37

4.2.1

全自由エネルギーの比較 . . . . . . . . . . . . . . . . . . . . . .

37

4.2.2

格子安定性と弾性係数 . . . . . . . . . . . . . . . . . . . . . . .

39

4.2.3

弾性係数成分と価電子密度の関係 . . . . . . . . . . . . . . . . .

41

4.2.4

O 原子溶解熱の検討 . . . . . . . . . . . . . . . . . . . . . . . .

45

Al および Ti 原子添加による系の安定性評価

475.1 解析方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

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

51

5.2.1

全自由エネルギーの比較 . . . . . . . . . . . . . . . . . . . . . .

51

5.2.2

格子不安定性解析による系の安定性評価 . . . . . . . . . . . . .

53

5.2.3

O 原子溶解熱が安定性におよぼす影響 . . . . . . . . . . . . . .

63

6

結論

68

A

関連学術講演

70

参考文献

75

謝辞

78

第1章 緒 論
軽水炉の高燃焼度化は，資源の有効利用や原子力エネルギーの高効率利用に欠かせない技術の一つであり，その最も有効な手段の一つとして，炉心被覆管の改良が進められている．現在の軽水炉炉心材にはジルコニウム合金やオーステナイト鋼が用いられている．ジルコニウム合金は中性子吸収性に優れているものの高温下での強度に課題があり，より一層の高燃焼度化は難しい．またオーステナイト鋼は加工性や水環境下における耐食性において優れた材料であり，高燃焼度化に向けた研究開発が進められている (1) が，高照射領域における照射脆化やスウェリング，照射誘起応力腐食割れなどの問題を抱えている．酸化物分散型強化鋼 (ODS 鋼) は，オーステナイト鋼に比べ耐照射性，耐ヘリウム脆化特性に優れていながら，高温での材料安定性にも優れているため，高温かつ高照射領域に耐える高燃焼度被覆管材料として着目されている．
また，現在の製鋼技術基盤の利用が可能であり，現在生産されている鉄鋼材料の中では総合的に判断すると最も優れた強度特性を持つことから，新たな構造用鉄鋼材料に発展する可能性も有している．ODS 鋼中の添加物として用いられる酸化イットリウム(Y2 O3 ) は主にジルコニアの安定剤や，化学蒸着法 (CVD) を用いたセラミックスのコー
ティング剤として知られており，ODS 鋼中に添加することで，高温でのクリープ特性の向上に寄与していることが実験的手法により明らかにされている (2) ．また，熱処理
後の ODS 鋼中には，ナノサイズの酸化物が観察され，Al や Ti に代表される第 3 の元素を伴い，YTiO3 ，YAlO3 のような複合酸化物として析出していることが知られている (3) ．しかしながら，製造プロセスのメカニカルアロイング中にとりだして観察すると酸化物の析出は認められず，酸化イットリウムは均一に分布をしているという報告もある (4) ．そのため，複雑な製造プロセスにおいて酸化物の分散や機械的特性向上のメカニズムの解明が求められている．
原子種と原子配置のみを必要情報とし，量子力学にのみ基づいて非経験的に材料物

1

性を予測可能な第一原理計算は，従来の実験による試行錯誤的な手法によらず，電子論に基づいて非経験的に材料特性を評価する手法として注目されており，実材料の特性を検討した例も多数報告されている (5)(6) ．イットリウムを対象とした第一原理解析
として，化合物の結晶構造や磁気的特性の予測 (7)(8) などがある．しかしながら膨大な計算時間とメモリを必要とする第一原理計算は，現時点では数十〜数百程度の原子数での解析が限界であり，変形時におけるイットリウムの役割などを直接扱うことは難しい．
少数の原子の解析でも，力学特性に重要な知見を与えるものに格子不安定性解析 (9)がある．格子不安定とは，外力下で変形している結晶格子が釣り合いを失い，外力の増加を必要とせずに変形が進行する状態を示している．本来，完全均一結晶の弾性安定限界を評価するために提案されたものであるが，多原子系における転位発生等の局
所変形開始の有力な基準となることが報告されている (10) ．格子不安定解析は少数の原子しか扱えないが，大変形下における電子状態を精密に評価できる第一原理計算に適した解析である．第一原理計算による不安定解析の例としては W,Mo,Nb,Fe といった

bcc 金属の理想強度評価 (12)(14) ，引張およびせん断下の鉄の理想強度 (15) ，Si と Al の[001] 引張下の格子不安定性 (11) ，fcc 金属 (Al,Cu,Ag) の安定性 (16) 等がなされてきた．
本論文では，酸化物分散強化のメカニズムに関する知見を得るべく，16 原子よりなる bcc-Fe 8 格子を基に Y,Al,Ti,O 原子を添加した系について第一原理計算による格子不安定性解析を行う．第 2 章では第一原理分子動力学法の基礎理論，および電子状態計算の高速化手法について説明する．第 3 章では，格子不安定性解析の概要と，弾性係数を用いた格子不安定性評価について説明する．第 4 章では，bcc-Fe に酸化イットリウムを添加した系において解析を行い，イットリウム原子および酸素原子の添加による特性向上を格子不安定性の観点から検討すると共に，全自由エネルギー，O 原子の溶解熱，異なる O 原子添加サイトの組み合わせによる弾性係数の違い等について検討する．第 5 章では，酸化イットリウムに第 3 の元素としてアルミニウムおよびチタンを添加した系において，第 4 章同様に格子不安定性の観点および，全自由エネルギー，

O 原子の溶解熱の観点から考察する．第 6 章では本研究で得られた結果を総括する．

2

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

2.1

断熱近似と平均場近似

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

3

に，通常，以下の 2 つの基本的な近似が導入される．

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

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

2.2

密度汎関数法Hohenberg と Kohn は，外場ポテンシャル v(r)(原子核からの電場) 中における多電子系 (N 電子系) の基底状態の全エネルギー Etot が電子密度 ρ(r) の汎関数として∫Etot [ρ] =+v(r)ρ(r) dr + T [ρ]12∫∫

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

Ts [ρ(r)] =

occ
∑
i

1
< ψi | − ∇2 |ψi >
2

(2.2)

と書くと，式 (2.1) は∫Etot [ρ] =1 ∫∫ ρ(r ′ )ρ(r) ′+dr dr + Exc [ρ]2|r ′ − r|(2.3)

Exc [ρ] = T [ρ] − Ts [ρ] + E[ρ]

(2.4)

∑occ

のように書ける．ここで，v(r)ρ(r) dr + Ts [ρ]iは電子が占有している準位についての和をとることを表す．Exc は一電子近似のもとでの交換相関エネルギーであり，電子間相互作用を考慮した電子の運動エネルギー T [ρ] から，相互作用のない電子の運動エネルギー Ts [ρ] を分離することによって，電子間の複雑な相互作用を全てこの項に押し込めている．
∫

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

1[− ∇2 + veﬀ (r)]ψi (r) = εi ψi (r)
2

(2.5)

ここで， veﬀ (r) は有効一電子ポテンシャルであり次式となる．
∫

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

(2.6)

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

ρ(r) =occ∑|ψi (r)|2

(2.7)

i

となる．
以上のようにして，多電子問題は式 (2.5)∼(2.7) を Self–Consistent に解く問題に帰
着される．

5

2.3

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

∫





εxc (r) ρ (r) dr
(2.8)
δExc [ρ (r)]dµxc (r) == εxc (r) + ρεxc (r) δρdρとして計算する．つまり，電子密度 ρ(r) の点 r における交換相関エネルギーを同じ

Exc [ρ] =電子密度の一様電子ガス中のそれで代用する．これを局所密度近似 (Local Density

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

εxc (r) = εx + εc

(2.9)

0.4582rs(2.10)εx (r) = −



 −

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)
ここで，

(

rs =

3 1
4π ρ

(2.11))13(2.12)

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

(2.13)

4µx (r) = εx3(2.14)µc (r) =

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

)]

rs
1.0529(2.15)

1)
+
1
+
(rs >
√

=
2

3(1 + 1.0529 rs + 0.3334rs )0.6668rs−0.0584 + 0.0311 ln rs − 0.0084rs + 0.00133rs ln rs (rs <
= 1)
6

となる．

2.4

逆格子空間

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

(n1 , n2 , n3 は整数)

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

b = 2π12
a1 · a2 × a3

a1 × a2 b3 = 2π(2.17)

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

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

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

G · R = 2π(m1 n1 + m2 n2 + m3 n3 )(2.19)

を満たす．結晶の並進対称性から，波動関数 ψ(r) と ψ(r + R) は同じ固有値をとる関数となり，

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

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

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

k=

h1
h2
h3
b1 + b2 + b3
n1
n2
n3
7

(h1 , h2 , h3 は整数)

(2.22)

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

2.5

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

∑

n
Ck+G
| k+G>(2.23)

1
exp[i(k + G) · r]Ω

(2.24)

Gと表す．ここで，|k + G >=であり (Ω は全結晶体積)，規格直交条件

1∫exp [−i (k + G) · r] exp [i (k + G′ ) · r] drΩ Ω1∫
exp [i (G′ − G) · r] dr=Ω Ω< k + G | k + G′ > == δGG′(2.25)

を満たす．式 (2.23) 中の Σ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.26)

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

は Brillouin ゾーン内の k 点についての和をとることを表す．

8

以上のように平面波を基底関数として波動関数を展開すると，Kohn–Sham 方程式(2.5) は次のように展開係数を固有ベクトルとする行列固有値問題となる．
∑

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

=⇒

∑
G′

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

(2.27)

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

1< k + G|f (r)|k + G > =f (r) exp[−iG · r] exp[iG′ · r]drΩ Ω
1∫=f (r) exp[−i(G − G′ ) · r]drΩ Ω′

= f (G − G′ )(2.28)

1
1< k + G | − ∇2 | k + G′ >= |k + G|2 δGG′
2
2

(2.29)を用いる．

(a) 運動エネルギーの項
運動エネルギーの項はとなる．
一方，式 (2.6) に示したように veﬀ は原子核からのクーロン相互作用項 (v)，電子間
クーロン相互作用項 (Vcoul )，交換相関項 (µxc ) からなる．平面波基底バンド計算では結晶結合に重要な役割を果たす価電子のバンド構造を効率的に計算するため，原子核からのクーロン項のかわりに内殻電子と原子核を正電荷をもったひとつのポテンシャ
ルとして扱う擬ポテンシャル法が用いられることが多い．擬ポテンシャル法を用いる
ことにより，膨大な平面波数を必要とする内殻電子の波動関数を直接扱うことなく価電子状態を正確に表すことができる (21)(22) ．擬ポテンシャルは 2.8 節で後述するよう
PPと，依存する非局所擬ポテに，電子の角運動量に依存しない局所擬ポテンシャル Vloc,lPPからなり，次式で表される．
ンシャル Vnloc
PPPP(r − Ra ) P̂l(r − Ra ) + Vnloc,lVlPP (r − Ra ) P̂l = Vloc

9(2.30)

ここで，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Ω Ω locPP= Vloc(G − G′ )(2.31)

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

∑∑
R

VaPP,loc (|r − r a − R|)(2.32)

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

∫

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

VaPP,loc (G) =

∫

= 2π
=

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

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

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

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

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

1 ∑exp[−i(G − G′ ) · r a ]VaPP,nloc (k + G, k + G′ )
Ωat a

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

VaPP,loc (k + G, k + G′ )10(2.34)

= 4π

∑

∫

PP,nlocVa,l(r)jl (|k + G|r)jl (|k + G′ |r)r2 dr(2l + 1)Pl (cos ω)l(2.35)となる (23) ．ここで，Pl は Legendre 多項式，jl は球 Bessel 関数であり，ω は k + G と

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

ρ(G) exp[iG · r]

(2.36)

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

(2.37)

ρ(r) =

G

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

ρ(G) exp[iG · r]

(2.38)

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

Vcoul (r) = 4π

∑ ρ(G)
G|G|2

exp[iG · r]

(2.39)

が得られる．これより，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′ ) ∫
′ 2

|G |

G′

= 4π

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

ρ(G)
|G|2

(2.40)

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

1
Vcoul (r) exp [−iG · r] exp [iG′ · r] dr< k + G |Vcoul (r) | k + G > =Ω Ω
1∫=
Vcoul (r) exp [−i (G − G′ ) · r] drΩ Ω′= Vcoul (G − G′ )11(2.41)

となる．

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

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

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′ =1PPPP|k + G|2 δGG′ + Vloc
(G − G′ ) + Vnloc(k + G, k + G′ )2
+Vcoul (G − G′ ) + µxc (G − G′ )(2.44)と逆空間での表式となる．

2.6

系のエネルギー

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

BZ ∑
occ
∑
k

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

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

]

+
−

1 ∑ ∑ a a′ ∑ erfc (|R + r a′ − r a | γ)Z Z
2 a a′ v v R|R + r a′ − r a |∑ Zva2 γ

√

a

Z 2π
2πZ 2
−
+ lim
π
2Ωat γ 2 G→0 Ωat |G|2

(2.46)である．
ここで∑

ρ(r) =

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

(2.47)

Gという関係を用いると
BZ
occ
∑
∑
∑
1∑
n
PP
Etot =
fk
fn|k + G|2 |Ck+G|2 + Ωat
Vloc(G)ρ(−G)
2 k
n
G
G

+

BZ
∑

fk

occ
∑
n

k

∑

fn

∑∑

GG

′

′
n∗
n
PP
Ck+G
Ck+G
′ Vnloc (k + G, k + G )∑1
Vcoul (G)ρ(−G) + Ωat+ Ωat
εxc (G)ρ(−G) + EEwald
2G
G(2.48)とフーリエ成分により表現できる．
PP
式 (2.33),(2.40) より，Vloc(G) と Vcoul (G) は G = 0 で発散するが，これらの発散成

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

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
∑
k

fk

occ
∑
n

fn

∑ ∑
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.49)′
は，式 (2.46) の第 5 項の発散成分を取り除いたものである．
ここで，EEwald

2.7

応力

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

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

Sa (G) = exp(−iG · r a )(2.50)

はスケーリングの元のもとで不変であるから，平均応力は∂K γ= −δαγ K β∂εαβ( K γ = (k + G)γ )∂Ωat= −δαβ Ωat∂εαβ(2.51)(2.52)という関係を用いることにより，

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

=−
−

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

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

+

+ δαβ

∑

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

G

+

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

2.8

擬ポテンシャル法

ブロッホの定理 (20) により，固体中の電子の波動関数は平面波基底により展開が可能である．しかし，平面波基底では原子核に強く引き付けられて局在している内殻電14子の波動関数や，価電子密度の著しい変動をを表現するには非常に多くの展開項数を要する．平面波数は解くべきハミルトニアンの次元数に比例し直接計算量に影響するので，これをできるかぎり少なくすることが望ましい．通常の固体材料では，内殻電子は原子核に強く引き付けられており，他の原子からの影響をほとんど受けず価電子がその特性を決定付けているといえるので，内殻電子と原子核をひとつのイオンと考え，原子間領域の価電子のみを取り扱うのが擬ポテンシャル法である．擬ポテンシャル
法は，その歴史の初期においては原子核付近で強い反発作用が現れたり，原子核領域
において真の波動関数と擬波動関数の 2 乗のノルムが一致していなかったりしたため，

self-consistent な計算には適用できなかった．そこで，Hamman らは，これらの問題を解決した HSC 型 (BHS 型) と呼ばれるノルム保存型擬ポテンシャルを開発した (21) ．しかし，第二周期元素や遷移金属では依然として非常に多くの平面波数が必要であった
ため，Troullier らはそれらの元素においても比較的少ない平面波数で扱える TM 型擬
ポテンシャルを開発した (22) ．また，Vanderbilt らはノルム保存条件をはずすことにより，さらに少ない平面波数で計算を行えるウルトラソフト型擬ポテンシャルを開発した (28) ．
本節では，まずノルム保存型擬ポテンシャルとして TM 型を説明する．その後，ノルム非保存型擬ポテンシャルとしてウルトラソフト型とそれを用いた場合の系のエネ
ルギー等について説明する．

2.8.1

TM 型擬ポテンシャル

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

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

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

15

(2.54)

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

2. 内殻領域で節を持たない擬波動関数 ψlPP (r) を次式のような解析関数形で表す．
{

ψlPP (r) =(r >ψlAE (r)= rcl )
lr exp [p(r)] (r <= rcl )

(2.55)

p (r) = c0 + c2 r2 + c4 r4 + c6 r6 + c8 r8 + c10 r10 + c12 r12(2.56)

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

′

′′

2

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

(r >= rcl )

(2.57)(r <= rcl )

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

∫ rcl
0

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

より，[∫ r

2c0 + ln0cl]r2(l+1) exp [2p(r) − 2c0 ] dr = ln

[∫ r
0

cl

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

P(rcl )p(rcl ) = ln
rcll+1P′ (rcl ) l + 1′
p (rcl ) =−P(rcl )rcl
2(l + 1) ′2
p (rcl ) − [p′ (rcl )]
rcl
2(l
+
1)
p′′′ (rcl ) = 2VlAE′ (rcl ) +p′ (rcl )rcl2
2(l + 1) ′′−p (rcl ) − 2p′ (rcl )p′′ (rcl )rclp′′ (rcl ) = 2VlAE (r) − 2εl −

16(2.60)(2.61)(2.62)(2.63)

4(l + 1) ′p (rcl )rcl2
2(l + 1) ′′′4(l + 1) ′′+p (rcl ) −p (rcl )2
rcl
rcl

p′′′′ (rcl ) = 2VlAE′′ (rcl ) −−2 [p′′ (rcl )] − 2p′ (rcl )p′′′ (rcl )2(2.64)

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

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

4. これらの非線形連立方程式を解く．まず c2 を仮定し，式 (2.65) から c4 を決める．
残りの 5 個の係数は式 (2.60)〜式 (2.64) の連立一次方程式であり，ガウス消去法
により求める．最後に求まった係数を用いて c2 が妥当であるか式 (2.58) により判断する．c2 の決定には bisection 法を用いる．

5. 以上により求まった擬ポテンシャルから，価電子による遮蔽効果を取り除くこと
により内殻電子を含めたイオンの裸のポテンシャルを得る．
PP
PP
PP
Vion,l(r) = Vscr,l(r) − Vcoul(r) − µPP
xc (r)(2.66)

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

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

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

擬ポテンシャルの KB 分離型表現
平面波展開による第一原理分子動力学法では，大きなハミルトニアン行列を繰り返し解く必要があるため，その繰り返しの中で変化しない量はメモリー上に記憶しておくことが高速化の基本となる．特に式 (2.35) の非局所項は，平面波の 2 乗のループを17

含んでおり計算時間がかかるとともに，記憶する量も平面波数の増加に対してその 2乗で増える．そのため，大規模な計算ではすぐにメモリー容量に破綻をきたす．そこ
で，非局所項に次式で表される KB 分離型表現 (29) を用いれば，平面波の 2 乗のループ
は 1 乗のループとなる．
PP
PP
(r)|
(r)ψlPP (r) >< ψlPP (r)Vnloc,l|Vnloc,l
KB
Vnloc,l (r) =P̂lPP< ψlPP (r)|Vnloc,l(r)|ψlPP (r) >(2.68)

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

×

ΩCl

{∫ ∞
0
l
∑

×

0

}
PP
ψlPP (r)Vnloc,l(r)jl (|k + G|r)r2 dr

PPψlPP (r)Vnloc,l(r)jl (|k + G′ |r)r2 dr}∗
Ylm (k + G)Ylm(k + G′ )(2.69)m=−lとなる．ここで，PP
Cl =< ψlPP (r)|Vnloc,l(r)|ψlPP (r) >PP
Vnloc,l(r) =∑∑PP,nlocVl,a(|r − ta − R|)(2.70)(2.71)aRである．したがって，PP,nloc
PP
PP
Cla =< ψla(r)|Vl,a(r)|ψla(r) >Ala (k + G) =∫ ∞0PP,nlocPPψla(r)Vl,a(r)jl (|k + G|r)r2 dr(2.72)(2.73)とおくと，(4π)2 ∑ ∑ 1
KB′< k + G|Vnloc,l (r)|k + G > =Ωat l a Cla×l∑{exp[−iG · r a ]Ala (k + G)Ylm (k + G)}

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

KB
n
< k + G|Vnloc,l(r)|k + G′ > Ck+G′

18

(2.74)

=(4π)2 ∑ ∑ 1Ωat l a Cla×l∑{exp[−iG · r a ]Ala (k + G)Ylm (k + G)}

m=−l

×


∑


G′

∗
n
′
′
′
Ck+G
′ exp[iG · r a ]Ala (k + G )Ylm (k + G )(2.75)となり，

AYlam (k + G) = Ala (k + G)Ylm (k + G)∗(2.76)をあらかじめ記憶しておけば計算が速くなる．また，この行列要素を計算した際に，

CAYnalkm (k + G) =

∑

n
∗
′
′
′
Ck+G
′ exp[iG · r a ]Ala (k + G )Ylm (k + G )(2.77)

Gを記憶しておけば後のエネルギーや原子に働く力の計算が高速化できる．

2.8.2

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

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

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

19

′
=
VNL

|χi >< χi |< χi |Ψi >(2.79)

′
ここで，VNL
は非局所ポテンシャルであり，関数 χi は

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

ただし|βi >= Σj (B −1 )ji |χj >Bij =< Ψi |χj > ,(2.82)

また，< Ψi |βj >= δij(2.83)である．擬ポテンシャルにノルムの保存条件を課さなかったことにより，内殻領域において電子密度がQij (r) = Φ∗i (r)Φj (r) − Ψ∗i (r)Ψj (r)(2.84)だけ不足している．また求められた波動関数も，ノルムが
∫

Qij =

r<rcl

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

S =1+

∑

Qij |βi >< βj |(2.86)

ijと定義すれば，規格直交条件が以下のように満足される．

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

∑

(Bij + εj Qij )|βi >< βi |

ij

となる．

20

(2.88)(2.89)

2.9

電子占有数

金属では Fermi エネルギー εF の近傍に多くのエネルギー準位が存在するため，整数の占有値では問題が生じる (30) ．たとえば時間とともに Fermi エネルギー近傍の２つの準位が交差してしまうと，電子密度が不連続に変化してしまう．このような問題を避けるために，Gaussian Broadening(31) という方法を用い，fn のかわりに非整数の占有
数 fi

[

(

1
εi − εFfi =1 − erf2
σ)](2.90)

を導入し．フェルミレベルに対して σ の幅で占有状態をぼかしてある程度の非占有状態も計算する．実際の数値計算では

2∑

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

S = −kB∑

Ef = Etot − T S(2.92){fi ln fi + (1 − fi ) ln(1 − fi )}(2.93)

iを考えなければならない．

2.10

FFT

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

∑

n
Ck+G

(2.94)nexp[iG · r]
Ck+G

(2.95)

ukn (G) =Gとおけば，フーリエ逆変換より

ukn (r) =∑G

21

となる．同様に

u∗kn (r) =

∑

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

(2.96)

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

(2.97)

Gであるから，

ukn (r)u∗kn (r) =

∑∑
G

G′

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

ρ(r) =

occ ∑
BZ
∑
n

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

(2.98)となり電子密度分布が得られる．すなわちハミルトニアンから求められる固有ベクトn
ル Ck+Gをフーリエ変換することにより，実空間の電子密度分布 ρ(r) を式 (2.26) に従って直接評価するより高速に計算できる．ρ(r) が求められれば交換相関エネルギー，交換相関ポテンシャルの実空間における値が得られ，フーリエ変換によって逆空間での値も求められる．このように実際の計算ではフーリエ変換を多用するため，一般に高速フーリエ変換 (Fast Fourier Transformation:FFT) のプログラムが用いられる．

2.11

電子系の最適化手法

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

共役勾配法の原理
共役勾配法は，一般には正定な係数行列をもつ連立１次方程式を最適化の考えに立って解くために，あるいは，多次元空間の２次関数 F (X) の最小化問題を解くために用いられる計算手法である．共役勾配法では，前者の問題は結局後者の問題に帰着され，

22適当な初期値 X0 から出発して順次修正を加えながら · · · ，Xm−1 ，Xm ，Xm+1 ，· · · と変化させて F (X) を最小にする X を探索する．
密度汎関数法に基づく電子状態計算では系の全エネルギー Etot は，電子密度すなわ
ち波動関数の汎関数で表され正しい波動関数によって最小化される．したがって，平面波基底の波動関数を用いた場合には，系の全エネルギーを最小にする係数ベクトル

C nk を規格直交条件のもとで求める計算を行えばよい．
（ここで，C nk は平面波展開係n
数 Ck+Gを成分に持つベクトルである．
）すなわち，′Etot= Etot −∑λmn (< Ψmk |Ψnk > −δmn )mn= Etot −

∑

λmn

(
∑

mn

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

(2.99)

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

∑∑
G

G′

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

]





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

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

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

方向) を求める．それをもとにして εkn を最小にするような新たな係数ベクトルを求める．以上の手順を εkn が収束するまで繰り返した後に，電子密度とハミルトニアンの更新を行い全エネルギーを計算する．

< 残差ベクトル >′Etot を εkn に置き換えることによって，式 (2.99) の Etot
は ε′kn に置き換わるとすると，残差ベクトルは式 (2.101) より次式で表される．
[
i
Rnk
=−

∂ε′kn
m∗
∂Ck+G

]i

(

ただし，式中の i は， i 回目のステップにおける[)i= − H − λin I · C nk

(2.102)という意味を表し，

]

n,i
i
i
i
C nk
= Ck+GH = [Hk+G,k+G′ ] , λin =< ΨnkĤ Ψnk>′ ,(2.103)である．これは，i のステップにおいて εkn を最小にする方向 (最急降下方向) を示すベクトルを表している．
iRnkには，最終的に得られる次のステップの波動関数 Ψi+1
nk が同じ k 点における n 以

外の全バンドの波動関数 Ψmk (m ̸= n) と直交するように，直交化処理が施される．
i′
i
Rnk
= Rnk
−

∑ (

)

i∗
i
C mk
· R ink C mk

(2.104)m̸=n<preconditioning>′i
残差ベクトル Rnk
に対して preconditioning という処理を施す．大きな逆格子ベクトルについては平面波の運動エネルギーが大きくなるが，このことが残差ベクトルに影響して収束性を悪化させる．preconditioning は，この問題を回避して収束を速めるた
i
めに行われる．preconditioning された残差ベクトルを Gnk
とすると
i
i′
Gnk
= K i · Rnk(2.105)と表される．ここで(27 + 18x + 12x2 + 8x3 )(27 + 18x + 12x2 + 8x3 + 16x4 )Ekin (G)x=i
Ekin

KGG′ = δGG′

24

(2.106)(2.107)

1
Ekin (G) = |k + G|2
2
1
i
i
i
Ekin
= < Ψnk| − ∇2 |Ψnk>
2

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

)

i′
i
i∗
i
i
Gnk
= Gnk
− C nk
· Gnk
C nk
−

∑ (

)

i∗
i
i
C mk
· Gnk
C mk

(2.110)m̸=nii
ここで，Gnk
は C nkと直交しなければならないことに注意が必要である．

< 探索方向 >探索方向は，次のようにして定められる．
′

i
i
i−1
F nk
= Gnk
+ γ i F nk

γi =






(2.111)i′∗
Gnk· R ink′i−1 ′ ∗
i−1 ′Gnk
· R nk 0(i > 1)(2.112)(i = 1)

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

)

i′
i
i∗
i
i
F nk
= F nk− C nk· F nk
C nk
i
D nk
=(

i′
F nk
i′∗
i′
F nk
· F nk

(2.113)(2.114))1
2< 新たな係数ベクトルの組み立て >新たな係数ベクトルの組立ては次のように行われる．
i
i
C i+1
nk = αC nk + βD nk(2.115)結合係数 α と β は，エネルギー期待値 εkn を最小化するように決定される．すなわち，

C ink ，D ink を基底とする 2 × 2 ハミルトニアン行列，[

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

]

[

=

ε11 ε12
∗
ε22
ε12

]

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

} 21

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

α= {

ε12

}1

ε12 ε∗12 + (ε11 − γ)2 2
ε11 − γβ = −{}1ε12 ε∗21 + (ε11 − γ)2 2(2.118)(2.119)以上の手順を εkn が収束するまで繰り返せばよい．計算の全体的な手順を以下に示す．

1. 係数ベクトル C nk の適当な初期値を，行列計算などによって，全 k 点の全状態について作成する．
i
i+1
から C nkを組立てる一連の計算を反復し，適当
2. 各 k 点の各状態について，C nk

な条件で打ち切る．打ち切り条件は，例えば，1 回のステップでの εkn の減少値
が，最初のステップでの減少値の 30%以下や一定値以下になることである．反復計算が打ち切られれば，同じ k 点における次の状態についての計算へと移る．

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

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

26

第 3 章 格子不安定性解析の概要
格子不安定とは，外力下で変形している結晶格子が釣り合いを失い，外力の増加を必要とせずに不安定に変形が進行する現象を指している．有限変形下の結晶の安定性
は，従来は結晶の変形をブラベー格子の変形で代表することによって系のエネルギーの変数を限定し，エネルギー関数の 2 階微分を解析的に求めることにより評価していた (37) ．一方，Wang らは，結晶の変形をひずみで代表させることによって，系の安定性を弾性係数・弾性剛性係数 (9) の正値性によって評価する手法を提案した (10) ．分子動力学シミュレーションによる検証の結果，原子の熱揺動の影響を含んだ結晶の安定性
が，系全体の弾性係数・弾性剛性係数で評価できることが示されている．この評価方
法は系のエネルギー関数の表式が求まっていない場合でも，数値的に弾性係数・弾性
剛性係数を求めれば安定性評価が可能であるため，第一原理解析でも適用可能である．
本章では，まず従来のエネルギー関数の 2 階微分に基づいて結晶の安定性を評価する手法を説明する．その後，結晶の熱力学関係式から応力と弾性係数の定義を示し，最後に無負荷平衡点で有効である弾性係数の正値性に基づく安定性評価について説明
する．

3.1

不安定条件

結晶の変形を理想化し，すべての結晶格子が外力を受けて均一に変形するものと仮定する．すると fcc を含む立方体格子の変形は図 3.1 に示すような 6 つの格子パラメー
タ a1 ∼ a6 で記述され，内部エネルギー U はこれらの関数 U (a1 , a2 , · · · , a6 ) ≡ U ({am })となる．ここで，本節では原子の運動は考慮しないため，U ≈ Etot である．このとき，{am } の変形状態下にある結晶の安定性は，以下のように微小変形増分 {∆am } による

27

a3
3
4

2
1

6

5

a2
1

Fig.3.1 Unit cell of fcc lattice

エネルギーの変化を考えることによって求められる (9)(37) ．状態 {am } 近傍での内部エネルギーの Taylor 級数展開はU ({am + ∆am }) = U ({am }) +6∑

Fm ∆am +

m=1

6 ∑
6
1 ∑
Amn ∆am ∆an + · · · (3.1)
2 m=1 m=1と表される．ただし，

Fm =

∂U
,
∂am {am }

Amn =

∂2U
∂am ∂an {am }(3.2)

であり，|{am } は状態 {am } における微係数を表す．3 次以上の高次項を省略すると次式のように変形できる．

[U ({am + ∆am }) − U ({am })] −6∑

Fm ∆am =

m=1

6 ∑
6
1 ∑
Amn ∆am ∆an
2 m=1 m=1

(3.3)

左辺第 1 項は系のエネルギー増加量，第 2 項は状態 {am } で周囲の結晶から受けている
力 Fm のもとで微小変形 ∆am をするときになされる仮想的な仕事であり，左辺全体はエネルギー消費量を表している．これが負になると，外力の増加を必要とせずに変形∆am が連続的に生じる不安定状態となる．これより，結晶の力学的安定性はヘッシア
ン [Amn ] の正値性に帰着される．

28

3.2

応力と弾性係数

熱力学の第 1 法則と第 2 法則から，

dU = T dS − dW

(3.4)

である (38) ．ここで，U は内部エネルギー，T は温度，S はエントロピ，dW は系が外界になす仕事である．外部応力 σ の負荷によって結晶が変形する際の dW を求めるた
め，結晶内の任意の点 X が応力の負荷によって X + ∆X に変化する均質一様な変形を考える．変形前の物体表面を S とし，その微小要素を dS とすると，dS において i方向に作用している力 fi は負荷応力 σij を用いて以下のように表せる．

fi = σij dSj(3.5)

X から X + ∆X への変位勾配テンソルを ∆u とすると，∆Xi = ∆uij Xj(3.6)

である．したがって，dS においてなされる仕事は∆W = −fi ∆Xi = −σij dSj ∆uik Xk(3.7)と表される．全仕事 dW は，Gauss の発散定理を用いて次のようになる．

dW = −∫Sσij ∆uik Xk dSj = −∫Vσij ∆uij dV = −σij ∆uij V (X)(3.8)

ここで，V (X) は初期状態 X における結晶の体積である．応力テンソル σij は対称テンソルであるため，式 (3.8) の dW には ∆uij の非対称成分は寄与しない．Lagrange のひずみテンソル1(uij + uji + uki ukj )
2

(3.9)

1
(∆uij + ∆uji ) ∼= ∆uij2(3.10)

ηij =の微小量を ∆uij に等しいとおく．

dηij =

29これより，式 (3.4) は次のようになる．

dU = T dS + V (X)σij dηij

(3.11)

dU = V (X)σij dηij

(3.12)

したがって，断熱過程ではとなり，基準配置における応力テンソルと弾性係数は，

1σij (X) =V (X)1Cijkl (X) =
V (X)

(

(

∂U
∂ηij

)

(3.13)
η′∂2U∂ηij ∂ηkl)(3.14)
η′となる．ここで，η ′ は ηij で偏微分する際に他のひずみ成分を固定することを意味する．これらの微係数を用いて，U を基準状態 X まわりのひずみ ηij について Taylor 展
開すると次式のようになる．

1
U (X, ηij ) = U (X) + V (X)σij ηij + V (X)Cijkl ηij ηkl + · · ·2(3.15)

Lagrange ひずみテンソルの対称性から，式 (3.13) の応力テンソルは対称テンソルである．また，式 (3.14) の弾性係数テンソルはさらにひずみの示数 ij と kl の交換対称性
から Voigt 対称性 (38) と呼ばれる次の対称性を持つ．

Cijkl = Cjikl = Cijlk = Cklij

3.3(3.16)

弾性係数による無負荷平衡点での格子不安定性評価

基準状態 X からの微小ひずみ ηij について，内部エネルギーの展開式である (3.15)の 3 次以上の高次項を省略して書き直すと以下のようになる．

1[U (X, ηij ) − U (X)] − V (X)σij ηij = V (X)Cijkl ηij ηkl
2(3.17)

左辺第 1 項がエネルギー変化，第 2 項が応力 σij が保持されたまま変形すると仮定した場合の外部負荷によってなされる仕事であり，無負荷平衡点での不安定性は弾性係30

数マトリクス Cijkl の正値性に帰着される．弾性係数が無負荷平衡点における応力–ひ∂σij
ずみの勾配を表すものであることを考えると，このクライテリオンは<0，すなわ∂ηkl
ち，変形に対する抵抗力を喪失する点を表すものと解釈できる．弾性係数 Cijkl の対称
部分
SYM
CIJ
=

)
1( T
CIJ + CIJ
2

(3.18)SYM
SYM
を Voigt 表記 (38) したもは Cijklの正値性が系の安定性を支配する (40) ．ここで CIJのであり，T は転置行列を表す．無負荷平衡点での系の安定性は次の 6 × 6 行列の正値性により評価される (37) ．

C11 C12 C13 0
00
C12 C22 C23 000
C13 C23 C33 0
0
0
SYM
CIJ
=
0
0
0 C44 00
0
0
0
0 C55 0
0
0
0
0
0 C66

(3.19)SYM|CIJ| <0 となる条件は，以下の 5 つのいずれかとなる．

C11 C12 C13
C12 C22 C23
C13 C23 C33
C11 C12C12 C22<0<0(3.20)(3.21)

C44 < 0(3.22)

C55 < 0(3.23)

C66 < 0(3.24)

第一式 (3.20) は，膨張による格子の崩壊，即ち体積弾性率が 0 になることを意味している (?) ．第二式 (3.21) は，横方向変形のバランスを表しており，横方向変形が等方変形から非等方変形に変形経路分岐が起こることを意味する．さらに，第三式 (3.22)〜
第五式 (3.24) はそれぞれの変形モードへのせん断不安定が生じることを表している．

31

第 4 章 酸化イットリウム含有 Fe の第一原理格子不安定解析
本章では，酸化イットリウムが Fe 中に固溶することによる特性向上を検討することを目的として，8 格子の bcc-Fe からなるスーパーセルを基に，Y 原子を 2 個，O 原子
を 3 個添加した系において第一原理解析を行い，構造安定性を全自由エネルギー，電子密度および格子不安定性の観点から検討する．Fe 中における Y 原子，O 原子の安定サイトの組み合わせが現時点では不明なため，添加サイトの違いによる構造安定性の変化，力学特性向上への寄与についても議論する．

4.1

解析方法

本研究での解析は全て Kresse らにより開発された平面波基底ウルトラソフト擬ポ
テンシャル法 (28) に基づく第一原理バンド計算コード VASP(42) (Vienna Ab-initio Sim-ulation Package) を用いて行った．交換相関項には局所密度近似 (19) （Local Density
Approximation, LDA）に勾配を考慮した一般化密度勾配近似 (43)（Generalized Gradient Approximation, GGA）を用いた．また収束計算には残差最小化手法 (44) (Residual
Minimization Method – Direct Inversion in the Iterative Subspace, RMM–DIIS) を採用している．カットオフエネルギー，バンド数，波動関数と電子密度の FFT メッシュ
は原子種，原子数と電子数から VASP が定める値 (45) を用いた．
解析は，図 4.1(a) に模式的に示すような bcc の単位セルを 2 × 2 × 2 並べたスーパー
セルを用いて行う．図には Fe 原子を示していないが立方体の頂点と中心に Fe 原子を配置している (bcc). このセルを基本として Y 原子を 2 個，O 原子を 3 個添加する．2
つの Y 原子は青球で示すように〈111〉方向に 1 つの Fe 原子をはさんだ位置に置換す32

る．O 原子については，図 4.2 に示すように bcc 格子の中心原子と頂点原子からなる八面体の中心位置が最も安定なサイトであることが報告されているので (3) ，Y-Y の周囲で考えられる添加サイトの組み合わせから図 4.1(b)〜(f) に示すようなモデル Octa1〜

Octa5 を考慮した．図では O 原子が存在する面を薄く着色して表示している．逆格子
空間の k 点メッシュは Monkhorst-Pack 法 (46) に基づき全ての系で 43 とした．その他の計算条件を表 4.1 にまとめて示す．セル内での原子移動は考えず，図 4.1 の位置に固定した状態で基準となる bcc 格子の寸法を等方的に拡大・縮小して電子状態計算を行い，bcc 格子長 a に対する全自由エネルギーの変化を調べた．そのとき平均応力が零となる格子長さを，系の基準格子長さ a0 として求めた．その後，格子長さ a0 の条件
でイオン固定の条件を外してセル内部での原子構造緩和を行い，得られた構造に対して格子不安定解析を行った．格子安定性の評価に必要な弾性係数は，図 4.3 に模式的に示すようにセルに ±0.05%の引張ひずみ，もしくは ±0.025%のせん断ひずみを与え，ひずみ摂動に対する応力の変化から数値的に求めた．

33

(a)Fe16(b)Octa1(c)Octa2(d)Octa3(e)Octa4(f)Octa5

Fig.4.1 Simulation models.

34

Fig.4.2 Schematic of oxygen doping site.

Table 4.1 Calculation conditions.
Number of atoms
Number of electrons
Number of bands
Cutoﬀ energy (eV)Number of k points

Fe14 Y2 O3 (Octa1) Fe14 Y2 O3 (Octa2) Fe14 Y2 O3 (Octa3)
19
19
19
136
136
136
106
106
106
495.0
495.0
495.0
3
3
4
4
43

Number of atoms
Number of electrons
Number of bands
Cutoﬀ energy (eV)Number of k points

Fe14 Y2 O3 (Octa4) Fe14 Y2 O3 (Octa5)
19
19
136
136
106
106
495.0
495.0
3
4
43

35

∆η11 = r0.01∆η23= r0.01

∆σ11
∆η11
∆σ22
C21 =
∆η11
∆σ33
C31 =
∆η11

C11 =

C44 =

∆σ23
∆η23

∆η22 = r0.01∆η31= r0.01

∆σ11
∆η22
∆σ22
C22 =
∆η22
∆σ33
C32 =
∆η22

C12 =

C55 =

∆σ31
∆η31

∆η33 = r0.01
∆σ11
∆η33
∆σ22
C23 =
∆η33
∆σ33
C33 =
∆η33

∆η12= r0.01

C13 =

C66 =

∆σ12
∆η12

z
y
x

Fig.4.3 Strain perturbation on supercell for evaluation of elastic coeﬃcients.

36

4.2

解析結果と考察4.2.1

全自由エネルギーの比較

格子長さに対する全自由エネルギーの変化を図 4.4 に示す．全自由エネルギーの絶対値は解析セル中の電子数によって変わるので，元素構成が異なる純 Fe については示
していない．また，Octa1〜Octa5 の各セルにおいて全自由エネルギーが最小値となる
格子長さ a0 を表 4.2 に示した．Y-O 原子の添加サイトの組み合わせにおいて，最も低いエネルギーを示したのは Octa5 であった．Octa3 は全域において高いエネルギーをとっており，エネルギー的には最も不利な構造であることが分かる．それ以外の構造を比較すると，格子長さ 6.5Å を境に Octa5 と Octa1 が最もエネルギーの低い構造となった．Octa2 においては 6.5Å 付近，Octa5 においては 6.4Å 付近において，エネル
ギー曲線が大きくぶれる現象が起き，その後のエネルギーが他の構造と比較して急激に増加するという現象が起きている．この原因として，O 原子の添加サイトの違いが考えられる．図 4.1 から分かるように，Octa2,Octa5 では Y 原子が所属する八面体にO 原子が存在しない．置換された Y 原子は他の原子から電子を補っているという報告がされており (47) ，Octa2,Octa5 では格子長さの変化によって Y 原子と O 原子間の電子のやりとりが変化しこのような不連続な挙動を示した可能性がある（図 4.5）．

Table 4.2 Lattice parameter for minimum total free energy in Octa-models.

Lattice parameter，a0

Octa1
6.612

Octa2
6.530

37

Octa3
6.699

Octa4
6.617

Octa5
6.488

Fig.4.4 Comparison of total free energy on Octa-model.

Fig.4.5 Comparison of Maximun of valence electron density on Octa-models.

38

4.2.2格子安定性と弾性係数

安定性判別式である弾性係数 CIJ の小行列式の値 (式 (3.20)〜(3.24) の左辺) を表 4.3に示す．今回解析を行った Fe 単元系・Octa1〜5 の全ての系において全ての値が正値を示し，力学的には存在し得る構造である事が確認できる．さらに詳細な機械的特性の検討として，表 4.4 に弾性係数 CIJ の各成分の値を示す．実材料が存在する Fe 単元
系について実験値と比較すると，実験値 (48) では C11 が 243(GPa),C12 が 138(GPa) であるのに対し，本解析による結果は C11 が 272(GPa)，C12 が 156(GPa) であり誤差はそれぞれ 11.9% ,11.3%である．Fe 単元系と置換が行われた系を比較すると，原子の添加が行われた系では弾性係数の全ての成分が低下している．Y 原子，O 原子の添加が行われた Octa1〜5 において比較すると，各成分で最も高い値をとったのは Octa1，もしくは Octa5 であった．これは前述の全自由エネルギーの比較より得られた結果と同様の傾向を示している．また，Octa2 モデルの C13 成分において値が負となっている
が，その値は小さくほぼ 0 とみなしてよいものと考える．定義した座標によってせん
断弾性係数に負の成分が生じることは可能であり，また小行列式の値は全て正値であるので，不安定とはならない．

39

Table 4.3 Minor determinant of CIJ on Octa-models.

Fe16
Octa1
Octa2
Octa3
Octa4
Octa5

Spinodal(GPa3 )7.94×106
4.77×105
5.25×105
3.69×105
5.09×105
6.93×105

Born(GPa2 )4.96×104
1.58×104
8.97×103
6.32×103
1.07×104
8.74×103

C44 (GPa) C55 (GPa) C66 (GPa)
56.0
56.0
54.0
32.0
24.0
26.0
18.0
34.0
40.0
28.0
22.0
10.0
36.0
30.0
6.0
42.0
40.0
42.0

Table 4.4 Components of CIJ (GPa) of Octa-models.

Fe16
Octa1
Octa2
Octa3
Octa4
Octa5

C11
272
85.0
69.0
43.0
66.0
111

C22
272
224
151
190
199
104

C33
278
57.0
80.0
77.0
65.0
121

C12
156
56.5
41.5
43.0
49.5
53.0

C13
159
44.0
-2.0
24.5
29.0
60.0

40

C23
159
57.0
48.5
50.5
49.0
55.5

C44
56.0
32.0
18.0
28.0
36.0
42.0

C55
56.0
24.0
34.0
22.0
30.0
40.0

C66
54.0
26.0
40.0
10.0
6.0
42.0

4.2.3弾性係数成分と価電子密度の関係

弾性係数が最も高い成分を持つ Octa1，Octa5 および，比較的低い値をとった Octa2について，図 4.6〜4.8 に O 原子の添加が行われた面での価電子密度を示す．各図中上部にはどの面の密度を見ているかを模式的に表示している．価電子密度が著しく低い箇所が見られるが，この位置に Y 原子が存在している．これは価電子数の違いによるもので，Fe の価電子数が 8, O 価電子数が 6 であるのに対し，Y の価電子数が 3 であることに起因する．Octa2 について，図 4.7 右および左において，Fe-O-Fe という非常に価電子密度の高い原子の並びが存在し，この原子群と他の原子との間には図中の点線部において結合力の低下が見られた．図 4.6 の Octa1, 図 4.8 の Octa5 においても同様に Fe-O-Fe 部において価電子密度が非常に高くなっているが，どちらの系においても原子間の結合力の低下は見られなかった．この原因として Octa1 においては，Fe-O-Fe結合部の他に Fe-O-Y 結合部のような Y 原子近傍に O 原子が存在することが挙げられる．Y 原子を置換したことによる価電子密度の減少を O 原子が補う事で結合力の低下を防いでいるものと考えられる．Octa5 については，Octa1 とは異なり Y 原子近傍にO 原子は添加されていないが，Fe-O-Fe 結合が [100][010][001] 方向にそれぞれ 1 つずつ存在し，各方向の変形抵抗を高めていると考えられる．

41

Fig.4.6 Distribution of valence electron density on oxide add plane (Octa1).

42

Fig.4.7 Distribution of valence electron density on oxide add plane (Octa2).

43

Fig.4.8 Distribution of valence electron density on oxide add plane (Octa5).

44

4.2.4O 原子溶解熱の検討

O 原子の溶解熱は，系の全自由エネルギー，O 原子単体での全自由エネルギー，対象となる O 原子を取り除いた系の全自由エネルギーより，以下のように定義される．

△H = Etotal − Eother − Eoxygen(4.1)

O 原子添加サイトを区別して算出したそれぞれの溶解熱を模式図と共に図 4.9 に示す．
Octa4(β) で最も低い-7.876(eV) となり，最大値をとった Octa3(α) の値である-2.319(eV)とは 3.4 倍の開きがあった．各サイトごとに検討を行うと，Octa1(γ),Octa3(α) および(β ),Octa4(γ ) のような Y 原子が所属する八面体にある O 原子添加サイトにおいては同モデルの他の添加サイトの値と比較して高い (絶対値の小さい) 値を取っている．構造間での違いに着目すると，全体的に Octa1,Octa4 の値が低く，逆に Octa2,Octa3 では今回の最小値-7.876(eV) より 50% 以上高い値を持つサイトが存在している．この原因として，高い値をとった Octa2,Octa3 では 3 つの O 原子添加サイト間の距離が近いことがある．両モデル共に (α)-(γ) 間の距離が bcc1 格子長さよりも短いため，O 原子が侵入しにくくなっているものと考えられる．一方，その他のモデルではこのような O
サイトの接近がないので，比較的大きな負の溶解熱エネルギーをもったものと考える．

45

α

β

-3.104(eV)

α
γ

β

γ

-2.884(eV)-2.956(eV)(a)Octa1(b)Octa2γ-4.610(eV)-4.430(eV)-2.463(eV)αβ-6.460(eV)-7.876(eV)-2.319(eV)(c)Octa3(d)Octa4γ-4.287(eV)βα-4.086(eV)-3.909(eV)(f)Octa5

Fig.4.9 Comparison of solution heat on Octa models.

46

第 5 章 Al および Ti 原子添加による系の安定性評価
前章では，bcc-Fe 中の Y 原子および O 原子の置換による系の安定性を検討した．一方，実際の ODS 鋼製造においては Y2 O3 は Al や Ti に代表される第 3 元素と共に添加され，鋼中において YAlO3 や YTiO3 といった複合酸化物として存在していることが確認されている (3) ．
本章では，前章で検討した Fe14 Y2 O3 を基に，Y 原子を Al，また Ti 原子に置換した系について，全自由エネルギー，格子不安定性の両者の観点より安定性を評価し，前章で検討した Fe14 Y2 O3 との比較を行う．同時に，O 原子の溶解熱が系の安定性におよぼす寄与についても検討する．

5.1

解析方法

前章で解析したスーパーセル (図 4.1) を基に，1 つの Y 原子を Al，または Ti に置換した系を解析対象とする（図 5.1，および 5.2）．図では赤球で Al 原子を，緑の球で Ti
原子を表示している．図において，手前下方の Y 原子を第 3 元素に置換した系を a モ
デル，奥側上方の Y 原子を置換したモデルを b モデルとした．Octa5 は，a,b の置換が等価なので区別していない．逆格子空間の k 点メッシュは全て 43 とした．その他の解析条件を表 5.1 にまとめて示す．
以上を用いて系の安定性，置換による特性向上を検討した．

47

(a)Al-octa1a(b)Al-octa1b(c)Al-octa2a(d)Al-octa2b(e)Al-octa3a(f)Al-octa3b(g)Al-octa4a(h)Al-octa4b(i)Al-octa5

Fig.5.1 Images of Al-add model.

48

(a)Ti-octa1a(b)Ti-octa1b(c)Ti-octa2a(d)Ti-octa2b(e)Ti-octa3a(f)Ti-octa3b(g)Ti-octa4a(h)Ti-octa4b(i)Ti-octa5

Fig.5.2 Images of Ti-add model.

49

Table 5.1 Calculation conditions for X-octa-models.
Al-octa
19
136
106
495.0

Number of atoms
Number of electrons
Number of bands
Cutoﬀ energy (eV)50

Ti-octa
19
138
106
495.0

5.2

解析結果と考察5.2.1

全自由エネルギーの比較

格子長さに対する系の全自由エネルギーの変化を図 5.3 に示す．前章と同様に，全自由エネルギーの絶対値は解析セル内の電子数によって変わる値であるため，構成元
素が同じ Al モデル，Ti モデル間での値を比較する．Al を置換した系について，最もエネルギーが低くなったのは octa1a と octa4b であった．a,b モデル間で比較すると，

octa1,octa4 モデルにおいてエネルギーの値に大きな差が生じ，octa2 モデルにおいても比較的小さい差が生じている．その他のモデルについては a,b 間でのエネルギー差はほとんど見られなかった．一方，Ti を置換した系においては octa1a モデルのエネ
ルギーが最小であった．Al を置換した系で最小の値をとった octa4b モデルは Ti を置換した場合にはエネルギーがわずかに上昇した．octa2a,b モデル間のエネルギー差も無くなっており，同一の構造でも置換元素の違いによって相対的なエネルギーの増減が起きることが確認された．Al の場合と同様に，octa1,octa4 モデルでは a,b モデル間において大きなエネルギーギャップを生じた．また，Al,Ti の置換はいずれにおいても，octa1b,3a,3b,4a は上側に位置する高いエネルギーを有するグループに属する．これら 4 つの系はいずれも Y 原子を頂点とする八面体サイトに O 原子が添加されている．

octa1a,octa4b については，その頂点の Y 原子を Al または Ti に置換したものである．
このことから，Al,Ti 原子と比較して Y 原子の近傍には O 原子が侵入しにくい可能性が示唆される．

51

(a)Al-model(b)Ti-model

Fig.5.3 Comparison of total free energy on Al or Ti substitutional models.

52

5.2.2格子不安定性解析による系の安定性評価

detCIJ の 5 つの小行列式の値を表 5.2 および 5.3 に示す．Al を置換した系では全ての値は正となり，力学的に存在し得ることが確認できた．Ti を置換した系では Ti-octa3b
モデルのみ Spinodal 条件式の値が負となり，存在し得ないことがわかる．Ti-octa3b は前節で，エネルギーの高いグループに属していたが，他の 3 つ (Ti-octa1b,Ti-octa3a,Ti-

octa4a) は格子不安定性の観点からは存在し得る系である．第 3 元素の添加による力学
特性変化を詳しく検討するため，弾性係数の各成分の値を表 5.4，および 5.5 に示す．

Ti-octa3b に非常に高い値を取る成分が存在するが，先述のように存在し得ない系であるので，以後の考察では割愛する．表 5.4，表 5.5，および前章にて示した表 4.4 の数値を基に，図 5.4〜5.8 に弾性係数成分を棒グラフで表示した．図 5.4 の Octa1 をみると，全体的に Ti を添加した系の弾性係数が低くなる一方，Al を添加した場合は 9 成
分中 8 成分の値が Y2 O3 のみの系を上回った．図 5.5 の Octa2 では Ti を添加した系のC22 ,C33 成分の値が大きく上昇している（緑色の棒グラフ）ものの，C55 ,C66 では添加前の 30% 以下に減少している．Octa3 においては，図 5.6 に薄い緑色で示す Ti-octa3aにおいて全ての成分の値が Ti 添加によって上昇した．Octa4 では図 5.7 に示すように

C11 ,C22 ,C33 の引張に関する成分が Al-octa4a の C22 成分のみを除いて，Al,Ti の添加によって上昇している．Octa5 では，図 5.8 に示すように Ti を添加すると C44 を除く 8
成分が上昇しており，また，Al,Ti 間で比較すると Al 添加より，Ti の添加がいずれも剛性が上昇している．弾性係数の増加率ではなく，第 3 元素添加前後の大小関係の変化に着目すると，octa1b では 9 成分中 4 成分しか添加前より上昇しなかったが，それ以外の系では 9 成分中 7 成分以上が添加によって弾性係数が上昇しており，定性的ではあるが Al,Ti の添加によって剛性がアップしている．
前節で Al,Ti 置換共に最も低いエネルギーをとった octa1a に着目すると, C11 ，C22
成分において，Y2 O3 のみの系が最大となった (図 5.4 中，薄い色の棒グラフ) ものの，

Al-octa1a,Ti-octa1a に分類すると，Al-octa1a では 7 成分が，Ti-octa1a では 5 成分がY2 O3 のみの系より上昇した．また Al,Ti 間を比較すると，8 成分において Al を添加したものが Ti を添加したものを上回った．この原因について検討すべく，図 5.9 に (01̄1)53

面の価電子密度分布を示した．図上部にはセル中の対象となる面を模式的に表示している．Y2 O3 のみの系と Ti を添加した系では，価電子密度の低い箇所が点線部のようにつながっており，結合が弱い部分が存在する (図 5.9(a),(c))．一方，Al を添加した系
では，図 5.9(b) で矢印で示すように結合が回復している．これは構造緩和の際に O 原子が移動し，価電子密度が上昇したためである．このような価電子密度の変化が Al 添
加時の弾性係数が Ti 添加時より高くなった原因と考えられる．Al 添加では octa1a と同程度の低いエネルギーとなった octa4b について，図 5.7 の弾性係数成分に再び着目すると，Al-octa4b，Ti-octa4b 共に 9 成分中 8 成分が Y2 O3 だけのものよりも上昇している．これらの系について，O 原子添加サイトと Y,Al,Ti 原子が存在する面の価電子密度分布を図 5.10 に示す．図中の矢印部にある密度の低い箇所が，第 3 元素を添加することで小さくなっていることが分かる．octa1a の場合と同様に，このような価電子密度の上昇が弾性係数の上昇に寄与していると考えられる．最後に，弾性係数成分の改善が少なかった octa1b について，価電子密度分布を図 5.11 に示す．図 5.11(a) の矢印部にて存在した，Fe-O-Fe 原子群とその隣の Fe との結合部が，第 3 元素を置換した
場合には存在しないことがわかる．このことが Y 原子のみ置換した系の弾性係数が高い理由であると考えられる．

54

Table 5.2 Minor determinant of CIJ on Al-octa-model.

Fe14 YAlO3 (Al-octa1a)
Fe14 YAlO3 (Al-octa1b)
Fe14 YAlO3 (Al-octa2a)
Fe14 YAlO3 (Al-octa2b)
Fe14 YAlO3 (Al-octa3a)
Fe14 YAlO3 (Al-octa3b)
Fe14 YAlO3 (Al-octa4a)
Fe14 YAlO3 (Al-octa4b)
Fe14 YAlO3 (Al-octa5)

Spinodal(GPa3 )4.58×105
9.55×105
4.61×105
4.24×105
2.83×105
3.67×105
6.11×105
8.65×105
6.55×105

Born(GPa2 )1.07×104
1.58×104
7.80×103
6.90×103
4.21×103
7.20×103
1.23×104
1.51×104
9.07×103

C44 (GPa) C55 (GPa) C66 (GPa)
44.0
36.0
42.0
40.0
34.0
28.0
38.0
40.0
32.0
38.0
30.0
46.0
36.0
32.0
22.0
42.0
34.0
28.0
34.0
24.0
20.0
36.0
30.0
34.0
36.0
34.0
34.0

Table 5.3 Minor determinant of CIJ on Ti-octa-model.

Fe14 YTiO3 (Ti-octa1a)
Fe14 YTiO3 (Ti-octa1b)
Fe14 YTiO3 (Ti-octa2a)
Fe14 YTiO3 (Ti-octa2b)
Fe14 YTiO3 (Ti-octa3a)
Fe14 YTiO3 (Ti-octa3b)
Fe14 YTiO3 (Ti-octa4a)
Fe14 YTiO3 (Ti-octa4b)
Fe14 YTiO3 (Ti-octa5)

Spinodal(GPa3 )5.22×105
1.89×105
3.57×106
3.10×106
8.54×105-1.66×109
4.81×105
1.11×106
1.14×10655

Born(GPa2 )1.30×104
8.59×103
2.10×104
1.83×104
1.21×104
8.57×105
1.20×104
1.71×104
1.30×104

C44 (GPa) C55 (GPa) C66 (GPa)
40.0
28.0
40.0
36.0
26.0
24.0
32.0
10.0
10.0
32.0
10.0
10.0
42.0
26.0
36.0
38.0
46.0
28.0
24.0
14.0
6.0
38.0
44.0
24.0
38.0
40.0
42.0

Table 5.4 Components of CIJ (GPa) of Al-octa-model.

Al-octa1a
Al-octa1b
Al-octa2a
Al-octa2b
Al-octa3a
Al-octa3b
Al-octa4aAl-octa4b
Al-octa5

C11
77.0
83.0
78.0
76.0
34.0
58.0
69.0
98.0
112

C22
203
229
134
132
162
169
221
194
107

C33
94.0
84.0
102
97.0
92.0
84.0
78.0
93.0
105

C12
70.5
56.5
51.5
56.0
36.0
51.0
54.0
62.5
54.0

C13
61.5
40.0
57.0
47.5
27.0
41.5
43.0
56.5
53.5

C23
72.5
56.0
48.5
58.0
49.0
57.0
51.0
58.0
50.0

C44
44.0
40.0
38.0
38.0
36.0
42.0
34.0
36.0
36.0

C55
36.0
34.0
40.0
30.0
32.0
34.0
24.0
30.0
34.0

C66
42.0
28.0
32.0
46.0
22.0
28.0
20.0
34.0
34.0

C55
28.0
26.0
10.0
10.0
26.0
46.0
14.0
44.0
40.0

C66
40.0
24.0
10.0
10.0
36.0
28.0
6.0
24.0
42.0

Table 5.5 Components of CIJ (GPa) of Ti-octa-model.

Ti-octa1a
Ti-octa1b
Ti-octa2a
Ti-octa2b
Ti-octa3a
Ti-octa3b
Ti-octa4a
Ti-octa4b
Ti-octa5

C11
79.0
48.0
87.0
82.0
83.0
830
69.0
92.0
132

C22
196
197
274
253
197
2358
226
218
131

C33
72.0
40.0
205
201
106
-1451
69.0
91.0
133

C12
49.5
29.5
53.0
49.5
65.0
1118
60.0
54.5
65.5

56

C13
47.5
28.5
49.5
46.5
49.0
-212
42.0
44.0
65.5

C23
53.5
31.5
71.5
62.0
69.5
387
60.5
55.0
67.5

C44
40.0
36.0
32.0
32.0
42.0
38.0
24.0
38.0
38.0

(GPa)

300
250
Octa modelAl-a modelAl-b model
Ti-a modelTi-b model

200
150
100
50
0

C11

C22

C33

C12

C13

C23

C44

C55

C66

Fig.5.4 Comparison of each component of CIJ on Octa1 model.

(GPa)

300
250
Octa modelAl-a modelAl-b model
Ti-a modelTi-b model

200
150
100
50
0

C11

C22

C33

C12

C13

C23

C44

C55

C66

Fig.5.5 Comparison of each component of CIJ on Octa2 model.

57

(GPa)

300
250
Octa modelAl-a modelAl-b model
Ti-a model

200
150
100
50
0

C11

C22

C33

C12

C13

C23

C44

C55

C66

Fig.5.6 Comparison of each component of CIJ on Octa3 model.

(GPa)

300
250
Octa modelAl-a modelAl-b model
Ti-a modelTi-b model

200
150
100
50
0

C11

C22

C33

C12

C13

C23

C44

C55

C66

Fig.5.7 Comparison of each component of CIJ on Octa4 model.

58

(GPa)

300
250
Octa model
Al modelTi model

200
150
100
50
0

C11

C22

C33

C12

C13

C23

C44

C55

C66

Fig.5.8 Comparison of each component of CIJ on Octa5 model.

59

(a)Octa1(b)Al-octa1a(c)Ti-octa1a

Fig.5.9 Distribution of valence electron density on (01̄1) plane (Octa1a).

60

(a)Octa4(b)Al-octa4b(c)Ti-octa4b

Fig.5.10 Distribution of valence electron density on O-add plane (Octa4b).

61

(a)Octa1(b)Al-octa1b(c)Ti-octa1b

Fig.5.11 Distribution of valence electron density on O-add plane (Octa1b).

62

5.2.3O 原子溶解熱が安定性におよぼす影響

O 原子の溶解熱を添加サイトごとに表 5.6, および 5.7 に示し，前章の結果 (図 4.9) とあわせてその大小関係をグラフ表示したものを図 5.12〜5.16 に示す．図中の右には，前章と同じく α,β,γ がそれぞれどの O 原子の添加サイトでの溶解熱であるかを模式的に示している．また棒グラフは，Y が Y2 O3 のみの系での溶解熱であるが, Al,Ti 置換については，図中 ⃝A に添加されているものを a-model, ⃝B に添加されたものを b-model として示している．前節同様，Octa5 では ⃝A ，⃝
B の区別がないので棒グラフに a,b-modelの区別はない．また，Octa3 では Ti の b-model が不安定であるので除いている．第 3元素の置換の有無に着目すると，Octa1,Octa4,Octa5 においては Al,Ti を添加すると負の溶解熱の絶対値が小さくなり，結合が弱くなっていることが示唆される．これらの系は前章において，溶解熱の絶対値が大きく，結合が強いとされたものであり，エネ
ルギー的に見ても，前章の 5 つの Octa モデルのうち低い値をとった 3 つである．一方，Octa2 においては Al を ⃝B サイトに置換した時，著しく負の溶解熱が大きくなっており，この構造の中では最も分離しにくいことが示唆される．Octa3 でも，Al,Ti 置換共に添加によって負の溶解熱が大きくなっていることが確認できる．5.2.1 節において全
自由エネルギーが低い値をとった octa1a，および octa4b に着目すると，図 5.12，5.15にみられるように，第 3 元素の添加によって溶解熱が正の方向に変化し，結合が弱くなっている．しかしながら，比較的高いエネルギーであった Octa2,Octa3 では，第 3元素の添加が行われることで酸素の負の溶解熱は各サイトで大きくなり，結合が強くなっている．すなわち，第 3 元素の添加によって，元々エネルギー的および溶解熱の観点からは結合が強かった構造での結合力は弱められるが，それ以外のサイトでの結合力が強くなっているので，酸素は Y 原子の近傍に様々な形態で結合できる可能性がある．

63

Table 5.6 Solution heat(eV) of Al-octa-model.
a-modelAl-octa1a
Al-octa2a
Al-octa3a
Al-octa4aAl-octa5

αsite-3.062-2.783-6.757-3.533-3.159βsite-3.062-3.094-6.061-2.975-3.159γsite-4.143-3.174-8.600-1.648-3.159b-modelAl-octa1b
Al-octa2b
Al-octa3b
Al-octa4b

αsite-3.745-5.047-6.605-2.975βsite-3.745-5.954-9.681-3.533γsite-1.917-5.711-9.291-4.104

Table 5.7 Solution heat(eV) of Ti-octa-model.
a-modelTi-octa1a
Ti-octa2a
Ti-octa3a
Ti-octa4a
Ti-octa5αsite-3.261-2.626-10.642-3.568-3.011βsite-3.261-2.877-8.064-3.568-3.011γsite-4.323-3.116-10.527-1.806-3.011b-modelTi-octa1b
Ti-octa2b
Ti-octa3b(Unstable)
Ti-octa4b

64

αsite
-3.7442.542-8.821-3.091βsite
γsite-3.744 -2.005-2.803 -3.033-11.758 -11.456-3.091 -4.116

Fig.5.12 Comparison of solution heat on Octa1 model.

Fig.5.13 Comparison of solution heat on Octa2 model.

65Fig.5.14 Comparison of solution heat on Octa3 model.

Y
Al

α

Ti

γ

Y
Al

β

Ti

α

Y
Al

β

γ

Ti

Fig.5.15 Comparison of solution heat on Octa4 model.

66

Y
Al

α

Ti
Y
Al

γ
β

Ti

α
β

Y
Al

γ

Ti

Fig.5.16 Comparison of solution heat on Octa5 model.

67

第6章

結論

本研究では，ODS 鋼における酸化物分散強化機構について電子レベルから非経験に知見を得ることを目的として，bcc-Fe に Y2 O3 ，YAlO3 ，YTiO3 を添加した系の第一原理計算を行い，その力学的安定性を弾性係数の正値性に基づく格子不安定性の観点から検討した．また，従来の構造安定性評価で用いられる自由エネルギーや溶解熱の大小についても考察した．得られた結果を以下に総括する．
第 2 章では，平面波基底での第一原理分子動力学法の基礎理論および電子状態計算の高速化手法について説明し，第 3 章では，格子安定性の概念と，弾性係数を用いた
格子不安定解析について述べた．
第 4 章では，酸化イットリウムが Fe 中に固溶することによる特性向上を, O 原子の添加サイトの違いから検討するために，bcc の単位セルを 2 × 2 × 2 並べたスーパーセルに Y 原子を 2 個，O 原子を 3 個添加した系の第一原理計算を行った．2 つの Y 原子
を bcc の第 3 近接位置に固定し，その結合の間の 3 つの O 原子の添加サイト (bcc の八面体中心) の組み合わせを 5 通り考慮して解析を行い，Y 原子と O 原子の間の電子のやり取りによって全自由エネルギーが増減すること，価電子密度の変化が弾性係数の大小に関係があること，などを明らかにした．弾性係数の行列式は全ての組み合わせ
で正の値となり，いずれも力学的には存在し得ることを示したが，弾性係数の値は Fe
単体よりもいずれも低下しており，Y2 O3 の添加が直接剛性向上に結びつくわけではな
いことが示された．また，O 原子の溶解熱は全自由エネルギーが低い系ほど大きな負の値を持ち，結合が強いことが分かった．
第 5 章では，固溶における第 3 元素の役割を検討するために，先の Y-Y の結合を

Y-Al または Al-Y のように Al および Ti に置換して全自由エネルギー，弾性係数，O 原
子の溶解熱を調べた．同じ価電子構造の系であっても Y,Al/Ti,O 原子の配置によって68

全自由エネルギーが変わるのは当然であるが，特に，Y 原子を頂点とする八面体サイ
トに O 原子を添加した系については全体的にエネルギーが高くなり，エネルギー的には不利な構造であることを示した．また，第 3 元素の添加による弾性係数の変化も構造によって様々であるため一概には言えないが，Al,Ti の添加によって Y2 O3 だけの場
合よりも弾性係数が増加し剛性が上昇するケースが多く見られた．また，O 原子の溶
解熱は，前章で全自由エネルギーが低く溶解熱が大きいことが示されたサイトでは添加によって溶解熱が小さくなったが，溶解熱が小さかったサイトでは逆に添加によって結合が強められることが分かった．このことは，第 3 元素の添加により，酸素は Y原子の近傍に様々な形態で結合できることを示している．
以上のように，Y,O,Al,Ti の添加では Fe 単元系を超える弾性係数をもつ系は今回はなく，単純に剛性アップするわけではないが，第 3 元素の添加によって O 原子が様々な酸化物構造をとり得ること，Y2 O3 のみの添加よりも剛性が上昇することなどを確認した．実際の ODS 鋼中におけるこれらの元素割合は今回の解析とは桁違いに小さく，
本来はもっと多数の原子の中で構造評価すべきであるが，第一原理計算で直接取り扱
うのは現実的には不可能である．今後，経験的な原子間ポテンシャルを用いた MC シ
ミュレーションなどに本解析結果が反映されることを期待する．

69

付録 A 関連学術講演
山本智，横川望，屋代如月，冨田佳宏”イットリウム酸化物添加鉄の第一原理格子不安定解析”
第 13 回分子動力学シンポジウム，鹿児島大学，(2008.5)

山本智，横川望，屋代如月，冨田佳宏”第一原理計算による Fe 中のイットリウム酸化物の格子不安定性解析”
M&M2008 材料力学カンファレンス，立命館大学，(2008.9)

山本智，横川望，屋代如月，冨田佳宏”鉄・イットリウム酸化物合金の第一原理格子不安定性解析”
日本機械学会第 21 回計算力学講演会，琉球大学，(2008.11)

山本智，横川望，屋代如月”第一原理計算による Y・Al 酸化物分散 Fe の構造安定性評価”
日本機械学会第 22 回計算力学講演会，金沢大学，(2009.10)

70

md

71

mm

72

cmd

73

cmd

74

参考文献
(1) H.Kishimoto, R.Kasada, O.Hashimoto, and A.Kimura Journal of Nuclear Materials, 386 (2009),533.
(2) Krajnikov, A.V., Demidik, A.N. , and H.M. Ortner, Materials Science and Engineering A, 234 (1997),357.
(3) Yong Jiang. John R, Smith. and G,R,Odette. Phys. Rev. B, 79 (2009),064103.
(4) Schaublin, R. Ramar, A. Baluc, N. V.de Castro, Monge, M.A. Leguey, T.
Schmid, N. and Bonjour, C., Journal of Nuclear Materials, 351 (2006),247.
(5) 尾方成信, 北川浩, 第一原理分子動力学計算による窒化アルミニウム／アルミニ
ウム界面の構造評価日本金属学会誌, 60, (1996), 1079.

(6) 北村隆行, 梅野宜崇, 西田昌弘, 原子鎖の強度に関する第一原理解析, 日本機械学会論文集 A 編, 67, (2001), 203.

(7) Stefano Curtarolo, Dane Morgan and Gerbrand Ceder Calphad, 29 (2005),163.
(8) Ching, W. -Y. and Yong-Nian Xu, Journal of Magnetism and Magnetic Materials, 209 (2000),28.
(9) Born, M., and Huang, K., Dynamical Theory of Crystal Lattices, (1954) Oxford
UP.
(10) Wang, J., Yip, S., Phillpot, S.R., and Wolf, D., Phys. Rev. B, 20 (1993), 4182.
(11) 屋代如月, 大穂正史, 山上勝也, 冨田佳宏, 材料, 52, (2002), 241.
(12) Šob M, Wang LG, Vitek V, Mat Sci. Eng. A, 234 (1997), 1075.
(13) Luo,W., Roundy,D., Cohen,M., Morris Tr,J.W., Phys. Rev. B, 66 (2002),
094110.
(14) Clatterbuck, D.M., Chrzan, D.C., Morris Jr, J.W., Phil. Mag, Lett, 82 (2002),
141.
75

(15) Clatterbuck, D.M., Chrzan, D.C., Morris Jr, J.W., Acta Materialia, 51(2003),2271.
(16) Černý,M. Pokluda, J. J.Alloys and Compounds, 378 (2004),159.
(17) Hohenberg, P., and Kohn,W., Phys. Rev., 136, (1964), 864.
(18) Kohn, W., and Sham,L.J., Phys. Rev., 140, (1965), 1133.
(19) Perdew, J.P., and Zunger,A., Phys. Rev.,B, 23, (1981), 5048.
(20) キッテル，固体物理学入門 上 第 7 版 (1998), 丸善.
(21) Bachelet,G.B., Hamann,D.R. and Schlüter,M., Phys. Rev. B, 26, (1982), 4199.
(22) Troullier,N. et al., Phys. Rev. B, 43, (1991), 1993.
(23) L.D. ランダゥ，E.M. リフシッツ著，佐々木健，好村磁洋訳，量子力学 1，(1992),
東京図書.

(24) 上田顕，コンピュータシミュレーション, (1990), 朝倉書店.
(25) 香山正憲，固体材料の電子状態の基礎, (1993), 材料学会勉強会資料.
(26) Nielsen,O.H., and Martin,R.M., Phys. Rev. B, 32, (1985), 3780.
(27) Dacosta,P.G., Nielsen,O.H., and Kunc,K., J. Phys. C : Solid State Phys., 19,(1986), 3163.
(28) Vanderbilt,D., Phys. Rev. B, 4, (1990), 7892.
(29) Kleinman,L., and Bylander,D.M., Phys. Rev. Lett., 48, (1982), 1425.
(30) 尾口多美夫，佐々木泰造，固体物理, 25, (1990), 857.
(31) Fu,C.L., and Ho,K.M., Phys. Rev. B, 28, (1983), 5480.
(32) Teter,M.P. et al., Phys. Rev. B, 40, (1989), 12255.
(33) King-Smith,R.D., and Vanderbilt,D., Phys. Rev. B, 49, (1994), 5828.
(34) Gillan,M.J., J. Phys. Condens. Matter., 1 (1989), 689.

76

(35) Stich,I. et al. , Phys. Rev. B, 39, (1989), 4997.
(36) Bylander,D.M. et al., Phys. Rev. B, 42, (1990), 1394.
(37) Milstein,F., Phys. Rev. B, 3 (1971), 1130.
(38) Wallace, D.C., Thermodynamics of Crystals, (1972), Wiley, Newyork.
(39) 屋代如月，学位論文，金属結晶の変形と破壊の分子動力学シミュレーション並びに局所格子不安定性解析，(1998).

(40) Wang,J., Li,J., Yip,S., Phillpot,S.R., and Wolf,D., Phys. Rev. B, 52 (1995),
12627.
(41) Li, J., and Wang, T., J. Phys. Condens. Matter 10, (1998), 9889.
(42) Kresse,G. and Furthmüller,J., Phys.Rev.B, 54, (1996)，11169 .
(43) Langreth,D.C., and Perdew,J.P., Phys. Rev. B, 21, (1980), 5469.
(44) Pulay,P., Chem.Phys.Lett., 73, (1980)，393 .
(45) Kresse,G., and Furthmüller,J., VASP the GUIDE, (2003).
(46) Monkhorst,H.J., and Pack,J.D., Phys. Rev. B, 13, (1976), 5188.
(47) Gopejenko,A., et al. , 25th Syposium on Fusion Technology, (2008).
(48) 金属データブック，日本機械学会編 (1998), 丸善.

77

謝辞
本研究を進めるにあたり，厚くご指導賜りました屋代如月准教授に心より感謝の意を表します．浅学な著者が論文執筆を完遂することができたのも先生に休日，平日を問わず親切丁寧にご指導いただたいた事によるものでした．阪上隆英教授，長谷部忠司准教授，冨田佳宏教授（現・福井工大）には研究内容について様々なご助言を頂き，研究者としての心構えを教えていただきました．心より感謝致します．研究室の環境づくりにご尽力くださり、よく声を掛けてくださいました古宇田由夫技官に感謝致します．
研究室配属時の何も分からない著者に対して良き指針となり温かく接して頂いた諸先輩方に感謝致します．また，第一原理グループとして常に著者に助言を下さり，時にお互いの意見をぶつけあい議論した横川望氏 (現・西日本旅客鉄道)，久馬雅彦氏（現・三菱重工業）に深くお礼申し上げます．お二人のおかげで研究活動をスムーズに行うことができたものだと思っております．3 年間の研究室生活において公私共に多くの時間をすごし，切磋琢磨しあいながら貴重な経験をしてきた同輩の大石直樹君，金谷
敬輔君，北村真留久君，坂本忠鷹君，富永智之君，平山敬史君，山本隆史君に感謝致します．さらに，有意義な議論を交わし，共に困難を分かち合ってきた ODS 班の睦門
賢憲君，田仲稔君，山口明宏君をはじめとする後輩の皆さんに感謝します．
最後に学生生活を送るにあたり経済的・精神的に支えてくださった両親，家族に厚くお礼を申し上げます．ありがとうございました．

78