修士論文

アモルファス合金設計のための
短距離秩序構造の第一原理格子不安定解析

指導教員：冨田 佳宏

060T371N 久馬 雅彦

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

Master thesis

Ab-initio Lattice Instability Analysis
on Short-Range Ordered Structure
for Amorphous Alloy Design

Masahiko KYUMA
February 2008
Department of Mechanical Engineering,
Graduate School of Science and Technology,
Kobe University,Kobe,Japan1

要約
アモルファス合金設計の指標の 1 つとすべく，過冷却液体状態の安定化のキーとなる 20 面体クラスターについて，第一原理計算により弾性係数 CIJ を求め，格子不安定性ならびに電子状態を検討した．
まず，結晶との違い，ならびに周囲原子の有無による変化を検討するものとして，

Cu，Ni，Al，Zr，Ti それぞれ単元素について孤立クラスター，クラスターが周期的に並んだ準結晶，完全結晶の解析を行った．いずれの元素も完全結晶では detCIJ は正値となり，大小関係は Ni ＞ Cu ＞ Ti ＞ Al ＞ Zr であった．クラスターとした場合は，孤立，準結晶のいずれも Cu，Ni が安定，Al，Zr，Ti は不安定となった．次に，孤立クラスターの中心原子を各元素で置換し，合金化による安定性変化を検討した．Cu，Niを含むクラスターは，中心原子が Zr の場合，または外殻原子が Al のときだけ不安定となり，他は安定であった．Al を外殻原子とするクラスターは全て不安定であったが，不安定な Ti クラスターの中心を Al にすると安定となった．最後に，クラスター中の元素割合を変えたときの安定性変化を調べるために，Cu を中心とし，外殻原子の Zr 原子の割合を変えた Cu-Zr クラスターの解析を行った．いずれのクラスターも detCIJ は正となり安定であったが，Zr の外殻原子を 1 つだけ Cu とし，Cu-Cu ボンドが中心−外殻原子に 1 つだけ存在するクラスターは detCIJ が他に比べ大きい値となった．

2

SUMMARY
Toward a new design for amorphous alloys, we have evaluated the elastic coeﬃcients, CIJ , of various icosahedral atomic clusters, which play an important role inthe stabilization of metallic supercooled liquid, and discussed the lattice stability andelectron structures based on the ab-initio molecular dynamics simulation. First, wehave analysed monatomic Cu, Ni, Al, Zr and Ti isolated icosahedrons, periodic array of icosahedrons and perfect lattices in order to disucuss the inherent diﬀerence
between perfect lattice and icosahedron or the eﬀect of surrounding atoms. All theelement show positive value and the magnitude relation is Ni ＞ Cu ＞ Ti ＞ Al ＞ Zr.
For icosahedrons, both isolated and periodic array are stable with the element of Cu or
Ni, while they are unstable with Al, Zr and Ti. Then, the center atom of the isolated
icosahedrons is substituted with other element to discuss the stability change by alloying. Icosahedrons which contain Cu or Ni become unstable only if they have the center
atom of Zr or they have the outer shell of Al. The icosahedrons cannot be stable if they
have outer shell of Al element, however the center Al atom stabilizes the unstable Ti
icosahedron. Finally, we have increased the number of Cu atoms on the outer shell of
the icosahedron with the Cu center and Zr outer shell, to discuss the stability changeby the ratio of element. All the icosahedrons show positive detCIJ in any Cu-Zr ratio,
especially the icosahedron which has only single Cu atoms on the outer shell and also
has single Cu-Cu bond between the Cu center shows drastic increase in detCIJ by one
order.

目次
1

緒論

1

2

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

3

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

3

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

4

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

5

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

6

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

7

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

12

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

13

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

14

2.8.1

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

15

2.8.2

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

18

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

20

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

20

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

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

26

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

26

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

28

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

29

3

i

目次

ii

単元系クラスターの安定性評価

31

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

31

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

37

4.2.1

孤立クラスターの安定性 . . . . . . . . . . . . . . . . . . . . . .

37

4.2.2

クラスター準結晶の安定性 . . . . . . . . . . . . . . . . . . . . .

40

4.2.3

原子間の結合距離 . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.2.4

結合中の電子密度 . . . . . . . . . . . . . . . . . . . . . . . . . .

46

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

48

中心置換二元系クラスターの安定性評価

49

5.1 解析方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

49

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

51

5.2.1

中心置換クラスターの安定性 . . . . . . . . . . . . . . . . . . .

51

5.2.2

主小行列式の値の変化 . . . . . . . . . . . . . . . . . . . . . . .

53

5.2.3

電子構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

55

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

61

元素割合を変えたクラスターの安定性評価

62

6.1 解析方法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

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

65

6.2.1

元素割合とクラスターの安定性 . . . . . . . . . . . . . . . . . .

65

6.2.2

電子構造 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

67

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

72

結論

73

参考文献

75

謝辞

80

4

5

6

7

第1章 緒 論
1960 年に発見されたアモルファス金属は高強度，高耐食性，超伝導性などの優れた物理的特性を持つため (1) ，当初新材料としての利用が期待されていた．しかしながら，溶融金属を急冷することで液体状態のランダム構造を凝固させて作製するため，バルク（塊）状のものを得るのは難しく，薄膜状のものに限られていた (2) ．一方，凝固過程において過冷却液体状態が安定化するために，比較的遅い冷却速度でもアモルファ
ス構造を保つ合金 (3) が発見された．これは従来のものと区別してバルク金属ガラスと呼ばれ，構造材料としての適用が期待されている (4) ．
アモルファス構造は，長周期的な構造を持たないランダム充填構造をしており，巨視的には均一である．しかし，微視的にはランダムではなく，X 線や中性子による構造解析 (5)(6)(7) により，特有の短距離・中距離秩序を持つことが知られており，局所的には不均一な構造をしている．また，RMC(逆モンテカルロ) シミュレーションにより金属ガラス中の局所構造を特定することが可能になっており (8) ，これらの局所構造がアモルファス構造の安定性にどのように寄与するのか，盛んに研究されている．才田らは，電子顕微鏡観察などにより，過冷却液体の安定性が異なる合金を比較し，安定性が高い合金にのみ内部に 20 面体クラスターが存在することを明らかにし，過冷却液体の安定化には 20 面体クラスターが特に強く関係していると指摘している (9)(10) ．
一方，これらの局所構造が安定であればあるほど，金属ガラスは共有結合的な性質となり，靭性が低下する．形成能と変形能という相反する性能を両立させた金属ガラ
スを試行錯誤的に見つけ出すことは，途方もない時間，労力そして費用を必要とするであろう．そこで，シミュレーションによる合金予測・設計が期待されている．特に第一原理計算は，密度汎関数理論に基づいて Kohn-Sham 方程式 (11) を解くことで，逐次電子状態を計算する手法であり，電子論に基づいて非経験的に材料特性を予測できる手法として注目されている．金属ガラス形成能を考える上で重要な溶解熱 (12) や，ガラス形成能のメカニズムを明らかにするべく，原子間ポテンシャルを用いた簡便な分1

子動力学シミュレーションでアモルファス構造を作成した後，局所構造を抽出してその電子構造を評価するといった検討 (13)(14) がされている．しかしながら，第一原理計
算では膨大な計算量を必要とするため，現時点ではアモルファス構造を直接扱うことは困難である．
本論文では，アモルファス構造に重要な 20 面体クラスターの安定性ならびに力学特性に焦点を当て，第一原理計算による格子不安定解析 (15)−(19) を行う．格子不安定解
析は，もともとは連続体の弾性限界を，内部の結晶格子のエネルギーバランスから求めるもので，Born の安定基準 (15) として知られているものである．第一原理計算では，fcc，bcc そして hcp 等の結晶構造でのエネルギー最小値の大小で安定相状態の議論が行われることが多いが，構造の壊れやすさとして，エネルギー曲線の 2 次勾配に相当
する格子不安定性の評価が必要と考える．計算量の問題から，現時点では周囲に原子が存在する中での 20 面体クラスターの解析を行うことは不可能であるが，結晶構造との違いを比較することで，様々な合金の 20 面体クラスターの力学特性に何らかの知見を与えるものと考える．
本論文は本章を含め 7 章で構成される．第 2 章では，第一原理分子動力学法の基礎理論，および，電子状態計算の高速化手法について説明する．第 3 章では，格子不安定解析の概要と，弾性係数を用いた格子不安定性評価について述べる．第 4 章では，Cu，

Ni，Al，Zr，Ti の単元素について，真空中の孤立クラスターとクラスター準結晶および完全結晶の解析を行い，単元系クラスターの基本的特性を調べる．先述のようにクラスター周囲の原子を考慮するのは困難であるので，孤立クラスターと周期的に並ぶ構造ではあるがまわりに他原子があるクラスター準結晶を比較し，その寄与を調べる．
第 5 章では，各単元系クラスターの中心原子を置換した二元系クラスターを解析し，合金化による安定性の変化を検討する．第 6 章では，元素割合を変えたときの安定性の変化を比較するために，実際に金属ガラスが得られている Zr70 Cu30 を想定し，Cu を中心原子として外殻原子の元素割合を変えたクラスターについて検討する．第 7 章では本研究で得られた結果を総括する．

2

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

2.1

断熱近似と平均場近似

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

3

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

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

2.2

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

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

Ts [ρ(r)] =

occ
∑
i

1
< ψi | − ∇2 |ψi >2と書くと，式 (2.1) は∫Etot [ρ] =v(r)ρ(r) dr + Ts [ρ]
4

(2.2)

1 ∫∫ ρ(r 0 )ρ(r) 0+
dr dr + Exc [ρ]2|r 0 − r|(2.3)

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

(2.4)

∑occ

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

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

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

(2.5)

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

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

(2.6)

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

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

(2.7)

i

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

2.3

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

∫






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

(2.8)

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

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

εxc (r) = εx + εc
0.4582εx (r) = −
rs

(2.9)
(2.10)

 −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
4µx (r) = εx3µc (r) =(2.13)(2.14)[
1

−0.1423√
1 + 1.0529 rs + 0.3334rs


(

)]

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

