SUMMARY
Recent rapid progress in computers has made it possible to elucidate the mechanicalproperties such as theoretical strength as well as static atomic structure by using
the ab-initio molecular dynamics. The theoretical strength, however, is evaluatedon the assumption that crystals deform in a deﬁnite deformation path because ofcomputational limitation, so that it is necessary to clarify the relationship between
the theoretical strength and bifurcation criteria to the other deformation path. The
lattice instability analysis based on the classical interatomic potentials have suggested
that the bifurcation to the anisotropic transverse contraction becomes more important
than the theoretical tensile strength in the uniaxial [001] tension.
In this study, the bifurcation point under tension is investigated by means of ab-initio
molecular dynamics. Si and Al single crystal are subjected to the uniaxial [001] tension
under isotropic and anisotropic transverse contractions. The diﬀerence in free energy
shows that the deformation path of isotropic Poisson contraction would bifurcate toanisotropic contraction at εzz = 0.093 and εzz = 0.051 for Si and Al, respectively, while
the theoretical strengths under isotropic contractions reach εzz = 0.25(Si) and εzz =0.18(Al). Then the lattice instability at these bifurcation points is investigated on thebasis of the positive deﬁniteness of the determinants of elastic stiﬀness coeﬃcients. It
is shown that the positiveness is violated after the bifurcation point and the instability
is caused by the negativeness of the minor determinants of the stiﬀness matrix, which
represents the compliance against the anisotropic transverse contractions.

要約
近年の計算機性能の飛躍的な向上によって，非経験的第一原理分子動力学法による変形解析が可能となり，理想強度等の機械的特性が数多く報告されている．しかしながら，扱える原子数の制限から変形経路を仮定した解析となるため，それより得られる理想強度と自由度が高い系において生じる他の変形経路への不安定分岐との関係は不明である．一方で，結晶の安定限界をエネルギーバランスより求める理想格子不安定性解析では，例えば引張変形において横方向の非等方収縮を考えた場合，あるひずみにおいて結晶が不安定になり，よりエネルギー的に低い変形経路への分岐を生じることが示されている．このように，理想格子不安定条件は均一変形における安定限界を求めるものであるが，より自由度の高い系における局所の不安定変形開始の指標となりうる可能性がある．
本論文では，理想格子不安定条件を第一原理分子動力学法により精密に評価することを目的とする．その第一歩として，引張変形下での転位発生クライテリアとして重要と報告されている，横方向非等方収縮への変形分岐条件，および，変形分岐前後での格子不安定性について，Si 単結晶ならびに Al 単結晶を解析対象として検討した．[001]
方向引張下において等方的な Poisson 収縮から非等方収縮の方がエネルギー，応力ともに低くなる変形経路への分岐ひずみは，Si は εzz > 0.093，Al は εzz > 0.051 である
ことが示された．これらの変形分岐は，等方 Poisson 収縮を仮定した変形経路で得られる理想引張強度 εzz = 0.25(Si) および εzz = 0.18(Al) よりはるかに低いひずみで生じる．また，各ひずみ下における弾性剛性係数行列式の正値性により，変形分岐前後の格子不安定性を調べた結果，Si,Al いずれも分岐臨界ひずみ以降は負の値となり不安定状態となっていることがわかった．また，この不安定は，弾性剛性係数行列において横方向の非等方変形に対する抵抗を表す部分が負となることによりもたらされることが明らかになった．

Lattice Instability Analysis Basedon ab initio Molecular Dynamics

February 2002
Division of Mechanical Engineering,
Graduate School of Science and Technology,
Kobe University, Kobe, Japan

Masashi OHO

修士論文
第一原理分子動力学法による理想格子不安定解析

平成 14 年 2 月
神戸大学大学院自然科学研究科
機械工学専攻005t358n
大穂正史

目次
第1章 緒 論1

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

4

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

4

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

5

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

7

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

8

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

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

15

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

16

2.8 電子占有数 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.9

FFT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

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

19

第 3 章 Si 単結晶の理想格子不安定解析

24

3.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

24

3.1.1

モデル及び予備解析 . . . . . . . . . . . . . . . . . . . . . . . .

24

3.1.2

引張シミュレーション . . . . . . . . . . . . . . . . . . . . . . .

27

3.2 シミュレーション結果および考察 . . . . . . . . . . . . . . . . . . . . .

28

3.2.1

理想引張強度と横方向変形分岐 . . . . . . . . . . . . . . . . . .

28

3.2.2

変形分岐後の電子構造変化 . . . . . . . . . . . . . . . . . . . . .

32

3.2.3

格子不安定性に関する考察 . . . . . . . . . . . . . . . . . . . . .

34

3.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

第 4 章 Al 結晶の理想格子不安定解析

38

i

目次

ii

4.1 シミュレーション条件 . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.2 シミュレーション結果及び考察 . . . . . . . . . . . . . . . . . . . . . .

41

4.2.1

理想引張強度と横方向変形分岐 . . . . . . . . . . . . . . . . . .

41

4.2.2

変形分岐後の電子構造変化 . . . . . . . . . . . . . . . . . . . . .

45

4.2.3

格子不安定に関する考察 . . . . . . . . . . . . . . . . . . . . . .

48

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

50

第5章 結 論

51

参考文献

53

第 A 章 原子単位系

56

第 B 章 関連講演論文

58

謝辞

59

第1章
緒論
電子デバイス分野等において，より一層の高密度化・高機能化を目的とした材料の微細化・無欠陥化が進められている．このように原子レベルまで精密に制御された微小材料では，原子配置のわずかな変化がその電子的特性に大きな影響を及ぼすため，原子レベルからの検討が必要となる．
近年の計算機能力の飛躍的な向上を背景に，分子動力学法をはじめとする原子シミュ
レーションによる変形破壊現象の解明が盛んに行われている (1) ．従来から行われている，原子間相互作用を簡略化したポテンシャル関数により表す古典的分子動力学法では，すでに 10 億もの原子の集団的挙動を追跡することが可能となっている (2)−(4) ．しかしながら，古典的分子動力学法では，平衡状態の弾性定数や格子定数にパラメータ
フィッティングされたポテンシャル関数を用いているために，フィッティング時に考慮した平衡状態から大きく変化するような場合（例えば，大変形時や界面など）における正当性は保証されない．このため，古典的分子動力学法では，それにより観察される現象のメカニズムは議論できても，例えば，臨界応力や界面エネルギー等については定量的な評価を行うことができない．
第一原理分子動力学法は，原子間相互作用の評価に簡略化されたポテンシャル関数
を用いずに，各瞬間における原子位置とその原子種のみを必要情報とし，量子力学に基づいて逐次電子状態を計算し原子に作用する力を精密に評価する手法である．実験から求めた諸物理量を用いないという意味で非経験的手法とも呼ばれる．ここ 20 年で飛躍的な高精度化・高効率化がなされ，単結晶 bulk 状態のみならず表面や粒界・界面

1などといった複雑な系の安定構造や機械的特性の評価への適用が可能となってきており，同手法による電子・原子レベルからの材料特性解明が活発に行われつつある．例
えば，香山らは SiC 粒界，SiC/Ti 界面の第一原理引張試験 (5) を行い，界面の理想引張強度を検討している．この理想強度は実際の材料と比較して二桁近く大きいが，これは破壊の起点となるクラック先端の局所的な臨界応力値に対応すると報告している．
尾方らは AlN/Al 界面の構造解析を行い (6) ，AlN/Al 界面の adhesion energy を計算している．このエネルギーが母材である Al の表面エネルギーと同等であったことから，

