要約ODS 鋼における，Y2 O3 添加による機械的特性向上について，電子論ならびに格子力
学的な観点から新たな知見を得るために，第一原理計算により酸化物を添加した bcc-Feの構造安定性評価，引張/せん断変形解析を行った．
まず bcc-Fe 単位格子を 3 × 3 × 3 並べたスーパーセルの 2 つの Fe を Y または Ti に置換した構造について，3 つの O の添加サイトの組み合わせを変えて，全自由エネルギー，
弾性特性，O の溶解熱から添加構造の安定性を評価した．今回解析した構造はすべて格子力学の点からは安定となったが，弾性異方性や O 溶解熱からは bcc 格子において原点に Fe，±1/2(a, a, a) に Y または Ti，1/2(−a, −a, 0), 1/2(−a, 0, a), 1/2(0, a, −a) に O を添加する構造が有利であった（a は bcc の格子定数）．これは Y-Fe-Y,(Ti) の中心軸に対して，O を含む 3 つの Fe 八面体の短軸方向がそれぞれ [100],[010],[001] に配向し，かつ中央の Fe で頂点が共有される構造である．次に上述の構造で Y2 O3 を添加した系およびFe 単元系に [001] 引張/[111] せん断解析を行い，価電子密度，格子不安定性の観点から理想強度を議論した．引張においては，スーパーセル中での原子構造緩和による各原子の結合間距離の変化も検討した．格子不安定性の観点からは理想引張強度は Fe 単元系
が ε = 0.03, σ = 3.8GPa，Y2 O3 を添加した系は ε = 0.11 ∼ 0.16, σ = 13.9 ∼ 16.8GPaとなった．先述のように八面体短軸方向の強い Fe-O-Fe 結合が [100],[010],[001] の各方向に存在するために，3 軸引張に近い ε′ = 0 の引張において高い抵抗を示したと考えられる．格子不安定性から算出した理想せん断強度は，Fe 単元系が ε = 0.25, σ = 10.1GPa，Y2 O3 を添加した系は ε = 0.30 または 0.29, σ = 8.7 または 10.0GPa となった．

2

Summary
For a new insight on the mechanical properties of ODS(Oxide Dispersion Strengthened) steels, we have evaluated the stable site for Y2 O3 or YTiO3 in bcc-Fe, ideal
tensile and shear strength of these Fe with oxides, based on the free energy, solution
heat and lattice stability by using the ab-initio DFT calculation.
First, we performed DFT calculations for bcc-Fe contains Y2 O3 or YTiO3 , by using
a supercell consist of 3×3×3 bcc unit lattices. Three models are considered for oxygensites around Y-Fe-Y or Y-Fe-Ti axis. It is revealed that the most stable structure is
the model that locates Fe atom at (0,0,0), Y or Ti atoms at ±1/2(a, a, a), O atoms
1/2(−a, −a, 0), 1/2(−a, 0, a), 1/2(0, a, −a), respectively, here a is the lattice parameterof bcc lattice. That is, the center axis of Y-Fe-Y or Y-Fe-Ti bond is surrounded by
the triangles of O atoms, and the three octahedrons of bcc-Fe containing O atoms ofwhich short axis is directed to [100],[010],[001], respectively, are shared at the centerFe atom of Y-Fe-Y or Y-Fe-Ti axis.
Second, we have performed the [100] tensile/[111] shear simulations on the Y2 O3 Fe models and pure Fe, and discussed the change in the valence electron density,
ideal strength based on the lattice instability criteria, etc. For the case of tensile
simulation, we have also evaluated the change in the bond length if we considered not
only the electron density calculation but also structual relaxation by atom motion.
The elastic limit or ideal tensile strength against [001] tension (ε′ = 0) is evaluatedas ε = 0.03, σ = 3.8GPa for Fe and ε = 0.11 ∼ 0.16, σ = 13.9 ∼ 16.8GPa forY2 O3 -Fe models. As mentioned above, the strong Fe-O-Fe bonds are oriented in the[100],[010],[001] directions in Y2 O3 -Fe model, so that these bonds show high resistanceagainst the tensile condition of ε′ = 0, which is close to the hydrostatic tension. The
elastic limit for shear is calculated as ε = 0.25, σ = 10.1GPa for Fe, ε = 0.30 or 0.29,
σ = 8.7 or 10.0GPa for Y2 O3 -Fe models.

修士論文

第一原理計算ならびに
格子不安定解析による酸化物分散強化メカニズムの基礎的研究

指導教員：屋代 如月

田仲 稔

2012 年 2 月
神戸大学大学院 工学研究科 博士課程前期課程 機械工学専攻

Fundamental Study on Oxide Dispersion
Strengthening Mechanism:
Ab-initio DFT Calculationand Lattice Instability Analysis

Minoru TANAKA
February 2012
Department of Mechanical Engineering,
Graduate School of Engineering,
Kobe University, Kobe, Japan

目次
1

緒論

1

2

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

4

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

4

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

5

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

7

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

8

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

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

13

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

14

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

15

2.8.1

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

16

2.8.2

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

20

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

22

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

22

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

23格子不安定性解析

28

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

28

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

30

3.3 応力-ひずみ関係と弾性剛性係数 . . . . . . . . . . . . . . . . . . . . . .

31

3.4 弾性剛性係数による格子不安定性評価

34

3

i

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

目次

4

ii

酸化物含有 Fe の安定性評価

36

4.1 解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

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

41

4.2.1

全自由エネルギーと弾性剛性係数行列式の比較. . . . . . . . .

41

4.2.2

弾性剛性係数行列の各成分の比較 . . . . . . . . . . . . . . . . .

44

4.2.3

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

46

4.2.4

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

48

4.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

Y2 O3 含有 Fe の引張シミュレーション

54

5.1 解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

54

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

55

5.2.1

応力-ひずみ関係と弾性剛性係数の行列式の変化 . . . . . . . . .

55

5.2.2

引張変形下の価電子密度分布の変化 . . . . . . . . . . . . . . . .

58

5.2.3

引張変形下の原子間結合距離の変化 . . . . . . . . . . . . . . . .

62

5.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

Y2 O3 含有 Fe のせん断シミュレーション

66

6.1 解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

66

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

69

6.2.1

応力-ひずみ関係と各方向の応力変化 . . . . . . . . . . . . . . .

69

6.2.2

弾性剛性係数の行列式の変化 . . . . . . . . . . . . . . . . . . .

71

6.2.3

せん断変形下の価電子密度分布の変化. . . . . . . . . . . . . .

73

6.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

7

結論

82

A

関連学術講演

85

5

6

参考文献

94

目次
謝辞

iii
100

第1章 緒 論
原子炉の高燃焼度化は限られた資源の有効利用や原子力エネルギーの高効率利用に欠かせない技術の一つである．高燃焼度化，安全性向上のための有効な手段の一つとして燃料被覆管材料の改良が進められている．燃料被覆管には優れた耐中性子照射特性，耐食性，高温強度などが求められ，これまで軽水炉炉心材料としては，ジルコニ
ウム合金やオーステナイト鋼が用いられている．ジルコニウム合金は中性子吸収性に優れているものの高温での強度に課題があり，より一層の高燃焼度化は難しい．また，
オーステナイト鋼も加工性に優れているが高温での寸法安定性やヘリウムおよび中性子照射による脆化や耐食性に課題を抱えている．酸化物分散型強化鋼（ODS 鋼:Oxide

Dispersion Strengthened）は，オーステナイト鋼に比べ耐照射性，耐ヘリウム脆化性に優れていながら，高温において材料組織が安定しているため，次世代原子炉の燃料被覆管材として最も有力な材料である．また，原子力システム用の材料としてばかりではなく，火力プラント用配管など同様の性能を要求される分野にも利用可能であり，高温強度と耐食性が要求される自動車用鋼板，鉛電池セルの隔壁材料，核融合炉のブ
ランケット配管としても期待され，実用化に向け溶接，量産技術の開発や耐照射性，耐食性の実証試験を進める段階である (1) ．

ODS 鋼は高温でも安定な Y2 O3 を添加し，メカニカルアロイングなどの製造プロセ
スにより，nm オーダーの微細な酸化物粒子を均一かつ緻密に分散させた材料であり，高温での高いクリープ強度 (2) や耐照射性を有する (3)(4) ．微小な酸化物が転位の移動を抑制することで ODS 鋼は高い強度を有していることが示唆されており (5) , これまで，

SEM や TEM，X 線回析といった手法により製造プロセス中の組織の観察が試みられている (6) ．ODS 鋼は熱間押出時に微細な酸化物が α 相から γ 相への変態を抑制し，組織中に酸化物が微細に分散している α 相が残留している (7)(8) ．熱処理後の ODS 鋼中では残留 α 相とマルテンサイト相の固さの異なる 2 相が存在し，残留 α 相が強化相と1

なり，複合材料のような効果を発揮することで，高温クリープ特性を発揮していると考えられており (9)(10) ，高温でのクリープ変形は比較的弱い粒界やパケット粒界で起こると考えられている (11) ．また中性子照射後でも Y 酸化物の分散状態，形状，組成は安定しており (12) ，照射条件での高い引張強度は YTi 酸化物の密度に起因していることが報告されている (13) ．Y の割合や Y2 O3 以外に添加する元素割合を変化させることで試行錯誤的に ODS 鋼の組織制御や力学特性の検討が行われている (14) -(16) ．Ti が添加されることで，YTi の複合酸化物を形成し (17) ，析出する酸化物が微細化し (18)(19) ，さらに熱処理後の残留 α 相の生成が容易になる (9) ．Al が添加されることにより高い耐食性を有するようになるが，酸化物粒子を粗大化，分散密度を低下させ，残留する α相を減少させる (20) ．現在は Hf や Zr を添加することで Al による酸化物の凝集を妨げ，結晶粒界に炭化物や酸化物を形成することで，粒界すべりを抑制し高温での強度を向上させた「スーパー ODS 鋼」が発明されるなど (21) ，多く研究により ODS 鋼についての知見が得られている．
計算材料科学の分野に目を向けると転位と析出物の相互作用というテーマで分子力学法 (Molecular dynamics：MD) や離散転位動力学 (Discrete dislocation dynamics：

DDD)(22)(23) シミュレーションにより bcc-Fe 中の銅析出物と，刃状およびらせん転位の相互作用についてのシミュレーションや (24)(25) , 析出物に刃状転位およびらせん転位を衝突させるシミュレーション (27) , ODS 鋼中における分散粒子をモデル化した試みが行われている (26) ．これまで著者らのグループでは原子レベルで Y や O による分散や力学特性への寄与について知見を得るために第一原理計算により Y2 O3 を添加した Fe の安定性ならびに力学特性を議論してきた (28) ．原子種と原子配置のみを必要情報とし，量子力学にのみ基づいて非経験的に材料物性を予測可能な第一原理計算は，従来の実験による試行錯誤的な手法によらず，電子論に基づいて非経験的に材料特性を評価する手法として注目されており，実材料の特性を検討した例も報告されている (29) -(31) ．

ODS 鋼を対象とした第一原理計算も行われており，Fe マトリクス中に Y,O,Ti を添加し，各原子間の結合エネルギーや Y2 O3 析出物の安定性をエネルギーにより評価している (32) -(34) ．しかしながら膨大な計算時間とメモリを必要とする第一原理計算は，現2

時点では数十〜数百程度の原子数での解析が限界であり，変形時におけるイットリウ
ムの役割などを直接扱うことは難しく，転位運動への影響等を直接評価することは現時点では不可能である．第一原理計算によって Y2 O3 析出物の安定性の評価も試みられているが，Y2 O3 の Bixbyite 構造を完全に再現するのではなく，構造の一部を模擬したモデルで解析が行われている (35) .
少数の原子の解析でも，力学特性に重要な知見を与えるものに格子不安定性解析 (36)がある．格子不安定とは，外力下で変形している結晶格子が釣り合いを失い，外力の増加を必要とせずに変形が進行する状態を示している．本来，完全均一結晶の弾性安定限界を評価するために提案されたものであるが，多原子系における転位発生等の局
所変形開始の有力な基準となることが報告されている (37) ．格子不安定解析は，少数の原子しか扱えないが大変形下における電子状態を精密に評価できる第一原理計算に適した解析である．第一原理計算による不安定解析の例としては W,Mo,Nb,Fe といった