となる．

2.4

逆格子空間

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

R = n1 a1 + n2 a2 + n3 a3

6(n1 , n2 , n3 は整数)

(2.16)

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

(2.17)と定義される．これらのベクトル 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 の定理 (22) より

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

k=

h1
h2
h3
b1 + b2 + b3
n1
n2
n3

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

G7(2.23)

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

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

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

ρ (r) =

BZ
occ ∑
∑
n

=

fn fk |Ψkn (r)|2

k

occ ∑
BZ
∑∑∑
G

G0

n

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

(2.26)

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

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

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

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

=⇒

∑
G0

n
n
Hk+G,k+G0 Ck+G
0 = εkn Ck+G

(2.27)

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

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

を用いる．

8

(2.28)

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

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

(b) 局所項
局所擬ポテンシャルの行列要素は，PP< k + G |Vloc(r) | k + G0 >∫1
V PP (r) exp [−i(k + G) · r] exp [i(k + G0 ) · r] dr=Ω Ω loc
PP(G − G0 )(2.31)= Vlocである．結晶全体の局所擬ポテンシャルは格子周期関数であり，周期セル内の原子 aからの距離 r に対する局所擬ポテンシャル VaPP,loc (r) を用いて
PP
Vloc
(r) =

∑∑
R

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

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

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

∫

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

VaPP,loc (G) =

∫

= 2π

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(r) により，
ンシャル Va,l1 ∑exp[−i(G − G0 ) · r a ]VaPP,nloc (k + G, k + G0 )
Ωat aPP= Vnloc(k + G, k + G0 )(2.34)

PP< k + G|Vnloc(r)|k + G0 > =VaPP,loc (k + G, k + G0 )
= 4π

∫

∑

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

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

ρ(r) =∑

ρ(G) exp[iG · r]

(2.36)

ρ(r) exp[−iG · r]

(2.37)

G

ρ(G) =

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

ρ(G) exp[iG · r]

(2.38)

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

Vcoul (r) = 4π

∑ ρ(G)
G|G|210

exp[iG · r]

(2.39)

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

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

Vcoul (G) == 4π∑ ρ(G0 ) ∫
0 2

|G |

G0

= 4π

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

ρ(G)
|G|2

(2.40)

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

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

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

(2.42)µxc (r) exp [−iG · r] dr(2.43)

Gµxc (G) =1Ω∫となる．したがってハミルトニアン行列要素は (2.41) 式と同様に

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

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

11(2.44)

2.6

系のエネルギー

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

BZ ∑
occ
∑
k

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

1 ∑ ∑ a a0 ∑4π
− |G|20
EEwald =Z v Zv
)]
exp
exp
[iG
·
(r
−
r
a
a
2
2 a a0
4γ 2
G6=0 Ωat |G|+−]1 ∑ ∑ a a0 ∑ erfc (|R + r a0 − r a | γ)Z Z
2 a a0 v v R|R + r a0 − r a |∑ Zva2 γ

√

a

π

−

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

(2.46)である．
ここで

ρ(r) =∑

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

(2.47)

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

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

fk

occ
∑
n

fn

∑∑

G G0

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

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

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

BZ
occ
¯
¯
∑
∑
∑
1∑
¯ n ¯2PPfnfk¯ + Ωat
Etot =|k + G|2 ¯Ck+G
Vloc(G) ρ (−G)
2 k
n
G=0
G6=0

12

+

BZ
∑

fk

occ
∑

fn

G=0 G0 =0

n

k

∑ ∑

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

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

2.7

応力

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

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

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

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

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

=−
−

BZ ∑
occ ∑
1 ∑
n
|2 (k + G)α (k + G)β
fk fn |Ck+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∗
Ck+G
fk fn Sa (G − G0 )Ck+G
0Ωat k n G G0 l a∂
1 PP,nloc
×(k + G, k + G0 )}{V∂εαβ Ωat a,l+13(2.51)(2.52)

+

1∑
2Gα Gβ− δαβ }
Vcoul (G)ρ(−G){
2 G
|G|2

+ δαβ

∑

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

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

2.8

擬ポテンシャル法

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

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

14

2.8.1TM 型擬ポテンシャル

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

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

]

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

(2.54)

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

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

ψlPP (r) =ψlAE (r)(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) が次式で表される．
PP
Vscr,l
(r) =


AE

 Vl (r)

0

00

2

l+1 0
p (r) + [p (r)] εl +
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

15

[∫ r
0

cl

|ψlAE (r)|2 r2](2.59)

(b) 式 (2.57) の 2 次微分までが rcl で連続である条件[]

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

(2.60)(2.61)

2(l + 1) 0
2
p (rcl ) − [p0 (rcl )]rcl
2(l + 1) 0
p (rcl )p000 (rcl ) = 2VlAE0 (rcl ) +rcl2
2(l + 1) 00−
p (rcl ) − 2p0 (rcl )p00 (rcl )rcl
4(l + 1) 0p0000 (rcl ) = 2VlAE00 (rcl ) −p (rcl )rcl2
2(l + 1) 0004(l + 1) 00
p (rcl ) −+p (rcl )2
rcl
rcl
p00 (rcl ) = 2VlAE (r) − 2εl −−2 [p00 (rcl )] − 2p0 (rcl )p000 (rcl )2(2.62)(2.63)(2.64)

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

PP
PP
(c) Vscr,l(r) の r = 0 における曲率が 0 である条件 (Vscr,l(r) = 0)c22 + 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
Vion,l(r) = Vion,loc
(r) +

∑
l

16

PP
Vnloc,l(r)P̂l(2.67)

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

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

これを用いると，行列要素の非局所項は，
KB< k + G|Vnloc,l(r)|k + G0 > =∑ (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 + G0 |r)r2 dr}∗
Ylm (k + G)Ylm(k + G0 )(2.69)m=−lとなる．ここで，PP
Cl =< ψlPP (r)|Vnloc,l(r)|ψlPP (r) >PP
Vnloc,l
(r) =

∑∑
R

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

m=−l∗(k + G0 )}× {exp[iG0 · r a ]Ala (k + G0 )Ylm

17

(2.74)

と書ける．平面波展開係数との積は，
∑
G0

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

m=−l


∑

×

G0




n
∗
0
0
0
Ck+G
(2.75)
0 exp[iG · r a ]Ala (k + G )Ylm (k + G )

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

CAYnalkm (k + G) =

∑

n
∗
0
0
0
Ck+G
0 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 とすると，擬波動関数の満たすべきシュレーディンガー方程式は以下のようになる．
0− εi )|Ψi >= 0 ,(T + Vloc + VNL

18

0
=
VNL

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

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

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

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

∑

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

ij

となる．

19

(2.88)(2.89)

2.9

電子占有数

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

[

(

εi − εF1
1 − erf
2
σ

fi =

)]

(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)n
Ck+G
exp[iG · r]

(2.95)

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

(2.96)

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

ukn (r) =∑Gとなる．同様にu∗kn (r) =

∑
G

20

であるから，ukn (r)u∗kn (r) =

∑∑
G

G0

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

(2.97)

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

ρ(r) =

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

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

面波基底の波動関数を用いた場合には，系の全エネルギーを最小にする係数ベクトル

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

∑

λmn

(
∑

mn

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

(2.99)

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

∑∑
G

G0

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

Rnk = −

0
∂Etot
m∗
∂Ck+G

]


∑
= −  (Hk+G,k+G0 − λnn ) C n

k+G0

G0


,

(G=G1 ,G2 ,···,Gmax )(2.101)以下に金属の電子状態計算において代表的な BKL 法 (34) について解説する．

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

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

∂ε0kn
m∗
∂Ck+G

]i

(

)

i
= − H − λin I · C nk

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

]

¯
¯

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

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

∑ (

)

i∗
i
C mk
· R ink C mk

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

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

)

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

∑ (
m6=n

23

)

i∗
i
i
C mk
· Gnk
C mk

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

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

i
i
i−1
F nk
= Gnk
+ γ i F nk

i0∗

Gnk
· R ink0
i−1 0 ∗
i−1 0 (i > 1)
γ i =  Gnk· R nk 0(i = 1)(2.111)(2.112)

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

)

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

i0
F nk

(2.113)(2.114)

)1

i0∗
i0 2
F nk
· F nk< 新たな係数ベクトルの組み立て >新たな係数ベクトルの組立ては次のように行われる．
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
∗
ε12
ε22

]

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

} 12

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

α= {

ε12

}1

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

24

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

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

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

25

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

3.1

不安定条件

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

26

a3
3
4

2
1

6

5

a2
1

Fig.3.1 Unit cell of fcc lattice

ネルギーの Taylor 級数展開は
6
∑

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

∂U ¯¯
,
¯
Fm =
∂am ¯{am }¯∂ 2 U ¯¯Amn =¯∂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 ] の正値性に帰着される．

27

3.2

応力と弾性係数

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

dU = T dS − dW

(3.4)

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

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

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

(3.11)

したがって，断熱過程では

dU = V (X)σij dηij

(3.12)となり，基準配置における応力テンソルと弾性係数は，

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

(

(

∂U
∂ηij

)

(3.13)
η0∂2U∂ηij ∂ηkl)(3.14)
η0となる．ここで，η 0 は η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 対称性 (19) と呼ばれる次の対称性を持つ．

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 が保持されたまま変形すると仮定した場合の外部負荷によってなされる仕事であり，無負荷平衡点での不安定性は弾性係数
マトリクス Cijkl の正値性に帰着される．

Cijkl は式 (3.16) の対称性を持つので，i,j の 6 つの組み合わせ (i,j)=(1,1)，(2,2)，(3,3)，(2,3)，(3,1)，(1,2) で Cijkl を CIJ で表すことができる．(Voigt 表記) したがって，無負荷平衡点での系の安定性は次の 6 × 6 行列の正値性により評価される (16) ．

29

¯
¯ C
¯ 11
¯
¯ C12
¯
¯
¯ C13
|CIJ | = ¯¯
¯ 0
¯
¯ 0
¯
¯
¯ 0

¯