AlN/Al 界面の結合力が母材のそれと変わらないことを示し，界面剥離が少ないとの実験事実を理論的に裏付けている．梅野ら (7) や楠ら (8) は，電子デバイス材料として一般的に用いられている無転位 Si 単結晶のせん断解析を行い，その理想せん断理想強度等を報告している．しかしながら，第一原理分子動力学法は電子状態計算を逐次行いながら原子運動を追跡するため，膨大な計算を必要とする．このため，古典的分子動力学法に比べると扱える原子数は極めて少なく，たかだか数十〜数百が限界となる．
それゆえ必然的に，理想強度等の解析では変形経路を仮定したものとなるため，自由度が高い実際の系において生じる他の変形経路への不安定分岐との関係を明らかにする必要がある．
変形経路の分岐について，従来より行われてきた格子不安定性解析 (9)−(16) が重要な指標を与える．格子不安定とは，結晶格子が変形に対する抵抗力を失うことにより，外から力を加えることなく連続的に変形が進行し，原子構造が元の構造を維持できなくなる状態を指す．その不安定条件を定量的に評価する手法は Born により体系化された(9)．この Born の不安定基準の考え方に基づき，理想結晶の不安定分岐点，および，その分岐経路を求める試みが 1970 年代に盛んに行われた (11)−(16) ．例えば，α 鉄の [001]方向の引張変形においては，等方的な Poisson 収縮の変形経路から，非等方収縮がエ
ネルギー的に低くなる変形経路への分岐を生じることが格子不安定性により示された(12)．この理想均一結晶の格子不安定性解析は Wang ら (18) によって拡張され，弾性剛性係数行列 (20) の正値条件を用いて結晶の不安定性を評価することにより，内部自由度を有する結晶に適用することが可能となった．また，屋代ら (21) は，自由度の高い

2

系において生じる局所変形（転位の発生）の開始が，局所の結晶格子が上述の横方向
非等方変形への不安定分岐を生じたことによって引き起こされていることを明らかにした．
上述のように，理想格子不安定解析は均一変形における安定限界を求めるものであるが，より自由度の高い系における局所の不安定変形開始の指標となる可能性がある．
しかしながら，これまでの格子不安定条件に関する検討は，基本的に原子間相互作用をポテンシャル関数で表したものであり，不安定分岐に対応する臨界ひずみ等，得られた数値は実際の材料のそれを表しているものではない．そこで本研究では，理想格
子不安定条件を第一原理分子動力学法により精密に評価することを目的とする．その第一歩として，引張変形下での転位発生クライテリアとして重要と考えられる，横方
向非等方収縮への変形分岐条件について考察する．本論文の各章の構成は以下のとおりである．
第 2 章では，第一原理分子動力学法の基礎理論について述べる．第 3 章では，[001]方向に引張を受ける Si 単結晶について，等方的な横方向 Poisson 収縮から，非等方収縮がエネルギー的に低い変形経路となる変形分岐点を第一原理分子動力学法により求める．また，その分岐前後の格子不安定性について，弾性剛性係数行列の正値性により議論する．第 4 章では，fcc 金属である Al 単結晶を対象に，第 3 章と同様に [001] 方
向引張時の横方向分岐について検討する．第 5 章では，本研究で得られた結果を総括
する．

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

Ts [ρ(r)] =

occ
∑
i

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

(2.2)

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

Exc [ρ] = T [ρ] − Ts [ρ] + E[ρ](2.4)

のように書ける．ここで，Exc は一電子近似のもとでの交換相関エネルギーであり，電子間相互作用を考慮した電子の運動エネルギー T [ρ] から，相互作用のない電子の運動
エネルギー Ts [ρ] を分離することによって，電子間の複雑な相互作用を全てこの項に押し込めている．
∫電子密度に関する拘束条件 ρ(r)dr = N のもとで式 (2.3) に変分原理を適用することにより，以下の一電子シュレディンガー方程式 (Kohn–Sham 方程式) が得られる (25) ．

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

となる．ここで，∑occ
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 の関数形 (26) を用いた．
ε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+(rs > 1)√2
3(1 + 1.0529 rs + 0.3334rs )0.6668rs−0.0584 + 0.0311 ln rs − 0.0084rs + 0.00133rs ln rs (rs < 1)
7

となる．

2.4

逆格子空間

結晶のように原子が周期的に配置された構造の第一原理計算には，逆格子空間が用いられる．実空間における格子点の位置ベクトル R が基本並進ベクトル a1 ,a2 ,a3 によって，(n1 , n2 , n3 は整数)

R = n1 a1 + n2 a2 + n3 a3(2.16)と表されるとすると，逆格子空間の基本並進ベクトル b1 ,b2 ,b3 は
a2 × a3
b1 = 2π
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 の定理 (34) より

ψ(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∫
f (r) exp[−iG · r] exp[iG′ · r]drΩ Ω
1∫=f (r) exp[−i(G − G′ ) · r]drΩ Ω

< k + G|f (r)|k + G′ > == f (G − G′ )(2.28)

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

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

表される．
pp
pp
Vlpp (r − Ra ) P̂l = Vloc(r − Ra ) + Vnloc,l(r − Ra ) P̂l

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

∑∑

Vapp,loc (|r − r a − R|)(2.32)raRppと表せることから，Vloc(G − G′ ) は以下より与えられる．
pp
Vloc(G − G′ ) =1 ∑exp[−i(G − G′ ) · r a ]Vapp,loc (G − G′ )Ωat aVapp,loc (G − G′ ) =∫Vapp,loc (r) exp(−i(G − G′ ) · r)dr

(2.33)(2.34)

ここで，Ωat は周期セルの体積，r a はセル内の原子 a の位置ベクトル，R はセルの位置ベクトルである．

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

pp
= Vnloc(k + G, k + G′ )Vapp,loc (k + G, k + G′ )= 4π∑

(2l + 1)Pl (cos ω)∫(2.35)

pp,nlocVa,l(r)jl (|k + G|r)jl (|k + G′ |r)r2 drl(2.36)
11

となる (32) ．ここで，Pl は Legendre 多項式，jl は球 Bessel 関数であり，ω は k + G と

k + G′ との間の角度である．
本研究で用いた BHS 擬ポテンシャルでは，

Zv
Vapp,loc (r) = −

[ 2
∑

r

pp,nloc
Vl,a
(r) =

3 (
∑

erf
ccore
i

[(√

αicore r)]](2.37)

i=1)[

Ai + r2 Ai+3 exp −αi r2](2.38)

i=1, αi は Bachelet らの論と定義されている (27) ．ここで，Zv は価電子の電荷，αicore , ccorei
文 (27) に一覧表となって示されており，Ai はそれらの関数となっている．図 2.1，図 2.2に Si, Al についての疑ポテンシャルを示す．

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

ρ(G) exp[iG · r]

(2.39)

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

(2.40)

ρ(r) =

G

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

ρ(G) exp[iG · r]

(2.41)

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

Vcoul (r) = 4π

∑ ρ(G)
G|G|2

exp[iG · r]

(2.42)

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

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

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

= 4π

′ 2

|G |

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

ρ(G)
|G|2

(2.43)

12

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

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

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

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

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′ )と逆空間での表式となる．

13(2.47)

42

l=0ion, a.u.

-2

V

0

-4

l =1
l =2
-6

V

-8

pp,loc

-10

0

1

2

3

4

5

6

r , a.u.
Fig.2.1 BHS type pseudo potential (Si)

2
l=0

-2

l =1

V

ion, a.u.

0

l =2
-4

V

pp,loc

-6
0

1

2

3

4

5

r , a.u.
Fig.2.2 BHS type pseudo potential (Al)

14

6

2.6

系のエネルギー

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

BZ ∑
occ
∑
k

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

1 ∑ ∑ a a′ ∑4π
− |G|2′)]
exp
EEwald =
Zv Zv
exp
[iG
·
(r
−
r
a
a
2
2 a a′
4γ 2G̸=0 Ωat |G|+−]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.49)である．
ここで

ρ(r) =∑

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

(2.50)

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

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

fk

occ
∑
n

fn

∑∑

G G′

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

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

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

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

15

+

BZ
∑
k

fk

occ
∑

fn

n

∑ ∑
G=0 G′ =0

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

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

2.7

応力

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

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

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

はスケーリングの元のもとで不変であるから，平均応力は∂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