bcc 金属の理想強度評価 (38)(39) ，Si と Al の [001] 引張下の格子不安定性 (40) ，fcc 金属(Al,Cu,Ag) の安定性 (41) 等がなされてきた．
本論文では酸化物分散強化メカニズムによる知見を得るべく，bcc-Fe を基に Y,O を添加した系に対して静力学的な変形シミュレーションならびに格子不安定解析を行う．
第 2 章では第一原理分子動力学法の基礎理論，および電子状態計算の高速化手法について説明する．第 3 章では，格子不安定性解析の概要と，弾性係数を用いた格子不安定性評価について説明する．第 4 章では，bcc-Fe に Y2 O3 を添加した系において解析を行い，Y 原子および O 原子の添加による特性向上を格子不安定性の観点から検討すると共に，全自由エネルギー，O 原子の溶解熱，異なる O 原子添加サイトの組み合わせによる弾性係数の違い等について検討する．第 5 章では，最もエネルギー的に安定なモデルに対して，[001] 方向引張シミュレーションを行い，弾性限界や価電子密度の変化等について検討する．第 6 章では，[111] 方向にせん断シミュレーションを行い，第 5 章同様に弾性剛性係数，価電子密度分布の観点から検討する．第 7 章では本研究で得られた結果を総括する．

3

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

2.1

断熱近似と平均場近似

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

4

に，通常，以下の 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)と表せることを明らかにした (48) ．右辺の各項はそれぞれ，原子核による電子のポテン
シャルエネルギー，相互作用する多電子系での電子の運動エネルギー，電子間クーロン相互作用エネルギー，他の全ての電子間多体相互作用を表す交換相関エネルギーである．この Etot を最小にする ρ(r) が基底状態での電子密度分布となる．
相互作用のない系での電子の状態を表す波動関数 (電子波動関数) を ψi とし，その運5動エネルギー 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 方程式) が得られる (49) ．

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 に解く問題に帰
着される．

6

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の関数形 (50) を示す．

ε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)
7

となる．

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 の定理 (51) より

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

k=

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