C12 C13 00
0 ¯¯¯
C22 C23 00
0 ¯¯¯
C23 C33 000 ¯¯00 C44 00 ¯¯¯
000 C55 0 ¯¯¯
000
0 C66 ¯

(3.18)|CIJ | <0 となる条件は，以下の 5 つのいずれかとなる．
¯
¯ C
¯ 11
¯
¯ C
¯ 12
¯
¯ C13

¯

C12 C13 ¯¯¯
C22 C23 ¯¯ < 0¯C23 C33 ¯

¯
¯ C
¯ 11
¯
¯ C12

(3.19)¯

C12 ¯¯¯<0C22 ¯(3.20)

C44 < 0(3.21)

C55 < 0(3.22)

C66 < 0(3.23)

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

30

第 4 章 単元系クラスターの安定性評価
本章では，単元素からなる 20 面体クラスターの構造安定性を検討することを目的として，Cu，Ni，Al，Zr，Ti それぞれについて，真空中孤立クラスター，クラスター準結晶，結晶単位格子の第一原理解析を行い，格子安定性および原子間の結合距離と電子密度について検討する．孤立クラスターと周囲に原子が存在しているクラスター準結晶を比較することでその差を検討するとともに，これらの単元系クラスターと金属単結晶を想定して解析した単位格子を比較することでクラスター構造と結晶状態の違いを議論する．

4.1

解析方法

本研究での解析は全て Kresse らにより開発された平面波基底ウルトラソフト擬ポ
テンシャル法 (30) に基づく第一原理バンド計算コード VASP(40) (Vienna Ab-initio Sim-ulation Package) を用いて行った．交換相関項には局所密度近似 (21) （Local Density
Approximation, LDA）に勾配を考慮した一般化密度勾配近似 (41)（Generalized Gradient Approximation, GGA）を用いた．また収束計算には残差最小化手法 (42) (Residual
Minimization Method – Direct Inversion in the Iterative Subspace, RMM–DIIS) を採用した．カットオフエネルギー，バンド数，波動関数と電子密度の FFT メッシュは原子種，原子数と電子数から VASP が定める値 (43) を用いた．
解析を行った単元系 20 面体クラスターを図 4.1(a) に示す．構成元素は Cu，Ni，Al，Zr，Ti であり，クラスター内での原子間距離はそれぞれの結晶状態での最近接原子間距離 (2.56 Å，2.49 Å，3.17 Å，2.86 Å，2.89 Å)(22) とし，クラスターを真空中の立方体中心に図 4.1(b) のように配置して解析セルとした．平面波基底のバンド計算では全31方向への周期性が要請されるので，セル辺長については，図 4.2(a) のように原子間距離の 4 倍と十分大きくした真空中孤立クラスターと，図 4.2(b) のように 3 倍として周囲クラスターとの距離を小さくしたクラスター準結晶の 2 通りを考えた．クラスター
準結晶は隣接クラスターとの最近接距離が原子間距離と等しく，クラスター間に結合が生じうる距離である．周囲に他の原子が存在するアモルファス中の 20 面体クラス
ターについて，直接第一原理計算を行うことは現時点では困難なため，孤立させたク
ラスターの解析と，自分自身が周期的に並ぶ構造ではあるが周囲に原子が存在するクラスターの解析を行うことでその差を議論することを目的としている．逆格子空間のk 点メッシュは Monkhorst-Pack 法 (44) に従い，33 の格子点をとっている．表 4.1 にそれぞれの計算条件をまとめて示す．
また，アモルファス金属中でみられるクラスター構造と，単結晶金属中の電子構造の違いを議論するために，単位格子による完全結晶の解析も行った．fcc の Cu，Ni，Alについては図 4.3(a) に示すような単位格子を解析セルとし，hcp の Zr，Ti は図 4.3(b)の着色部分を解析セルとする．完全結晶の解析条件を表 4.2 に示す．

VASP による電子状態及び原子構造の緩和計算によって，無負荷平衡状態の原子配置を求めた後，図 4.4 に模式的に示すようにセルに ±1.0% の引張ひずみ，もしくはせん断ひずみを与え，原子構造緩和をせずに電子状態の収束計算を行った．このとき，クラスターおよび単位格子はセルの変形に対応したアフィン変形を与えている．各ひず
み摂動による応力変化から弾性係数 CIJ を数値的に評価し，その行列式 detCIJ の正値性よりそれぞれの安定性を評価した．

32

Length of cell edgezyx(a)Monatomic icosahedral structure.

(b)Simulation cell.

Fig.4.1 Simulation model.

Table 4.1 Calculation conditions for monatomic icosahedrons.
Number of atoms
Number of electrons
Distance between atoms (Å)
Length of cell edge (Å)
Number of bands
Cutoﬀ energy (eV)Number of k points

Cu13
13
143
2.56
10.24, 7.68
98
292.2
33

Ni13
13
130
2.49
9.96, 7.47
78
302.2
33

Number of atoms
Number of electrons
Distance between atoms (Å)
Length of cell edge (Å)
Number of bands
Cutoﬀ energy (eV)Number of k points

Zr13
13
52
3.17
12.68, 9.51
52
187.4
33

Ti13
13
52
2.89
11.56, 8.67
52
226.5
33

33

Al13
13
39
2.86
11.44, 8.58
46
161.5
33

(a)Isolated icosahedrons.

Possibility of bond(b)Quasi-crystal of icosahedrons.
Fig.4.2 Isolated icosahedron and quasi-crystal of icosahedrons by periodic
boundary condition.

34

z
y
x

(a)fcc (Cu, Ni, Al).

(b)hcp (Zr, Ti).

Fig.4.3 Unit lattices of fcc and hcp structures.

Table 4.2 Calculation conditions for perfect lattices.

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

Cu
fcc
4
44
30
292.2
153

Ni
fcc
4
40
28
302.2
153

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

Zr
hcp
2
8
8
187.4
153

Ti
hcp
2
8
8
226.5
153

35

Al
fcc
4
12
14
161.5
153

∆η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.4 Strain perturbation on supercell for evaluation of elastic coeﬃcients.

36

4.2

解析結果と考察4.2.1

孤立クラスターの安定性

孤立クラスターおよび完全結晶における弾性係数の行列式 detCIJ の値を図 4.5 に示
す．上段に孤立クラスターの値を，下段に完全結晶の値を示している．正値を右側，負値を左側に示し，絶対値の大きさを対数として棒グラフで表している．孤立クラスターの detCIJ は，Cu，Ni が正値となって安定であるが，Al，Zr，Ti は負値となりいずれも不安定となる．完全結晶における detCIJ はいずれの元素も正値であり，安定と判定された Cu，Ni クラスターの値と比較していずれも桁違いに大きい．安定性解析においては，行列式の値の大小関係の物理的意味はなく，正値性のみが意味を有するが，原子間ポテンシャルを用いたシミュレーション (45) では合金系における行列式の値の大小関係が引張変形時の最大応力に対応していることが報告されており，行列式の大小関係と機械的特性に関係がある可能性も示唆されている．Cu，Ni 孤立クラスターの

detCIJ の値を比較すると，Ni が Cu よりも 1 桁程度大きく，この傾向は完全結晶でも同じとなっている．
表 4.3 に孤立クラスターと完全結晶における弾性係数の各成分 C11 ，C22 ，C33 ，C12 ，

C23 ，C31 ，C44 ，C55 ，C66 の値を示す．孤立クラスターおよび hcp では対称性がないので全ての成分を表示しているが，fcc 完全結晶においては C11 =C22 =C33 ，C23 =C31 =C12 ，

C44 =C55 =C66 なのでまとめて 1 つの値で示している．孤立クラスターのいずれの成分も完全結晶に比べて大きく減少しており，変形抵抗が極めて小さいことがわかる．

detCIJ が正値であった Cu，Ni の孤立クラスターはいずれの成分も正値となっており，
ひずみ摂動に対してクラスターは変形抵抗を有している．一方，detCIJ が負となった

Al，Zr，Ti クラスターは負の弾性係数が垂直成分 C11 〜C33 およびせん断成分 C44 〜
C66 に生じており，ひずみ摂動でクラスターが崩壊する「不安定な」構造であること
がわかる．

37

Absolute value of determinant of elastic coefficients, detCIJ (GPa6).
109

106

103

100 100

103

Isolated icosahedron

StableUnstable

Cu
Ni
Al
Zr
Ti

106

109

1012

Positive value

Cu13
Ni13
Al13
Zr13
Ti13

Negative value

1012

Perfect latticeFig.4.5 Determinant of elastic coeﬃcients of isolated monatomic icosahedrons
and monatomic perfect lattices.

38

Table 4.3 Components of CIJ (GPa) of isolated monatomic icosahedrons andmonatomic perfect lattices.

Cu13

C11
15.7

Cu
Ni13

C33
14.4

C12
11.3

C13
11.0

187.5
21.2

Ni
Al13

C22
15.3

21.0

Al

1.0

C44
1.7

122.7
22.2

17.1

17.3

259.7
0.9

C23
11.2

10.8

11.8

128.6

C66
2.2

90.7
17.1

2.5

170.5
−0.2

C55
1.7

2.2

2.7

102.1
11.7

−5.0

46.4

−5.0

−5.0

53.9

Zr13

−14.8

−14.7

−14.0

24.0

23.2

23.2

−18.9

−18.9

−18.9

Zr

152.8

152.7

169.6

56.4

70.7

70.1

18.3

18.0

48.3

Ti13

−22.1

−21.6

−19.6

26.7

25.6

25.3

−24.2

−23.0

−23.5

Ti

181.1

185.1

166.3

53.2

67.8

73.3

37.9

38.9

64.8

39

4.2.2クラスター準結晶の安定性