1 ∑∑∂V pp,loc (G)Sa (G){ a 22Gα Gβ + Vapp,loc (G)δαβ }ρ(−G)Ωat G a∂(G )BZ ∑occ ∑ ∑ ∑ ∑
1 ∑
n
n∗
Ck+G
fk fn S(G − G′ )Ck+G′Ωat k n G G′ l a∂
1 pp,nloc
×{Va,l(k + G, k + G′ )}∂εαβ Ωat+16(2.54)(2.55)

+

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

+ δαβ

∑

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

G

1 ∂EEwaldZ ∑αa+− δαβ 2
Ωat ∂εαβΩat a(2.56)と表すことができる．ここで疑ポテンシャルに BHS 型を用いる場合，局所項は24πZa ∑∂Vapp,loc (G)
G2
G2
core
=
exp(−
(1
+
)c
)
G4 i=14αicore i
4αicore∂(G2 )(2.57)と表され，非局所項は4π
1 pp,nloc
Va,l(K, K ′ ) =(2l + 1)Pl (cos ω)Fa,l (K, K ′ )
Ωat
Ωat

′

Fa,l (K, K ) =∫ ∞0pp,nloc
jl (Kr)jl (K ′ r)Va,l(r)r2 dr(2.58)(2.59)

により，∂1 pp,nloc{V(k + G, k + G′ )}∂εαβ Ωat a,l
4π∂Fa,l (K, K ′ )′=(2l + 1)[−δαβ Pl (cos ω)Fa,l (K, K ) + Pl (cos ω)
Ωat
∂εαβ
′
′
′
′
K αK β + K αK β ′K αK β K αK β
+ {cos(+}Pl (cos ω)Fa,l (K, K ′ )] (2.60)′ 2 )−
2
|K|
|K |
KK ′
と表される．

2.8

電子占有数

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

fi =1
1 + exp{(εi − εF )/kB T }17(2.61)

を導入する．kB はボルツマン定数，T は電子系の温度である．実際の数値計算では

2∑

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

S = −kB∑

Ef = Etot − T S(2.63){fi ln fi + (1 − fi ) ln(1 − fi )}(2.64)

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

2.9

FFT

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

∑

n
Ck+G

(2.65)n
Ck+G
exp[iG · r]

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

(2.67)

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

(2.68)

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

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

ukn (r)u∗kn (r) =

∑∑
G

G′

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

ρ(r) =

occ ∑
BZ
∑
n

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

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

18

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

2.10

電子系の最適化手法

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

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

G

19

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

(2.70)

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

∑∑
G

G′

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

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

]





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

(2.72)

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

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

< 残差ベクトル >残差ベクトルは前述のとおり次式で表される．ここで ε′kn は，式 (2.70) の Etot を εknに置き換えたものである．また式中の i は， i 回目のステップにおける
を表す．

[
i
Rnk
=−

∂ε′kn
m∗
∂Ck+G

]i

(

)

i
= − H − λin I · C nk

20

という意味

(2.73)

ただし，
[

]

n,i
i
i
i
C nk
= Ck+GH = [Hk+G,k+G′ ] , λin =< ΨnkĤ Ψnk>′ ,(2.74)

である．これは，i のステップにおいて εkn を最小にする方向 (最急降下方向) を示すベクトルを表している．
iには，最終的に得られる次のステップの波動関数 Ψi+1
Rnk
nk が同じ k 点における n 以

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

∑ (

)

i
i∗
C mk
· R ink C mk

(2.75)

m̸=n<preconditioning>′i
残差ベクトル Rnk
に対して preconditioning という処理を施す．大きな逆格子ベクトルについては平面波の運動エネルギーが大きくなるが，このことが残差ベクトルに影響して収束性を悪化させる．preconditioning は，この問題を回避して収束を速めるた
i
めに行われる．preconditioning された残差ベクトルを Gnk
とすると
i
i′
Gnk
= K i · Rnk(2.76)と表される．ここで(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>
= < Ψnk
Ekin
2
KGG′ = δGG′

(2.77)(2.78)(2.79)(2.80)

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

)

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

∑ (

)

i∗
i
i
C mk
· Gnk
C mk

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

< 探索方向 >21(2.81)

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

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

γi =






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

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

(

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

i′
F nk

(2.84)(2.85)

)1

i′∗
i′ 2
F nk
· F nk< 新たな係数ベクトルの組み立て >新たな係数ベクトルの組立ては次のように行われる．
i
i
C i+1
nk = αC nk + βD nk(2.86)

結合係数 α と β は，エネルギー期待値 ε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.87)

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

} 21

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

α= {

ε12

}1

ε12 ε∗12 + (ε11 − γ)2 2
ε11 − γβ = −{}1ε12 ε∗21 + (ε11 − γ)2 2(2.89)(2.90)

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

22

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

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

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

23

第3章

Si 単結晶の理想格子不安定解析
古典的ポテンシャル関数を用いて行われた理想格子不安定性解析では，[001] 方向に単軸引張を受ける結晶において最初に表れる不安定条件は，横方向の非等方収縮に対する変形分岐点に対応することが示されている (12) ．また，多原子系の分子動力学シ
ミュレーションにより，転位の発生は局所の結晶格子が横方向の非等方変形を開始す
ることによりもたらされていることが報告されている (21) ．本章では，この横方向変形
分岐に対する格子不安定条件を，第一原理分子動力学法により精密に評価することを目的とし，[001] 方向に単純引張を受ける Si 単結晶のシミュレーションを行う．

3.1

シミュレーション条件

3.1.1

モデル及び予備解析

図 3.1 に示すダイヤモンド構造 Si 単結晶単位格子（原子数 8）を解析対象とし，k 点のサンプリング数は全ブリルアンゾーン内に対して 8 点，平面波のカットオフエネル
ギーは 8.0a.u. とした．また，擬ポテンシャルには BHS 型擬ポテンシャル (6) を，交換
相関エネルギーの表式には Perdew と Zunger による関数形 (26) を使用し，単位系として Hartree 単位系 (付録 A 参照) を用いた．計算に用いた条件を表 3.1 に示す．

まず，予備解析として以下の二つの解析を行った．

(1) 0K の格子定数 (0.5429nm) において，カットオフエネルギー Ecut を変化させた24

時の全エネルギー Etot の変化を計算し，平面波の打ち切りエネルギーを決定する．

(2) 打ち切りエネルギーを (1) で決定した値に固定し，全原子をダイヤモンド構造の格子位置に配置した状態で格子定数 a を変化させ，格子の膨張・収縮に対する系の自由エネルギーの変化を計算する．
解析 (1) で得られた，カットオフエネルギー Ecut と全エネルギー Etot との関係を図3.2 に示す．カットオフエネルギーが大きいほど基底状態に近づき全エネルギーは低い値となる．一方，カットオフエネルギーを大きくとると平面波の数が増加する．それに伴い，記憶すべき量が平面波の 2 乗に，ハミルトニアンの対角化に要する時間は 3乗に比例するため，解析精度と計算時間とのバランスをとることが重要である．本解析では，Ecut = 8.0 a.u. で十分収束したとしてこの値を採用する．
解析 (2) で得られた，格子定数と全エネルギーとの関係を図 3.3 に示す．a = 0.5374nmのときに全エネルギーが極小値をとり，誤差は約 1.01% にとどまっている．

[001]

z
x [100]y [010]
Fig.3.1 Simulation cell

Table 3.1 Calculation condition for analysis
Number of atoms
Lattice constantCutoﬀ energy
Number of k points
Boundary condition
Pseudopotential

25

80.5374 nm
8 a.u.
2×2×2
Supercell
BHS typeTotal energy Etot , a.u.

-30.6

-30.9

-31.2

-31.5

-31.8

0

2

4

6

8

10

Cutoff energy Ecut , a.u.

12

Fig.3.2 Relationships between Cutoﬀ energy and Total energy

Total energy Etot , a.u.

-31.7324-31.7326-31.7328-31.733-31.7332
0.536