(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 点についての和をとることを表す．

9

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

10

(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′ )11(2.34)

= 4π

∑

∫

PP,nlocVa,l(r)jl (|k + G|r)jl (|k + G′ |r)r2 dr(2l + 1)Pl (cos ω)l(2.35)となる (54) ．ここで，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′ )12(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 の方法 (55) によって表したもので，[

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
13

]

+
−

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 の発散項とうまく打ち消し合うため，次式のように表すことができる (56) ．

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 というスケーリングを適用し，それを対応するひずみテンソルの成分で14微分することによって得られる (57)(58) ．Ω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

擬ポテンシャル法

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

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

16

(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 −17(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 乗のループを

18

含んでおり計算時間がかかるとともに，記憶する量も平面波数の増加に対してその 2乗で増える．そのため，大規模な計算ではすぐにメモリー容量に破綻をきたす．そこ
で，非局所項に次式で表される KB 分離型表現 (60) を用いれば，平面波の 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′

19

(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

20

′
=
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

となる．

21

(2.88)(2.89)

2.9

電子占有数

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

22となる．同様に

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)) の対角化を繰り返す方法であるが，この方法では対象とする系によっては多大な計算労力を必要とする．そこで，近年電子状態計算を効率的に行う方法が開発された (63)−(66) ．本節では
共役勾配法についてその概要を示す．

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

23適当な初期値 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 法 (63) について解説する．

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

方向) を求める．それをもとにして ε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′

25

(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
26

} 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 を行う．

27

第3章

格子不安定性解析

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

3.1

不安定条件

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

28

a3
a3
a2
a1

a4
a5 a6

a2
a1

Fig.3.1 Unit cell of fcc lattice.

{am } の変形状態下にある結晶の安定性は，以下のように微小変形増分 {∆am } によるエネルギーの変化を考えることによって求められる (36)(68) ．状態 {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 ] の正値性に帰着される．

29

3.2

応力と弾性係数

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

dU = T dS − dW

(3.4)

である (69) ．ここで，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 =

30

これより，式 (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 対称性 (69) と呼ばれる次の対称性を持つ．

Cijkl = Cjikl = Cijlk = Cklij

3.3(3.16)

応力-ひずみ関係と弾性剛性係数

一般に弾性係数は応力–ひずみ関係の勾配と考えられているが，前節で示したとおり，結晶の熱力学関係式からは，弾性係数は基準配置 X におけるエネルギーの 2 次導関数と定義される．応力をひずみと関係づけるためには，2 つの基準配置 X と X̄ における応力 σij と σ̄ij の関係を，X と X̄ 間のひずみに対応させる必要がある．
状態 X̄ から状態 X への変形勾配テンソルを Jij とすると，

Jij = Jji =
31

∂Xi
∂ X̄j

(3.17)

である．対応する Lagrange ひずみテンソル nij は，

nij =

1
(Jki Jkj − δij )2(3.18)と表される．基準状態 X̄ から任意の状態 x への Lagrange ひずみテンソルを η̄ij ，基準状態 X から x へのそれを ηij とする．(3.18) と同様に考えると，η̄ij と ηij に関して，

1
1 ∂xk ∂xk
η̄ij + δij =2
2 ∂ X̄i ∂ X̄j
1 ∂xk ∂Xr ∂xk ∂Xs=
2 ∂Xr ∂ X̄i ∂Xs ∂ X̄j()
1= ηrs + δrs Jri Jsj
2(3.19)

が導ける．これを ηrs で微分すると次式が得られる．
(∂ η̄ij∂ηrs)= Jri Jsj

(3.20)

η′

基準状態 X における応力 σij は，

1
σij =
V
(

=

(

V̄
V

∂U
∂ηij

)

1
V̄

)
(

η′

∂U
∂ η̄kl) (η̄ ′∂ η̄kl∂ηij)(3.21)
η′

ここで，V̄ −1 (∂U/∂ η̄kl )η̄′ は，基準状態 X からの値であり，σ̄kl ではない．基準状態 X̄
から X へのひずみが nij であることより，U は次式のように Taylor 展開できる．
(∂ Ū
U = Ū +∂nij)

nij + · · ·(3.22)

n′

右辺の微係数は基準状態 X̄ からのひずみに対するものであることを考慮すると，

1
V̄

(

∂U
∂ η̄kl)= σ̄kl + C̄klmn nmn + · · ·(3.23)

η̄ ′である．式 (3.20)(3.23) を用いて，(3.21) は以下のように変形できる．
(

σij =

V̄
V

)

{

}

Jik Jjl σ¯kl + C̄klmn nmn + · · ·32(3.24)

一方，状態 X における結晶の 3 つのベクトル a, b, c からなるテンソルを h とし，状態 X̄ のそれを h̄ とすると，
[ ]

[

]

[

V̄ = det h̄ = det J −1 h = det J −1
ij V]

(3.25)

Jij−1 は式 (3.18) の逆を考えると ηij について以下のように展開することができる．
Jij−1 = δij − nij + · · ·(3.26)[]V̄= det Jij−1 = 1 − nii + · · ·V(3.27)

これより，また，Jij も Jij−1 と同様に ηij について以下のように展開できる．

Jij = δij + nij + · · ·(3.28)

式 (3.27) と式 (3.28) から，式 (3.24) のひずみ nij についての一次展開は[

σij = σ̄ij + nkl −σ̄ij δkl + σ̄il δjk + σ̄jk δil + C̄ijkl](3.29)となる．上式が基準状態 X̄ からのひずみに対する応力の関係である．ひずみの対称性を考慮し，

1
Bijkl = (σil δjk + σjl δik + σik δjl + σjk δil − 2σij δkl ) + Cijkl
2(3.30)

なる弾性剛性係数 (elastic stiﬀness coeﬃcient) を導入すると，応力とひずみの関係は以下のようになる．

σij (x) = σij (X) + Bijkl ηkl + · · ·(3.31)

これより基準配置 X からのひずみ ηij に対する応力変化は次式で与えられる．
(

∂σij
∂ηij

)

= Bijkl(3.32)

η′式 (3.30) より，無負荷平衡点における Bijkl は Cijkl に一致する．これゆえ，σij =0 における応力–ひずみの勾配は，弾性係数と一致するが，σij ̸=0 の点では Bijkl となる．

33

3.4

弾性剛性係数による格子不安定性評価

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

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

左辺第 1 項がエネルギー変化，第 2 項が応力 σij が保持されたまま変形すると仮定した場合の外部負荷によってなされる仕事であり，結晶の安定性は弾性係数マトリクス

Cijkl の正値性に帰着される．ただし実際には Cijkl による不安定条件は，式 (3.3) のヘッ
シアンによる不安定条件には無負荷平衡点でのみ有効であり，変形が大きくなると差が生じることが示されている (70) ．その理由については議論がなされているところであるが，式 (3.3) と式 (3.33) を比較した時，状態 X におけるエネルギー展開時に体積一定という拘束条件があるのが一因と思われる．そこで，Wang らは，弾性係数の代わりに弾性剛性係数 Bijkl のヘッシアンの正値性で不安定性を評価することを提案している (37)(71) ．弾性剛性係数が応力-ひずみの勾配を表すものであることを考えると，この∂σij< 0，すなわち，変形に対する抵抗力を喪失する点を表すもの
クライテリオンは∂ηijと解釈でき，弾性剛性係数の対称部分

1 T
SY M
BIJ
= (BIJ+ BIJ )
2

(3.34)の正値性が安定性を支配する (71) ．ここで，Cartesian 座標における指標 i〜l は 1 から3 までの値をとるが，対称性を考えると独立な成分は 21 であり，Bijkl は(i, j)=(1,1),(2,2),(3,3),(2,3),(3,1),(1,2) の 6 つの組み合わせを用いて BIJ（I,J=1〜6）のように表記される（Voigt 表示）(69) ．また，T は転置行列を表す．
本研究で対象とする [001] 方向の変形では，系の安定性は次の 6 × 6 行列の正値性により評価される (68) ．

34

B11 B12 B13 000
B12 B22 B23 000
B13 B23 B33 0
0
0
SY M
BIJ
=
0
0
0 B44 00000
0 B55 0
0
0
0
0
0 B66

(3.35)|BIJ | <0 となる条件は，以下の 5 つのいずれかとなる．
B11 B12 B13
B12 B22 B23
B13 B23 B33
B11 B12B12 B22<0<0(3.36)(3.37)

B44 < 0(3.38)

B55 < 0(3.39)

B66 < 0(3.40)

第一式 (3.36) は，膨張による格子の崩壊，即ち体積弾性率が 0 になることを意味している (72) ．本論文では以降これを spinodal 条件と呼び，それによって生じる不安定を

spinodal 不安定と呼ぶ．第二式 (3.37) は，横方向変形のバランスを表しており，横方向変形が等方変形から非等方変形に変形経路分岐が起こることを意味する．これを同様に以降 Born 条件 (Born 不安定) と呼ぶ．さらに，第三式 (3.38)〜第五式 (3.40) はそ
れぞれの変形モードへのせん断不安定が生じることを表している．同様にそれぞれ B44
条件 (B44 不安定)，B55 条件 (B55 不安定)，B66 条件 (B66 不安定) と呼ぶ．

35

第4章
酸化物含有 Fe の安定性評価
これまでの検討では，2 × 2 × 2 の 8 格子の bcc-Fe からなるスーパーセルを基に，Y原子を 2 個，O 原子を 3 個添加した系と，さらに 1 つの Y 原子の代わりに Ti 原子を置換した系において第一原理解析を行い，構造安定性を格子不安定性の観点から議論してきた (7) ．本章では，3 × 3 × 3 の 27 格子に解析セルを拡大し，構造安定性を全自由エネルギー，電子密度および格子不安定性の観点から検討する．Fe 中における Y，O，Ti 原子の安定サイトの組み合わせが現時点では不明なため，添加サイトの違いによる構造安定性の変化，力学特性向上への寄与についても議論する．

4.1

解析条件

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

セルに，図中の大球で示すように Y 原子を 2 個置換，小球で示すように O 原子を 3 個添加し解析を行った．また，比較のためにこれまで行ってきた 8 格子の bcc-Fe を基本とするモデル (Small モデル) を図 4.1(d)∼(g) に示す．2 つの Y 原子は [111] 方向に 1 つの Fe 原子を挟んだ位置に置換し，O 原子については，図 4.2 に示すように bcc 格子の中心原子と頂点原子からなる八面体の中心が最も安定なサイトであることが報告されているので (34) ，Y-Y の周囲で 3 通りの添加サイトを考慮してきた（O 原子が存在する面を薄く着色している．
）. ただし，2 × 2 × 2 のスーパーセルでは Y 原子が連続する構造であるため今回は 27 格子に拡大し，Y 原子が連続しない条件で解析を行った．
図 4.1(e),(f) の S-octa1 と S-octa2 モデルは 8 格子モデルで Y,O を添加したときに最も自由エネルギーが低くなったモデルである．Small モデル (2 × 2 × 2 モデル) は 16 個の Fe 原子，Large モデル (3 × 3 × 3 モデル) は 54 個の Fe 原子で構成されるスーパー
セルに Y2 O3 をそれぞれ添加しており，スーパーセル内における Y2 O3 の含有割合はSmall モデル >Large モデルとなっている．さらに図 4.3(a)∼(e) に示すように octa1,2
ならびに S-octa1,2 の一つの Y 原子を Ti 原子に置換した系についても解析した（それぞれ Tiocta1A,1B,2 ならびに S-Tiocta1A,1B,2 と称する）．octa2 と S-octa2 については対称性があるため，図中の大球の上下どちらを Ti に置換しても等価なので 1 つのモデル（Tiocta2,S-Tiocta2）としている．逆格子空間の k 点メッシュは Monkhorst-Pack
法 (77) に基づき Fe16 モデルでは 113 , その以外モデルで 43 とした．その他の計算条件
を表 4.1 にまとめて示す．セル内での原子核の移動は考えず，図 4.1,4.3 の位置に固定した状態で基準となる bcc 格子の寸法を等方的に拡大，縮小して電子状態計算を行い，bcc 格子長 3a に対する全自由エネルギーの変化を調べた．全自由エネルギーが最小となる格子長さを系の基準格子長さ 3a0 とし，ひずみ摂動を与えて系の断性剛性係数を求めて格子不安定解析を行った．弾性剛性係数 Bij は，図 4.4 に模式的に示すようにセルに ±0.05 の引張/圧縮ひずみ，もしくは ±0.025 のせん断ひずみを与え，ひずみ摂動に対する応力変化から数値的に求めた．

37

3a
a
3a

a

z
y
x

a

3a
(a) octa1

Y(b) octa2(c) octa3O(d) Fe16(e) S-octa1(f) S-octa2(g) S-octa3

Fig.4.1 Simulation models for Y2 O3 in Fe.

Short axis

Fig.4.2 Schematic of oxygen doping site.

38

3a
a
3a

a

a

z
y
x

3a
(a) Tiocta1A(b) Tiocta1B(b) Tiocta2(c) S-Tiocta1A(d) S-Tiocta1B(e) S-Tiocta2

Y
TO

Fig.4.3 Simulation models for YTiO3 in Fe.

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

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

Fe14
16
96
296
113

Fe14 Y2 O3
19
106
495
43
Fe16 YTiO3
19
107
495
43

39

Fe52 YTiO3
57
335
495
43

Fe52 Y2 O3
57
334
495
43

∆η1= r5.0010 -4
∆σ1
∆η1
∆σ2
B12 =
∆η1
∆σ3
B13 =
∆η1

∆η4= r2.5010-4∆η5= r2.5010

-4

B11 =

B44 =

∆η2= r5.0010

∆σ4
∆η4

-4

∆σ1
∆η2
∆σ2
B22 =
∆η2
∆σ3
B23 =
∆η2

B21 =

B55 =

∆σ5
∆η5

∆η6= r2.5010 -4∆η3 = r5.0010 -4
∆σ1
∆η3
∆σ2
B32 =
∆η3
∆σ3
B33 =
∆η3

B31 =

B66 =

∆σ6
∆η6

z
y
x

Fig.4.4 Strain perturbation on supercell for evaluation of elastic stiﬀness coeﬃcients.

40

4.2

解析結果と考察4.2.1

全自由エネルギーと弾性剛性係数行列式の比較

表 4.2 に各モデルで最小のエネルギー E0 と，その時の格子長さ a0 ，ひずみ摂動により得られた弾性剛性係数の行列式 detBij の値を示す．また全自由エネルギーは系を構成す
る原子数 N によって変わるので一原子当りに換算した値 E0 /N を示した．8 格子を基本とした，Small モデルの結果も合わせて示している．a0 の大きさを比較すると，Fe 単元
系に Y2 O3 を添加することにより格子定数が増加しており，Small モデル (2 × 2 × 2 モデル) は Large モデル (3×3×3 モデル) に比べ，Y2 O3 の含有割合が高いため全体的に a0 が大きくなっている．ひずみ摂動により得られた detBij はいずれのモデルも正値を示し，今回解析したモデルは力学的には存在し得る構造であることが確認できる．Large モデ
ルと Fe16 の detBij の大きさを比較すると Fe16>octa2>octa1>octa3 となり，Y2 O3 を添加することで無負荷状態での弾性剛性係数の行列式の値は低下した．また Y2 O3 を同じ組み合わせのサイトに添加したモデルの Small モデルと Large モデルを比較するとLarge モデル >Small モデルとなりスーパーセル内の Y2 O3 の含有割合が高いほど低い
弾性剛性係数の行列式の値を示し, Y2 O3 を Fe に添加することが原子レベルで直接的に剛性の向上に寄与しないことがわかる．原子一個当たりの全自由エネルギーを比較すると S-octa3>S-octa1>S-octa2>octa3>octa1>octa2>Fe16 となり Y2 O3 の割合が高い

Small モデルの方が全体的に全自由エネルギーが高くなる．図 4.3 で示したように Y2 O3を添加したいずれのモデルも，添加する 3 個の O 原子は bcc 格子の頂点と中心原子から構成される八面体の中心サイトに添加されている．Fe のみで構成される八面体の中心に 3 個の O 原子が添加されている S-octa2 と octa2 は，それぞれ Small モデルと Large
モデルの中で最小の全自由エネルギーを示し，有利な構造であった．表 4.3 に Ti を添加したモデルの結果を示す．各モデルの格子定数 a0 は Ti を置換する前のモデルに比べ小さくなっている．Ti を置換したいずれのモデルも Y2 O3 を添加したモデルと同様に，

detBij が正値を示し，力学的に存在し得ることが確認できた．Large モデルと Fe16 の
detBij の大きさを比較すると Fe16>Tiocta1A>Tiocta1B>Tiocta2 となり，やはり酸化41

物の添加によって弾性剛性係数の行列式の値は低下している．各モデルの Small モデルと Large モデルの値を比較すると Y2 O3 モデルと同様に Large モデル >Small モデルとなりスーパーセル内の YTiO3 の含有割合が高いほど低い弾性剛性係数を示した．先のY2 O3 の場合と同様，Y,Ti の複合酸化物は直接的に Fe の剛性向上には寄与しない．octa1の構造で Y を Ti に置換すると detBij の値は増加したが (Tiocta1A,1B>octa1)，octa2では減少している (Tiocta2<octa2)．Ti を置換する前の detBij の値は octa2>octa1 であるが，Ti 添加系では Tiocta1A>Tiocta2>Tiocta1B となり Ti 近傍に O 原子を置換したモデルが最大の detBij の値を示し，O 原子は Y に比べ Ti と強く結び付き，剛性も増加することが推測される．各モデルの原子一個当りの全自由エネルギーを比較すると S-Tiocta1B>S-Tiocta2>S-Tiocta1A>Tiocta1B>Tiocta1A>Tiocta2>Fe16 となる．

Tiocta1A は Ti 原子近傍に，Tiocta1B は Y 原子近傍に O 原子が存在しており O 原子が Y 近傍に置換されるより Ti 近傍に置換される方が高い剛性を示し，エネルギー的にも安定な構造になることが示された．これは O 原子は Fe 中で Y よりも Ti と結びつ
きやすく，Y,Ti と複合酸化物を形成することで高い剛性で安定な酸化物になることが示唆している．また，YTiO3 添加モデルの中では八面体サイトに Y,Ti が存在しない

Tiocta2 が最小の値を示している．これは酸化物添加モデルの中でも最小であり，今回解析したモデルの中で Tiocta2 が最もエネルギー的に有利な構造である．

42

Table 4.2 Determinant of elastic coeﬃcients and total free energy(Y2 O3 models).

Fe16
octa1
octa2
octa3
S-octa1
S-octa2
S-octa3a0 (Å)
2.856
3.020
2.999
3.040
3.306
3.344
3.242

detBij (GPa6 )1.27×1012
5.03×1011
9.54×1011
8.78×1010
5.64×109
5.20×109
2.30×1010

E0 (eV)-132.9-441.5-446.3-436.8-136.8-138.9-134.9

E0 /N (eV)
-8.31
-7.75
-7.83
-7.66
-7.20
-7.31
-7.10

Table 4.3 Determinant of elastic coeﬃcients and total free energy(YTiO3 models).

Fe16
Tiocta1A
Tiocta1B
Tiocta2
S-Tiocta1A
S-Tiocta1B
S-Tiocta2

a0 (Å)
2.856
2.990
3.005
2.986
3.232
3.271
3.221

detBij (GPa6 )1.27×1012
1.07×1012
9.77×1011
4.91×1011
2.15×109
1.67×109
5.05×101043

E0 (eV)-132.9-449.4-444.7-449.6-141.0-138.1-140.3

E0 /N (eV)
-8.31
-7.88
-7.80
-7.89
-7.42
-7.27
-7.38

4.2.2弾性剛性係数行列の各成分の比較

図 4.5 に Y2 O3 添加モデルの弾性剛性係数行列の各成分 Bij の値を示す．ただし，2×2×2 モデルの Fe 単元系 (Fe16 モデル) で評価した B11 ∼ B66 の値で無次元化している．Small モデルの各成分の値を比較すると，S-octa1,3 モデルの O 原子と Fe 原子が隣接する [001] 方向の弾性剛性係数 B33 が特に大きな値を示しているが，その他の弾性剛性係数成分の値は Fe16 より著しく低下している．S-octa2 モデルは, 垂直成分である B11 , B22 , B33 は比較的 Fe16 に近い値を示すが，S-octa1,3 モデルと同様にその他の成分は低い値を示す．Large モデルの値を比較すると，Small モデルと同様に octa1,3の B33 成分は比較的高い値を示している．またその他の弾性剛性係数成分が Fe 単元系(Fe16) の値に近づき，他方向への影響が減少している．octa2 モデルも S-octa2 に比べ，すべての成分が Fe16 モデルの値に近づき，Y2 O3 を添加したモデルのなかで最も弾性剛性係数各成分の異方性が小さくなっており，弾性剛性係数の各成分は Fe16 モデルの同じ方向の成分の 92% ∼ 102% の範囲に収まった．図 4.6 に Ti 添加モデルの弾性剛性
係数行列各成分を示す．こちらも 2×2×2 モデルの Fe 単元系 (Fe16 モデル) で評価した

B11 ∼ B66 の値で無次元化している．Y2 O3 モデルと同様に，Large モデルを Small モ
デルと比較すると異方性が減少し，octa2 モデルと同様に八面体サイトに Y,Ti が存在しない Tiocta2 が Ti 添加モデルの中で最も異方向性が小さい．図 4.7 に Large モデルの弾性剛性係数行列の各成分 Bij の値を，Ti 添加前後の値を並べて再掲する．Ti 添加
前と比べて Tiocta1A は 9 成分中 8 成分で，Tiocta1B モデルは 9 成分中全成分で Bij の値が増加したのに対し，Tiocta2 モデルは 9 成分中 7 成分の Bij の値が減少した．

44

2.5

Fe
2

Normarized elastic stiffness, Bij/Bij

Normarized elastic stiffness, Bij/Bij

Fe

2.5

S-octa1
S-octa2
S-octa3

1.5

1

0.5

0
B11

B22

B33

B12

B13

B23

B44

B55

octa1
octa2
octa3

2

1.5

1

0.5

0

B66

B11

B22

(a) 2×2×2model (Small models)

B33

B12

B13

B23

B44

B55

B66

(b) 3×3×3model (Large models)Fig.4.5 Comparison of each component of Bij (Y2 O3 models)．
2.5

Fe
2

Normarized elastic stiffness, Bij/Bij

Normarized elastic stiffness, Bij/Bij

Fe

2.5

S-Tiocta1A
S-Tiocta1B

1.5

S-Tiocta2

1

0.5

Tiocta1A

2

Tiocta1B
Tiocta2

1.5

1

0.5

0

0
B11

B22

B33

B12

B13

B23

B44

B55

B66

B11

B22

(a) 2×2×2model (Small models)

B33

B12

B13

B23

B44

B55

B66

(b) 3×3×3model (Large models)Fig.4.6 Comparison of each component of Bij (YTiO3 models)．
2.5

Fe

octa1
Tiocta1A
Tiocta1B

2

Normarized elastic stiffness, Bij/Bij

Normarized elastic stiffness, Bij/Bij

Fe

2.5

1.5

1

0.5

0

octa12
Tiocta2

2

1.5

1

0.5

0
B11

B22

B33

B12

B13

B23

B44

B55

B66

(a) octa1 and Tiocta1A,1B models.

B11

B22

B33

B12

B13

B23

B44

(b) octa2 and Tiocta2 models.

Fig.4.7 Comparison of each component of Bij (Large models)．

45

B55

B66

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

図 4.8，4.9 に YTiO3 添加モデル（Large モデル）の中で比較的エネルギーの低かった Tiocta1A および Tiocta2 の価電子密度分布を示す．ここで等高線は 0.07(eV/Å 3 )ごとに表示し，図の上部は表示している面を模式的に示している．価電子密度が著しく低い個所が見られるが，この位置に Y または Ti が存在している．これは価電子数の違いによるもので，Y および Ti の価電子が Fe に比べ少ないことに起因する．Y の価電子が少ないことに起因する結合の低下や，酸化物を添加することでスーパーセル
が大きくなり，Fe 原子間の結合が低下することにより，前節で示したように Y2 O3 やYTiO3 を添加したモデルの弾性剛性係数行列式の値が低い値を示したと考えられる．
特に Small モデルでは，このような結合の低下が周期境界条件下で顕著になることで，無負荷状態の弾性剛性係数行列式の値が Large モデル >Small モデルとなったと考えられる．図 4.8(a)(b)，図 4.9(a)(c) では垂直方向の Fe-O-Fe という非常に価電子密度の高い原子の並びが存在しており，さらに円 (A) で示すように垂直方向の Fe-Fe 間の価電子密度も増加している．octa1,3 モデルの弾性剛性係数行列の B33 以外の成分が，Fe単元系モデルに比べ小さな値を示したのにもかかわらず，B33 のみ特に大きな値を示したのはこの Fe-O-Fe 及び Fe-Fe の価電子密度の高い並びが [001] 方向に 3 箇所存在し変形抵抗を高めているためであると考えられる．一方，垂直方向の弾性剛性係数成分

B11 , B22 , B33 が比較的高い値を示した Tiocta2 は，Fe-O-Fe の結合が [100][010][001] 方
向に 1 つずつ存在するため異方性が小さくなったことが分かる．また Tiocta1A の (1̄10)面（図 4.8(c)）に注目すると，Fe 原子間に円 (B) のような価電子密度の谷が存在する
が，O 原子の存在により円 (C) で示す Fe 原子間には価電子密度の谷は見られない．O原子はこのように Fe 原子間の価電子密度を補っており，周囲の Fe 原子と強く結合していると考えられる．一方，Tiocta2 の O 原子が添加された面 (図 4.9(b)) では，O 原子近傍に価電子が集中するために，周囲の Fe 原子は O 原子や他の O 原子近傍の Fe 原子と強く結合し（円 (D)），O 原子から離れた Fe 原子との間で破線で示すような価電子密度の谷が生じている．

46

z
y
x

Fe

Fe

A

B

A

Ti

Fe

Fe

O

O

Fe

Fe

B

7.4

C
B
C

Y

0 (eV/Å3)
0

(a)

(b)

(c)

Fig.4.8 Distribution of valence electron density(Tiocta1A)．

z
y
x

D Fe

Fe

Fe

A
7.4

Ti

D

O

D

A

O
Fe

Fe

O

Fe

(a)

Y
0

Fe

(b)

(c)

Fig.4.9 Distribution of valence electron density(Tiocta2)．

47(eV/Å3)

4.2.4O 原子溶解熱の検討

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

△H = Etotal − Eother − Eoxygen(4.1)

それぞれの O 原子添加サイトについて算出した Y2 O3 添加モデルの溶解熱 △H を模式図と共に図 4.10 に示す．Small モデルの結果も合わせて示している．Large モデルの溶解熱の値は Small モデルに比べ，全体的に小さくなっている．表 4.2 で示したように，Large モデルのユニットセルの辺長 a0 は，Small モデルに比べ小さくなっており，

Fe 原子間が狭い．このため Large モデルは Small モデルに比べ Fe 原子間に O 原子が侵入しにくくなり，O 原子の溶解熱が全体的に小さくなったと考えられる．Large モデ
ルの溶解熱を比較すると，八面体サイトに Y 原子が存在する octa1(α) と octa3(α), (β)
が正の値を示し，熱力学的にはこのサイトに存在しえないことを示している．八面体
がすべて Fe 原子で構成されるサイトの値はすべて負の値を示し，O 原子が Fe 原子のみで構成される八面体中心サイトに侵入しやすいことを示唆している．octa2(α)(β)(γ)の溶解熱はいずれも負で同程度の値を示している．octa1(β)(γ) や octa3(γ) のサイト
はそれより大きい値をとっているが，これは octa1 ならびに octa3 は原子サイズの大きな Y 原子の近傍に O 原子が侵入しているため，octa2 に比べ格子定数 a0 が大きくなっており (表 4.2)，Fe 八面体の中心サイトが拡げられたためと考えられる．
図 4.11 に TiYO3 添加モデル（Large モデル）の溶解熱を模式図と共に示す．Ti 置換
前のモデルの結果も合わせて示している．Ti 置換前のモデルの溶解熱と比べると，Y原子を Ti 原子に置換したことですべてのモデルで格子の大きさ a0 が小さくなっている
ので O 原子が侵入しにくくなり，Tiocta1A(α) を除き全体的に溶解熱が小さくなっている．Ti 添加前と同様に，八面体サイトに Y 原子が存在する Tiocta1B(α) は正の値を示し，O 原子が侵入しえないことが示唆される．一方，octa1(α) は正の値を示していたが，Y 原子を Ti 原子に置換することにより溶解熱が負の値を示し（Tiocta1A(α)），

Tiocta1A は存在可能な添加サイトの組み合わせであることが示された．Tiocta1A と48

Tiocta1B の全自由エネルギーを比較しても Tiocta1B>Tiocta1A となっており，O 原
子は Y より Ti 原子近傍に侵入しやすく安定な配置であることが示唆される．先述した
ように O 原子近傍の Y 原子を Ti 原子に置換することによる detBij の値の増加や，異方性の減少も見られ，Ti 原子を添加することで O 原子は Y 原子近傍で様々なサイトの組み合わせをとることが可能になることが推測される．しかし，Y 原子や Ti 原子が八面体の頂点となる octa1(α) は八面体サイトが Fe のみで構成される (β), (γ) に比べ，溶解熱が小さく，O 原子は Fe のみで構成される八面体サイトの中心に存在する方が安定であることが示唆される．octa2 の構造では (α), (β), (γ) いずれも八面体を Fe 原子が構成しており，Ti 添加によって格子長さが小さくなるため値は Small モデルに比べ減少しているものの，いずれも負の値をとりかつサイトによる差がほとんどない．先の
全自由エネルギーの大きさ，ならびに Bij の等方性や O 原子溶解熱がほぼ等価である
ことからも octa2 および Tiocta2 が最も存在する可能性が高いと結論づけられる．

49

3.604(eV)-4.690(eV)-7.016(eV)-7.016(eV)-1.167(eV)-1.167(eV)S-octa1octa1-0.684(eV)-3.908(eV)-0.683(eV)-4.287(eV)-4.085(eV)-0.684(eV)S-octa2octa2

3.069(eV)-2.613(eV)-1.529(eV)-4.995(eV)-2.613(eV)S-octa3

3.069(eV)octa3

Fig.4.10 Comparison of solution heat(eV)(Y2 O3 models).

50

-0.289(eV)-0.627(eV)-0.627(eV)
3.604(eV)-1.167(eV)Tiocta1A-1.167(eV)3.952(eV)-0.911(eV)octa1-0.911(eV)Tiocta1B-0.684(eV)-0.391(eV)-0.683(eV)-0.684(eV)-0.393(eV)-0.393(eV)octa2Tiocta2

Fig.4.11 Comparison of solution heat(eV)(YTiO3 models).

51

4.3

結言

酸化物を Fe 中に固溶することによる特性向上を，添加サイトの違いを考慮し検討するために，bcc-Fe 単位セルを 3 × 3 × 3 並べたスーパーセルに Y 原子を 2 個，O 原子を3 個添加した系，ならびに Y 原子 1 個の代わりに Ti 原子を 1 個置換した系について第一原理計算による全エネルギーの評価，弾性剛性係数の正値性による格子安定性，O原子の溶解熱などを検討した．得られた結果を以下に示す．

(1) 先行研究の 2 × 2 × 2 モデルで安定であることが示されていた添加サイトの組み合わせで解析したところ，弾性剛性係数の行列式（detBij ）の値は 2 × 2 × 2 のそれに比べて低下するものの，いずれも正の値をとり「力学的に」存在し得る構造であることが示された．

(2) 全自由エネルギーの大きさを比較すると，Y2 O3 , YTiO3 のいずれを添加してもFe 単元系よりエネルギーが大きくなるが，3 個の O 原子を Fe のみで構成される八面体の中心サイトに添加した構造が Y2 O3 , YTiO3 を添加した系では低い全自由エネルギーを示し，エネルギー的には有利な構造であることが示された．

(3) 同じ O 原子の添加構造で比較すると，2 × 2 × 2 モデルと 3 × 3 × 3 モデルの detBij
は Fe 単元系 > 3 × 3 × 3 モデル > 2 × 2 × 2 モデルとなり，酸化物を添加することが直接 bcc-Fe の剛性向上に寄与しないことが示された．

(4) 弾性剛性係数 Bij の各成分を比較したところ，Ti や Y が八面体頂点を構成するサイトに O 原子を添加した構造では異方性が極めて強くなり，酸化物の長手方向に対応する [001] 方向の弾性剛性係数 B33 が Fe 単元系より著しく高く，他は著
しく低くなった．八面体が Fe のみで構成されるサイトに O 原子を添加したモデルは Y2 O3 , YTiO3 いずれにおいても異方性が小さくなり，Y2 O3 を添加したモデルの各成分の値は Fe 単元系の 92% ∼ 102% の範囲となった．

(5) YTiO3 を添加した系において価電子密度分布を調べた結果，O 原子は八面体サ

52

イトを構成する隣接 Fe 原子と強く結合し，さらに次の隣接する Fe 原子との間にて価電子密度を増加させることが分かった．

(6) 強い異方性を示した構造は，3 つの八面体の短軸方向が [001] であり，かつ八面体の辺が共有されるように隣接した構造のため，(5) の O 原子による結合が強められ強い異方性を生じたことが明らかとなった．

(7) 異方性が最も小さくなった構造は，O 原子を含む八面体の短軸方向が [100],[010],[001]となっており，各隣接八面体は 1 つの Fe 頂点でのみで接する構造である．このため O 原子による結合が等方的になったものと考えられる．

(8) (7) の構造の O 原子の溶解熱はすべて同程度の負の値を示したが，(6) のように隣接する八面体が同じ方向に有する構造で Y 原子を頂点とする八面体中心の O
原子溶解熱は正の値を示しそのサイトに O 原子は侵入しえないことが示された．

53

第5章

Y2O3 含有 Fe の引張シミュレーション
本章では，第 4 章で最もエネルギー的に安定な構造であった octa2 ならびに S-octa2モデルの引張シミュレーションを行い，弾性剛性係数，電子状態等について検討する．
また引張過程の各ひずみにおいて第 4 章と同様に弾性剛性係数の正値性に基づく格子不安定解析を行い，Y2 O3 が理想強度に与える影響について検討する．

5.1

解析条件

前章と同様に，Kresse らにより開発された平面波基底ウルトラソフト擬ポテンシャ
ル法 (59) に基づく第一原理バンド計算コード VASP(73) (Vienna Ab-initio Simulation

Package) を用いて行った．k 点サンプリング，交換相関表現ならびに電子状態の収束計算は前章と同様である．解析には，16 原子で構成される bcc-Fe モデル (図 5.1(a)) な
らびに，第 4 章の Y2 O3 添加モデルにおいて 3 × 3 × 3 モデルと 2 × 2 × 2 モデルで最もエネルギーの低い値を示した，octa2 と S-octa2 を対象とした（図 5.1(b)(c))．原子核固定の条件のもと，まず [001] 方向に引張ひずみ増分 △ε33 = 0.002 ∼ 0.1 を与え，電子状態計算を行った．このひずみ増加/電子状態計算を繰り返し，静力学的に応力-ひずみ曲線を算出した．横方向は拘束 (ε11 = ε22 = 0) とした．また各ひずみ下において，前章で述べたように微小ひずみ摂動を与え，応力変化から弾性剛性係数を数値的に算出し，系の安定性を評価した．octa2 モデルのみ，ひずみ △ε33 = 0.05 ごとにセル内での原子構造緩和計算を行い，octa2 モデルで構造緩和を考慮した状態 (octa2(relaxed))についても検討した．

54

Y
O

z
y
x
(a) Fe16(b) octa2(c) S-octa2

Fig.5.1 Simulation models.

5.2

解析結果と考察5.2.1

応力-ひずみ関係と弾性剛性係数の行列式の変化[001] 引張変形で得られた各モデルの応力-ひずみ関係を図 5.2 に示す．図中の矢印は応力ピークのひずみを示している．octa2 は無負荷状態で緩和すると全方向に 5.1(GPa)の引張ひずみの状態となったため，圧縮変形を与え構造緩和計算を行い，応力が 0 となるひずみを算出し，octa2(relaxed) の引張の基準とした．引張初期の勾配はほぼ同じに見えるが，ひずみ ε33 = 0 でのひずみ摂動から得られる引張方向の弾性剛性係数 B33 はFe16 (232GPa)>octa2 (213GPa)>octa2(relaxed) (164GPa)>S-octa2 (151GPa) でありわずかに差がある．Fe16 モデルの応力はひずみ ε33 = 0.02 まではほぼ線形に上昇したが，ε33 = 0.02 以降で鈍化し非線形性が強く現れ，ひずみ ε33 = 0.03 付近で変曲点をとる．その後再び応力は上昇し，ε33 = 0.4 でピーク応力 σ33 = 24.1(GPa) を示した．一方，

Y2 O3 を添加したモデルでは Fe16 のような引張初期の変曲点は存在せず，応力はひずみ
ε33 = 0.05 までほぼ線形に上昇する．応力-ひずみ曲線のピーク応力 σ33 の大きさは Fe16
(24.1GPa)>octa2 (22.8GPa)>S-octa2 (18.8GPa)>octa2(relaxed) (17.8GPa)，ピーク
ひずみ ε33 は Fe16 (ε33 = 0.4)>octa2,S-octa2 (ε33 = 0.3)>octa2(relaxed) (ε33 = 0.18)
であり，Fe の初期の変曲点を無視すればいずれも Fe 単元系に Y2 O3 を添加することで低55

下している．図 5.3 に弾性剛性係数の行列式 detBij の値の変化を示す．図の左は detBijの変化をそのまま示したもの，右は引張前の各系における detBij の値で割り，無次元化したものである．detBij の絶対値では Y2 O3 を 2×2×2bcc-Fe に添加した系 (S-octa2，赤線) のオーダーが他と異なるため，左の図では常に 0 に見えるが初期値で無次元化す
ると 1 から減少している．矢印は初めて detBij が負になるひずみを示している．格子
安定性の観点からは，Fe16 はやはり最初の変曲点が弾性限界であり，Y2 O3 を添加した他の系に比べてはるかに小さいひずみである．格子不安定ひずみは octa2(relaxed)(ε33 = 0.16)>octa2 (ε33 = 0.13)>S-octa2 (ε33 = 0.11)>Fe16 (ε33 = 0.03) となり，
Y2 O3 を添加したすべてのモデルで弾性限界が著しく増加した．図 5.4 に各モデルの応力-ひずみ曲線を再掲し，弾性限界を矢印で示した．格子不安定解析による弾性限界での応力（理想引張強度）は octa2(relaxed) (16.8GPa)>octa2 (16.5GPa)>S-octa2(13.9GPa)>Fe16 (3.8GPa) となった．
30

25

,GPa

20

Stress

15

Fe16
10

octa2
S-octa25octa2(relaxed)

0
0

0.1

0.2
0.3
Tensile strain

0.4

Fig.5.2 Stress-strain curves.

56

0.5

12

(x10 )
1.2

Magnitude of 6x6 determinant (normalized)

1.2

Fe16
1

octa2
S-octa2

0.8

octa2(relaxed)
0.6
0.4
0.2
0

1
0.8
0.6
0.4
0.2
0
-0.2

-0.2
0.05

0.1

0.15

0.2

0

0.05

Fig.5.3 Changes in the determinants of Bij .

30

25

20

15

Fe16
10

octa2
S-octa25octa2(relaxed)

0
0

0.1

Tensile strain

Tensile strain

,GPa

0

StressMagnitude of 6x6 determinant, GPa 6

1.4

0.1

0.2
0.3
Tensile strain

0.4

0.5

Fig.5.4 Stability limits on stress-strain curves.

57

0.15

0.2

5.2.2

引張変形下の価電子密度分布の変化

図 5.5，5.6 に octa2，S-octa2 の引張変形下の価電子密度分布を示す．図の上部には表示している面を模式的に示しており，価電子密度は 0.072(eV/Å 3 ) ごとに等高線で示している．octa2，S-octa2 のひずみ ε33 = 0 の価電子密度分布では Y,O 原子が添加されることで [001] 方向の Fe-O-Fe の間で高い価電子密度となっており，また円 (A) で示すよ
うに Y,O 周囲で部分的に価電子が高くなることで Fe-Fe 間の価電子が補われている．引張変形が進行しても，Fe-O-Fe の間や円 (A) 部分の価電子密度は高い値を示し，[001] 方
向で強く結合していると推測される．octa2 および S-octa2 は (100),(010),(001) 面にも同様に O 原子を添加したため，この Fe-O-Fe 結合は [100],[010] 方向にも存在している．
つまり，octa2 および S-octa2 の弾性限界が Fe 単元系（Fe16）に比べ高い値を示した理由として，この Fe-O-Fe の結合が [100],[010],[001] 方向に存在し，3 軸引張に近い ε′ = 0の引張で結合を維持したためと考えられる．また図 5.5 の octa2 の ε33 = 0.13 の価電子密度分布に注目すると矢印 (B) で示す部分に価電子密度の谷が生じ，ε33 = 0.18 ではさらに拡大しているが，これにより Y や O 原子近傍に価電子が集中し，O 原子から離れた

Fe-Fe 原子間から価電子密度が低下していくことが示唆される．図 5.7 に octa2(relaxed)の価電子密度変化を示す．ひずみ ε33 = 0 の価電子密度分布では矢印で示す Fe-Fe-O-Fe
原子によって高価電子密度の帯が [001] 方向に存在する．この Fe-Fe-O-Fe の帯は octa2では見られなかったが，これは緩和により (001) 面に侵入した O 原子が [001] 方向に隣接する Fe 原子間を広げ，広げられた Fe 原子がさらに [001] 方向に隣接する Fe 原子と接近することで Fe-Fe 間の価電子密度が上昇したためである．octa2(relaxed) が octa2に比べ高い弾性限界を示し，今回解析したモデルの中で最大となったのは，このように引張方向に配向した Fe-Fe-O-Fe 結合に起因すると考えられる．O 原子近傍の Fe(a)と Fe(c) に注目すると 引張変形下で O 原子に相対的に接近することで，Fe(a) と Fe(b)間で高価電子密度の帯にくびれが発生し，不安定後のひずみ ε33 = 0.18 では破線で示
すように Fe(a),(c) と Fe(b),(d) の間で価電子密度の谷が生まれており，さらに円 (C) に示す O-Fe(c) 間に価電子が局在している．

58

Fe
(A)
Y

z

7.4

Fe

y

O
Fe

x
O

0

Y

(B)

(B)

(A)

0.03

0
(eV/Å3)

0.08

0.13

Fig.5.5 Changes in the distribution of valence electron density(octa2)．

59

0.18

(A)
7.4

Fe
z

O
y
Fe
x

(A)
O

0

Y

0
(eV/Å3)

(A)

(A)
0.03

0.08

0.13

Fig.5.6 Changes in the distribution of valence electron density(S-octa2)．

60

0.18

Fe(b)

Fe(d)
7.4

Fe(a)

z
y

Fe(c)

O
Fe

x
O

0
(eV/Å3)

0

Y

(C)

0.03

0.08

0.13

0.18

Fig.5.7 Changes in the distribution of valence electron density(octa2(relaxed))．

61

5.2.3

引張変形下の原子間結合距離の変化

図 5.8 に引張変形時の Fe,Y,O の結合距離変化を示す．図中の破線は各モデルの弾性限界ひずみである．いずれの距離もひずみ ε33 = 0 の値で割り無次元化している．
図 (f) に模式的に示すように，3 つの O-O 結合，6 つの Y-O 結合，2 つの Y-Fe 結合の平均を算出するとともに，図中 O(α) の引張方向 ([001] 方向) に存在する Fe 原子について, O-Fe ならびに Fe-Fe 結合の距離変化を算出した．octa2(relaxed) の Y-Fe 間の距離（図 5.8(a)）は octa2 とほぼ同じ割合で離れていくが，O-O 間の距離（図 5.8(b)）はひずみ ε33 = 0.08 までセルが変形しても変化せず，さらにそれ以降では O-O 結合は

ε33 = 0 の時より収縮している．この O-O 結合の収縮によって Y-O 結合 (図 5.8(c)) も
ひずみ ε33 = 0.08 以降に減少している．すなわち引張による Y-Fe 結合はアフィン変形に近いが，原子構造の緩和により Y-Fe 軸を囲む O-O 結合が収縮して Fe 付近に集まり，

Y-Fe 間の結合の低下を補っていることが示唆される．図 5.8(d) の Fe-O の距離をみると，octa2(relaxed) では octa2 より増加割合が小さく，引張方向の Fe-O 結合が変形しにくく，強く結合していることがわかる．逆に，引張方向の Fe-Fe 結合は octa2 に比べ大きく延伸されている．すなわち，Fe 原子が O 原子に引き寄せられたため，隣接する

Fe-Fe 間がより大きく引きのばされる．

62

1.25

1.25

Normalized distance, d/d0

Normalized distance, d/d0

1.2

octa2
1.15

octa2(relaxed)

1.1
1.05
1
0.95

1.2
1.15
1.1
1.05
1
0.95

0.9

0.9
0

0.05

0.1
Tensile strain

0.15

0.2

0

0.05

0.15

0.2

(b) Distance between O-O.

1.25

1.25

1.2

1.2

Normalized distance, d/d0

Normalized distance, d/d0(a) Distance between Y-Fe.

0.1
Tensile strain

1.15
1.1
1.05
1
0.95

1.15
1.1
1.05
1
0.95

0.9

0.9
0

0.05

0.1
Tensile strain

0.15

0.2

0

0.05

0.1
Tensile strain

0.15

0.2

(d) Distance between Fe-O.

(c) Distance between Y-O.
1.25

Normalized distance, d/d0

1.2
1.15
1.1

Y

1.05

O(

O

1

z

0.95

O
y

0.05

0.1
Tensile strain

0.15

0.2

Fe

Y

0.9
0

)

x

(e) Distance between Fe-Fe.

(f) Schematic of bonding.

Fig.5.8 Changes in bond length under tension(octa2 and octa2(relaxed))．

63

5.3

結言

エネルギー，弾性特性，溶解熱いずれの点でも有利であった Y2 O3 添加モデル (octa2,

S-octa2) について，第一原理計算による引張シミュレーション（横ひずみ ε′ = 0）を行った．比較のため Fe 単元系での引張も 2 × 2 × 2 の bcc セルで行った．いずれも原子移動を考えない静力学解析であるが，3 × 3 × 3bcc-Fe に Y2 O3 を添加した octa2 モ
デルについては原子核構造緩和を考慮したシミュレーションも実施した．得られた結果を以下に示す．

(1) Fe 単元系の応力-ひずみ曲線には引張初期に変曲点が現れた. 応力-ひずみ曲線のピーク応力 σ33 の値は Fe 単元系 (24.1GPa)>octa2 (22.8GPa)>S-octa2 (18.8GPa)>octa2(緩和) (17.8GPa) となり，ピークひずみ ε33 の値は Fe 単元系 (ε33 = 0.4)>
octa2,S-octa2 (ε33 = 0.3)>octa2(緩和) (ε33 = 0.18) となった．
(2) 各ひずみ下での弾性剛性係数を評価し，detBij が初めて負になる弾性限界 ε33 の値は，octa2(緩和) (ε33 = 0.16)>octa2 (ε33 = 0.13)>S-octa2 (ε33 = 0.11)>Fe
単元系 (ε33 = 0.03) となった．またその時の応力 (理想引張強度)σ33 は octa2(緩和) (16.8GPa)>octa2 (16.5GPa)>S-octa2 (13.9GPa)>Fe 単元系 (3.8GPa) となり，格子安定性の観点からは bcc-Fe 中に Y2 O3 を添加することにより弾性限界が増加し，理想引張強度も増加した．

(3) 引張変形下の価電子密度分布を検討した結果，O 原子を含む八面体の短軸方向の Fe-O-Fe 結合の価電子密度が高いことがわかった．octa2 構造ではこの結合が[100],[010],[001] のいずれの方向にも存在しており，3 軸引張に近い ε′ = 0 の引張で弾性限界を向上させたものと考えられる．

(4) 構造緩和した octa2 と静力学的にアフィン変形させた octa2 それぞれについて，
dz2 原子軌道のような形状の Y-(FeO3 )-Y 構造の O-O,Y-O,Y-Fe の結合距離を調べた結果，中心軸の Y-Fe-Y はほぼアフィン変形しているが，中心 Fe を囲む Oの三角形リングは収縮することがわかった．

64

(5) (001) 面上の O 原子とその上の Fe 原子（酸化物の周辺原子）の距離変化からは，第一近接の Fe 原子を強く引き寄せるため，引張によって [001] 方向の第二近接のFe-Fe 間がアフィン変形より大きく延伸されることが示された．

65

第6章

Y2O3 含有 Fe のせん断シミュレーション
本章では Fe 単元系ならびに Y2 O3 を添加したモデルに対して [111] 方向のせん断シ
ミュレーションを行い，前章と同様に格子安定性，電子状態等について検討する．Y2 O3の添加サイトの組み合わせは第 4 章で最もエネルギー的に安定な構造であった S-octa2や octa2 と等価なものとする．またせん断過程の各ひずみにおいて弾性剛性係数の正値性に基づく格子不安定解析を行い，Y2 O3 がせん断方向の理想強度に与える影響について検討する．

6.1

解析条件

前章と同様に，Kresse らにより開発された平面波基底ウルトラソフト擬ポテンシャ
ル法 (59) に基づく第一原理バンド計算コード VASP(73) (Vienna Ab-initio Simulation

Package) を用いて行った．交換相関表現ならびに電子状態の収束計算は前章と同様である．解析には図 6.1 で模式的に示すように bcc-Fe に対して x, y, z 軸を [111], [112̄], [1̄00]
方向にとり，6 原子の bcc-Fe で構成されるスーパーセル (Fe6 モデル) を基本とした. 図では [111] 方向に 2 個の基本となるスーパーセルを並べて表示し，bcc 格子の [1̄10] 面に着色している．さらにこの Fe6 モデルを 4 × 1 × 2 または 4 × 1 × 1 並べたスーパー
セルに Y を 2 個置換，O 原子を 3 個添加し, それぞれ octa2 モデルと S-octa2 モデルとした（図 6.2）．その他の計算条件を表 6.1 に示す．まずセル内での原子核の移動は考えず，図の位置に固定した状態で基準となるスーパーセルの寸法を等方的に拡大，縮小して電子状態計算を行い，bcc 格子長さ a に対する全自由エネルギーの変化を調べた．その時全自由エネルギーが最小となる格子長さを系の基準格子長さ a0 とした．次66

に原子核固定の条件のもと，[111] 方向にせん断ひずみ △ϵ13 = 0.01 ∼ 0.03 を与え電子状態計算を行った．ひずみ増加/電子状態計算を繰り返し，静力学的に応力-ひずみ曲線を算出した．y 方向と z 方向は拘束 ε22 = ε33 = 0 とした．また各ひずみ下において，

x, y, z 軸を [100],[010],[001] 方向を基準にとり，微小ひずみ摂動を第 4 章と同様に与え，応力変化から弾性剛性係数を数値的に算出し，系の安定性を評価した．

[010]

2.44949 x aaz [110][001]

ay [112]a1.414214 x ax [111]0.866025 x a[100]

Fig.6.1 Simulation model(Fe6 model).

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

Fe6 model
636296
22 × 10 × 14

67

octa2 model
51
298
495
4×4×4

S-octa2 model
27
154
495
4×4×4

3.4641 x a

3.4641 x az [110]Oy [112]Y

2.828428 x a

O

Y

Y

Y

O

x [111]OO

2.44949 x a(a) octa2 model

2.44949 x a(b) S-octa2 modelz [110]y [112]x [111](c) Front view of octa2 and S-octa2.

z [110]y [112]x [111](d) Side view of octa2 and S-octa2.

y [112]z [110]1.414214 x aOx [111](e) Top view of octa2 and S-octa2.

Fig.6.2 Simulation models(octa2 and S-octa2 models).
68

6.2

解析結果と考察6.2.1

応力-ひずみ関係と各方向の応力変化[111] 方向のせん断変形で得られた各モデルの応力-ひずみ関係を図 6.3 に示す．図中の矢印は応力ピークのひずみを示している．せん断初期（ε13 = 0 ∼ 0.03）の応力-ひずみ勾配は Fe6 (51.1GPa)>octa2 (49.8GPa)>S-octa2 (41.6GPa) となっており，
引張シミュレーションと同様に Fe 単元系に Y2 O3 を添加することで剛性は低下している．変形前の各モデルの bcc 格子のセル辺長 a0 は Fe6 (2.857Å ), octa2 (3.015Å ),S-octa2 (3.162Å ) である．したがって，Y2 O3 を添加することですべり面の間隔が広がったこともせん断方向の剛性が低下した一因と考えられる．いずれのモデルもひずみ

ε13 = 0.13 までほぼ線形的に応力は上昇し，Fe 単元系モデルにおいても引張変形に見られた変形初期の変曲点は存在しない．応力-ひずみ曲線のピーク応力 σ13 の大きさは

Fe6 (10.48GPa)>octa2 (10.35GPa)>S-octa2 (9.31GPa) であり，Fe 単元系に Y2 O3 を添加することにより最大応力で見た弾性限界は低下するが，ピークひずみ ε13 は S-octa2(ε13 = 0.39)>octa2 (ε13 = 0.35)>Fe6 (ε13 = 0.30) となり，ひずみでみると Y2 O3 の添加で限界は向上する．図 6.4 に模式的に示すように，各モデルの [111] 方向に隣接する

Fe 原子間距離を l0 ，ピーク応力時の [111] 方向への原子の移動量 △l とすると，Fe6 モ
デルは △l/l0 = 0.49 となり，[111] 方向に原子間距離をほぼ半分シフトさせた位置でピーク応力を示している．一方，octa2，S-octa2 はそれぞれ △l/l0 = 0.57, 0.63 と半分
以上シフトした点でピークを生じる．x, y, z 軸を [100],[010],[001] 方向に取り直し，せん断変形下の各方向の応力の変化を図 6.5 に示す．せん断変形によって bcc セルが非等方的に変形するため，各方向の応力状態が大きく異なる．いずれのモデルでも [111]方向のせん断変形によって bcc セルが [010] 方向へ延伸しているため，せん断初期において [010] 方向の応力は上昇しているが，[100] 方向へは収縮しているため圧縮の応力状態となり減少していく．ひずみ ε[111] = 0.1 までの値を比較すると，どちらの成分でも Fe6 の変化が最も大きく，[100] と [010] 方向の応力差が最も大きい．

69

12

Fe6
octa2
10

S-octa2, GPa

8Stress

6

4

2

0
0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

0.4

0.45

0.5

Tensile strain

Fig.6.3 Stress-strain curves.

l
C

A

A
C
z [110]y [112]Bz [110]x [111]Dy [112]x [111]

Fig.6.4 Deﬁnition of atom shift △l.

70

B

l0

D

10

10

5

5

,GPa

15

0

Stress,GPa
Stress

15

-5

0
-5

-10

-10

-15

-15

-20

-20
0

0.05

0.1

0.15

0.25

0.3

Shear strain

0.2

[111]

0.35

0.4

0.45

0.5

0

(a) Fe6 model

0.05

0.1

0.15

0.25

0.3

Shear strain

0.2

[111]

0.35

0.4

0.45

0.5

(b) octa2 model15
xx [100]

10yy [010]zz [001]

Stress,GPa5yz [011]xz [101]0xy [110]

-5
-10
-15
-20
0

0.05

0.1

0.15

0.25

0.3

Shear strain

0.2

[111]

0.35

0.4

0.45

0.5

(c) S-octa2 model

Fig.6.5 Stress-strain curves in x[100] y[010] z[001] orientation.

6.2.2

弾性剛性係数の行列式の変化

図 6.6 に弾性剛性係数の行列式 detBij の値の変化を示す．図左は detBij の値をそのまま，右は引張前の detBij の値で割り，無次元化して示している．矢印は初めて

detBij が負になるひずみを示している．格子不安定ひずみは S-octa2 (ε13 = 0.30)>octa2(ε13 = 0.29)>Fe6 (ε13 = 0.25) となり，せん断に対しては Y2 O3 の含有割合の高いモデ
ルの弾性限界が高い値を示したが，引張時のような Fe 単元系との大きな差はない．図6.7 に各モデルの応力-ひずみ曲線を再掲し，弾性限界を矢印で示した．格子不安定解
析による弾性限界での応力（理想せん断強度）は Fe6 (10.14GPa)>octa2 (10.04GPa)>

S-octa2 (8.56GPa) となり，応力を基準とすれば Fe 単元系と Y2 O3 の含有割合が小さ
い octa2 はほぼ同じで，Y2 O3 の含有割合が高い S-octa2 が 1.6GPa 低下する．

71

12

(x10 )
1.5

Magnitude of 6x6 determinant (normalized)

Fe6
octa2

0.6

S-octa2

0.4

0.2

0

-0.2

1

0.5

0

-0.5

-1

-0.4
0.05

0.1

0.15

0.2

0.25

0.3

0

0.35

0.05

0.1

Shear strain

0.15

Fig.6.6 Changes in the determinants of Bij .

12

Fe6
octa2
10

S-octa2
8

6

4

2

0
0

0.05

0.1

0.2

Shear strain

, GPa

0

StressMagnitude of 6x6 determinant, GPa

0.8

0.15

0.2

0.25

0.3

0.35

0.4

0.45

Tensile strain

Fig.6.7 Stability limits on stress-strain curves.

72

0.5

0.25

0.3

0.35

6.2.3

せん断変形下の価電子密度分布の変化

図 6.8∼6.12 に Fe 単元系モデル（Fe6）と Y2 O3 を添加したモデル（S-octa2）のせん
断変形下の価電子密度分布の変化を示す．図右上には表示している面を模式的に示した．なお，octa2 および S-octa2 の Y2 O3 を添加するサイトの組み合わせは同じである
ため，応力-ひずみの関係から，Y2 O3 添加による効果が大きく現れた S-octa2 の価電子密度分布のみ示し，octa2 は省略した．Fe6 モデルの価電子密度分布および模式図はS-octa2 との比較のため [111] 方向にスーパーセルを 4 つ並べた状態で示している．Fe6の図 6.8 に示した面では，せん断変形が進行するとともに矢印で示す価電子密度の谷の部分に等高線が増えていることから，価電子密度が減少しすべり面間の結合が低下していると考えられる．また図 6.9 に示す面をみると，せん断変形下で第一近接の関係にある Fe(a)-Fe(b) 間の距離が離れていき，第二近接の関係にある Fe(a)-Fe(c) が近づ
く．ひずみ ε13 = 0.2 および格子不安定となるひずみ ε33 = 0.25 では Fe(a)-Fe(b) 間の価電子密度が低下し，Fe(a)-Fe(c) 間の価電子密度と同程度になっている．そして，不安定後のひずみ ε13 = 0.3 では Fe(a)-Fe(c) 間が高い価電子密度を示している．つまり第一近接の関係にある Fe 原子間の結合が低下し，原子間の結合が変化する過程で不安定となっている．一方，S-octa2 の O 原子を添加した面（図 6.10）をみると，せん断変形下で第二近接の関係にある Fe(e) と Fe(f) の距離は短くなっていくが，挟まれた O
原子が Fe の接近に対して抵抗になっていると考えられる．また，ひずみ ε13 = 0 で O
原子近傍の Fe(d) 原子は第一近接の Fe(e) 原子よりも O 原子との間で高い価電子密度を示し，せん断変形が進行しても高い値を維持している．O 原子は Fe 原子の八面体
中心（(001) 面の面心位置）にあり，本来第二近接にある [001] 方向の Fe-Fe 結合に強い Fe-O-Fe 結合を形成しており，第一近接→第二近接への結合の切り替えの大きな抵抗となる．O 原子を添加した他の面（図 6.11）では，O 原子と隣接する Fe 原子と第一
隣接 Fe 原子との間の価電子密度 （円 A）はその他の Fe-Fe 原子間（円 B）より高い．
また，不安定となる前後では，図 6.11 に矢印で示すように O 原子が存在しないすべり
面の間で等高線が増えており，価電子が減少していることがわかる．Y 原子が添加された面（図 6.12）では，Y によってすべり面の間の価電子が低く，すべり面の間の結73

合が低下している．また，Fe 単元系同様に第一近接→第二近接への結合の切り替えが生じている．ひずみ ε13 = 0.1 の図中に矢印で示した部分では，Y によって価電子密度が低下している箇所がせん断変形下で Fe 間を移動していることから，Y 原子によってFe-Fe 間の結合の低下および B13 が低下していると推測できる．

74

z [110]
8.4

y [112]x [111]0(eV/Å3)00.3 (Stress peak point)

0.1

0.35

0.2

0.39

0.25 (Stability limit)

0.4

Fig.6.8 Changes in the distribution of valence electron density(Fe6)．

75

z [110]
8.4

y [112]x [111]0(eV/Å3)

Fe(c)

Fe(b)

Fe(a)
0

0.3 (Stress peak point)

0.1

0.35

0.2

0.39

0.25 (Stability limit)

0.4

Fig.6.9 Changes in the distribution of valence electron density(Fe6)．

76

O
O

z [110]8.4

y [112]O

YYOx [111]0(eV/Å3)

Fe(e)
O
Fe(d)

Fe(f)

0

0.3 (Stability limit)

0.1

0.35

0.2

0.39 (Stress peak point)

0.4

0.25

Fig.6.10 Changes in the distribution of valence electron density(S-octa2)．

77

O
O

z [110]8.4

y [112]O

YYOx [111]0(eV/Å3)OBA
0

0.3 (Stability limit)

0.1

0.35

0.2

0.39 (Stress peak point)

0.25

0.4

Fig.6.11 Changes in the distribution of valence electron density(S-octa2)．

78

O
O

z [110]8.4

y [112]O

YYOx [111]0(eV/Å3)YY

0

0.3 (Stability limit)

0.1

0.35

0.2

0.39 (Stress peak point)

0.25

0.4

Fig.6.12 Changes in the distribution of valence electron density(S-octa2)．

79

6.3

結言

第 4 章でエネルギー，弾性特性，溶解熱いずれの点でも有利であった octa2 構造と同じ添加サイトに Y2 O3 を添加したモデル（octa2，S-octa2）について，x, y, z 軸を

[111], [112̄], [1̄10] 方向にとり，第一原理計算による [111] 方向せん断シミュレーションを行った．比較のため Fe 単元系 (Fe6 モデル) に対しても同様の解析を行った．得られた結果を以下に示す．

(1) せん断初期の応力-ひずみ勾配は Fe6 (51.1Gpa)>octa2 (49.8GPa)>S-octa2 (41.6GPa)
であり，Y2 O3 を添加することで剛性は低下した．Y2 O3 を添加することですべ
り面の間隔が広がっており，これも剛性低下の一因と考えられる．

(2) 応力-ひずみ曲線のピーク応力 σ13 の大きさは Fe6 (10.48GPa)>octa2 (10.35GPa)>
S-octa2 (9.31GPa) となり，ピーク応力による理想強度は Y2 O3 を添加することにより低下する．ただし，ピークひずみ ε13 は S-octa2 (ε13 = 0.39)>octa2(ε13 = 0.35)>Fe6 (ε13 = 0.30) であり，ひずみでみると Y2 O3 の添加で限界ひずみは向上する．

(3) [111] 方向に隣接する Fe 原子間距離を l0 ，ピーク応力時の [111] 方向への原子の移動量 △l として評価したところ，Fe6 モデルは △l/l0 = 0.49 となり，ピークは原子間距離をほぼ半分シフトさせた位置であった．一方，octa2，S-octa2 はそれ
ぞれ △l/l0 = 0.57, 0.63 と半分以上移動した後にピークを生じていた．

(4) detBij が初めて負になる格子不安定ひずみは S-octa2 (ε13 = 0.30)>octa2 (ε13 =0.29)>Fe6 (ε13 = 0.25) となり，(2) のピーク応力に対する限界ひずみよりそれぞれ低下するものの大小関係は変わらない．また，引張のような Fe 単元系との大きな差はなかった．

(5) 格子不安定点での応力は Fe6 (10.14GPa)>octa2 (10.04GPa)>S-octa2 (8.56GPa)となり，Y2 O3 の含有割合が高い S-octa2 が他よりも小さくなる．

80

(6) せん断変形下の価電子密度変化を調べたところ，Fe 単元系（Fe6）において，格子不安定となるひずみ ε33 = 0.25 では第一近接の関係にある原子間の価電子密度が低下し，第二近接の関係にある原子間の価電子密度と同程度になった．不安
定後のひずみ ε13 = 0.3 では第二近接の関係にあった原子間の価電子密度が高くなっていることから，この点ですべり面間の結合の切断・再結合を生じていることがわかる．

(7) Y2 O3 を添加した系では，O 原子は Fe 原子の八面体中心（(001) 面の面心位置）にあり，本来第二近接である [001] 方向の Fe-Fe 結合に強い Fe-O-Fe 結合を形成する．このため [111] 方向のせん断による第一近接→第二近接への結合の切り替
えの大きな抵抗となる．一方，Y を含むすべり面では価電子密度が疎であり，格子不安定となる点ではこの面で結合の切り替わりを生じていた．

81

第7章

結論

本研究では，ODS 鋼の機械的特性，強化メカニズムの解明につながる新たな知見を得ることを目的として，bcc-Fe に Y2 O3 または YTiO3 を添加した系の第一原理計算を行い，その構造安定性を自由エネルギーや溶解熱の大小，弾性剛性係数の正値性に基づく格子不安定性の観点から検討した．また，Y2 O3 を添加したモデルに対して，引張並びにせん断変形シミュレーションを行い，Y2 O3 が理想強度に与える影響を検討した．得られた結果を以下に統括する．
第 2 章では，平面波基底での第一原理分子動力学法の基礎理論および電子状態計算の高速化手法について説明し，第 3 章では，格子安定性の概念と，弾性剛性係数を用いた格子不安定解析について述べた．
第 4 章では，Y2 O3 の bcc-Fe 中における安定構造や，力学特性の変化を O 原子の添
加サイトの違いから検討するために，Fe の bcc 単位格子を 3 × 3 × 3 並べたスーパーセルに Y 原子を 2 個，O 原子を 3 個添加した系，ならびに Y 原子 1 個の代わりに Ti 原子
を 1 個置換した系について第一原理計算を行った．今回考慮した構造の弾性剛性係数の行列式（detBij ）の値は，いずれも正の値をとり「力学的に」存在し得る構造であることが示唆された．価電子密度分布をみると O 原子は八面体サイトを構成する隣接

Fe 原子と強く結合し，さらに次に隣接する Fe 原子との間にも価電子密度を増加させていた．3 個の O 原子を Fe のみで構成される八面体の中心サイトに添加した構造は，八面体の短軸方向が [100],[010],[001] であり，各隣接八面体が 1 つの Fe 頂点でのみで接するため，この O 原子による結合が等方的になり，最も弾性剛性係数の異方性が小さく，O 原子の溶解熱もすべて同程度の負の値を示した．Y2 O3 , YTiO3 のいずれを添加しても Fe 単元系より全自由エネルギーは大きくなるものの，この構造の全自由エネ
ルギーが低い値を示しており，エネルギー，弾性特性，溶解熱いずれの点でも有利で82

あることが明らかになった．
第 5 章では，エネルギー，弾性特性，溶解熱いずれの点でも有利であった上記の添加
サイトの組み合わせに Y2 O3 を添加したモデルおよび Fe 単元系について，第一原理計
算による引張シミュレーション（横ひずみ ε′ = 0）を行った．ひずみ ε33 = 0 におけるヤング率は，Y2 O3 を添加することで低下し，応力-ひずみ曲線においてもピーク応力 σ33
およびピークひずみ ε33 は Fe 単元系モデルが最大の値を示した．ただし，Fe 単元系は引張初期で変曲点が生じており，実際にはその点で detBij が負になっていた．したがって，格子不安定性の観点からは理想引張強度は Fe 単元系が ε33 = 0.03, σ33 = 3.8GPa，

Y2 O3 を添加した系は ε33 = 0.11 ∼ 0.16, σ33 = 13.9 ∼ 16.8GPa （モデルの大小，構造緩和の有無による差）となった．Y2 O3 を添加した系の引張変形下の価電子密度分布を検討した結果，O 原子を含む八面体の短軸方向の Fe-O-Fe 結合の価電子密度が高く，この結合が [100],[010],[001] のいずれの方向にも存在するために，3 軸引張に近い ε′ = 0の引張で弾性限界が向上したことが示唆された．また，構造緩和も行った系において各結合の引張による距離変化を調べたところ，dz2 原子軌道のような形状の Y-(FeO3 )-Y
構造において，中心軸の Y-Fe-Y はほぼアフィン変形しているが，中心 Fe を囲む O の三角形リングは収縮することがわかった．
第 6 章では，x, y, z 軸を [111], [112̄], [1̄10] 方向とするスーパーセルで，Fe 単元系ならびに第 5 章と等価なサイトに Y2 O3 を添加した系について，[111] 方向にせん断を与えるシミュレーションを行った．
せん断初期の応力-ひずみ勾配は，Y2 O3 の含有割合が高い程低下した．応力-ひずみ曲線のピーク点は，Fe 単元系が ε13 = 0.3, σ13 = 10.5GPa，Y2 O3 を添加した系は

ε13 = 0.35 または 0.39, σ13 = 9.3 または 10.4GPa （モデルの大きさによる差）となり，
Y2 O3 を添加することでピーク応力は低下するが，ピークひずみは向上した．[111] 方
向に隣接する Fe 原子間距離を l0 ，[111] 方向への原子の移動量を △l とすると Fe6 モデ
ルは △l/l0 = 0.49 でピーク応力を示しており，[111] 方向に半分シフトさせた点が応力ピークに相当する．一方，Y2 O3 を添加したモデルは △l/l0 = 0.57 または 0.63 と，半分以上シフトした点でピークを示した．格子不安定となるひずみで弾性限界を評価す83ると，いずれも応力ピークより低ひずみとなったが Fe 単元系と Y2 O3 を添加した系の大小関係は変わらず，理想強度の関係も変わらなかった．せん断変形下の価電子密度分布を調べた結果，Fe 単元系では [111] 方向のせん断によって，bcc 格子の第一近接原子が離れ，第二近接原子との原子間結合の切り替えが生じる点で格子不安定となっていた．Y2 O3 を添加した系では O 原子が bcc 格子の (001) の面心位置にあるため，[001]方向の Fe-Fe 間（第二近接の関係にある）に強い Fe-O-Fe 結合を形成し，先の Fe 単元系で観察した結合の切り替えの抵抗となっていた．一方，Y が存在する原子面では価電子密度が疎になっており，格子不安定となる点では，この面で第一近接から第二近接への Fe-Fe 結合の切り替えを生じていた．
以上のように，本研究では bcc-Fe 中に Y2 O3 を添加し，安定な構造や力学特性を様々な観点から評価した．実際の ODS 鋼中では酸化物は数 nm の大きさで析出しており，添加元素の割合も今回の解析に比べ桁違いに小さい．したがって，より多数の原子の中で Y2 O3 の析出構造（Bixbyite 構造）を再現した評価が望ましいが，現時点では第一原理計算で直接取り扱うことは不可能である．今後，経験的なポテンシャルを用いた MC および MD シミュレーションなどに本解析結果が反映されることを期待する．

84

付録 A 関連学術講演
田仲稔，睦門賢憲，屋代如月”第一原理計算による酸化物含有 Fe の引張シミュレーションと格子不安定解析”
第 15 回分子動力学シンポジウム，札幌コンベンションセンター，(2010.6)

田仲稔，古賀健治，世木隆，屋代如月”第一原理格子不安定解析による酸化物含有 Fe の構造安定性ならびに理想引張強度評
価”
日本機械学会第 23 回計算力学講演会，北見工業大学，(2010.9)

田仲稔, 古賀健治, 世木隆, 奥田隆成, 屋代如月

”Fe マトリックス中における酸化イットリウムの安定サイト探索：第一原理格子不安
定解析”
第 1 回マルチスケールマテリアルモデリングシンポジウム （第 16 回分子動力学シンポジウム），大阪大学コンベンションセンター，(2011.5)

田仲稔，世木隆，古賀健治，屋代如月”Fe マトリックス中における Y・Ti 酸化物の安定サイト探索：第一原理格子不安定解
析”
日本機械学会第 24 回計算力学講演会（CMD2011），岡山大学 津島キャンパス，(2011.10)

85

md

86

md

87

cmd

88

cmd

89

md

90

md

91

cmd

92

cmd

93

参考文献
(1) 井上賢紀, “酸化物分散型フェライト鋼の開発で理想の高速炉燃料被覆管材料を実現”, JAEA ニュース, 27, (2011), 2.

(2) D.K. Mukhopadhyay, F.H. Froes, D.S. Gelles, “Development of oxide dispersion
strengthened ferritic steels for fusion”, Journal of Nuclear Materials, 258-263,(1998), 1209-1215.
(3) P. Pareige a, M.K. Miller, R.E. Stoller, D.T. Hoelzer, E. Cadel, B. Radiguet,“Stability of nanometer-sized oxide clusters in mechanically-alloyed steel under ion-induced displacement cascade damage conditions”, Journal of Nuclear
Materials, 360, (2007), 136-142.
(4) Y. Yano, R. Ogawa, S. Yamashita, S. Ohtsuka, T. Kaito, N. Akasaka, M. Inoue,
T.Yoshitake, K. Tanaka, “Eﬀects of neutron irradiation on tensile propertiesof oxide dispersion strengthened (ODS) steel claddings”, Journal of Nuclear
Materials, 419, (2022), 305-309.
(5) S. Ukai, M. Fujiwara, “Perspective of ODS alloys application in nuclear environments”, Journal of Nuclear Materials, 307-311, (2002), 749-757.
(6) M. Brocq, B. Radiguet, J.M. Le Breton, F. Cuvilly, P. Pareige, F. Legendr,“Nanoscale characterisation and clustering mechanism in an Fe-Y2O3 model
ODS alloy processed by reactive ball milling and annealing”, Acta Materialia,
58, (2010), 1806-1814.
(7) M. Yamamoto, S. Ukai, S. Hayashi, T. Kaito, S. Ohtsuka, “Formation of residual
ferrite in 9Cr-ODS ferritic steels”, Materials Science and Engineering, A-527,(2010), 4418-4423.
(8) S. Kim, S. Ohtsuka, T. Kaito, S. Yamashit, M. Inoue, T. Asayama, T. Shobu,“Formation of nano-size oxide particles and d-ferrite at elevated temperature in
9Cr-ODS steel”, Journal of Nuclear Materials, 417, (2011), 209-212.

94

(9) H. Sakasegawa, “Radiation Eﬀect on Oxide-Dispersion-Strengthened Ferritic/Martensitic Steel under Inﬂuence of Applied Stress and Improvement of
Shape Stability under Irradiation”, JAEA-Research, 53, (2007).
(10) H. Sakasegawa, S. Ukai, M. Tamura, S. Ohtsuka, H. Tanigawa, H. Ogiwara,
A. Kohyama, M. Fujiwara, “Stability of Y-Ti complex oxides in Fe-16Cr-0.1Ti
ODS ferritic steel before and after heavy-ion irradiation”, Journal of Nuclear
Materials, 373, (2008), 82-89.
(11) S. Ukai, S. Ohtsuka, T. Kaito, H. Sakasegawa, N. Chikata, S. Hayashi, S.
Ohnuki, “High-temperature strength characterization of advanced 9Cr-ODS ferritic steels”, Materials Science and Engineering, 510-511, (2009), 115-200.
(12) H. Kishimoto, R. Kasada, O. Hashitomi, A. Kimura, “Stability of Y-Ti complex
oxides in Fe-16Cr-0.1Ti ODS ferritic steel before and after heavy-ion irradiation”, Journal of Nuclear Materials, 386-388, (2009), 533-536.
(13) D.A. McClintock, M.A.Sokolov, D.T.Hoelzer, R.K. Nanstad, “Mechanical properties of irradiated ODS-EUROFER and nanocluster strengthened 14YWT”,
Journal of Nuclear Materials, 392, (2009), 353-359.
(14) R. Schaeublin, T. Leguey, P. Spatig, N. Baluc, M. Victoria, “Microstructure
and mechanical properties of two ODS ferritic/martensitic steels”, Journal of
Nuclear Materials, 307-311, (2002), 778-782.
(15) R.L. Klueh, J.P. Shingledecker, R.W. Swindeman, D.T. Hoelzer, “Oxide
dispersion-strengthened steels: A comparison of some commercial and experimental alloys”, Journal of Nuclear Materials, 341, (2005), 103-114.
(16) S. Ohtsuka, S. Ukai, M. Fujiwara, T. Kaito, T. Narita, “Improvement of 9CrODS martensitic steel properties by controlling excess oxygen and titanium
contents”, Journal of Nuclear Materials, 329-333, (2004), 372-376.
(17) M. Ohnuma, J. Suzuki, S. Ohtsuka, S.W. Kim, T. Kaito, M. Inoue, H. Kitazawa“A new method for the quantitative analysis of the scale and composition ofnanosized oxide in 9Cr-ODS steel”, Acta Materialia, 57, (2009), 5571-5581.

95

(18) H. Kishimoto, M.J. Alinger, G.R. Odette, T. Yamamoto, “TEM examinationof microstructural evolution during processing of 14CrYWTi nanostructured
ferritic alloys”, Journal of Nuclear Materials, 329-333, (2004), 369-371.
(19) M.J. Alinger, G.R. Odette, D.T. Hoelzer, “On the role of alloy composition andprocessing parameters in nanocluster formation and dispersion strengthening in
nanostuctured ferritic alloys”, Acta Materialia, 57, (2009), 392-406.
(20) S. Ohtsuka, T. Kaito, M. Inoue, T. Asayama, S.W. Kim, S. Ukai, T. Narita,
H. Sakasegawa, “Eﬀects of aluminum on high-temperature strength of 9Cr-ODS
steel”, Journal of Nuclear Materials, 386-388, (2009), 479-482.
(21) “スーパー ODS 鋼”, 公開特許公報, 2010-65302, (2010).
(22) H.M. Zibib, M. Rhee, J.P. Hirth, “On Plastic Deformation and The Dynamics
of 3D Dislocations”, Int. J. Mech. Sci., 40–2–3 (1998), 113–127.
(23) V. Giessen, A. Needleman, Discrete Dislocation Plasticity, Modeling Simulation
Mater. Sci. Eng., 3 (1995), 689–735.
(24) 赤星 保浩, 福田 忠生, 原田 昭治, 学術講演会講演論文集, 50, (2001), 277-278.
(25) Y. Aoki, A. Takahashi, M. Kikuchi, JSME annual meeting 2006(6), (2006),
57-58.
(26) B. Bako, D. Weygand, M. Samaras, J. Chen, M.A. Pouchon, P. Gumbsch, W.
Hoﬀelner, Philosophical Magazine, 87 (2007),3645-3656.
(27) 睦門賢憲, “分子動力学法による酸化物分散強化メカニズムの検討：転位のカッ
ティング抵抗の評価”, 神戸大学修士論文, (2011).

(28) 山本智, “第一原理格子不安定解析による酸化物分散強化機構の基礎的研究”, 神戸大学修士論文, (2010).

(29) S. Ogata, J. Li, N. Hirosaki, Y. Shibutani, S. Yip, “Ideal shear strain of metals
and ceramics”, Physical Review, B 70, (2004), 104.
(30) 尾方成信, 北川浩, “第一原理分子動力学計算による窒化アルミニウム／アルミ
ニウム界面の構造評価”, 日本金属学会誌, 60, (1996), 1079.

96

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

(32) A. Gopejenko, Y. Zhukovskii, P.V. Vladimirovb, E.A. Kotomin, A. Moslang,“Modeling of yttrium, oxygen atoms and vacancies in γ-iron lattice”, Journalof Nuclear Materials, 416, (2011), 40-44.
(33) A. Gopejenko, Y.F. Zhukovskii, P.V. Vladimirov, E.A. Kotomin, A. Moslang,“Ab initio simulation of yttrium oxide nanocluster formation on fcc Fe lattice”,
Journal of Nuclear Materials, 406, (2010), 345-350.
(34) Y. Jiang, J.R. Smith, G.R. Odette, “Formation of Y-Ti-O nanoclusters in nanostructured ferritic alloys: A ﬁrst-principles study”, PHYSICAL REVIEW, B 79,(2009), 064103.
(35) L.Barnard, G.R.Odette, I. Szlufarska, D. Morgan, “An ab initio study of Ti-YO nanocluster energetics in nanostructured ferritic alloys”, Acta Materialia, 60,(2012), 935-947.
(36) M. Born, Huang, Dynamical Theory of Crystal Lattices, (1954) Oxford UP.
(37) J. Wang, S. Yip, S.R. Phillpot, D. Wolf, Phys. Rev. B, 20 (1993), 4182.
(38) M. Šob, L.G. Wang, V. Vitek, Mat Sci. Eng. A, 234 (1997), 1075.
(39) D.M. Clatterbuck, D.C. Chrzan, J.W. Morris, Phil. Mag, Lett, 82 (2002), 141.
(40) 屋代如月, 大穂正史, 山上勝也, 冨田佳宏, 材料, 52, (2002), 241.
(41) M. Černý, J. Pokluda, J.Alloys and Compounds, 378 (2004),1(42) A.V. Krajnikov, A.N. Demidik, H.M. Ortner, Materials Science and EngineeringA, 234 (1997),357.
(43) R. Schaublin, A. Ramar, N. Baluc, V.de Castro, M.A. Monge, T. Leguey, N.
Schmid, C. Bonjour, Journal of Nuclear Materials, 351 (2006),247.
(44) S. Curtarolo, D. Morgan, G. Ceder, Calphad, 29 (2005),163.
(45) W.Y. Ching, N.X. Yong, Journal of Magnetism and Magnetic Materials, 209(2000),28.
97

(46) W. Luo, D. Roundy, M. Cohen, J.W. Morris Tr, Phys. Rev. B, 66 (2002), 094110.
(47) D.M. Clatterbuck, D.C. Chrzan, J.W. Morris Jr, Acta Materialia, 51(2003),2271.
59.
(48) P. Hohenberg, W. Kohn, Phys. Rev., 136, (1964), 864.
(49) W. Kohn, L.J. Sham, Phys. Rev., 140, (1965), 1133.
(50) J.P. Perdew, A. Zunger, Phys. Rev.,B, 23, (1981), 5048.
(51) キッテル，固体物理学入門 上 第 7 版 (1998), 丸善.
(52) G.B. Bachelet, D.R. Hamann, M. Schlüter, Phys. Rev. B, 26, (1982), 4199.
(53) N. Troullier Phys. Rev. B, 43, (1991), 1993.
(54) L.D. ランダゥ，E.M. リフシッツ著，佐々木健，好村磁洋訳，量子力学 1，(1992),
東京図書.

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

98

(65) M.J. Gillan, J. Phys. Condens. Matter., 1 (1989), 689.
(66) I. Stich, Phys. Rev. B, 39, (1989), 4997.
(67) D.M. Bylander, Phys. Rev. B, 42, (1990), 1394.
(68) F. Milstein, Phys. Rev. B, 3 (1971), 1130.
(69) D.C. Wallace, Thermodynamics of Crystals, (1972), Wiley, Newyork.
(70) 屋代如月, 学位論文，“金属結晶の変形と破壊の分子動力学シミュレーション並びに局所格子不安定性解析”,(1998).

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

99

謝辞
本研究を進めるにあたり，厚くご指導賜りました屋代如月准教授に心より感謝の意を表します．浅学非才な著者が研究を続け，論文執筆を完遂することができたのも先生に親切丁寧にご指導いただたいた事によるものでした．研究活動に対して幅広い視点
からの御助言により，著者の狭い見解を広げて下さった田中克志教授に感謝致します．
本研究を進める上で貴重な御助言・御協力をいただいた古賀健治様（株式会社コベルコ科研）をはじめ，株式会社コベルコ科研の方々に厚く御礼申し上げます．研究室配属時の何も分からない著者に対して良き指針となり温かく接して頂いた諸先輩方に感謝致します．また，第一原理計算グループとして丁寧に研究についてご助言頂いた山
本智氏 (現・川崎重工業株式会社)，ODS 研究グループとして常に助言を下さり時にお互いの意見を議論した睦門賢憲氏（現・株式会社神戸製鋼所）に深くお礼申し上げます．お二人のおかげで研究活動をスムーズに行うことができたものと思っております．

3 年間の研究室生活において公私共に多くの時間をすごし，切磋琢磨しあいながら貴重な経験をしてきた同輩の谷口誠典君，山口明宏君，小楠貴治君，工藤修平君，永井
和広君，中山隆史君に感謝致します．最後に学生生活において，著者が経済的・精神的困難に陥ることのないように積極的に支援してくれた両親，家族に最大限の感謝と敬意を表してここに記します．ありがとうございました．

2012 年 2 月

100

田仲 稔