クラスター準結晶の弾性係数の行列式 detCIJ の値を，前節で示した孤立クラスターのそれとあわせて図 4.6 に示す．Zr 準結晶は detCIJ ＞ 0 と正値を示したが，安定性判
別式 (3.19)〜(3.23) の 2 個が負となった結果，負×負で正となっているので不安定として示した．Zr 以外ではクラスター準結晶の detCIJ の正負は孤立クラスターにおける傾向と同じである．安定クラスター Cu，Ni の detCIJ の値の大小関係も Ni ＞ Cu で一致している．Cu，Ni 準結晶クラスターの detCIJ の値は，孤立クラスターのそれと完全結晶の値の中間程度となった．
表 4.4 に (a) 孤立クラスターと (b) クラスター準結晶の弾性係数の各成分 CIJ の値を示す．周囲クラスターとの最近接原子間距離が x, y, z 方向で同一ではないので，準結晶では異方性がより強くなっているが，Cu，Ni ではほとんどが孤立クラスターよりも変形抵抗が増加している．不安定な Al，Zr，Ti クラスターは，準結晶とすることで垂直成分 C11 〜C33 およびせん断成分 C44 〜C66 が大幅に増加し，正値となっている．しかしながら，いずれのクラスター準結晶も負の成分が生じており，安定に存在しえな
い構造であることがうかがえる．

40

Absolute value of determinant of elastic coefficients, detCIJ (GPa6).
109

106

103

100 100

103

Isolated icosahedron

StableUnstable

Cu13
Ni13
Al13
Zr13
Ti13

106

109

1012

Positive value

Cu13
Ni13
Al13
Zr13
Ti13

Negative value

10

12

Quasi-crystal of icosahedrons
Fig.4.6 Determinant of elastic coeﬃcients of isolated monatomic icosahedrons
and monatomic quasi-crystals of icosahedrons.

41

Table 4.4 Components of CIJ (GPa). (a)Isolated monatomic icosahedrons,(b)Monatomic quasi-crystals of icosahedrons.

(a)

C11
15.7

C22
15.3

C33
14.4

C12
11.3

C13
11.0

C23
11.2

C44
1.7

C55
1.7

C66
2.2

(b)

29.0

27.6

30.2

12.0

18.1

18.9

11.0

9.3

8.0

(a)

21.2

21.0

22.2

17.1

17.3

17.1

2.5

2.2

2.7

(b)

36.4

52.8

63.8

20.0

16.0

27.6

20.6

11.3

19.9

(a)

0.9

1.0

−0.2

10.8

11.8

11.7

−5.0

−5.0

−5.0

(b)

16.0

18.5

41.0

21.5

−1.88

−1.35

16.0

17.4

17.0

(a)

−14.8

−14.7

−14.0

24.0

23.2

23.2

−18.9

−18.9

−18.9

(b)

21.1

15.1

30.2

23.9

13.1

15.3

4.0

2.7

−6.6

(a)

−22.1

−21.6

−19.6

26.7

25.6

25.3

−24.2

−23.0

−23.5

(b)

25.5

33.1

26.3

12.7

20.0

11.4

5.1

−1.1

7.65

Cu

Ni

Al

Zr

Ti

42

4.2.3原子間の結合距離

図 4.7 に孤立クラスター構造内の結合距離平均を示す．各元素ごとに中心原子と外殻原子間での結合と外殻原子どうしの間での結合の距離をそれぞれ示しており，完全結晶の第一近接の原子間結合長に対する比率をグラフの横に併記した．
孤立クラスターの中心原子−外殻原子間の結合距離は完全結晶の原子間距離よりも小さくなっており，孤立クラスターは完全結晶よりも原子が密に詰まった状態である．
また，外殻原子どうしの間の結合距離はいずれの元素でもこれより増加し，完全結晶の結合長に近い．detCIJ が正であった Cu，Ni の結合距離は中心原子−外殻原子間で結晶状態結合長の 94 ％程度，外殻原子どうしの間で 99 ％程度である．負の CIJ 成分
個数が Zr，Ti に比べて少ない Al のクラスターにおいても同等の値となっている．これらのクラスターは図 4.8(a) に示した等電子密度面からわかるように，2 原子間での直線的な結合を形成している．一方，多くの CIJ 成分が負となっていた Zr，Ti の結合距離は，完全結晶での結合長に対する比率が Cu，Ni，Al クラスターと比べて小さくなっている．その原子間結合は，図 4.8(b) に示すように外殻原子 3 つと中心原子の四面体の中点に価電子を有する結合となっている．
図 4.9 にクラスター準結晶内の結合距離平均を示す．図 4.7 と同様に，2 種類の結合の結合距離と完全結晶の結合長に対する比率を示している．元素ごとに 2 種類の結合の値を比較すると，クラスター準結晶でも孤立クラスターと同様に外殻原子どうしの間の結合のほうが大きくなっている．また，Zr，Ti クラスターより Cu，Ni，Al クラスターのほうが結晶状態結合長に対する比率が大きくなっているのも孤立クラスターと同じ傾向である．孤立クラスターとクラスター準結晶の結合距離を比較すると，クラスター準結晶のほうがいずれの元素でも大きくなっている．完全結晶の結合長に対する比率も増加しており，結合距離が結晶状態に近づいているが，Cu，Ni，Al では完
全結晶の結合長よりも大きくなっていることがわかる．

43

Distance, (Έ).
0

0.5

1

1.5

2

2.5

3

Percentage against
nearest neighbor distance
3.5 in perfect crystal.

Cu

Center - outer
Outer - outer

94.5
99.3

Ni

Center - outer
Outer - outer

93.6
98.4

Al

Center - outer
Outer - outer

93.4
98.2

Zr

Center - outer
Outer - outer

89.3
94.0

Ti

Center - outer
Outer - outer

89.6
94.1

Fig.4.7 Bond lengths in isolated icosahedrons, between the center and outer
atoms and between two outer atoms. The nearest neighbor bonds areshown for perfect lattices.

(a)Cu, Ni, Al cluster(b)Zr, Ti cluster

Fig.4.8 Equivalent surface for valence electrons.

44Distance, (Έ).
0

0.5

1

1.5

2

2.5

3

3.5

Percentage against
nearest neighbor distancein perfect crystal.

Cu

Center - outer
Outer - outer

96.4
101.4

Ni

Center - outer
Outer - outer

97.9
103.2

Al

Center - outer
Outer - outer

95.2
100.2

Zr

Center - outer
Outer - outer

91.1
95.8

Ti

Center - outer
Outer - outer

91.4
96.1

Fig.4.9 Bond lengths in quasi-crystals of icosahedrons, between the centerand outer atoms and between two outer atoms. The nearest neighborbonds are shown for perfect lattices.

45

4.2.4結合中の電子密度

図 4.10 に結合中の最大の電子密度の値を示す．図 4.10(a) に中心原子と外殻原子の間，図 4.10(b) に外殻原子どうしの間，図 4.10(c) に周期境界条件による隣接クラスターとの間の結合中電子密度をそれぞれ表している．図中では□で孤立クラスターを，△でクラスター準結晶を示しており，×は完全結晶での第一近接の原子間結合における値である．
図 4.10(a), (b) に示す孤立クラスター内の結合中の電子密度 (□) は，ほとんどが結晶での電子密度 (×) に近い値となっているが，安定な Cu，Ni の中心原子−外殻原子間結合では結晶より高い値をとっており，結合力が強いことが示唆される．ただし，隣接クラスターとの間の電子密度がほとんど 0 である部分が，ひずみ摂動を与えると容易に変形するため，完全結晶よりも変形抵抗が小さくなっていると考える．
一方，クラスター準結晶の結合中の電子密度 (△) は，中心原子−外殻原子間，2 外殻原子間の両方において孤立クラスター (□)，完全結晶 (×) より低い値となっている．
これは，図 4.10(c) に示したように隣接クラスター間に電子密度が生じたため，クラス
ター中の電子密度が減少したものである．このように，クラスター準結晶の個々のボンド結合力は弱まっているが，準結晶としての構造により剛性は孤立クラスターより

C11 〜C33 の垂直成分，C44 〜C66 のせん断成分が上昇する．ただし，対称性がないため
C12 ，C23 ，
C31 の各成分が低下する．

46

Electron density, e/Έ3

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

Cu

Ni

Al

Zr

0

Ti

(a)center - outer bond

Cu

Ni

Al

Zr

Ti

(b)outer - outer bondElectron density, e/Έ3

0.4

0.3

Isolated Quasi-crystal
center-outer
outer-outercluster-cluster

0.2

Perfect lattice

0.1

0

Cu

Ni

Al

Zr

Ti

(c)cluster - cluster
Fig.4.10 Maximum valence electron density between atoms and clusters.

47

4.3

結言

単元素からなる 20 面体クラスターの構造安定性を検討するために，Cu，Ni，Al，Zr，

Ti の真空中孤立クラスター，クラスター準結晶，完全結晶についての第一原理解析を行い，弾性係数の行列式 detCIJ の値，原子間の結合距離，結合中の電子密度について検討した．得られた結果をまとめて以下に示す．

1. 孤立クラスター，クラスター準結晶いずれにおいても Cu，Ni は弾性係数 CIJ の行列式が正値をとり安定であった．完全結晶ではいずれの元素も正値をとるが，Al，Zr，Ti クラスターは負値となり不安定であることが示された．
2. Cu，Ni の detCIJ の大小関係は，孤立クラスター＜クラスター準結晶＜完全結晶でいずれも Cu ＜ Ni となることが示された．

3. 安定な Cu，Ni は孤立クラスター，クラスター準結晶いずれにおいても CIJ の成分が全て正である．一方，不安定な Al，Zr，Ti クラスターでは CIJ に負の成分が存在した．

4. クラスター内の原子間結合距離を調べた結果，安定／不安定，孤立／準結晶クラスターのいずれにおいても，外殻原子間の距離よりも，中心原子−外殻原子の距離が短い．また，不安定なクラスターは完全結晶の最近接原子間距離からの差が大きい．

5. 結合の電子密度分布から，Cu，Ni，Al クラスターは 2 原子間で直線的な結合をしているのに対し，Zr，Ti クラスターは 3 つの外殻原子と中心原子の四面体中
心に電子を有する結合を示した．

6. detCIJ ＞ 0 である Cu，Ni 孤立クラスターは中心原子と外殻原子の間の結合中に電子を多く局在させており，結合が強くなっていることを示唆した．