0.537

0.538

0.539

Lattice constant a , nm

0.54

Fig.3.3 Relationships between lattice constant and Total energy

26

3.1.2

引張シミュレーション

予備解析 (1)(2) において求めた，全エネルギーが最小となる格子定数 0.5374nm で表されるスーパーセル (図 3.1) を無負荷基準状態とした．その後，種々の引張ひずみ

εzz をスーパーセルの z 軸方向に ∆εzz = 0.005(または 0.001) 刻みで与えて構造緩和を行った．ここで，横方向の収縮に関して次の 2 通りの解析を行った．
解析 I：等方収縮

εxx = εyy の等方収縮条件で横方向圧縮ひずみを与え，それぞれの原子に働くヘルマ
ン–ファインマン力に基づき内部緩和を行う．圧縮ひずみは 5 × 10−4 刻みで変化させて計算を行い，与えた引張ひずみ εzz に対して全エネルギーが最小となるスーパーセルを決定する (図 3.4 手順 (a))．
解析 II：非等方収縮
解析 I で求めたスーパーセルを基準とし，セルの x 方向長さ Lx ，y 方向長さ Ly のア
スペクト比 Ly /Lx を体積一定の条件下で 1.00 から 1.012 まで 0.001 刻みで変化させる．
解析 I と同様に内部構造緩和を行い，各アスペクト比のセルにおける原子構造と全エネルギーを求める (図 3.4 手順 (b))．
(a) Deformation path of isotropic
Poisson contraction

Total Energy Etot

Bifurcation pointzyx(b) Deformation path ofanisotropic contraction

Tensile strain εzz

Fig.3.4 Schematic of two deformation path under tension

27

3.2

シミュレーション結果および考察

3.2.1

理想引張強度と横方向変形分岐

解析 I（等方収縮) より得られた全エネルギー – 引張ひずみ関係，および，引張応力– 引張ひずみ関係を図 3.5，図 3.6 に示す．応力は文献 (48) に従い評価した．ひずみの増加とともに応力は単調に増加し，εzz = 0.25 でピーク応力 18.73GPa を示した．従って，等方収縮を仮定した解析で得られる理想引張強度は εzz = 0.25，σzz = 18.73GPaとなる．

εzz = 0.093〜0.099 において，等方収縮下の安定構造からアスペクト比 Ly /Lx を変えて計算した全エネルギーの変化を図 3.7 に示す．図中の破線は解析 I の等方収縮の場
合 (Ly /Lx = 1 に相当) の全エネルギーである．εzz = 0.093 以下の引張ひずみにおいては，図 3.7(a) に示すように非等方横ひずみ下ではエネルギーが上昇し，等方収縮する
変形経路が最もエネルギーが低い値をとっていたのに対し，εzz = 0.094 において等方収縮の場合よりも非等方収縮の場合にエネルギーが低くなる経路が現れる (図 3.7(b)，Ly /Lx = 1.001 の点)．以降は引張ひずみが大きくなるとともに最もエネルギーの低い点がアスペクト比の大きな側へとシフトした．等方収縮より低いエネルギーとなるこ
れらの変形経路を全エネルギー–引張ひずみ，および，引張応力–引張ひずみ関係に示すと図 3.8，図 3.9 のようになる．これより，εzz = 0.094 の点から等方収縮の場合に比べてエネルギー，応力ともに低くなる変形経路への分岐を生じることがわかる．また，
この変形分岐を安定限界としたときの理想引張強度は εzz = 0.093, σzz = 10.73GPa となる．
変形分岐の臨界ひずみ εzz = 0.093 における横方向ひずみは εxx = εyy = −0.033 であり，これより ⟨ 011 ⟩ 方向の分解せん断ひずみは εzx = 0.295 と計算される．また，分解せん断応力は σzx = 10.35GPa となる．(111) 面上の理想均一せん断を仮定した解析
では，垂直方向 (⟨ 111 ⟩, ⟨ 110 ⟩, ⟨ 112 ⟩) 応力を零に制御したときの理想せん断強度がεzx = 0.30，約 10GPa となることが報告されている (7) ．したがって，(111) 面の理想せん断強度と，[001] 引張下で生じる横方向非等方変形分岐には密接な関係があることが28示唆される．

Total energy Etot , a.u.

-31.62-31.64-31.66-31.68-31.7-31.72-31.74

0

0.04

0.08

0.12

0.16

0.2

0.24

0.28

Tensile strain εzz
Fig.3.5 Total energy–Strain curve under tension

Tensile stress σzz , GPa

20
18
16
14
12
10
8
6
4
2
0

0.04

0.08

0.12

0.16

0.2

0.24

0.28

Tensile strain εzz
Fig.3.6 Stress–Strain curve under tension

29

Lx

Total energy Etot , a.u.

-31.7132

Ly

[001]

z

-31.7134

Energy underisotropic contraction-31.7138-31.714-31.7142-31.7144

x [100]

Energy underanisotropic contraction-31.7136

y [010]1

1.002 1.004 1.006 1.008

1.01

1.012

Aspect ratio Ly/Lx-31.7128-31.7128-31.713-31.713

Total energy Etot , a.u.

Total energy Etot , a.u.

(a) εzz = 0.093-31.7132-31.7134-31.7136-31.7138-31.71411.002 1.004 1.006 1.008

1.01-31.7136-31.71381

1.002 1.004 1.006 1.008

1.01

Aspect ratio Ly/Lx

Aspect ratio Ly/Lx(b) εzz = 0.094(c) εzz = 0.095

1.012-31.7118

Total energy Etot , a.u.

Total energy Etot , a.u.

-31.7134-31.7141.012-31.7124-31.7126-31.7128-31.713-31.7132-31.7134-31.7136-31.7132-31.712-31.7122-31.7124-31.7126-31.7128-31.71311.002 1.004 1.006 1.008

1.01

1.012

1

1.002

1.004

1.006

1.008

1.01

Aspect ratio Ly/Lx

Aspect ratio Ly/Lx(d) εzz = 0.097(e) εzz = 0.099

1.012

Fig.3.7 Relationships between total energy and aspect ratio of transverse
strain under anisotropic contraction30

Total energy Etot , a.u.

-31.707-31.71-31.713

Isotropic contraction-31.716

Anisotropic contraction-31.719
0.08

0.09

0.1

0.11

Tensile strain εzz
Fig.3.8 Total energy–strain curve under tension

Tensile stress σzz , GPa

12

11

10
Isotropic contraction
Anisotropic contraction
9
0.08

0.09

0.1

Tensile strain εzz
Fig.3.9 Stress–strain curve under tension

31

0.11

3.2.2

変形分岐後の電子構造変化

εzz = 0.099 における，等方収縮および，非等方収縮それぞれの (110) 面の価電子密度分布を図 3.10 に示す．内殻電子を表示していないため，原子核近傍では密度が低くなっている．等方収縮下の変形経路では対称性を有しているため A∼D のボンドは全て同じ強さであるのに対し，非等方収縮下では B,D のボンドが強められ，A,C のボンドが弱くなっている．(1̄10) 面価電子密度分布の図 3.11 では，(110) 面と同様に，非等方収縮下では B,D のボンドが強まり，A,C のボンドが弱まっている．以上のボンドの変化を図 3.12 に模式的に示した．太線が強まったボンドを，波線が弱まったボンドを示している．すなわち，引張下で広がっていたボンド間隔が，太線部分では戻る側に変化し，波線部分ではより広がる方向へ変化していることを意味している．これは，図3.12 の右図に矢印で示したような部分転位を生じる原子移動に対応するものと考えら
れる．

[001]

z

A

x [100]y [010]BCD(a) εzz = 0.099 (Isotropic)(b) εzz = 0.099 (Anisotropic)

Fig.3.10 Distribution of valence electron density on (110) plane32

[001]zx [100]y [010]

B

C

A

D

(a) εzz = 0.099 (Isotropic)(b) εzz = 0.099 (Anisotropic)

Fig.3.11 Distribution of valence electron density on (1̄10) plane

A

A

A
A

A
A

C

C

[001]

A

z

A

A

x [100]AA

y [010]

Fig.3.12 Schematic of change in bond strength33

A

3.2.3格子不安定性に関する考察

均一変形する結晶を格子パラメータで表すことにより，系のエネルギー関数の支配変数を減らして 2 次導関数を導出し，その正値性によって系の安定性を評価する従来の理想格子不安定解析 (11) は，ひずみで表した系全体の変形を，系の弾性剛性係数行列の正値性に基づいて評価するように拡張された (17) ．弾性剛性係数は有限変形下の結晶の変形抵抗を表す物理量であり，応力と弾性係数により次式のように定義される．
(

Bijkl =

∂σij
∂εkl

)

= Cijkl + (σil δjk + σjl δik + σik δjl + σjk δil − 2σij δkl )/2(3.1)

ここで εkl は無負荷平衡状態からのひずみであり，δij はクロネッカーのデルタである．
sym
Bijkl の対象部分 Bijkl= 21 (Bijkl + Blkji ) の正値性が系の安定性を支配する (19) ．本解

析で対象とした横方向等方圧縮条件下の [001] 方向引張では，系の安定性は次の 6 × 6行列の正値性により評価される．


B11 B12 B13 000 B12 B11 B13
0
0
0 




 B13 B13 B33

0
0
0


detBIJ = det 00 B44 0
0 
 0



 0
0
0
0 B44 0 


0
0
0
0
0 B66

(3.2)

sym
ここで BIJ は Bijkl
を Voigt 表記 (20) したものである．

変形分岐前後の格子不安定性を評価するため，解析 I の等方収縮の εzz = 0.090 ∼0.094 の各平衡状態における弾性係数を，微小ひずみ摂動を与えたときの応力勾配から数値的に評価した（模式図 3.13）．得られた弾性係数，および，式 (3.1) より計算される弾性剛性係数の値を，分岐点前後のひずみ εzz = 0.093 と εzz = 0.094 の場合について表 3.2 に示す．また，εzz = 0.090 ∼ 0.094 における detBIJ の変化を図 3.14 に示す．

detBIJ は εzz の増加によって減少し，εzz = 0.093 では 0 に非常に近い値をとっているが正値であり安定である．一方，引張ひずみ 0.094 で detBIJ が負となった．この不安定は，図 3.15 に示すように行列式 (3.2) の中の横方向変形のバランスを表す主小行列
式

B11
B12

B12
B11

の値が負となることによりもたらされていた．このように，格子不安34

定性の観点からも，εzz > 0.093 が横方向変形に対する分岐クライテリアであることが
示された．

∆ε33=

}0.0005∆ε11=

}0.0005∆ε12=0.00035

y
x
C33 =

∆σ33
∆ε33

C31 = C13 =

C11 =
∆σ11
∆ε33

∆σ11
∆ε11

C12 = C21 =

∆σ22
∆ε11

C44 =

∆σ23
∆ε23

C55 =

∆σ13
∆ε13

C66 =

∆σ12
∆ε12

Fig.3.13 Calculating procedure of elastic constant

Table 3.2 Elastic constants and elastic stiﬀness coeﬃcientsεzz = 0.093 εzz = 0.094
εzz = 0.093 εzz = 0.094
C11 (GPa)
196.4
195.2
B11 (GPa)
196.4
195.2
C12 (GPa)
195.2
197.5
B12 (GPa)
195.2
197.5
C33 (GPa)
223.5
221.0
B33 (GPa)
234.3
231.8
C13 (GPa)
173.0
164.0
B13 (GPa)
162.2
153.2
C44 (GPa)29.5
28.9
B44 (GPa)29.5
28.9
C66 (GPa)28.9
34.1
B66 (GPa)
39.8
39.6

35

6

Magnitude of 6x6 determinant, detBIJ, GPa

0. 4

x1012
x1012
0.0 4

0. 3

0.02

0.0 0

0. 2
-0.02
0.091

0.092

0.093

0.094

0. 1

0.0
0.09

0.091

0.092

0.093

0.094

Tensile strain εzz
Fig.3.14 Relationship between tensile strain and crystal stability under [001]tension (isotropic contraction)

x104

Magnitude of 2x2 determinant,2
detBIJ (I, J = 1,2), GPa

1.50
1.25
1.00
0.75
0.50
0.25
0.00
-0.25
0.09

0.091

0.092

0.093

0.094

Tensile strain εzz
22
under [001] tension− B12
Fig.3.15 Relationship between tensile strain and B11(isotropic contraction)

36

3.3

結言

[001] 方向に単軸引張を受ける Si 単結晶について，横方向等方収縮とした場合と非等方収縮とした場合の系のエネルギー，および，電子状態を第一原理分子動力学法により検討した．得られた結果を要約して以下に示す．

(1) 引張に対して横方向等方収縮を仮定した場合と，非等方収縮とした場合の系のエ
ネルギーの変化より，εzz = 0.094 以降の引張ひずみにおいては，非等方収縮の方がエネルギー的に低くなる変形分岐が存在することが示された．

(2) (1) の変形分岐時の引張ひずみ，および，横方向 Poisson 収縮ひずみから計算される ⟨011⟩ 方向分解せん断ひずみは εzx = 0.295 であった．これは，第一原理分子動力学法によって求められた ⟨011⟩ 方向の理想せん断強度 εzx = 0.30(7) に近い．

(3) 非等方収縮変形時の電子状態変化から，分岐後の変形経路では部分転位に対応す
る方向へ原子が移動するように，ボンドの強弱が変化していることが示された．
これと (2) より，本解析で対象とした引張時の横方向変形分岐は，転位発生に対する理想せん断強度と密接に関係していることが示唆される．

(4) 各引張ひずみ下における弾性剛性係数を数値的に求め，その正値性に基づき変形
分岐前後の格子不安定性を評価した．その結果，εzz = 0.094 以降は弾性剛性係数の行列式の値が負となり，不安定状態となっていることが示された．また，こ22が負となること− B12の不安定は，横方向変形に対する抵抗を表す小行列式 B11によりもたらされることが判明した．このように，格子不安定性の観点からも，

εzz > 0.093 が横方向変形に対する分岐クライテリアであることが示された．

37

第4章

Al 結晶の理想格子不安定解析
本章では，前章と同様に，fcc 金属である Al 単結晶を対象として，[001] 引張変形下における変形経路の分岐とその分岐をもたらす電子構造の変化，ならびに格子不安定性について検討する．

4.1

シミュレーション条件

図 4.1 に示す Al 単結晶単位格子（原子数 4）を解析対象とした．k 点のサンプリン
グ点数は Monkhorst と Pack のサンプリング手法 (29) に従い全ブリルアンゾーン内に対して 27 点とし，擬ポテンシャルには BHS 型擬ポテンシャル (27) を，交換相関エネル
ギーの表式には Perdew と Zunger による関数形 (26) を使用した．

前章と同様に，予備解析として次の解析を行った．

(1) 0K の格子定数 (0.402nm) において，カットオフエネルギー Ecut を変化させた時の自由エネルギー Efree の変化を計算し，平面波の打ち切りエネルギーを決定する．

(2) 打ち切りエネルギーを (1) より求めた値に固定し，全原子を fcc 構造の格子位置に配置して格子定数 a を変化させ，格子の膨張・収縮に対する系の自由エネルギーの変化を計算する．

解析 (1) より求めた自由エネルギー–カットオフエネルギー関係を図 4.2 に示す．カッ38

トオフエネルギーの増加に伴って自由エネルギーは減少し，Ecut >= 10a.u. で十分収束した．これより，カットオフエネルギーを 10a.u. とした．
解析 (2) より求めた自由エネルギー–格子定数の関係を図 4.3 に示す．自由エネルギーが最小となるのは a =0.397nm のときであり，Al の格子定数の実験値からのずれは1.14%となった．
予備解析により決定したカットオフエネルギー，および，基準格子定数を用いて前章と同様の引張シミュレーション (解析 I および解析 II) を行った．ただし，等方収縮の解析 I における横ひずみの決定は，系の全エネルギーではなく垂直応力 σxx = σyy が0 となるセルを探索することにより行った．