7. クラスター準結晶では個々のボンド結合力は孤立クラスターより弱まっているが，準結晶の構造により孤立クラスターより剛性が上昇することを示唆した．

48

第 5 章 中心置換二元系クラスターの安定性評価
前章では，単元素からなる孤立クラスター，クラスター準結晶，完全結晶の違いを検討した．一方，実験的に発見されているアモルファス構造は全て二元素以上の多元系であり (46) ，20 面体クラスターもアモルファス構造中では通常多元素で構成されているものと考えられる．
本章では，前章で検討した単元系クラスターの中心原子を他元素で置換した二元系
孤立クラスターについて，格子安定性と電子密度分布を調べる．前章で検討した単元系孤立クラスターと中心原子置換クラスターを比較し，合金化による安定性の変化を
原子間の電子構造から詳細に議論する．

5.1

解析方法

前章で解析した，Cu，Ni，Al，Zr，Ti の 5 つの元素の組み合わせで中心原子のみ置
換した二元系クラスターを解析対象とする．解析条件を表 5.1 にまとめて示す．逆格子空間の k 点メッシュは全て 33 の格子点をとっている．初期構造作成時のクラスター
内での原子間距離は，構成原子 2 種類のうち原子半径が大きいほうの最近接原子間距
離の値を用いた．前章と同様に，辺長が原子間距離の 4 倍の立方体セル中心に配置した孤立クラスターについて，原子構造緩和，摂動計算を行って，弾性係数 CIJ を求め，その行列式 detCIJ の値からクラスターの格子安定性を検討した．

49

Table 5.1 Calculation conditions for A1 B12 binary icosahedrons with the center
atom A.
Number of electrons
Distance between atoms (Å)
Length of the cell edge(Å)
Number of bands
Cutoﬀ energy (eV)Ni1 Cu12
142
2.56
10.24
97
302.2

Al1 Cu12
136
3.17
11.44
94
292.2

Zr1 Cu12
135
2.86
12.68
94
292.2

Ti1 Cu12
136
2.89
11.56
94
292.2

Number of electrons
Distance between atoms (Å)
Length of the cell edge(Å)
Number of bands
Cutoﬀ energy (eV)Cu1 Ni12
131
2.56
10.24
92
302.2

Al1 Ni12
123
2.86
11.44
88
302.2

Zr1 Ni12
124
3.17
12.68
88
302.2

Ti1 Ni12
124
2.89
11.56
88
302.2

Number of electrons
Distance between atoms (Å)
Length of the cell edge(Å)
Number of bands
Cutoﬀ energy (eV)Cu1 Al12
47
2.86
11.44
50
292.2

Ni1 Al12
46
2.86
11.44
49
302.2

Zr1 Al12
40
3.17
12.68
46
187.4

Ti1 Al12
40
2.89
11.56
46
226.5

Number of electrons
Distance between atoms (Å)
Length of the cell edge(Å)
Number of bands
Cutoﬀ energy (eV)Cu1 Zr12
59
3.17
12.68
56
292.2

Ni1 Zr12
68
3.17
12.68
60
302.2

Al1 Zr12
51
3.17
12.68
52
187.4

Ti1 Zr12
52
3.17
12.68
52
226.5

Number of electrons
Distance between atoms (Å)
Length of the cell edge(Å)
Number of bands
Cutoﬀ energy (eV)

Cu1 Ti12
59
2.89
11.56
56
292.2

Ni1 Ti12
58
2.89
11.56
55
302.2

Al1 Ti12
51
2.89
11.56
52
226.5

Zr1 Ti12
52
3.17
12.68
52
226.5

50

5.2

解析結果と考察

5.2.1

中心置換クラスターの安定性

弾性係数の行列式 detCIJ の絶対値を，図 4.5 と同様に正負で左右にわけて図 5.1 に示す．値は上から外殻原子ごとにまとめて示しており，それぞれの一番下に前章で解析した単元系孤立クラスターの値を示している．Zr1 Cu12 ，Cu1 Al12 は，前章で Zr 準結晶について述べたように，偶数個の不安定条件が重なって負×負で detCIJ ＞ 0 となっているので不安定として示した．また，Zr1 Ni12 の値は −1.24 × 10−4 GPa6 とほとんど0 であるため棒グラフではなく図中に値を記した．なお，Zr1 Cu12 と Zr1 Ni12 は原子構
造緩和を行うと，クラスターの膨張がとどまることなく進行し，無負荷状態を求めることができなかったので，原子構造緩和なしの摂動計算で求めた detCIJ の値を示している．これらはいずれも不安定と判定されるが，外力なしに膨張し続けたことからこの判定結果は妥当である．
前章の解析で detCIJ が正値をとった Cu，Ni 単元系クラスターの中心原子を他元素で置換したクラスターは，中心原子を Zr とした系以外でいずれも detCIJ が正値をとった．一方，不安定な Zr，Ti 単元系クラスターの中心原子を Cu，Ni 原子で置換すると，detCIJ が正となっている．また，Al 原子も Ti 単元系クラスターを安定化しているが，外殻原子が Al であるクラスターは全て不安定となっている．安定なクラスターの detCIJ の大小関係をみると，Cu，Ni を外殻原子とするクラスターは，Al1 Cu12 を除いていずれも中心原子を他の元素に置換することで単元系より高い値となっている．
また，Cu，Ni を中心とし，Zr，Ti を外殻とするクラスターも，Cu，Ni 単元系より高くなっている．

51

Absolute value of determinant of elastic coefficients, detCIJ (GPa6).
108

106

StableUnstable

104

102

100

100

102

104

106

108

Ni1Cu12
Al1Cu12Zr1Cu12
Ti1Cu12
Cu13

Cu1Al12
Ni1Al12Zr1Al12
Ti1Al12Al13

Cu1Zr12Ni1Zr12Al1Zr12Ti1Zr12
Zr13

Cu1Ti12Ni1Ti12
Al1Ti12
Zr1Ti12
Ti13

Fig.5.1 Determinant of elastic coeﬃcients of A1 B12 binary icosahedrons.

52

Positive value

Negative value

Cu1Ni12
Al1Ni12
㧙1.2410-4 (unstable) Zr1Ni12
Ti1Ni12Ni13

5.2.2

主小行列式の値の変化

合金化による安定化について，より詳しく検討するために，detCIJ の主小行列式の値 (式 (3.19)〜(3.23) の左辺) を表 5.2 に示す．先述のようにこれらの条件式の値のいずれかが負となると，系は不安定と判定されるが，その条件によって不安定モードが異なる．(a) の Spinodal 条件はダイラタンシー (膨張) に対するものであり，先述の膨張する Zr1 Cu12 ，Zr1 Ni12 が実際にこの条件が負となっていることがわかる．Born 条件は引っ張りに対して横方向の Poisson 収縮の等方性に関するもので，C44 〜C66 はせん断変形に対する抵抗である．
単元系で不安定であった Al，Zr，Ti はいずれも Spinodal 条件以外負となっている．
これらのクラスターは，中心原子の置換によって Born 条件とせん断条件 3 方向の値が単元系から増加しているが，Al が外殻原子であるクラスターは合金化してもせん断条
件が全て負のため不安定のままである．一方，Zr，Ti 単元系の中心原子を Cu，Ni とすると，Born 条件，せん断条件 3 方向が大きく増加して正となり，安定となっている．

53

Table 5.2 Minor determinants of detCIJ . (a)Spinodal stability, (b)Born stability, (c)C44 stability, (d)C55 stability, (e)C66 stability.

Ni1 Cu12
Al1 Cu12
Zr1 Cu12
Ti1 Cu12
Cu13
Cu1 Ni12
Al1 Ni12
Zr1 Ni12
Ti1 Ni12
Ni13
Cu1 Al12
Ni1 Al12
Zr1 Al12
Ti1 Al12
Al13
Cu1 Zr12
Ni1 Zr12
Al1 Zr12
Ti1 Zr12
Zr13
Cu1 Ti12
Ni1 Ti12
Al1 Ti12
Zr1 Ti12
Ti13(a)Spinodal(GPa3 )
1.82E+03
3.15E+02−2.30E+01
7.27E+02
5.81E+02
3.39E+03
2.75E+03−1.24E-01
1.64E+03
1.03E+03−4.94E+00
3.53E+02
2.71E+02
1.44E+02
2.73E+03
1.28E+03
2.02E+03
6.81E+01
7.46E+03
4.67E+04
2.40E+03
1.78E+03
9.00E+02
7.09E+03
6.73E+04(b)Born(GPa2 )
1.80E+02
7.16E+01
8.20E+00
9.97E+01
1.11E+02
2.66E+02
2.29E+02−5.73E-01
1.72E+02
1.52E+02−5.39E+00−6.98E+01
4.94E+01−6.41E+01−1.16E+02
1.56E+02
1.90E+02−3.77E+01−2.31E+02−3.60E+02
2.10E+02
1.83E+02
1.47E+02−1.54E+02−2.39E+0254(c)C44
(GPa)
3.6
1.5
−1.6
2.8
1.7
3.4
4.6
0.1
3.6
2.5
−1.2
−2.7
−2.1
−1.9
−5.0
3.9
4.4
−0.7
−9.7
−18.9
4.5
5.2
3.5
−9.6
−24.2

(d)C55
(GPa)
4.1
1.6
−1.6
2.7
1.7
3.6
4.6
0.1
3.7
2.2
−1.2
−2.6
−2.2
−1.8
−4.9
3.9
4.4
−1.0
−10.0
−18.9
4.4
3.6
0.8
−9.6
−23.0

(e)C66
(GPa)
4.5
1.8
−1.6
2.8
2.2
3.8
4.6
0.1
3.6
2.7
−0.5
−2.9
−2.2
−3.6
−4.9
3.9
4.4
−1.4
−9.9
−18.9
4.4
3.1
3.4
−9.8
−23.5

5.2.3電子構造

図 5.2 に Cu を外殻原子とするクラスターの断面の電子密度分布を示す．等高線間隔
は 0.02 e/Å 3 としている．Cu 単元系より大きな detCIJ を示した Ni1 Cu12 ，Ti1 Cu12 は，Cu 単元系のそれと似ているが，電子密度の拡がりが小さくなっており，クラスター内部に電子が局在化している．Ti1 Cu12 は特にその傾向が顕著であるが，Ti と Cu 間の密度は低くなっているため，Ni1 Cu12 より detCIJ は小さくなったものと考える．Al1 Cu12も電子密度の拡がりは Cu13 に比べて小さくなっているが，中心原子と外殻原子の間の結合が切れており，このために単元系よりも変形抵抗が減少したと考えられる．また，中心原子を Zr で置換すると，全ての原子間の電子密度が著しく低くなっており，原子間の結合が切れて不安定となったものと考えられる．
図 5.3 は Ni を外殻原子とするクラスターの電子密度分布である．やはり，中心原子を置換すると，電子雲の拡がりが小さくなり，クラスター内の電子の局在化を生じている．Al で置換したときに Al への電子の局在がなくなっているのは図 5.2 の Al1 Cu12と同じである．しかし Al1 Ni12 の detCIJ は Ni 単元系よりは大きい値を示した．図 5.2と図 5.3 の Al で置換したクラスターでの分布を比べると，Ni まわりに局在した電子は円ではなく中心方向にわずかに配向している．不安定な Zr1 Ni12 クラスターはやはり原子間の結合が消滅している．
いずれも不安定である外殻原子が Al のクラスターは，図 5.4 に示すように中心原子を置換しても各原子間に結合を生じていない．一方，図 5.5 および 5.6 に示した Zr，Ti
を外殻原子とするクラスターは，Cu，Ni，Al と異なりほぼ円形の分布を示す．断面において外殻原子間の結合が存在すれば図 5.2〜5.4 のような分布になるので，Zr，Ti ク
ラスターでは外殻原子間の結合が極めて弱いことが示唆される．
図 5.5，5.6 中の丸で囲んだ部分は図 4.8(b) で述べた外殻原子 3 つと中心原子の四面体の中点に価電子を有する結合である．単元系のこの部分の電子密度は勾配がなく低い値であるが，中心原子を置換して安定となったクラスターは明確な中心原子−外殻原子間の強い結合をしている．結合部の電子密度は Cu，Ni が外殻原子である安定クラスター中の直線的結合と比べて低い値であるが，クラスターの detCIJ は近いオーダーである．これは四面体中心に電子を有する結合は 4 原子の間で形成されているため，2原子間での直線的結合よりも変形抵抗を増加させやすい結合形状であるからと考えら

55

れる．

StableCut surface
3
0.4 (e/Έ)0

Ni1Cu12StableAl1Cu12Zr1Cu12StableStableTi1Cu12

Cu13

Fig.5.2 Distribution of valence electron density in A1 Cu12 .

56

StableCut surface
3
0.4 (e/Έ)0

Cu1Ni12StableAl1Ni12Zr1Ni12StableStableTi1Ni12

Ni13

Fig.5.3 Distribution of valence electron density in A1 Ni12 .

57

Cu1Al12

Ni1Al12

Cut surface
3
0.4 (e/Έ)0

Zr1Al12

Ti1Al12

Al13

Fig.5.4 Distribution of valence electron density in A1 Al12 .

58

StableStableNi1Zr12

Cu1Zr12

Cut surface
3
0.4 (e/Έ)0Al1Zr12Ti1Zr12

Zr13

Fig.5.5 Distribution of valence electron density in A1 Zr12 .

59

StableStableCu1Ti12Ni1Ti12

Al1Ti12

Zr1Ti12StableCut surface
3
0.4 (e/Έ)

0

Ti13

Fig.5.6 Distribution of valence electron density in A1 Ti12 .

60

5.3

結言

Cu，Ni，Al，Zr，Ti 単元系クラスターの中心原子を他元素で置換した二元系クラス
ターについて，合金化による単元系孤立クラスターとの格子安定性の差を比較するた
めに，弾性係数の行列式と小行列式の値，および電子構造変化について調べた．得られた結果をまとめて以下に示す．

1. 安定な Cu，Ni 単元系の中心原子を置換すると，Zr で置換した系以外でいずれもdetCIJ が正値となった．さらに detCIJ の値は Al1 Cu12 を除いていずれも単元系
より大きな値を示した．

2. 不安定な Zr，Ti 単元系の中心原子を Cu，Ni で置換すると，いずれも detCIJ が正値となり安定化した．Ti 単元系については Al 原子で置換した系も detCIJ が正となった．また，これらの detCIJ の値は Cu，Ni 単元系より大きくなった．

3. 外殻原子が Al のクラスターは全て detCIJ が負となり不安定であった．
4. 中心原子を Zr で置換して不安定となった Zr1 Cu12 ，Zr1 Ni12 クラスターは，クラスター内部の電子密度が著しく低くなり，原子間の結合が切れていた．

5. 合金化によって Cu，Ni 単元系クラスターより detCIJ が大きくなった系は，Cu，
Ni を外殻原子とする場合は主としてクラスター内部への電子の局在化を生じていた．Cu，Ni を中心として Zr，Ti を外殻原子としたクラスターは，外殻原子 3
つと中心原子の四面体の中点に電子を有する結合が強くなっていた．

61

第 6 章 元素割合を変えた
クラスターの安定性評価
前章では，中心置換二元系クラスターの安定性を調べ，単元系との比較，元素の組
み合わせによる変化を原子間の電子構造から議論した．本章では，クラスター中の元素割合を変えた二元系クラスターの安定性変化を調べるために，実際に金属ガラスが構成される Zr70 Cu30 (9) を想定して，Cu 単元系の外殻原子を全部または一部 Zr に置換した二元系クラスター Cux Zr13−x について，これまでと同様に安定性変化ならびに電子構造変化を議論する．

6.1

解析方法

解析対象の 20 面体クラスター Cux Zr13−x (x = 1〜12) を図 6.1 に示す．いずれも中心原子を Cu とし，周囲の外殻原子の元素割合を変化させたものである．外殻原子配置の自由度も考慮すべきであるが，ここでは簡単のために，それぞれの元素を上下に偏析させた．解析条件を表 6.1 に示す．クラスター内での原子間距離は Zr の結晶状態
での最近接原子間距離 3.17 Åとした．前章までと同様，辺長が原子間距離の 4 倍の立
方体セル中心に配置して孤立クラスターとして，原子構造緩和および摂動計算を行い弾性係数 CIJ ならびにその行列式 detCIJ を評価した．

62

Table 6.1 Calculation conditions for Cux Zr13−x icosahedrons.

Number of electrons
Distance between atoms (Å)
Length of cell edge (Å)
Number of bands
Cutoﬀ energy (eV)Number of k points

63

Cux Zr13−x
7x + 52
3.17
12.68
7x + 782
292.2433

x=1

Cu1 Zr12

x=2

Cu2 Zr11

x=3

Cu3 Zr10

x=4

Cu4 Zr9

x=5

Cu5 Zr8

x=6

Cu6 Zr7x=7

Cu7 Zr6

x=8

Cu8 Zr5x=9

Cu9 Zr4

x = 10

Cu10 Zr3

x = 11

Cu11 Zr2

x = 12

Cu12 Zr1

Cu atom

Zr atom

Fig.6.1 Cux Zr13−x (x = 1〜12) icosahedrons.
64

6.2

解析結果と考察

6.2.1

元素割合とクラスターの安定性

Cux Zr13−x の弾性係数の行列式 detCIJ の値を図 6.2 に示す．図は横軸に Cu 原子の数，縦軸に行列式の値をとっている．中心原子はいずれも Cu である．左端の x = 1 は前章で示した，中心が Cu で外殻原子が全て Zr の Cu1 Zr12 クラスターである．そのクラスターの外殻原子をひとつだけ Cu に置換すると，値は著しく増加し 1 オーダー高い値となる (4.16×105 GPa6 )．この Cu2 Zr11 は 1 つだけ Cu-Cu ボンドを有する．外殻原子に Cu を 2 つ以上 (x ≧ 3) 有する系では detCIJ の値は大きく低下し，Cu1 Zr12 よりも低い値となる．クラスター外殻に Zr 三角格子が 1 つだけ存在する Cu10 Zr3 でやや detCIJの値は大きくなるが，Zr をひとつだけ外殻原子に含む Cu12 Zr1 クラスターの detCIJ の値は非常に小さくなり，Cu 単元系クラスターに近いオーダーとなる．
図 6.3 に Cux Zr13−x の弾性係数の各成分の変化を示す．最も大きな detCIJ の値を示
した x=2 の系では，垂直成分 C11 〜C33 ，せん断成分 C44 〜C66 のいずれも他のクラス
ターより大きくなっており，また異方性も小さい．他のクラスターについても，全体的な傾向は C11 〜C33 ，C44 〜C66 とも，図 6.2 の detCIJ のそれと対応している．一方，

x=7 のクラスターでは C33 は大きな値を示しているが，C11 ，C22 は低く，異方性が強いため detCIJ の値は低くなっている．

65

Determinant of elastic coefficinents, detCIJ x105 (GPa6)

5

4

3

2

1

0

1

2

3

4

5

6

7

8

9

10 11 12

Number of Cu atoms, x

Fig.6.2 Determinant of elastic coeﬃcients of Cux Zr13−x icosahedrons.

Compornents of CIJ (GPa)

20

16

12

8

4

0

1

2

3

4

5

6

7

8

9

10 11 12

Number of Cu atoms, x
C12
C13
C23

C11
C22
C33

C44
C55
C66

Fig.6.3 Components of CIJ of Cux Zr13−x icosahedrons.

666.2.2電子構造