[001]

z

y [010]x [100]
Fig.4.1 Simulation cell

Table 4.1 Calculation condition for analysis
Number of atoms
Lattice constant a
Cutoﬀ energy Ecut
Number of k points
Mesh fort FFT
Pseudo potential
Boundary condition

4
0.397 nm
10 a.u.
3×3×3
32 × 32 × 32
BHS type
Supercell39

Free energy Efree , a.u.

-8.15
-8.2
-8.25
-8.3
-8.35
-8.4
-8.45

0

2

4

6

8

10

12

Cutoff energy Ecut , a.u.
Fig.4.2 Free energy – Cutoﬀ energy curveFree energy Efree , a.u.

-8.3808-8.381-8.3812-8.3814-8.3816-8.3818-8.3820.39 0.392 0.394 0.396 0.398

0.4

0.402 0.404

Lattice constant a , nm
Fig.4.3 Free energy – Lattice constants curve

40

4.2

シミュレーション結果及び考察

4.2.1

理想引張強度と横方向変形分岐

解析 I(等方収縮) より求めた自由エネルギー–引張ひずみ関係，および，引張応力–
引張ひずみ関係を図 4.4，図 4.5 に示す．ひずみの増加とともに応力は単調に増加し，εzz = 0.18 でピーク応力 23.48GPa を示した．従って，等方的収縮を仮定した本解析で得られる Al の理想引張強度は εzz = 0.18，σzz = 23.48GPa となる．
等方収縮の安定構造からアスペクト比 Ly /Lx を変えて計算した自由エネルギーの変化（εzz = 0.050〜0.055）を図 4.6 に示す．図中の破線は，解析 I の，σxx = σyy = 0 の等方収縮の場合の自由エネルギー (Ly /Lx = 1 に相当) を示し，実線が各アスペクト比における自由エネルギーを示している．εzz = 0.051 までは，図 4.6(a)(b) に示すようにアスペクト比を変化させるとエネルギーが上昇し，等方収縮する変形経路が最もエ
ネルギーが低い値をとっていたのに対し，εzz = 0.052 において等方収縮の場合よりもエネルギーが低くなる経路が現れる (図 4.6(c)，Ly /Lx = 1.014 の点)．以降は引張ひずみが大きくなるとともに最もエネルギーの低い点がアスペクト比の大きな側へとシフトした．図 4.7，図 4.8 にこれらのアスペクト比の横方向ひずみにおけるエネルギーおよび引張応力を，それぞれ等方収縮の場合の値とあわせて示す．εzz = 0.052 の点から等方収縮の場合に比べてエネルギー，応力ともに低くなる変形経路への分岐を生じている．また，この不安定分岐を安定限界としたときの理想引張強度は εzz = 0.052，

σzz = 5.65GPa となる．

41

Free energy Efree , a.u.

-8.36

-8.37

-8.38

-8.39

-8.4

-8.41

0

0.04

0.08

0.12

0.16

0.2

0.24

Tensile strain εzz
Fig.4.4 Free energy–Strain curve under tension

Tensile stress σzz , GPa

28
24
20
16
12
8
4
0
0

0.04

0.08

0.12

0.16

0.2

Tensile strain εzz
Fig.4.5 Stress–strain curve under tension

42

0.24

Lx

Free energy Efree , a.u.

-8.4025

Ly

[001]

z

Energy underanisotropic contraction-8.40254

Energy underisotropic contraction-8.40256-8.40258

y [010]x [100]-8.40252

1

1.01

1.02

1.03

1.04

Aspect ratio Ly/Lx(a) εzz = 0.050-8.40238

Free energy Efree , a.u.

Free energy Efree , a.u.

-8.40248-8.4025-8.40252-8.40254-8.40256

1

1.01

1.02

1.03

-8.402457-8.402458-8.40242-8.402459
1

1.01

1.02

1.03

1.04

-8.40244

1

1.01

1.02

1.03

Aspect ratio Ly/Lx

Aspect ratio Ly/Lx(b) εzz = 0.051(c) εzz = 0.052

1.04-8.40225

Free energy Efree , a.u.

Free energy Efree , a.u.

-8.4024-8.40246

1.04-8.40238-8.4024-8.40242-8.40244-8.40246-8.402456

1

1.01

1.02

1.03

-8.40227-8.40229-8.40231-8.40233

1.04

Aspect ratio Ly/Lx(d) εzz = 0.053

1

1.01

1.02

1.03

Aspect ratio Ly/Lx(e) εzz = 0.055

Fig.4.6 Relationships between free energy and aspect ratio of transverse
strain under anisotropic contraction
43

1.04

Free energy Efree , a.u.

-8.401-8.4015-8.402-8.4025

Isotropic
Anisotropic-8.403
0.045

0.05

0.055

0.06

0.065

Tensile strain εzz
Fig.4.7 Free energy–Strain curve under tension

Tensile stress σzz , GPa

8

7

6

5

Isotropic
Anisotropic

4
0.045

0.05

0.055

0.06

Tensile strain εzz
Fig.4.8 Stress–Strain curve under tension

44

0.065

4.2.2

変形分岐後の電子構造変化

変形分岐後の εzz = 0.052 における (111) 面の価電子密度分布の変化を図 4.9 に示す．等方収縮の変形経路では，図 4.9(a) に示すように，最近接原子間の電子密度が高い部分（図中破線）のうち，A 部は B,C 部に比べ電子密度が高く，B,C 部は等しい電子構造となっている．これは，(001) 面上では，等方収縮により原子間距離が減少し，(100),(010) 面上では延伸しているためである．横方向変形のアスペクト比 Ly /Lx を増加（x 方向に収縮，y 方向に延伸）させる非等方収縮変形下では，図 4.9(b) に示すように，(001) 面上の A 部はほとんど変化せず，B 部で電子密度が減少し，C 部で増加する
ような構造に変化した．図 4.10 の (111̄) 面の価電子密度分布では B 部で増加，C 部で減少し，図 4.11 の (11̄1) 面では B 部で増加，C 部で減少している．以上の電子状態の変化を図 4.12 に模式的に示した．太線が電子密度が増加した最近接原子間，波線が減少したそれを表している．なお，二重線は Poisson 収縮により原子間隔が減少している (001) 面上のものであり，引張によって延伸した太線と波線に比べ大きな変化はな
い．引張下で延伸していた太線の最近接原子間の距離が戻る方向に，波線部分はより広がる方向へ変化することにより，図中矢印で示したように部分転位に対応する原子移動を生じることが示唆される．

45

[001]zx [100]y [010]ABC(a) εzz = 0.052 (Isotropic)(b) εzz = 0.052 (Anisotropic)

Fig.4.9 Distribution of valence electron density on (111) plane

[001]

z
x [100]y [010]B(a) εzz = 0.052 (Isotropic)

C(b) εzz = 0.052 (Anisotropic)

Fig.4.10 Distribution of valence electron density on (111̄) plane46

[001]zx [100]y [010]CB(a) εzz = 0.052 (Isotropic)(b) εzz = 0.052 (Anisotropic)

Fig.4.11 Distribution of valence electron density on (11̄1) plane

A

A

A
A

A

A

A
C

[001]

z

A
C
A

A

x [100]Ay [010]

Fig.4.12 Schematic of change in valence electron density

47

A4.2.3格子不安定に関する考察