図 6.4 および 6.5 に Cux Zr13−x クラスターの断面の電子密度分布を示す．等高線間隔
は 0.02 e/Å 3 で記しており，中心原子 Cu −外殻原子 Zr 間の等高線には後の考察のためにその値を記した．赤く高密度になった円の中心が Cu 原子であり，原子核近傍が低い密度になっている原子が Zr である．detCIJ が著しく大きな値を示した Cu2 Zr11 では Cu-Cu ボンドができたため，図 6.4(b) に示すようにクラスター内部での電子密度が増加しており，中心原子− Zr 間の電子局在も Cu1 Zr12 より強くなっている．一方，外殻の Cu 原子が増えると，電子が Cu 原子に引き寄せられてクラスター上部に移動して
おり，クラスター下部にある Zr と中心原子との間の結合部分の電子密度は Cu5 Zr8 を除いて著しく低下している．Cu が 8 つ以上の系では，クラスター下部の Cu 原子まわりの電子密度によって，Zr まわりの電子密度は再び増加している．

x=10 までの Cux Zr13−x クラスターについて，図 6.6 および 6.7 に図中右下に示す位
置の Zr 三角格子を作る面の電子密度分布を示す．Zr 三角格子では図 4.8(b) に示すような四面体中心に電子を有する結合が形成されており，これは前章で述べたように二元系クラスター内で電子密度が高くなった場合，直線的結合よりも変形抵抗を増加させると考えられる．Cu2 Zr11 は三角格子内部での電子密度が Cu1 Zr12 より増加しており，中心原子− Zr 間と外殻部の三角格子に電子が多く局在するという 2 つの要因が重なって detCIJ は急増したものと考えられる．x=3〜6 では三角格子内部の電子密度の値は同程度であり，中心原子− Zr 間の電子密度の低下によって detCIJ が減少したことが示唆される．一方，x=7〜10 では Cu を増やすと三角格子内部の電子密度の高い部分が広がっており，三角格子の結合力は強くなっていることが示唆され，そのために detCIJ が再び増加したものと考えられる．

67

0.24
0.22
0.24

0.22
0.20

3
0.4 (e/Έ)(a) x=1(b) x=2

0

0.20

0.24

0.20

0.24

0.20

Cut surface(c) x=3(d) x=4

0.14
0.22
0.18

0.20

0.24

0.22

(e) x=5(f) x=6

Fig.6.4 Distribution of valence electron density in Cux Zr13−x (x = 1〜6).

68

0.16

0.18

0.18

0.24

0.20

3
0.4 (e/Έ)(a) x=7(b) x=8

0

0.20

0.20
0.20

0.22

(c) x=9(d) x=10

Cut surface

0.24

0.22

0.20
0.22

(e) x=11(f) x=12

Fig.6.5 Distribution of valence electron density in Cux Zr13−x (x = 7〜12).

69

0.20

0.22

0.20

(a) x=1(b) x=2

0.22

0.22

0.22

(c) x=3(d) x=43
0.4 (e/Έ)

0.22

Zr
Zr
Zr

0

(e) x=5

Fig.6.6 Distribution of valence electron density of the triangle lattice in
Cux Zr13−x (x = 1〜5).

70

0.24
0.22

(a) x=6(b) x=7

0.26

0.28

(c) x=8(d) x=93
0.4 (e/Έ)

Zr

0.28

Zr
Zr

0

(e) x=10

Fig.6.7 Distribution of valence electron density of the triangle lattice in
Cux Zr13−x (x = 6〜10).

71

6.3

結言

クラスター中の元素割合を変えたときの安定性変化を調べるために，Zr70 Cu30 (9) の金属ガラスを想定した解析として，Cu 単元系の外殻原子を全部または一部 Zr に置換した二元系クラスター Cux Zr13−x の第一原理格子不安定解析を行い，弾性係数の行列式 detCIJ の値，弾性係数の各成分および電子構造の変化について検討した．得られた
結果をまとめて以下に示す．

1. Cu を中心とした二元系クラスター Cux Zr13−x の detCIJ はいずれも正値であったが，外殻原子を 1 つだけ Cu とし，中心原子と 1 つだけ Cu-Cu ボンドを有する

Cu2 Zr11 の値は他と比べて一桁大きい値となった．
2. 外殻原子の Cu 原子を 2 つ以上とすると，detCIJ は中心原子 Cu，外殻全てが Zrのクラスターより低い値となった．x=3〜10 の範囲では x=7 近傍で detCIJ の値
が極小値をとり，x=10 で比較的大きな値をとった．

3. detCIJ の値が極めて大きい Cu2 Zr11 は，弾性係数の垂直成分 C11 〜C33 およびせん断成分 C44 〜C66 いずれも他のクラスターに比べて大きくなり，かつ異方性が小さかった．また，元素割合を変えたときの弾性係数の垂直成分 C11 〜C33 とせん断成分 C44 〜C66 の変化の全体的な傾向は，detCIJ のそれと対応していた．

4. Cu2 Zr11 は，中心原子と Zr の間，外殻の Zr 原子が作る三角格子間のいずれにおいても，Cu1 Zr12 よりも高い電子密度を示し，Cu を増やすと x=7 までは中心原子− Zr 間の電子密度が低下し，その後 x=10 までは Zr 三角格子の電子密度が上昇していることを明らかにした．これらの電子密度変化が detCIJ の変化に対応している可能性を示唆した．

72

第7章 結 論
アモルファス合金設計の指標の 1 つとすべく，アモルファス構造の安定化に重要な

20 面体クラスターの格子安定性を，平面波基底ウルトラソフト擬ポテンシャル法に基づく第一原理解析により評価した．本研究で得られた結果を以下に総括する．
第 2 章では，第一原理分子動力学法の基礎理論および電子状態計算の高速化手法について説明し，第 3 章では，格子不安定解析の概要ならびに弾性係数を用いた安定性
評価について述べた．
第 4 章では，クラスターと結晶の安定性の違い，およびクラスターの周囲の原子の有無による安定性変化を明らかにすべく，Cu，Ni，Al，Zr，Ti の各単元素について，
孤立クラスター，クラスター準結晶，そして完全結晶の単位格子の解析を行った．ク
ラスターはいずれも Cu，Ni が弾性係数 CIJ の成分が全て正となり，その行列式も正値をとり安定であった．Al，Zr，Ti クラスターは CIJ に負の成分を生じる不安定な構造であった．完全結晶はいずれの元素も正値をとった．Cu，Ni の detCIJ の大小関係
は，孤立クラスター＜クラスター準結晶＜完全結晶でいずれも Cu ＜ Ni となった．ク
ラスター内の原子間結合距離はいずれにおいても外殻原子間の距離よりも中心原子−
外殻原子の距離が短く，不安定なクラスターは完全結晶の最近接原子間距離からの差が大きくなっていた．結合の電子密度分布をみると，Cu，Ni，Al クラスターは 2 原子間で直線的な結合をしているのに対し，Zr，Ti クラスターは 3 つの外殻原子と中心原子の四面体中心に電子を有する結合をしていた．また，Cu，Ni 孤立クラスターは中心原子−外殻原子間に結晶の最近接結合より電子を多く有しており，クラスター準結晶
では個々のボンド結合力は孤立クラスターより弱いが，周期的な構造となることにより孤立クラスターと比べて剛性が上昇していることを示した．
第 5 章では，合金化による安定性の変化を検討することを目的として，Cu，Ni，Al，Zr，Ti 単元系の中心原子をそれぞれ別の元素で置換したクラスターについて検討した．

73

Cu，Ni を含むクラスターは，それらが外殻原子である場合は中心原子が Zr のとき，
Cu，Ni が中心原子の場合は外殻原子が Al のときだけ不安定となった．Cu，Ni を含む
その他の安定なクラスターは，Al を中心とし，Cu を外殻とするクラスター以外，全て Cu，Ni 単元系より値が大きくなった．これらのクラスターは，Cu，Ni を外殻原子とするものはクラスター内部へ電子がより局在し，中心原子を Cu，Ni とするものは 3つの外殻原子と中心原子の作る四面体の中心に電子を有する結合となっていた．外殻原子を Al とするクラスターは全て detCIJ が負となり不安定であったが，不安定な Ti
クラスターの中心原子を Al としたものは安定となった．中心が Zr で外殻が Cu またはNi のクラスターは，クラスター内部の電子密度が著しく低くなり結合が切れていた．
第 6 章では，クラスター中の元素割合を変えたときの安定性変化を調べるために，実在する Zr70 Cu30 の金属ガラスを想定して，Cu 単元系の外殻原子を全部または一部 Zrに置換した二元系クラスター Cux Zr13−x を解析した．これらのクラスターはいずれもdetCIJ ＞ 0 で安定であった．Cu を中心として Zr を外殻とするクラスターの外殻原子を 1 つだけ Cu で置換したクラスターは前者に比べて detCIJ が一桁大きい値となった．
外殻の Cu 原子を 2 つ以上とすると，detCIJ の値は低下し，Cu 単元系のクラスターに近づいた．
本研究では 5 種類の元素について単元系，二元系クラスターの格子不安定性を評価し，その変化をもたらす要因について，クラスター構造内の電子状態変化から検討した．detCIJ の大小関係の物理的意味については議論の余地があるが，CIJ の増加に対応するものと考えれば，クラスターの力学的特性を表す一つの基準と考えることができる．本論文では 5 つの元素しか対象としなかったが，周期表にしたがって様々な元素のクラスターの格子不安定解析を行い，さらに 3 元系以上についても検討すること
で，将来アモルファス合金設計の指針になるデータベースの構築につながるものと期
待される．

74

参考文献
(1) 増本健, 深道和明, ”アモルファス合金”, (1982), アグネ.
(2) 井上 明久, 新しい精密工学用材料としてのバルク金属ガラス, 精密工学会誌, 70,
9, (2004), 1154-1158.
(3) A. Inoue, K. Ohtera, K. Kita, and T. Masumoto, New amorphous Mg-Ce-Ni
alloys with high strength and good ductility, Jpn. J. Appl. Phys., 27 (1998),
2248-2251.
(4) A. Inoue, Stabilization of Metallic Supercooled Liquid and Bulk Amorphous
Alloys, Acta Mater., 48, (2000), 279-306.
(5) T. Nasu, M. Sakurai, T. Kamiyama, T. Usuki, O. Uemura, K. Tokumitsu andT. Yamasaki, Structural Comparison of M-W (M = Fe, Ni) Alloys Produced byElectrodeposition and Mechanical alloying, Mater. Sci. Eng. A, 375, (2004),163170.
(6) E. Matsubara, M. Sakurai, T. Nakamura, M. Imafuku, S. Sato, J. Saida and
A. Inoue, Environmental structural studies in amorphous and quasicrystalline
Zr70 Al6 Ni10 Pt14 alloys, Scripta Materialia, 44, (2001), 2297-2301.
(7) E. Matsubara, T. Ichitsubo, J. Saida, S. Kohara and H. Ohsumi, Structural
study of Zr-based metallic glasses, Journal of Alloys and Compounds, 434,
(2007), 119-120.
(8) W. K. Luo, H. W. Sheng, F. M. Alamgir, J. M. Bai, J. H. He and E. Ma,
Icosahedral Short-Range Order in Amorphous Alloys, Phys. Rev. Lett., 92,(2004),145502.
(9) J. Saida, M. Imafuku, S. Sato, E. Matsubara, A. Inoue, Local structure in
Quasicrystal-forming Zr-based metallic glasses correlated with a stability of thesupercooled liquid state, Journal of Non-Crystalline Solids, 353, (2007), 37043708.

75

(10) J. Saida, M. Imafuku, S. Sato, T. Sanada, E. Matsubara, A. Inoue, Correlation
between local structure and stability of supercooled liquid state in Zr-based
metallic glasses, Mater. Sci. Eng. A, 449-451, (2007), 90-94.
(11) W. Kohn, and L. J. Sham, Self-Consistent Equations Including Exchange andCorrelation Eﬀects, Phys. Rev., 140, (1965), 1133-1138.
(12) 上杉徳照, 東健司, ガラス系性能予測のための溶解熱と寸法因子の第一原理データベースの構築, 日本機械学会第 18 回計算力学講演会講演論文集, (2005), 471-472.

(13) G. Kresse and J. Hafner, Ab initio molecular-dynamics simulation of the liquidmetal amorphous-semiconductor transition in germanium, Phys. Rev. B, 49,(1994), 14251-14269.
(14) X. D. Wang, S. Yin, Q. P. Cao and J. Z. Jiangb, Atomic structure of binary
Cu64.5 Zr35.5 bulk metallic glass, APPLIED PHYSICS LETTERS, 92, (2008),
011902-011905.
(15) M. Born and K. Huang, Dynamical Theory of Crystal Lattices, (1954) Oxford
UP.
(16) F. Milstein, Indirect Exchange Coupling of Magnetic Moments in Rare-Earth
Metals, Phys. Rev. B, 3, (1971), 1125-1133.
(17) J. Wang and S. Yip, Crystal instabilities at ﬁnite strain, Phys. Rev. Lett, 71,(1993), 4182-4185.
(18) J. Wang, S. Yip, S. R. Phillpot and D. Wolf, Mechanical instabilities of homogeneous crystals, Phys. Rev. B, 52, (1995), 12627-12635.
(19) D. C. Wallace, Thermodynamics of Crystals, (1972), Wiley, Newyork.
(20) P. Hohenberg and W. Kohn, Inhomogeneous Electron Gas, Phys. Rev., 136,
(1964), 864-871.
(21) J. P. Perdew and A. Zunger, Self-interaction correction to density-functional
approximations for many-electron systems, Phys. Rev.,B, 23, (1981), 5048-5079.
(22) C., Kittel. 著，宇野良清ら訳，キッテル固体物理学入門 上 第 7 版，(1998), 丸善.
76

(23) G. B. Bachelet, D. R. Hamann and M. Schlüter, Pseudopotentials that work:
From H to Pu, Phys. Rev. B, 26, (1982), 4199-4228.
(24) N. Troullier and Jose Luriaas Martins, Eﬃcient pseudopotentials for plane-wave
calculations, Phys. Rev. B, 43, (1991), 1993-2006.
(25) L.D. ランダゥ，E.M. リフシッツ著，佐々木健，好村磁洋訳，量子力学 1，(1992),
東京図書．

(26) 上田顕，コンピュータシミュレーション， (1990)， 朝倉書店．
(27) 香山正憲，固体材料の電子状態の基礎，(1993)， 材料学会勉強会資料．
(28) O. H. Nielsen and R. M. Martin, Quantum-mechanical theory of stress and force,
Phys. Rev. B, 32, (1985), 3780-3791.
(29) P. G. Dacosta, O. H. Nielsen and K. Kunc, Stress theorem in the determinationof static equilibrium by the density functional method, J. Phys. C : Solid State
Phys.,, 19, (1986), 3163-3172.
(30) D. Vanderbilt, Soft self-consistent pseudopotentials in a generalized eigenvalue
formalism, Phys. Rev. B, 41, (1990), 7892-7895.
(31) L. Kleinman and D. M. Bylander, Eﬃcacious Form for Model Pseudopotentials,
Phys. Rev. Lett., 48, (1982), 1425-1428.
(32) 尾口多美夫，佐々木泰造，固体物理， 25， (1990)， 857．
(33) C. L. Fu and K. M. Ho, First-principles calculation of the equilibrium groundstate properties of transition metals: Applications to Nb and Mo, Phys. Rev. B,
28, (1983), 5480-5486.
(34) M. P. Teter, M. C. Payne and D. C. Allan, Solution of Schrodinger s equation
for large systems, Phys. Rev. B, 40, (1989), 12255-12263.
(35) R. D. King-Smith and D. Vanderbilt, First-principles study of crystalline silica,
Phys. Rev. B, 49, (1994), 5828-5831.

77

(36) M. D. Gillan, J. G. Hamilton, Dehydrogenation of the ring-opened metathesis polymers of bicyclo[2.2.1]hepta-2,5-diene: a novel approach to conducting
elastomers, J. Phys. Condens. Matter 1, (1989), 23-26.
(37) I. Stich, R. Car, M. Parrinello and S. Baroni, Conjugate gradient minimizationof the energy functional: A new method for electronic structure calculation,
Phys. Rev. B, 39, (1989), 4997-5004.
(38) D. M. Bylander, Leonard Kleinman and Seongbok Lee, Self-consistent calculations of the energy bands and bonding properties of B12 C3 , Phys. Rev. B, 42,(1990), 1394-1403.
(39) W. Li and T. Wang, Ab initio investigation of the elasticity and stability ofaluminium, J. Phys.: Condens. Matter, 10, (1998), 9889-9904.
(40) G.Kresse and J.Furthmüller, Eﬃcient iterative schemes for ab initio total-energy
calculations using a plane-wave basis set, Phys.Rev.B, 54, (1996), 11169-11186.
(41) D. C. Langreth and J. P. Perdew, Theory of nonuniform electronic systems. I.
Analysis of the gradient approximation and a generalization that works, Phys.
Rev. B, 21, (1980), 5469-5493.
(42) P.Pulay, Convergence acceleration of iterative sequences. the case of scf iteration,
Chem.Phys.Lett., 73, (1980)，393-398.
(43) Kresse,G., and Furthmüller,J., VASP the GUIDE, (2003).
(44) H. J. Monkhorst and J. D. Pack, Special points for Brillouin-zone integrations,Phys. Rev. B, 13, (1976), 5188-5192.
(45) 屋代如月, 西村正臣, 冨田佳宏, 不均一系の大域的・局所的安定性：ナノ多結晶体とアモルファス金属の局所格子不安定解析, 日本機械学会論文集 A 編, 73, (2007),

216-223.
(46) 井上明久, 金属ガラスの発展の経緯, 金属, 75 No.1, (2005), 7-11.
(47) 尾方成信, 材料特性評価のための非経験的分子動力学法の基礎的研究, 大阪大学
博士論文.

78

関連講演論文(1) 久馬雅彦，屋代如月，冨田佳宏，“Cu-Zr 二十面体クラスターの第一原理格子不安定解析”, 日本機械学会第 19 回計算力学講演会 講演論文集, (2006), 639-640.

(2) 久馬雅彦，屋代如月，冨田佳宏，“正 20 面体短距離秩序クラスターの中心原子の安定性変化：第一原理格子不安定解析”, 日本機械学会第 20 回計算力学講演会 講
演論文集, (2007), .265-266.

79謝辞
本研究は，著者が学部四回生時に研究室に配属されてから三年間で冨田佳宏教授，ならびに屋代如月准教授のご指導の下で行ったものであります．研究を遂行していくにあたって，両先生の終始丁寧かつ適切な御教示によって浅学な著者は的確に導いていただき，本論文を完成させることができたものであります．心より深謝し，御礼申し上げます．また，研究活動へのお二人の情熱は，工学の世界に入ってまだ日も浅い著者にとって，畏敬の存在であるとともに大きな心の支えになったことは言うまでもありません．研究活動をこのような素晴らしい両先生の下で推進し，工学を学べたのは
著者の何よりの幸福であります．
同研究室の長谷部忠司准教授と，高木知弘准教授（京都工芸繊維大学）には研究会等で貴重なご助言を頂いたり，相談にのっていただきました．厚く御礼申し上げます．
また，研究生活についての種種の相談に乗っていただき，研究室の環境作りに尽力くださいました古宇田由夫技官に心より感謝いたします．
さらに，第一原理グループの先輩として，研究室配属当初から研究手法について直接ご指導していただいた久保圭佑氏 (現 富士ゼロックス株式会社) に感謝いたします．
ならびに格子不安定グループの先輩である西村正臣氏には様々なご助言を頂きました，
お礼申し上げます．本研究はお二人の先輩方のおかげで一層スムーズに進めることができたものであります．また，有益な議論および情報交換をした同じ第一原理グループの横川望君，山本智君および，格子不安定グループの樋口昌宏君に感謝いたします．
また，研究室でお互い励まし，助け合いながら研究生活をすごした同輩の朝西誠君，古賀裕君，田中克昌君に感謝いたします．
最後に 6 年間の学生生活で経済的，精神的に積極的に支えていただきました父・久
馬武彦，母・久馬和代，祖父・久馬義人をはじめとする家族に最大限の感謝と敬意を表します．ありがとうございました．
平成 20 年 2 月
久馬雅彦

80