前章 3.3.3 節と同様に，系の安定性を 6 × 6 の弾性剛性係数行列の正値性により評価する．変形分岐前後の格子不安定性を評価するため，解析 I の等方収縮の εzz = 0.035 ∼0.075 の各平衡状態における弾性係数および弾性剛性係数を計算した．表 4.2 に得られた弾性係数および弾性剛性係数の値を，変形分岐前後の εzz = 0.051 と εzz = 0.052 の
場合について示す．また，εzz = 0.035 ∼ 0.075 における detBIJ の変化を図 4.13 に示
す．detBIJ は εzz の増加によって減少し，εzz <= 0.051 では detBIJ は正値をとり安定であるが，εzz = 0.052 となった時に detBIJ が負となる．この不安定は図 4.14 に示すように，行列式 (3.2) の中の横方向変形のバランスを表す主小行列式

B11
B12

B12
B11

の値

が負となることによりもたらされる．このように，fcc 金属の Al についても，格子不安定性の観点から εzz > 0.051 が横方向変形に対する分岐クライテリアであることが示される．

Table 4.2 Elastic constants and elastic stiﬀness coeﬃcientsεzz = 0.051 εzz = 0.052
εzz = 0.051 εzz = 0.052
C11 (GPa)
225.6
225.4
B11 (GPa)
225.6
225.4
C12 (GPa)
224.3
225.5
B12 (GPa)
224.3
225.5
C33 (GPa)
235.7
235.5
B33 (GPa)
242.2
242.2
C13 (GPa)
117.9
115.7
B13 (GPa)
111.4
109.0
C44 (GPa)42.2
42.2
B44 (GPa)42.2
42.2
C66 (GPa)24.8
34.6
B66 (GPa)
31.6
31.3

48

0.30
Magnitudes of 6x6 determinant,detBIJ, GPa6

Magnitudes of 6x6 determinant, detBIJ, GPa6x 1012

0.25
0.20
0.15
0.10

0.0375x1012

0.02500.01250.00-0.0125
0.046

0.05

0.048 0.05 0.052 0.054
Tensile strain εzz

0.00
-0.05
0.03

0.04

0.05

0.06

0.07

0.08

Tensile strain εzz

Fig.4.13 Relationship between tensile strain and crystal instability under [001]tension (isotropic contraction)
x104

Magnitude of 2x2 determinant,detBIJ (I, J = 1,2), GPa2

0.4
0.3
0.2
0.1
0.0
-0.1
-0.2
0.046

0.048

0.05

0.052

0.054

Tensile strain εzz
22
under [001] tension− B12
Fig.4.14 Relationship between tensile strain and B11(isotropic contraction)

49

4.3

結言

fcc 金属である Al 単結晶について，前章と同様に [001] 方向引張変形下での横方向非等方収縮への変形分岐条件，および，変形分岐前後の格子不安定性について検討した．
得られた結果は以下のとおりである．

(1) 横方向の等方収縮を仮定して得られる理想引張強度 εzz = 0.18 に対し，εzz = 0.051
を臨界ひずみとする非等方収縮がエネルギー的に低くなる変形経路への分岐が示された．

(2) 電子構造の変化から，(1) の非等方収縮変形によって，部分転位に対応する方向への原子の移動が生じることが示された．

(3) 変形分岐ひずみ εzz = 0.052 以降は弾性剛性係数行列の行列式が負となる不安定状態であることが示された．前章同様，この不安定は横方向変形に対する抵抗を22表す小行列式 B11− B12が負となることによってもたらされていた．このように，

fcc 金属である Al についても，格子不安定性の観点から εzz > 0.051 が横方向変形に対する分岐クライテリアであることが示された．

50

第5章
結論
第一原理分子動力学法による変形破壊解析は，膨大な計算を必要とするため，必然的に変形経路を仮定したものとなる．このため，それより得られる理想強度と自由度の高い実際の系で生じる他の変形経路への分岐との関係を明らかにする必要がある．
本研究では，格子不安定条件に基づく変形分岐解析を第一原理分子動力学法により非経験的に行うことを目的とする．その第一歩として，引張変形下での転位発生クライテリアとして重要と報告されている，横方向非等方収縮への変形分岐について検討した．以下に本論文で得られた結果を総括する．
第 3 章では，[001] 方向 (z 軸方向) に単軸引張を受ける Si 単結晶について検討した．
等方的な Poisson 収縮を仮定して得られる理想引張強度は εzz = 0.25 であったのに対し，非等方収縮を考慮した場合，εzz = 0.094 以降の引張ひずみにおいて，非等方収縮の方がエネルギー，応力ともに低くなる変形経路への分岐が存在することがわかった．
分岐臨界ひずみ εzz = 0.093 における ⟨110⟩ 方向の分解せん断ひずみは εzx = 0.295 であり，垂直方向応力 0 の条件下における (111) 面の理想せん断強度 εzx = 0.30 に近い．
また，非等方収縮変形下の電子構造変化からも，部分転位に対応する原子の移動を示
唆する Si 間のボンドの変化が認められた．これらより，本解析における引張時の横方向変形分岐は，転位発生に対する理想せん断強度と密接に関係していることが示唆された．変形分岐に関して格子不安定性の観点から考察するため，各ひずみ下の弾性剛性係数を数値的に評価したところ，εzz = 0.094 以降は弾性剛性係数の行列式が負となる不安定状態になっていることが示された．この不安定は，横方向変形に対する抵抗

51

2
2
を表す小行列式 B11− B12が負となることによってもたらされていた．このように格子不安定性の観点からも εzz > 0.093 が [001] 引張を受ける Si の横方向変形に対する分岐クライテリアであることが示された．
第 4 章では，[001] 方向に単軸引張を受ける Al 単結晶について，第 3 章と同様の検討を行った．その結果，等方 Poisson 収縮下の理想引張強度 εzz = 0.18 に対して，εzz = 0.052以降の引張ひずみで非等方収縮への変形経路の分岐を生じることが示された．Si と同様，電子構造変化から非等方収縮により部分転位に対応する原子移動を生じることが示された．また，弾性剛性係数の正値性に基づく格子不安定性評価からも εzz = 0.052
以降が不安定であることが示され，fcc 金属である Al 単結晶においても同手法が不安定変形開始の指標となりうることが示された．

52

参考文献
(1) Kitagawa, Aihara, Kawazoe 編, Mesoscopic Dynamics of Fracture, (1998)(2) Abrahan, F. F., Materials Science for the 21st Century; The Society of Meterials
Science, A (2001), 195–202.
(3) Vashishta, P., Bachilechner, M. E., Nakano, A., Campbell, T. J., Kalia, R.
K., Kodiyalam, S., Ogata, S. and Shimojo, F., Materials Science for the 21stCentury; The Society of Meterials Science, A (2001), 213–219.
(4) Li, J., Cai, W., Chang, J. and Yip, S. , Materials Science for the 21st Century;
The Society of Meterials Science, A (2001), 220–233.
(5) 香山正憲, 田中真梧, 材研連講演論文集, (2000), 342–343.
(6) 尾方成信, 北川浩, 金属学会誌, 60–11, (1996), 1079–1083.
(7) 梅野宜崇, 北村隆行, 日本計算工学会第 6 回講演論文集, 6–2, (2001), 571.
(8) 楠幸久, 渋谷陽二, 冨田佳宏, 神戸大学修士論文, (1999)(9) Born, M and Huang, K., Dynamical Theory of Crystal Lattices, (1954), Oxford
UP.
(10) Huang, K., Proc. R. Soc. London, 204A (1950), 178.
(11) Milstein, F., Phys. Rev. B, 3–4B (1971), 1130.
(12) Milstein, F. and Hill, R. , J. Mechanics and Physics of Solids, 25 (1977), 457.
(13) Milstein, F. and Hill, R. , J. Mechanics and Physics of Solids, 26 (1978), 213.
(14) Milstein, F. and Farber, B. , Phys. Rev.Lett., 44–4, (1980), 277.
(15) Hill, R., Mathematical Proceedings of the Cambridge Philosophical Society, 77,(1975), 225.
(16) Milstein, F. and Hill, R. , J. Mechanics and Physics of Solids, 23 (1975), 21.
53

(17) Wang, F. et al., Phys.Rev.lett., 45–14, (1980), 1196.
(18) Wang, J., S.Tip, S.R.Phillpot and Wolf, D., Phys. Rev.Lett., 71–25 (1993), 4182.
(19) Wang, J., Li, S.Tip, S.R.Phillpot and Wolf, D., Phys. Rev., 52–17B (1995),
12627.
(20) Wallace, D. C., Thermodynamics of Crystals, (1972), Wiley, Newyork(21) 屋代如月, 冨田佳宏, 日本機械学会第 12 回計算力学講演会, (1999)(22) 屋代如月, 北村隆行, 大谷隆一, 京都大学博士論文, (1998)(23) 屋代如月, 北村隆行, 大谷隆一, 機論, 63–615, (1997), 2380–2386.
(24) Hohenberg,P. and Kohn,W., Phys. Rev.,, 136–3B, (1964), 864–871.
(25) Kohn,W. and Sham,L.J., Phys. Rev.,, 140–4A, (1965), 1133–1138.
(26) Perdew,J.P. and Zunger,A., Phys. Rev.,B, 23–10, (1981), 5048–5079.
(27) Bachelet, G. B., Hamann, D. R. and Schlüter, M., Phys. Rev.,B, 26–8, (1982),
4199–4228.
(28) 上田顕, コンピュータシミュレーション, (1990), 朝倉書店.
(29) Monkhorst,H.J. and Pack,J.D., Phys. Rev.,B, 13–12, (1976), 5188–5192.
(30) Chadi,D.J.and Cohen,M.L., Phys.Rev.B., 8-12(1973), 5747-5753.
(31) 香山正憲, 固体材料の電子状態の基礎, (1993), 材料学会勉強会資料.
(32) L.D. ランダゥ , E.M. リフシッツ著, 佐々木健, 好村磁洋訳, 量子力学 1, (1992),
東京図書.

(33) 尾口多美夫, 佐々木泰造, 固体物理, 25–11, (1990), 857–865.
(34) キッテル, 固体物理学入門 上 第 7 版 (1998), 丸善(35) Car,R. and Parrinello,M., Phys. Rev.Lett., 55, (1985), 2471–2474.
(36) Teter,M.P. et al., Phys. Rev.B, 40, (1989), 12255. ; Payne,M.C. et al.,Rev.Mod.Phys., 64, (1992), 1045.
54

(37) Bylander,D.M. et al., Phys. Rev.B, 42, (1990), 1394. ; Kresse,G. and Hafner,J,
Phys. Rev.B, 49, (1994), 14251.
(38) King-Smith,R.D. and Vanderbilt,D. , Phys. Rev.B, 49, (1994), 5828.
(39) Gillan,M.J. , J. Phys. Condens. Matter 1, (1989), 689.
(40) Stich,I. et al. , Phys. Rev.B, 39, (1989), 4997.
(41) 吉田貞史 電子技術総合研究所彙報, 62–10,11, (1999),137–147(42) Rignanese,G.M. , Ghosez,Ph. , Charlier,J.C. , Michenaud,J.P. and Gonze,X.,
Phys. Rev.B, 52, (1995), 8160.
(43) Bernasconi,M. , Chiarotti,G.L. , Focher,P. , Scandolo,S. , Tosatti,E. and Parrinello,M. , J.Phys.Chem.Solids 56, (1995), 501.
(44) Caperley,D.M. and Alder,B.J., Phys. Rev.Lett., 45–7, (1980), 566–569.
(45) Kerker,G.P., Phys. Rev.,B, 23–6, (1981), 3082–3084.
(46) Pui,K.L. and Marbin L.C., Phys. Rev.,B, 24–8, (1981), 4224–4227.
(47) Sasaki,T. , Rappe,A.M. and Louie,S.G, Phys. Rev.,B, 52–17, (1995), 12760–
12765.
(48) Nielsen,O.H.,Martin,R.M., Phys. Rev.,B, 32–6, (1985), 3780–3791.
(49) Dacosta,P.G.,Nielsen,O.H.and Kunc,K., J. Phys. C : Solid State Phys.,, 19,(1986), 3163–3172.

55

付録 A
原子単位系
従来の単位系を用いると，電子の質量は me = 9.1095 × 10−31 kg，原子の大きさの目安となるボーア半径は a0 = 5.2918 × 10−11 m となり，原子または分子を扱う場合に大変不便である．このため，原子単位 (atomic unit; a.u.) とよばれる単位系が用いられる．今，水素原子に対する Schrödinger 方程式は，
}

{

e2
h̄2 2
Ψ = EΨ
∇ −−2me
4πε0 r(A.1)となる．ここで h̄ は，Plank の定数 h を 2π で割ったものであり，me は電子の質量，−e は電子の電荷である．この方程式を無次元化するために，
x, y, z, −→ λx′ , λy ′ , λz ′(A.2)

とすると，式 (A.1) は，
{

}

h̄2
e2
′2
−
∇
−
Ψ′ = EΨ′
2me λ24πε0 λr′(A.3)

運動エネルギ，およびポテンシャルエネルギ演算子の前の係数は，λ を，

h̄2
e2
=
= Ea
m e λ24πε0 λ

(A.4)

を満たすように置くと，くくり出すことができる．ここで Ea はハートリー（hatree）と呼ばれるエネルギの原子単位である．また，式 (A.4) を λ について解くと，

λ=4πε0 h̄2= a0
m e e256(A.5)

となり，λ は Bohr 半径 a0 に等しく，長さの原子単位である．式 (A.3)，式 (A.4) に
より，

{

}

1
1
− ∇′2 − ′ Ψ′ = EΨ′
2
r

Ea

(A.6)

となるので，E ′ = E/Ea とすると，無次元の方程式，
}

{

1
1
− ∇′2 − ′ Ψ′ = E ′ Ψ′
2
r

(A.7)

を得ることができ，これが，原子単位系での Schrödinger 方程式である．表 A に原子単位系と SI 単位系との関係を示す．
物理量
長さ
質量
時間
エネルギー
電荷
力
応力

a.u. (Hartree)
1.0
1.0
1.0
1.0
1.0
1.0
1.0

SI 単位
5.2918 × 10−11 m
9.1095 × 10−31 kg
2.4189 × 10−17 sec
4.3598 × 10−18 J
1.6022 × 10−19 C
8.2388 × 10−8 N
2.9417 × 1013 Pa

また，SI 単位系ではないが，よく用いられる物理量として，長さの単位では，1[a.u.] =0.52981[A]，エネルギの単位では，1[a.u.] = 27.211[eV] などがある．

57

付録 B
関連講演論文(1) 大穂正史，屋代如月，冨田佳宏，第一原理計算によるシリコンの理想格子不安定解析，日本機械学会第 14 回計算力学講演会講演論文集，(2001)(2) 大穂正史，屋代如月，冨田佳宏，第一原理計算による Al 単結晶の理想格子不安
定解析，日本機械学会関西支部第 77 期定時総会講演会講演論文集，(2002)，講
演予定

58

謝辞
本研究を進めていくにあたって終始懇切丁寧な御指導を賜った冨田佳宏教授に深く感謝致します．また，副査として貴重な御助言をいただいた大前伸夫教授，田中拓助
教授，研究全般にわたり直接御指導頂いた屋代如月助手，研究内容について様々な御助言を頂いた安達泰治助教授，研究室の環境作りにご尽力くださった古宇田由夫技官に深く感謝致します．
さらに，研究室に配属された当初から，倦むことなく基礎から御指導くださった楠幸久氏（現・(株) 住友電気工業）をはじめとする諸先輩方に感謝の意を示すともに，共に励ましあい助け合いながら研究生活を過ごした，同輩の河野雄二君，佐藤克也君，澤田幸秀君，内藤正登君，ならびに，同じ研究分野で有益な議論を交わした山上勝也君ら講座内諸氏に心から御礼申し上げます．
最後に，いつ如何なる時も常に支え励まして下さった家族に最大限の謝意を示します．ありがとうごさいました．

59