要約Ni 基超合金の γ/γ ′ 界面の力学特性解明につながる新たな知見を得ることを目的として，Ni および Ni3 Al それぞれ単結晶について，[001] 方向引張 · 圧縮，静水圧引張 ·圧縮を第一原理分子動力学法により行い，応力–ひずみ関係を求めるとともに，各ひず
み下の格子不安定性について，弾性剛性係数の正値性から議論した．

[001] 方向の引張では，Ni，Ni3 Al のいずれもひずみ 0.3 近傍で応力–ひずみ関係にピークが現れるが，それよりはるかに小さいひずみ (0.1 程度) で格子不安定となることが示された．その不安定は Born 条件によりもたらされており，横方向の非等方変形への分岐点であることがわかった．格子不安定となるひずみを弾性限界とした場合のNi と Ni3 Al の理想引張強度はそれぞれ ε33 =0.100，σ33 =13.5GPa，および，ε33 =0.080，
σ33 =10.8GPa となる．なお，横方向の Poisson 収縮を拘束した場合，格子不安定となるひずみは高ひずみ側に，応力–ひずみのピークは低ひずみ側にシフトするが，Ni とNi3 Al の強度の大小関係は変わらない．
[001] 方向の圧縮では，横方向の応力を零とした場合は fcc または L12 格子中の bct 構造が bcc，B2 構造となる Bain の関係により応力–ひずみ曲線が 3 次関数的なものとなった．応力–ひずみ曲線の極小点（圧縮を正としたときのピーク）は，格子不安定となるひずみと一致した．その不安定は Spinodal 条件によりもたらされており，弾性剛性係数の変化からは結晶の相変態に対するものであることが示唆される．これより得られる理
想圧縮強度は Ni が ε33 =–0.100，σ33 =–5.44GPa，Ni3 Al が ε33 =–0.120，σ33 =–5.94GPaとなり，上記の引張と強度関係が逆転する．一方，横方向の Poisson 収縮を拘束した場
合は応力–ひずみ曲線は単調減少し極小値を持たなかったが，Ni が ε33 =–0.150，Ni3 Alが ε33 =–0.140 で Spinodal 不安定となった．
静水圧引張および圧縮下では，応力–ひずみ関係は横方向変形を拘束した [001] 引張および [001] 圧縮と似た傾向を示した．しかしながら，引張において格子不安定となる
ひずみは応力–ひずみのピークに一致し，かつ Born 条件ではなく Spinodal 条件によるものであった．弾性剛性係数の変化から，本不安定条件は膨張に対する不安定ではなく，相変態に対するものであることが示唆される．本条件による理想静水圧引張強度は Ni が ε33 =0.150，σ33 =30.0GPa，Ni3 Al が ε33 =0.150，σ33 =28.3GPa と評価される．
なお，圧縮側は格子不安定となるひずみは存在せず，応力–ひずみ関係も極値を持たない．したがって，完全に等方な圧縮条件下では結晶の変形抵抗は増加し，系は変形の増加とともにより安定となる．ただし，横方向の Poisson 収縮を拘束した [001] 圧縮で示したように，実際には等方性が崩れて格子不安定となるひずみが現れる．

SUMMARY
Ab initio molecular dynamics study for two single crystals, Ni and Ni3 Al, subjectedto the [001] tension/compression and hydrostatic tension/compression have been done.
The “ideal strength” is discussed from the view point of the lattice instability evaluated by the positive deﬁniteness of the elastic stiﬀness matrix as well as the peakof the stress–strain curve. The stress–strain curve shows a peak near at the strain
ε33 =0.3 for both of Ni and Ni3 Al under [001] tension. The crystals, however, become
unstable at far smaller strain about ε33 =0.1. The instability is due to the “Born condition” which leads the deformation bifurcation to the anisotropic transverse contraction.
The ideal strengths are evaluated as ε33 =0.100，σ33 =13.5GPa for Ni and ε33 =0.080，
σ33 =10.8GPa for Ni3 Al based on the lattice instability criteria. It is also suggested
that the criteria for lattice instability shifts to the higher strain while that of peakdecreases when the deformation is restricted in the transverse direction. However, therelationship of the strength between Ni and Ni3 Al is unaltered. In the case of [001] compression with stress-free condition in the transverse direction, the stress–strain curvesof both Ni and Ni3 Al become cubic line because of the Bain’s crystallographic relationship. The minimum of the curve, or the peak of compressive stress, coincides with thecritical strain for lattice instability. The instability is caused by the “Spinodal condition” which leads the crystallographic transformation. Thus the ideal strength for the[001] compression are evaluated as ε33 =–0.100, σ33 =–5.44GPa for Ni and ε33 =–0.120,
σ33 =–5.94GPa for Ni3 Al, reversing the strength relationship. On the other hand, thestress monotonically decreases without minimum when there is the deformation constraint in the transverse direction. The Spinodal condition, however, decides the ideal
strengths of [001] compression with transverse constraint as ε33 =–0.140, σ33 =–49.5GPafor Ni and ε33 =–0.130, σ33 =–39.9GPa for Ni3 Al. The stress–strain curves under thehydrostatic tension/compression resemble that of the [001] tension/compression withthe transverse constraint. The critical strain for lattice instability in tension, however, coincides with that of the peak of the curve. The Spinodal condition causes theinstability corresponding to the onset of crystallographic transformation. The ideal
strengths for the hydrostatic tension are evaluated as ε33 =0.150, σ33 =30.0GPa for Ni
and ε33 =0.150, σ33 =28.3GPa for Ni3 Al. There is no peak nor lattice instability strain
under the hydrostatic compression.

修士論文
第一原理分子動力学法を用いた

Ni 及び Ni3Al の理想格子不安定性解析

平成 16 年 2 月
神戸大学大学院 自然科学研究科
機械工学専攻023T420N
山上勝也

Ab initio Molecular Dynamics Studyon Lattice Instability of Ni and Ni3Al

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

Katsuya YAMAGAMI目次
第1章 緒 論1

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

5

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

6

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

6

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

8

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

9

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

10

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

14

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

15

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

16

2.8.1

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

17

2.8.2

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

22

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

23

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

24

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

25

第 3 章 理想格子不安定性解析の概要

30

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

30

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

31

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

33

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

36. . . . . . . . . . . . . . . . . .

第 4 章 平衡状態における物性評価

38

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

i

38

目次

ii

4.1.1

モデルおよび予備解析 . . . . . . . . . . . . . . . . . . . . . . .

38

4.1.2

平衡状態における格子定数，弾性定数の評価 . . . . . . . . . . .

41

4.2 解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

4.2.1

カットオフエネルギーの決定 . . . . . . . . . . . . . . . . . . .

42

4.2.2

k 点数の決定 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

43

4.2.3

平衡状態における格子定数および弾性定数 . . . . . . . . . . . .

44

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

46

第 5 章 [001] 方向単軸引張変形下の理想格子不安定性解析

475.1 解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

5.2 解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

5.2.1

自由エネルギー–引張ひずみ，引張応力–引張ひずみ関係. . . .

505.2.2

格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

52

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

56

第 6 章 [001] 方向単軸圧縮変形下の理想格子不安定性解析

58

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

58

6.2 解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

6.2.1

自由エネルギー–圧縮ひずみ，圧縮応力–圧縮ひずみ関係. . . .

596.2.2

格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

62

6.2.3

圧縮変形における結晶構造関係 . . . . . . . . . . . . . . . . . .

66

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

68

第 7 章 静水圧変形下の理想格子不安定性解析

70

7.1 解析条件 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

70

7.2 解析結果及び考察 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

71

目次

iii

7.2.1

自由エネルギー–ひずみ関係，応力–ひずみ関係. . . . . . . . .

71

7.2.2

格子不安定性 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

7.3 結言 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

第8章 結 論

77

参考文献

80

第 A 章 原子単位系

84

第 B 章 関連発表論文 · 関連講演論文

86

謝辞

87

第1章
緒論
理想強度とは，転位や粒界といった欠陥を一切含まない理想的な均一結晶に対して変形を加えた場合に，弾性限界を超えて塑性変形が開始する最小のひずみまたは応力である (1) ．材料の無欠陥化を進めていった時に到達する限界強度であるとも言える．
電子デバイス分野などでは，電子特性の改善の観点から転位等の格子欠陥が極力排除されている．無転位のシリコンウェハから界面転位が発生する臨界応力は，理想強度に近いとの報告もある (2) ．このような欠陥生成によるデバイスの信頼性低下予測のみ
ならず，その材料が本来持つ強度の指標として，理想強度は材料設計開発に重要な意
味を持つ．
ここ十数年の飛躍的な計算機能力の向上を背景に，古典的分子動力学法や第一原理
分子動力学法等の計算機シミュレーションによる原子レベルからの材料設計への期待が高まっている．原子間相互作用を簡便なポテンシャル関数で表す古典的分子動力学
法では，大規模なものでは数億オーダーの原子数で材料の変形 · 破壊現象をシミュレー
トすることもできる (3)−(5) ．その反面，用いているポテンシャル関数が平衡状態の材料
物性を表現するようにフィッティングされているため，平衡状態から大きく異なる状態 (例えば粒界 · 界面部や，大変形時) における正当性は保証されない．一方，第一原
理分子動力学法は，原子間相互作用を仮定することなく，原子種と原子配置のみを必要情報とし，量子力学に基づいて原子間相互作用を厳密に評価できる手法である．そのため，電子状態が単結晶 bulk 状態のそれと大きく異なっていると考えられる粒界や異種材料界面の特性評価も行われている．例えば，香山らは SiC 粒界，SiC/Ti 界面の1

第一原理引張試験を行い (6) ，界面の理想引張強度を検討している．この理想強度は実際の材料と比較して二桁近く大きいが，これは破壊の起点となるクラック先端の局所的な臨界応力値に対応すると報告している．尾方らは AlN/Al 界面の第一原理引張試
験を行い (7)(8) ，破壊が界面に隣接する Al 層間より起こることや，そのときの引張応
力が 18GPa 程度であることなどを報告している．しかしながら，第一原理分子動力学
法では，逐次電子状態を計算しながら原子運動を追跡するために計算量は膨大であり，これらの理想強度は少数の原子で解析した応力–ひずみ曲線のピーク値を用いて評価されたものである．そのため，周期境界条件に強く影響を受け，より実際に近い多数原子の解析において生じる他の変形経路への不安定分岐を表現できておらず，冒頭で述べたような理想強度をより高強度側に過大評価する (9) ．
少数の原子の解析でも，巨視的な材料強度に重要な知見を与えるものに格子不安定性解析がある (10)−(13) ．格子不安定とは，外力下で変形している結晶格子が釣り合いを失い，外力の増加を必要とせずに変形が不安定に進行する現象を指している．結晶の安定限界を評価する試みは Born によって行われ，「Born の不安定基準」を呼ばれる，弾性係数マトリックスの正値条件が提案された (10) ．一方，Born の基準は無負荷平衡点におけるものであるとして，有限変形下の結晶安定性を評価する試みが，単純な原子間ポテンシャルを仮定した解析により 1970 年代に活発になされた (14)−(16) ．近年
では，少数の原子しか扱えないが，大変形時にも原子間相互作用を精密に評価できる
第一原理分子動力学法を用いた格子不安定性解析が試みられている (17)−(19) ．例えば，Milstein らは fcc または bcc 構造のアルカリ金属に単軸引張変形を与える解析を行い，格子不安定となる臨界ひずみ以降では変形経路の分岐が起こり，bcc↔fcc 相変態現象が生じることを報告している (17) ．Krenn らは，W の理想せん断強度を評価し，ウィ
スカーのナノインデンテーション試験より得られた強度とそれに対応がみられることから，理想格子不安定性解析が実材料の理想強度評価法として有効であることを示した (18) ．著者らは，[001] 方向単軸引張を受ける Si および Al 単結晶の格子不安定性解析を行い (19) ，応力–ひずみ関係がピークを示すひずみよりよりもはるかに小さなひずみ
で系は Born 条件による不安定点に達していることを明らかにした．弾性剛性係数の変2化から，Born 条件の不安定は横方向の非等方変形への経路分岐に対応するものである
ことを示し，ピークより評価した Si と Al の硬度等の大小関係は実験値と矛盾するが，Born 条件を弾性限界とした場合にはその矛盾が解消されることを示した．また，電子状態の詳細な観察から，格子の横方向崩壊により転位が発生する可能性があることを示唆した．
本研究では，Ni と Ni3 Al 単結晶について，種々の変形下の格子不安定性を第一原理
分子動力学法により検討する．Ni と Ni3 Al を解析対象としたのは，Ni 基耐熱超合金のγ/γ ′ 界面構造を想定している．実用高温構造用材料として用いられている Ni 基超合
金は, 母相である fcc 構造の Ni(γ) 相中に L12 構造の Ni3 Al(γ ′ ) 相がサブミクロンオー
ダーで格子状に析出した微視構造を有しており，無数に存在する Ni/Ni3 Al 界面が材料の力学特性に重要な役割を果たしていると考えられている．このため，γ ′ 相の体積分率 (20) ，γ ′ 相形状および大きさ (21) ，格子ミスフィット (22) 等に関する実験的検討や，古典的分子動力学法等による界面近傍の転位挙動に関するシミュレーションによる検討(23)(24)などが数多くなされているが，微視的な界面そのものの構造については不明な点が多く残されている．
界面近傍の力学特性の解析を精密に行うには，異種原子間の相互作用を電子状態から計算する第一原理分子動力学法を用いるべきである．しかしながら，界面の変形挙動を直接モデル化するには多数の原子が必要であり，第一原理分子動力学法では計算量が膨大となって現実的ではない．そこで本研究では，Ni 単結晶および Ni3 Al 単結晶の格子不安定条件を第一原理的に評価することにより，Ni 基超合金の力学特性解明のための新しい知見を得ることを目的とする．本論文は以下のように構成される．
第 2 章では，第一原理分子動力学法の基礎理論，および，電子状態計算の高速化手法について説明する．第 3 章では，格子の安定性の概念を説明し，弾性剛性係数を用いた格子不安定性解析の方法について述べる．第 4 章では，Ni および Ni3 Al それぞれの単結晶について，第一原理計算から平衡状態の物性を評価し，それを実験値と比較することによって本研究で得られる解析結果の精度を確認する．第 5 章では，[001] 方
向引張を受けるそれぞれの単結晶の理想強度を求めるとともに，格子不安定となる臨

3

界ひずみを明らかにする．このとき，横方向の変形に等方的な Poisson 収縮を仮定した場合と，横方向変形を拘束した場合を考慮し，それらの変形拘束が理想強度ならびに格子不安定性に及ぼす影響についても検討する．第 6 章では，[001] 方向圧縮を受けるそれぞれの単結晶について，第 5 章と同様に，変形拘束が理想強度や格子不安定性に及ぼす影響について検討する．第 7 章では，静水圧引張または圧縮のひずみを与えた場合の理想強度，格子不安定性について検討する．第 8 章では，本研究で得られた結果を総括する．

4

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

5

2.1

断熱近似と平均場近似

通常，我々が扱う系は多数の原子核と電子からなる集合体である．そして電子間，原子核間，および電子と原子核との間の相互作用は多体問題であり，一般的に解くことができない．このような複雑な問題を実際に解くことが可能な問題へと帰着するために，通常，以下の 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)と表せることを明らかにした (25) ．右辺の各項はそれぞれ，原子核による電子のポテン
シャルエネルギー，相互作用する多電子系での電子の運動エネルギー，電子間クーロ

6

ン相互作用エネルギー，他の全ての電子間多体相互作用を表す交換相関エネルギーである．この Etot を最小にする ρ(r) が基底状態での電子密度分布となる．
相互作用のない系での電子の状態を表す波動関数 (電子波動関数) を ψi とし，その運動エネルギー Ts を

Ts [ρ(r)] =

occ
∑
i

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

(2.2)

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

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

(2.4)

∑occ

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

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

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

(2.5)

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

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

(2.6)

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

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

(2.7)

i

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

7

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

εxc (r) = εx + εc

(2.9)

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



 −

0.1423(rs >√= 1)
εc (r) =  1 + 1.0529 rs + 0.3334rs
−0.0480 + 0.0311 ln rs − 0.0116rs + 0.0020rs ln rs (rs <
= 1)
ここで，

(

rs =

3 1
4π ρ

(2.11))13(2.12)

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

(2.13)

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

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

)]

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

となる．

2.4

逆格子空間

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

(n1 , n2 , n3 は整数)

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

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

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

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

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

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

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

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

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

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

k=

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

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

10

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

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

=⇒

∑
G′

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

(2.27)

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

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

= f (G − G′ )(2.28)

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

(2.29)を用いる．

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

11(2.30)

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

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

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

∑∑
R

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

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

∫

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

VaPP,loc (G) =

∫

= 2π
=

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

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

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

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

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

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

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

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

= 4π

∑

∫

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

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

ρ(G) exp[iG · r]

(2.36)

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

(2.37)

ρ(r) =

G

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

ρ(G) exp[iG · r]

(2.38)

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

Vcoul (r) = 4π

∑ ρ(G)
G|G|2

exp[iG · r]

(2.39)

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

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

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

|G |

G′

= 4π

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

ρ(G)
|G|2

(2.40)

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

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

(2.41)

となる．

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

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

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

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

2.6

系のエネルギー

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

BZ ∑
occ
∑
k

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

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

]

+
−

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

√

a

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

(2.46)である．
ここで∑

ρ(r) =

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

(2.47)

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

+

BZ
∑

fk

occ
∑
n

k

∑

fn

∑∑

GG

′

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

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

Etot =

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

+

BZ
∑
k

fk

occ
∑
n

fn

∑ ∑
G=0 G′ =0

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

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

2.7

応力

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

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

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

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

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

=−
−

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

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

+

+ δαβ

∑

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

G

+

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

2.8

擬ポテンシャル法

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

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

2.8.1

TM 型擬ポテンシャル

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

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

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

17

(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) が次式で表される．

AE Vl (r)PP
Vscr,l(r) = 

′

′′

2

(r >
= rcl )

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

∫ rcl
0

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

より，[∫ r

2c0 + ln

]

cl

r
0

2(l+1)

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

[∫ r
0

cl

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

P(rcl )p(rcl ) = ln
rcll+1P′ (rcl ) l + 1
p′ (rcl ) =−P(rcl )rcl
2(l + 1) ′2
p (rcl ) − [p′ (rcl )]rcl
2(l + 1) ′p′′′ (rcl ) = 2VlAE′ (rcl ) +p (rcl )rcl2
2(l + 1) ′′−p (rcl ) − 2p′ (rcl )p′′ (rcl )rcl
4(l + 1) ′p′′′′ (rcl ) = 2VlAE′′ (rcl ) −p (rcl )rcl2
4(l + 1) ′′2(l + 1) ′′′+p (rcl ) −p (rcl )2rclrcl
p′′ (rcl ) = 2VlAE (r) − 2εl −−2 [p′′ (rcl )] − 2p′ (rcl )p′′′ (rcl )218(2.60)(2.61)(2.62)(2.63)(2.64)

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

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

∑

PP
Vnloc,l(r)P̂ll

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

19(2.67)

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

PP
PP
(r)|
(r)ψlPP (r) >< ψlPP (r)Vnloc,l|Vnloc,l
P̂l
PP
PP
PP
< ψl (r)|Vnloc,l (r)|ψl (r) >(2.68)

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

×

ΩCl

{∫ ∞
0
l
∑

×

0

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

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

∑∑
R

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

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

20

(2.74)

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

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

m=−l

×


∑


G′

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

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

CAYnalkm (k + G) =

∑

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

21(2.77)

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

′
VNL
=

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

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

(2.81)

ただし

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

また，< Ψi |βj >= δij
22(2.83)

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

Qij =

r<rcl

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

S =1+

∑

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

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

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

∑

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

ij

となる．

2.9

電子占有数

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

[

(

1
εi − εFfi =1 − erf2σ23)](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)′n∗
Ck+G′ exp[−iG · r]

(2.96)

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

(2.97)

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

24

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

2.11

電子系の最適化手法

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

共役勾配法の原理
共役勾配法は，一般には正定な係数行列をもつ連立１次方程式を最適化の考えに立って解くために，あるいは，多次元空間の２次関数 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

25

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

(2.99)

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

∑∑
G

G′

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

Rnk = −

′
∂Etot
m∗
∂Ck+G

]

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

k+G′

G′


,

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

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

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

∂ε′kn
m∗
∂Ck+G

]i

(

)

i
= − H − λin I · C nk

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

[

]

n,i
i
i
i
C nk
= Ck+GH = [Hk+G,k+G′ ] , λin =< ΨnkĤ Ψnk>′ ,(2.103)である．これは，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.104)m̸=n<preconditioning>′i
残差ベクトル Rnk
に対して preconditioning という処理を施す．大きな逆格子ベクトルについては平面波の運動エネルギーが大きくなるが，このことが残差ベクトルに影響して収束性を悪化させる．preconditioning は，この問題を回避して収束を速めるた
i
めに行われる．preconditioning された残差ベクトルを Gnk
とすると
i
i′
Gnk
= K i · Rnk(2.105)と表される．ここで(27 + 18x + 12x2 + 8x3 )(27 + 18x + 12x2 + 8x3 + 16x4 )Ekin (G)
x=
i
Ekin
1
Ekin (G) = |k + G|2
2
1
i
i
i
| − ∇2 |Ψnk>
= < Ψnk
Ekin
2
KGG′ = δGG′

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

)

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

∑ (

)

i∗
i
i
C mk
· Gnk
C mk

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

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

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

γi =






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

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

(

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

i′
F nk

(2.113)(2.114)

)1

i′∗
i′ 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 )2
ε11 + ε22−
+ ε12 ε∗12
γ=
2
4

} 21

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

α= {

ε12

}1

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

28

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

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

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

29

第3章
理想格子不安定性解析の概要
格子不安定とは，外力下で変形している結晶格子が釣り合いを失い，外力の増加を必要とせずに不安定に変形が進行する現象を指している．有限変形下の結晶の安定性
は，従来は結晶の変形をブラベー格子の変形で代表することによって系のエネルギーの変数を限定し，エネルギー関数の 2 階微分を解析的に求めることにより評価していた(11)．一方，Wang らは，結晶の変形をひずみで代表させることによって，系の安定性

を弾性剛性係数 (10) の正値性によって評価する手法を提案した (12) ．分子動力学シミュ
レーションによる検証の結果，原子の熱揺動の影響を含んだ結晶の安定性が，系全体の弾性剛性係数で評価できることが示されている．弾性剛性係数による評価は，系のエネルギー関数の表式が求まっていない場合でも，数値的に弾性剛性係数を求めれば安定性評価が可能であるため，第一原理解析でも適用可能である．
本章では，まず従来のエネルギー関数の 2 階微分に基づいて結晶の安定性を評価する手法を説明する．その後，結晶の熱力学関係式から応力と弾性係数 (10) の定義を示
し，非線形弾性変形における応力とひずみの関係を表す弾性剛性係数について説明する．最後に，弾性剛性係数の正値性に基づく安定性評価について説明する．

3.1

不安定条件

結晶の変形を理想化し，すべての結晶格子が外力を受けて均一に変形するものと仮定する．すると fcc を含む立方体格子の変形は図 3.1 に示すような 6 つの格子パラメー

30

タ a1 ∼ a6 で記述され，内部エネルギー U はこれらの関数 U (a1 , a2 , · · · , a6 ) ≡ U ({am })となる．ここで，本節では原子の運動は考慮しないため，U ≈ Etot である．このとき，{am } の変形状態下にある結晶の安定性は，以下のように微小変形増分 {∆am } によるエネルギーの変化を考えることによって求められる (10)(11) ．状態 {am } 近傍での内部エネルギーの Taylor 級数展開はU ({am + ∆am }) = U ({am }) +6∑

Fm ∆am +

m=1

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

Fm =

∂U
,
∂am {am }

Amn =

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

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

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

Fm ∆am =

m=1

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

(3.3)

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

3.2

応力と弾性係数

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

dU = T dS − dW

(3.4)

である (46) ．ここで，U は内部エネルギー，T は温度，S はエントロピ，dW は系が外界になす仕事である．外部応力 σ の負荷によって結晶が変形する際の dW を求めるた
め，結晶内の任意の点 X が応力の負荷によって X + ∆X に変化する均質一様な変形

31

a3
3
4

2
1

5

6

a2
1

Fig.3.1 Unit cell of fcc latticeを考える．変形前の物体表面を S とし，その微小要素を dS とすると，dS において i方向に作用している力 fi は負荷応力 σij を用いて以下のように表せる．

fi = σij dSij

(3.5)

X から X + ∆X への変位勾配テンソルを ∆u とすると，∆Xi = ∆uij Xij(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 のひずみテンソルηij =1(uij + uji + uki ukj )
2
32

(3.9)

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

dηij =

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

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

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

(3.11)

dU = V (X)σij dηij

(3.12)

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

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

(

(

∂U
∂ηij

)

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

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

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

Cijkl = Cjikl = Cijlk = Cklij

3.3(3.16)

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

一般に弾性係数は応力–ひずみ関係の勾配と考えられているが，前節で示したとおり，結晶の熱力学関係式からは，弾性係数は基準配置 X におけるエネルギーの 2 次導33

関数と定義される．応力をひずみと関係づけるためには，2 つの基準配置 X と X̄ における応力 σij と σ̄ij の関係を，X と X̄ 間のひずみに対応させる必要がある．
状態 X̄ から状態 X への変形勾配テンソルを Jij とすると，

∂Xi
∂ X̄j

Jij = Jji =(3.17)

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

nij =

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

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

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

(3.20)

η′

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

1
σij =
V
(

=

(

V̄
V

∂U
∂ηij

)

1
V̄

)
(

η′

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

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

nij + · · ·(3.22)

n′

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

1
V̄

(

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

34

(3.23)

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

σij =

V̄
V

)

{

}

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

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

[

]

[

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

(3.25)

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

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

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

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

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

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

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

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

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

∂σij
∂ηij

)

= Bijkl(3.32)

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

35

3.4

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

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

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

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

)
1( T
BIJ + BIJ
2

(3.34)

SYM
SYMの正値性が系の安定性を支配する (13) ．ここで BIJ
は Bijkl
を Voigt 表記 (46) したものである．また，T は転置行列を表す．
本研究で対象とする [001] 方向の変形では，系の安定性は次の 6 × 6 行列の正値性により評価される (11) ．

B11 B12 B13 000
B12 B11 B13 000
B13 B13 B33 0
0
0
SYM
BIJ
=
0
0
0 B44 00000
0 B44 0
0
0
0
0
0 B66
36

(3.35)SYM|BIJ| <0 となる条件は，以下の 4 つのいずれかとなる．
2
B33 (B11 + B12 ) − 2B13<0(3.36)

2
2
B11
− B12
<0

(3.37)

B44 < 0(3.38)

B66 < 0(3.39)

第一式 (3.36) は，体積弾性率が 0 になることを意味しており (48) ，本論文では以降これを Spinodal 条件と呼び，それによって生じる不安定を Spinodal 不安定と呼ぶ．第二式 (3.37) は，横方向変形のバランスを表しており，横方向変形が等方変形から非等方変形に変形経路分岐が起こることを意味する．これを同様に以降 Born 条件 (Born 不安定) と呼ぶ．さらに，第三式 (3.38)，第四式 (3.39) はそれぞれの変形モードへのせん
断不安定が生じることを表しており，同様にそれぞれ B44 条件 (B44 不安定)，B66 条件(B66 不安定) と呼ぶ．
また，本研究では静水圧変形についても検討するが，その場合独立な BIJ が減少し，系の安定性は次の 6 × 6 行列の正値性により評価される (11) ．

B11 B12 B12 000
B12 B11 B12 000
B12 B12 B11 0
0
0
SYM
BIJ
=
0
0
0 B44 00000
0 B44 0
0
0
0
0
0 B44

(3.40)

SYM| <0 となる条件は，以下の 3 つのいずれかとなる．
|BIJ

B11 + 2B12 < 0(3.41)

2
2
<0
− B12
B11

(3.42)

B44 < 0(3.43)

第一式 (3.41) が Spinodal 条件，第二式 (3.42) が Born 条件，第三式 (3.43) は B44 条
件である．

37

第4章
平衡状態における物性評価
第一原理バンド計算では，理論上，実験により得られた値を一切必要としないが，計算量を減らすために第 2 章で述べたような近似をいくつか導入している．よって，得られる解析結果の信頼性 (精度) を確認しておくことが重要である．その方法としては，
実在の物質の解析の場合，解析により得られる平衡格子定数，体積弾性率等の物性値を，実験により得られた値と比較することによって確認するのが一般的である (8) ．
本章では，fcc 構造の Ni 単結晶，および，二元系金属間化合物である L12 構造の Ni3 Al
単結晶の解析に必要な平面波のカットオフエネルギーの値および k 点数を，系のエネ
ルギーの収束性から決定する．次に，それらの値を用いて計算される平衡状態での Ni
および Ni3 Al の物性値を理論値と比較することにより，本研究で得られる解析結果の精度を確認する．

4.1

解析条件

4.1.1

モデルおよび予備解析

計算は，Kresse らにより開発された VASP(Vienna Ab-initio Simulation Packcage)(49)を用いた．Ni および Al の擬ポテンシャルには Ultrasoft 型 (31) を用い，k 点のサンプリ
ングは Monkhorst-Pack 法 (50) に従った．交換相関項の表現には一般化密度勾配近似 (51)(Generalized Gradient Approximation，GGA) を用い，残差最小化手法 (52) (Residual
Minimization Method – Direct Inversion in the Iterative Subspace，RMM–DIIS) を用
38

いて電子状態の収束計算を行った.
図 4.1(a)(b) に示す fcc 構造の Ni 単位格子 (原子数 4)，および，L12 構造の Ni3 Al 単位
格子 (原子数 4) をスーパーセルとし，格子定数 a を理論値 (Ni : 3.52 Å，Ni3 Al : 3.567

Å ) (53) に固定して次の二つの解析を行った．
(1) 平面波のカットオフエネルギー Ecut を決定するため，k 点数を 133 個に固定してEcut を変化させ，系の自由エネルギー Efree の変化を求める解析(2)k 点のサンプリング点数を決定するため，カットオフエネルギーを (1) より決定した Ecut に固定して k 点数を変化させ，系の自由エネルギー Efree の変化を求める解析

Ni，Ni3 Al それぞれの解析 (1),(2) における計算条件を表 4.1，4.2 に示す．

39

a

a

[001]
[010]
[100]

a

a
a

(a)

Ni unit cell

(b)

Ni

Ni3Al

unit cellAl

Fig.4.1 Unit lattices

Table 4.1 Conditions for ab-initio calculation of Ni
Number of atoms
Lattice constant (Å)
Cutoﬀ energy (eV)Number of k points
Number of bands
Boundary conditioncase(1)4case(2)4
3.52

150.0〜400.0
133

302
33 ∼ 153

28
Supercell

Table 4.2 Conditions for ab-initio calculation of Ni3 Al
Number of atoms
Lattice constant (Å)
Cutoﬀ energy (eV)Number of k points
Number of bands
Boundary conditioncase(1)4case(2)4
3.567

150.0〜400.0
133
25Supercell

40302
3 ∼ 1533

4.1.2平衡状態における格子定数，弾性定数の評価

前節の予備解析で得た Ecut と k 点数を用いて，無負荷平衡状態における Ni および

Ni3 Al の格子定数，弾性定数を以下の手順で求めた．まず，格子定数 a を変化させて，系に働く平均応力 σm を評価し，応力が零となる平衡格子定数 a0 を求めた．その後，図4.2 に示すように，平衡格子点で微小ひずみ摂動 ∆ηij (=0.01 あるいは 0.005) を与えたときの応力変化 ∆σij から弾性係数 C11 ，C12 ，C44 を数値的に求めた．さらに，得られた弾性係数を用いて体積弾性率 B0 を次式で評価した (12)(28) ．

B0 =

11=

C11 + 2C12
3

(4.1)

∆η

0.01

C

11 =

12=0.005

∆σ
∆η
∆σ
∆η

11

11

[001]
[010]
[100]

C

12 =

22

11

Fig.4.2 Schematic of evaluation of elastic coeﬃcients

41

C

44

∆σ
∆η

12

=

12

-12

Ni

Free energy

free

, eV

Ni3Al

-16

-20

-24
100

200

300

Cutoff energy

400

500

E , eVcut

Fig.4.3 Relationships between cutoﬀ energy and free energy4.2

解析結果及び考察

4.2.1

カットオフエネルギーの決定

カットオフエネルギー Ecut を変化させたときの系の自由エネルギー Efree の変化を
図 4.3 に示す．図より，Ni，Ni3 Al ともに 242eV 以上のカットオフエネルギーでは傾きがほぼ零となって収束している．カットオフエネルギーが大きいほど基底状態に近づいて自由エネルギーは低い値となり，電子状態をより正確に求めることできる．一方，解析に必要な平面波数 Npw は
3(2Ecut ) 2
Ω
6π 2

(4.2)

と見積もられ，平面波数の増加に対し記憶すべき量は 2 乗に，計算時間は 3 乗に比例する (Ω はセルの体積)．よって，解析精度と計算量とのバランスをとることが重要である．一般に，カットオフエネルギーを増加させたときに自由エネルギーの変化が

1mRy/atom(13.6meV/atom) 以下となれば十分収束したみなす (54) が，Ni，Ni3 Al ともに Ecut >=242eV で本条件を満たしている．ただし，次章以降では大きな変形を扱うため，Ni，Ni3 Al とも Ecut は 302eV とした．

42

-21.7

Free energy

E

free

, eV

-21.6

-21.8

-21.9

Ni
Ni3Al
-22
3

3

3

6

3

9

Number of

3

12

3

15

pointsFig.4.4 Relationships between number of k points and free energy

4.2.2

k 点数の決定

Ecut を 302eV に固定し，k 点数を変化させた時の系の自由エネルギー Efree の変化を3
図 4.4 に示す．Ni，Ni3 Al ともに k<73 以下では Efree が大きく変化しているが，k>=7で収束する傾向がある．金属系ではフェルミエネルギー近傍における占有数が非整数となるため (8) に，そのフェルミ面を正確に表現するには多数の k 点が必要となる．し
たがって，k 点数が多いほど電子状態をより正確に求めることできるが，やはり k 点
数の増加に対し記憶すべき量および計算時間は比例して増加するため，カットオフエ
ネルギーの決定と同様に計算精度と計算量のバランスをとることが重要である．一般に，金属で 1000 個/atom 以上，遷移金属では 5000 個/atom 以上の k 点をとれば，k点数を増加させたときの自由エネルギーの変化が 1mRy/atom(13.6meV/atom) 以下と3なるが (55) ，本解析では k>=5 で自由エネルギーの変化は 1mRy/atom(13.6meV/atom)

以下となっており，k=53 でも十分な k 点数であることがわかる．ただし，Ecut と同様，本研究では次節以降大きな変形を取り扱うため，Ni，Ni3 Al とも k 点数は十分大きな 133 個とした．

43

4.2.3

平衡状態における格子定数および弾性定数

Ni，Ni3 Al それぞれについて，格子定数 a を変化させたときの平均応力 σm の変化を
図 4.5(a)(b) に示す．平均応力 σm が 0 となる平衡格子定数 a0 は，Ni が 3.5238Å，Ni3 Al
が 3.5748Å となった．Ni，Ni3 Al の格子定数の理論値 (53) 3.52Å，3.567Å とのずれはそ
れぞれ 0.11%，0.22%と極めて小さい．交換相関項の表現に GGA を用いた場合，平衡格子定数の誤差は 1%以内となる (56) が，本解析は十分その範囲に収まっている．
本解析により求めた，平衡格子点における弾性定数 C11 ，C12 ，C44 および，体積弾性率 B0 の値を表 4.3 に示す．表にはそれぞれの理論値 (28)(53) をあわせて示している．格子定数に比べ，弾性係数の誤差はやや大きくなる．GGA を用いた場合，体積弾性率は10%程度の誤差を生じるとされている (56) ．本解析でも，文献 (28) と (53) の Ni，Ni3 Alの C11 ，C12 から求めた B0 の理論値とのずれはそれぞれ 7.98%，1.64%であり，誤差の範囲内となっている．ただし，文献 (53) で報告されている B0 の値 (181GPa) に対しては，ずれは 12.1%となった．これは，体積弾性率の評価方法に起因するものと思われ，本解析で用いた Ecut および k 点数は，Ni，Ni3 Al の解析に十分な精度を有している．

44

Mean stress

Mean stress

σm, GPa

5

, GPa

5

0

-5
3.49

3.5

3.51

3.52

3.53

Lattice constant
(a)

ð

3.54

a,

3.55

0

-5
3.56

3.54

3.55

3.56

3.57

3.58

Lattice constant
(b)

Ni

a,

ð

3.59

Ni3Al

Fig.4.5 Change in mean stress against lattice expansion/contraction

Table 4.3 Material constants of Ni and Ni3 Al

a0
C11
C12
C44
B0

Ni
calc. theo.(53)(Å)
3.5238
3.52
(GPa)
248
247
(GPa)
180
147
(GPa)
100
125
(GPa)
203
181

45

(28)

theo.
3.52
261
151
132
(188)

Ni3 Alcalc. theo.(53)
3.5748
3.567
235
230
162
159
117
131
186
(183)

3.6

4.3

結言

本章では，まず Ni および Ni3 Al 単結晶について，それぞれの解析に必要なカットオ
フエネルギー，k 点数を決定する予備解析を行った．それらを決定した後，それぞれの平衡状態における物性を評価することにより，本研究で得られる解析結果の精度について検討した．得られた結果を要約して以下に示す．

(1) k 点数を固定し，カットオフエネルギー Ecut を変化させる解析を行い，系の自由エネルギー Efree が収束するカットオフエネルギーを決定した．Ni，Ni3 Al ともに

Ecut >=242eV で自由エネルギーは十分収束することが示されたが，次章以降の有限変形領域の解析でも精度を保てるよう，Ni，Ni3 Al のカットオフエネルギーとして 302eV を決定した．

(2) カットオフエネルギーを固定し，k 点数を変化させる解析を行い，系の自由エネ
ルギー Efree が収束する k 点数を決定した．Ni，Ni3 Al ともに k=53 でも自由エネ
ルギーが十分収束することが示された．(1) と同様，大変形領域の解析のために，

Ni，Ni3 Al とも k 点数は十分大きな 133 個に決定した．

(3) Ecut =302eV，k=133 個の条件で，無負荷平衡状態における Ni および Ni3 Al の格子定数 a0 を求めた．平衡格子定数 a0 は，Ni が 3.5238Å，Ni3 Al が 3.5748Å となり，理論値からのずれはそれぞれ 0.11%，0.22%と極めて小さく，交換相関項に GGAを用いた場合の誤差範囲内に十分収まった．

(4) Ni，Ni3 Al それぞれの平衡格子点で弾性係数 C11 ，C12 ，C44 を数値的に求め，それらを用いて体積弾性率 B0 を評価した．Ni および Ni3 Al それぞれの体積弾性率は，203GPa，186GPa となり，理論値とのずれはそれぞれ 7.98%，1.64%であり，誤差範囲内に収まるため，(1)(2) から決定した Ecut ，k 点数は，Ni，Ni3 Al の解析に十分な精度を有していることが示された．

46

第5章

[001] 方向単軸引張変形下の理想格子不安定性解析
本章では，単位格子の変形で表される Ni および Ni3 Al の理想均一結晶について，[001]方向に引張変形を与える解析を行い，理想引張強度を求めるとともに，各ひずみ下の弾性剛性係数を数値的に評価し，格子不安定となる臨界ひずみを明らかにする．このとき，横方向の変形に等方的な Poisson 収縮を仮定した場合と，横方向の変形を拘束した場合の 2 通りについて考慮し，これらの変形拘束が Ni，Ni3 Al の [001] 方向理想引
張強度ならびに格子不安定性に及ぼす影響についても検討する．

5.1

解析条件

前章と同様，Ni および Ni3 Al の単位格子 (図 4.1) を用いて解析を行った．前章で求めた平衡格子定数 a0 を基準状態とし，[001] 方向に ∆ε33 =0.03(または 0.01) 刻みで引張ひずみを与え，第一原理分子動力学計算を行い電子状態の収束計算 · 内部構造緩和を行った．カットオフエネルギー Ecut ，k 点数の条件は前章と同じとした．Poisson 収
縮に関する変形拘束を考慮し，横方向のひずみは以下の 2 通りに制御した．

Case(a) 引張ひずみ増分 ∆ε33 を与えた後で横方向ひずみ ε11 =ε22 を ∆ε=0.0001 刻み
で変化させ，横方向応力 σ11 =σ22 が零となるように横ひずみ ε11 =ε22 を制御 (模式図 5.1(a))．

47

Case(b) 横方向変形を拘束し，ε11 =ε22 =0 に固定 (模式図 5.1(b))
Case(a),(b) それぞれの変形経路における各ひずみ下の安定構造を求めた後，図 5.2に模式的に示すように，微小ひずみ摂動 ∆ηij (=0.01 あるいは 0.005) を与え，応力変化∆σij から [001] 方向変形における安定性評価に必要な 7 つの独立な弾性剛性係数 (B11 ，
B12 ，B13 ，B33 ，B31 ，B44 ，B66 ) を，数値的に得た．得られた弾性剛性係数を用いて式 (3.35) を評価し，その正値性から各ひずみ下における系の安定性を評価した．

48

(a)

σ =σ
11

22

= 0

(b)

ε =ε
11

22

= 0

Fig.5.1 Simulation condition in lateral direction under [001] tension

η

11=

}

∆η

0.01

33=

}

0.01

∆σ
∆η

11

B

11 =

11

33 =

B

12 =

B

∆σ
∆η

33

B

33

21

∆σ
∆η

22

=

B

13

[001]

∆σ
∆η

11

11

B

=

B

=

31

33

∆σ
∆η

33

=

11

[010]
[100]

η

∆η

13=0.005B

44=

12=0.005

66

B

55

∆σ
∆η

12

12

∆σ
∆η

13

=

13

Fig.5.2 Schematic of evaluation of elastic stiﬀness coeﬃcients under [001]deformation

49

-19

30

σ33 , GPa

40

-20

Stressfree , eVFree energy

-18

-21

20

10

σ33)
σ33)

Ni

Ni (

Ni3Al

Ni3Al (

-22
0

0.1

0.2

0.3

Applied strain

0.4

0.5

ε33

0

0.1

0.2

0.3

Applied strain(a) Free energy - strain curves

0.4

0.5

ε33

(b) Stress - strain curves

Fig.5.3 Relationships between free energy/stress and applied strain under[001] tension (Case(a):stress free condition)

5.2

解析結果及び考察

5.2.1

自由エネルギー–引張ひずみ，引張応力–引張ひずみ関係

Case(a) : 横方向応力 0Case(a) の解析により得られた Ni および Ni3 Al の自由エネルギー–引張ひずみ関係，
および，引張応力–引張ひずみ関係をそれぞれ図 5.3(a)(b) に示す．Ni，Ni3 Al ともに，
系の自由エネルギーは 2 次関数的に増加した．応力–ひずみ関係は，ε33 =0.030 程度までは Ni，Ni3 Al ともにほとんど同一であり，わずかに下に凸となる非線形が現れている．平衡状態近傍で応力–ひずみ関係が下に凸となる特徴は，fcc 格子の [001] 方向変形
時に現れる特徴である (57) ．その後，Ni，Ni3 Al ともに ε33 =0.030∼0.120 程度まではほぼ線形に応力は増加するが，ε33 =0.060 以降，同じひずみに対し Ni3 Al はわずかに低い応力を示し，応力–ひずみ曲線はわずかに異なったものとなった．この応力–ひずみ関係のピーク点を理想引張強度とする場合，Ni は ε33 =0.390，σ33 =35.3GPa，Ni3 Al は

ε33 =0.360，σ33 =34.0GPa と評価される．

50

40

σ33 , GPa

-19

-20

StressFree energy

free , eV

-18

-21

30

20

σ33)
σ33)
Ni (σ11)
Ni Al (σ11)
Ni (

10

Ni3Al (

Ni
Ni3Al

3

-22
0

0.1

0.2

0.3

Applied strain

0.4

0.5

ε33

0

0.1

0.2

0.3

Applied strain(a) Free energy - strain curves

0.4

0.5

ε33

(b) Stress - strain curves

Fig.5.4 Relationships between free energy/stress and applied strain under[001] tension (Case(b):deformation constraint condition)

Case(b) : 横方向変形拘束
Case(b) の解析により得られた Ni および Ni3 Al の自由エネルギー–引張ひずみ関係，
および，引張応力–引張ひずみ関係をそれぞれ図 5.4(a)(b) に示す．なお，図 5.4(b) には
横方向応力 σ11 (=σ22 ) も合わせて示した．Case(a) と同様，Ni，Ni3 Al の系の自由エネ
ルギーはいずれも 2 次関数的に増加した．また，Case(b) では [001] 方向の引張だけでなく，[100],[010] 方向にも引張を受ける 3 軸応力状態にあるため，同じ引張ひずみでの値を Case(a) と比較すると自由エネルギーは高い値となっている．一方，応力–ひずみ関係は，Ni，Ni3 Al とも引張初期から非線形性が現れており，ε33 =0.030 以降で Ni3 Al
が低応力を示すため，わずかにずれる．この応力–ひずみ関係のピークより評価される Ni，Ni3 Al の理想引張強度は，Ni が ε33 =0.330，σ33 =35.1GPa，Ni3 Al が ε33 =0.330，

σ33 =33.6GPa となった．Case(a) と同様，Ni が Ni3 Al よりも高い降伏応力となったが，降伏ひずみで見ると同一となる．

51

5.2.2格子不安定性Case(a) : 横方向応力 0系の安定性を評価する 4 つの式 (3.36)∼(3.39) の値の変化を図 5.5(a)∼(d) に示す．Ni
では，ε33 =0.110 で式 (3.37)(Born 条件)，ε33 =0.240 で式 (3.38)(B44 条件)，ε33 =0.390
で式 (3.36)(Spinodal 条件) がそれぞれ負となった．式 (3.39) の B66 条件は解析した引張ひずみの範囲では負になることはない．Spinodal 条件のひずみ ε33 =0.390 は応力–ひずみ関係のピーク (理想引張強度) に一致する．表 5.1 に式 (3.36)∼(3.39) の値が負になるひずみ前後の弾性剛性係数を示したが，ε33 =0.390 の点では引張方向に対する変形抵抗 B33 が負になっており，これが Spinodal 不安定をもたらしている．すなわち，先の応力–ひずみ関係のピーク点より得られる理想引張強度は，(001) 原子面内のへき開に対するものである．Born 条件が負となる ε33 =0.110 では，表 5.1 に示したように

B11 と B12 の大小関係が逆転している．この不安定により，横方向変形のバランスが崩れ，等方変形からよりエネルギーの低い安定な変形経路である非等方変形へと変形
経路分岐が生じることが示されている (19) ．B44 が負となる ε33 =0.240 では，(100) 面[001] 方向やそれらに等価な方向へのせん断不安定変形を生じるものと考えられる．ただし，ε33 =0.240 の B44 条件の不安定も，ε33 =0.390 の Spinodal 条件の不安定も，それより小さいひずみ ε33 =0.110 の Born 不安定で生じる変形分岐を拘束した場合の第二，第三の不安定条件である．

Ni3 Al は，ε33 =0.090，0.300，0.360 でそれぞれ Born，B44 ，Spinodal 条件が負となる．Ni と同様に B66 条件は解析した引張ひずみの範囲では負になることはない．Born
不安定と Spinodal 不安定のひずみは Ni よりも小さくなったが，B44 に対するそれは Ni
よりも大きくなった．Spinodal 条件が理想引張強度と一致し，B33 が負となっている
ことも表 5.1 からわかる．最初に現れる Born 不安定条件の不安定ひずみを結晶の弾性限界と考えると，降伏に対する理想強度は Ni が ε33 =0.100，σ33 =13.5GPa，Ni3 Al がε33 =0.080，σ33 =10.8GPa となり，応力–ひずみ関係のピークより評価した理想引張強度よりずっと低いひずみ · 応力となる．ただし，Ni と Ni3 Al の強度の大小関係は変わ
らない．

52

~

2

10

2
11 -

33

0
Ni
Ni3Al

0

0.1

0.2

0.3

Applied strain

(a)

0.4

0

-200
Ni
Ni3Al

-400

0.5

0

ε33

0.1

2
13

B (B +B )-2B
33

11

12

0.2

0.3

Applied strain

(b)

150

2
11 -

B

0.4

0.5

ε33

2
12

B

B , GPa

200

100

150

Magnitude of66

44

B , GPa

2

200

B
Magnitude of

200

11

12

400

-200

Magnitude of

10

B

2
13 ,

2
12 ,

600

B (B +B )-2B

Magnitude of

400

2
GPa

2
GPa

800

~

50

0
Ni

100

50

0
Ni

Ni3Al

Ni3Al

-50

-50
0

0.1

0.2

0.3

Applied strain

(c)

0.4

0.5

ε33

0

0.1

0.2

0.3

Applied strain

B

(d)

44

0.4

0.5

ε33

B

66

Fig.5.5 Changes in the minor determinants of detBIJ under [001] tension(Case(a):stress free condition)

Table 5.1 Elastic stiﬀness coeﬃcients under [001] tension(Case(a):stress free condition)
ε33
B11 (GPa)
B12 (GPa)
B33 (GPa)
B13 (GPa)ε33
B11 (GPa)
B12 (GPa)ε33
B44 (GPa)

Ni
0.360 0.390
206.8 212.1
255.3 248.7
4.7-3.3
20.4 14.8
0.100 0.110
239.7 239.2
237.3 239.6
0.210 0.240
5.2
-9.2

53

Ni3 Al0.330 0.360
181.1 179.5
239.8 234.8
13.0
-0.6
23.1 17.5
0.080 0.090
212.0 219.1
207.6 225.5
0.270 0.300
7.6
-0.2

Case(b) : 横方向変形拘束
Case(b) における式 (3.36)∼(3.39) の値の変化を図 5.6(a)〜(d) に示す．Ni，Ni3 Al とも，Case(a) と同様に引張ひずみの増加に対して Born 条件，B44 条件，Spinodal 条件の順で負になる．B66 条件は，やはり解析した引張ひずみの範囲では負になることはない．
Spinodal
B44
それぞれの条件が負となるひずみは，Ni が εBorn=0.360，33 =0.190，ε33 =0.270，ε33
Spinodal
B44
Ni3 Al が εBorn=0.360 となった．表 5.2 に不安定前後の弾
33 =0.130，ε33 =0.360，ε33性剛性係数の値を示した．B33 はひずみ εSpinodal=0.330 でほぼ零であり，Spinodal 不安
33
定は Case(a) と同様に応力–ひずみ関係のピークに一致すると考えられる．また，Ni，Ni3 Al とも εBorn
は Case(a) に比べてそれぞれ 0.080，0.040 大ひずみ側にシフトしてい33
る．一方，εSpinodal
は Ni3 Al が Case(a) と同じひずみ，Ni が 0.030 小ひずみ側にシフト33している．このため，横方向変形を拘束した場合，第一の不安定である Born 不安定と，第三の不安定である Spinodal 不安定の臨界ひずみの差が小さくなり，へき開を生じる可能性がある (47) ．最初に系が不安定となる εBorn
を弾性限界とした場合の理想強33
度は Ni が ε33 =0.180，σ33 =28.9GPa，Ni3 Al が ε33 =0.120，σ33 =21.3GPa となり，横方向拘束の場合も応力–ひずみ関係のピークより低いひずみ · 応力となる．

54

~

2

10

200

2
11 -

B
Magnitude of

200

11

12

400

33

0
Ni
Ni3Al

-200
0

0.1

0.2

0.3

Applied strain

(a)

0.4

0

-200
Ni
Ni3Al

-400

0.5

0

ε33

0.1

2
13

B (B +B )-2B
33

11

12

0.2

0.3

Applied strain
2
11 -

0.4

0.5

ε33

2
12

B

B

0.2

0.3

(b)
200

GPa

150

100

150

66,

44,

GPa

2

10

B

2
13 ,

2
12 ,

600

B (B +B )-2B

Magnitude of

400

2
GPa

2
GPa

800

~

Magnitude ofMagnitude of

100
50

0
Ni

50

0
Ni

Ni3Al

Ni3Al

-50

-50
0

0.1

0.2

0.3

0.4

0.5

0

0.1

Applied strain

Applied strain

(c)

(d)

44

0.4

0.5

66

Fig.5.6 Changes in the minor determinants of detBIJ under [001] tension(Case(b):deformation constraint condition)

Table 5.2 Elastic stiﬀness coeﬃcients under [001] tension(Case(b):deformation constraint condition)
ε33
B11 (GPa)
B12 (GPa)
B33 (GPa)
B13 (GPa)ε33
B11 (GPa)
B12 (GPa)ε33
B44 (GPa)Ni
0.330 0.360
90.3 91.3
104.6 104.7
0.7
-10.7
5.0
-1.5
0.180 0.190
131.9 127.7
130.3 128.8
0.240 0.270
3.4
-2.4

55

Ni3 Al0.330 0.360
86.1 78.4
102.6 97.7
0.8
-12.0
4.2
3.1
0.120 0.130
148.8 137.6
147.5 139.8
0.330 0.360
1.4
-2.0

5.3

結言

[001] 方向に引張変形を受ける Ni および Ni3 Al 単結晶について，等方的な Poisson
収縮を仮定した横方向応力 0 の場合 (Case(a))，および，横方向変形を拘束した場合(Case(b))，の 2 通りの解析を行った．それぞれの変形経路における応力–ひずみ曲線のピークから理想引張強度を評価するとともに，各ひずみ下の弾性剛性係数を求めて格子不安定性を評価した．得られた結果を要約して以下に示す．

(1) 応力–ひずみ曲線のピークより得られる理想引張強度は，Case(a) では Ni が ε33 =0.390，
σ33 =35.3GPa，Ni3 Al が ε33 =0.360，σ33 =34.0GPa となった．Case(b) では Ni がε33 =0.330，σ33 =35.1GPa，Ni3 Al が ε33 =0.330，σ33 =33.6GPa となり，Ni>Ni3 Alの大小関係は変わらないが降伏ひずみは同一となった．

(2) 各ひずみ下における系の安定性を評価した結果，Case(a),(b) とも Born 条件，B44
条件，Spinodal 条件の順に系が不安定分岐点にあることがわかった．

(3) Spinodal 不安定となるひずみは Case(a),(b) とも (1) の応力–ひずみ曲線のピークに一致し，かつ (001) 面のへき開方向の弾性剛性係数 B33 が負となることによりもたらされていた．

(4) 最初に現れる Born 不安定では，横方向変形に関する弾性剛性係数 B11 ，B12 の大小関係が逆転し，非等方変形への変形分岐を生じることが示唆される．

を結晶の弾性限界とした場合(5) 最初に現れる Born 不安定条件の臨界ひずみ εBorn
33の理想強度は，Case(a) では Ni が ε33 =0.100，σ33 =13.5GPa，Ni3 Al が ε33 =0.080，

σ33 =10.8GPa，Case(b) では Ni が ε33 =0.180，σ33 =28.9GPa，Ni3 Al が ε33 =0.120，
σ33 =21.3GPa となる．Ni，Ni3 Al の強度の大小関係は変わらないが，(1) の強度よ
りはるかに小さい応力 · ひずみとなる．

56

(6) 横方向の変形拘束は最初の Born 不安定を大ひずみ側にシフトさせ，Spinodal 不
安定を低ひずみ側にシフトさせる．これより，弾性限界である Born 不安定とへき開条件である Spinodal 不安定の差が縮まり，へき開を生じやすくなることが示唆
される．

57

第6章

[001] 方向単軸圧縮変形下の理想格子不安定性解析
6.1

解析条件

前章と同様，Ni および Ni3 Al の単位格子 (図 4.1) を用いて解析を行った．第 4 章で求めた平衡格子定数 a0 を基準状態とし，[001] 方向に ∆ε33 =–0.03(または–0.01) 刻みで
圧縮ひずみを与え，第一原理分子動力学計算を行い電子状態の収束計算 · 内部構造緩和を行った．カットオフエネルギー Ecut ，k 点数の条件は第 4 章と同じとした．ここで，横方向の膨張に関する変形拘束を考慮し，前章と同様に横方向のひずみは以下の2 通りに制御した．
Case(a) 圧縮ひずみ増分 ∆ε33 を与えた後，横方向ひずみ ε11 =ε22 を ∆ε=0.0001 刻み
で変化させ，横方向応力 σ11 =σ22 が零となるように制御 (模式図 6.1(a))．

Case(b) 横方向変形を拘束し，ε11 =ε22 =0 に固定 (模式図 6.1(b))
Case(a),(b) それぞれの変形下における安定構造を求めた後，前章と同様に 7 つの独立な弾性剛性係数 (B11 ，B12 ，B13 ，B33 ，B31 ，B44 ，B66 ) を数値的に求め，その正値性から系の安定性を評価した．

58

(a)

σ =σ
11

22

= 0

(b)

ε =ε
11

22

= 0

Fig.6.1 Simulation condition in lateral direction under [001] compression

6.2

解析結果及び考察

6.2.1

自由エネルギー–圧縮ひずみ，圧縮応力–圧縮ひずみ関係

Case(a) : 横方向応力 0Case(a) の解析により得られた Ni および Ni3 Al の自由エネルギー–圧縮ひずみ関係，
および，圧縮応力–圧縮ひずみ関係をそれぞれ図 6.2(a)(b) に示す．前章の引張の場合と異なり，Ni，Ni3 Al ともに圧縮ひずみ ε33 =–0.200 近傍でエネルギー–ひずみ曲線に変曲点が見られ，応力–ひずみ関係に極小および極大点が現れている．ε33 =–0.060 程度までの圧縮では Ni，Ni3 Al の応力–ひずみ関係はほぼ同一である．その後 Ni は ε33 =–0.100，

Ni3 Al は ε33 =–0.120 で極小点を示した．極小点における応力は Ni が σ33 =–5.44GPa，
Ni3 Al が σ33 =–5.94GPa となった．理想圧縮強度を示すこれらの極小点の後，圧縮ひず
みの増加とともに応力は再び増加し極大を示す．これについては，圧縮における結晶
構造関係 (Bain の関係) の観点から 6.2.3 節で述べる．

59

-21

5
Ni (

σ33 , GPa

free , eV-21.5

StressFree energy

σ33)
σ33)

Ni3Al (

0

-5

Ni
Ni3Al

-22

-10

-0.5

-0.4

-0.3

-0.2

Applied strain

-0.1

0

-0.5

ε33

-0.4

-0.3

-0.2

Applied strain(a) Free energy - strain curves

-0.1

0

ε33

(b) Stress - strain curves

Fig.6.2 Relationships between free energy/stress and applied strain under[001] compression (Case(a):stress free condition)

-5

100

σ33)
σ33)
Ni (σ11)
Ni Al (σ11)

0σ33 , GPa

-10

-15

StressFree energy

free , eVNi (Ni3Al (

-20

3

-100

-100

-200

-300

-400

Ni

-250

Ni3Al

-0.26

-0.34

-25

-500

-0.5

-0.4

-0.3

-0.2

Applied strain

-0.1

0

-0.5

ε33

-0.4

-0.3

-0.2

Applied strain(a) Free energy - strain curves

-0.1

ε33

(b) Stress - strain curves

Fig.6.3 Relationships between free energy/stress and applied strain under[001] compression (Case(b):deformation constraint condition)60

0Case(b) : 横方向変形拘束
Case(b) の解析により得られた Ni および Ni3 Al の自由エネルギー–圧縮ひずみ関係，
および，圧縮応力–圧縮ひずみ関係をそれぞれ図 6.3(a)(b) に示す．図 6.3(b) には横方
向応力 σ11 (=σ22 ) も合わせて示している．系の自由エネルギーは，Ni，Ni3 Al ともに前章の引張と同様に 2 次関数的に増加した．しかしながら，応力は上に凸の 2 次関数的な曲線となり，極大 · 極小を示さなかった．したがって，Case(b) では Ni，Ni3 Al ともに
応力–ひずみ関係からは [001] 方向の理想圧縮強度は定義できない．なお，図 6.3(b) 中に拡大して示したように，Ni は ε33 =–0.300 のとき，σ33 と σ11 がほぼ等しくなる．一方，Ni3 Al も ε33 =–0.300 で σ33 と σ11 が最も近づくが，Ni と異なり一致することはな
い．これらについては，結晶構造の Bain の関係から 6.2.3 節で後述する．

61

6.2.2格子不安定性Case(a) : 横方向応力 0
図 6.4(a)∼(d) は，系の安定性を評価する 4 つの式 (3.36)∼(3.39) の値の変化を示したものである．Ni，Ni3 Al ともに最初に Spinodal 条件，次に B66 条件が負となった．Born
条件，B44 条件は，ε33 >=–0.420 の範囲において負にならず，これらの不安定モードが発生することはない．Spinodal，B66 条件が負となるひずみは，Ni が εSpinodal=–0.110，33
Spinodal
66
66=–0.120，εB
εB
33 =–0.180 であり，Ni3 Al が ε33
33 =–0.210 となった．Spinodal 条件の不安定は，前章の引張と同様応力–ひずみ関係の極小点 (圧縮を正とした場合の応力–
ひずみ関係のピーク) にほぼ一致にした．しかしながら，表 6.1 に示した弾性剛性係数の前後で B33 は負になっの変化からわかるように，[001] 方向引張時と異なり，εSpinodal
33ていない．したがって，圧縮変形における Spinodal 不安定は，fcc の構造を保ったまま(001) 原子面間の抵抗を喪失して崩壊するのではなく，別の結晶構造変化への不安定と考えられる．なお，B66 条件の不安定は，(100) 面 [010] 方向やそれに等価な方向のせん断変形への抵抗が失われることを表す．fcc 構造のアルカリ金属では B66 条件が負になることにより，せん断変形への経路分岐が生じるとの報告がある (17) ．
等方的な Poisson 膨張下の [001] 圧縮では，最初に現れる Spinodal 不安定が応力–ひずみ曲線の極小点 (圧縮を正としたときのピーク) に一致したことから，格子不安定の観点からも理想圧縮強度は同じとなる．すなわち，横方向変形を許容した [001] 圧縮の理想強度は，Ni と Ni3 Al の強度の大小関係が前章の [001] 引張におけるそれと逆転
する．

62

~

~

2

10

GPa

2

2
12 ,

GPa

2
11 -

B
Magnitude of

33

Ni
Ni3Al

-500
-0.5

-0.4

-0.3

-0.2

Applied strain

(a)

-0.1

0

0
Ni
Ni3Al

-100

ε33

-0.4

2
13

B (B +B )-2B
33

11

12

-0.3

-0.2

Applied strain
2
11 -

-0.1

0

ε33

2
12

B

B

-0.3

-0.2

(b)
200

B , GPa

100

66

400

Magnitude of44

B , GPa

100

-0.5

600

Magnitude of

200

B

2
13 ,

0

11

12

B (B +B )-2B

Magnitude of

300

2

500

3

10

200

0

Ni

0

-100

-200
Ni

Ni3Al

Ni3Al

-200
-0.5

-300
-0.4

-0.3

-0.2

Applied strain

(c)

-0.1

0

-0.5

ε33

-0.4

Applied strain

B

(d)

44

-0.1

ε33

B

66

Fig.6.4 Changes in the minor determinants of detBIJ under [001] compression(Case(a):stress free condition)

Table 6.1 Elastic stiﬀness coeﬃcients under [001] compression(Case(a):stress free condition)
ε33
B11 (GPa)
B12 (GPa)
B33 (GPa)
B13 (GPa)ε33
B66 (GPa)

Ni
Ni3 Al-0.100 -0.110 -0.110 -0.120
290.0 292.1 269.9 264.5
105.2
96.2
85.0
68.4
254.5 251.7 216.3 219.6
223.5 228.8 195.2 196.4-0.150 -0.180 -0.180 -0.210
9.6
-8.0
10.2
-15.6

63

0

Case(b) : 横方向変形拘束
Case(b) における式 (3.36)∼(3.39) の値の変化を図 6.5(a)∼(d) に示す．Case(a) と同様，Spinodal 条件，B66 条件の順でそれぞれの値が負となった．そのひずみは，Ni がSpinodal
66
66
εSpinodal=–0.150，εB=–0.140，εB
33
33 =–0.300，Ni3 Al が ε33
33 =–0.300 である．Born

条件および B44 条件は ε33 >=–0.420 の範囲において負にならなかった．表 6.2 に Spinodal
不安定，B66 不安定前後における弾性剛性係数の値を示す．Case(a) と同様，εSpinodal
33の前後で B33 は負になっておらず，結晶構造変化開始への分岐点と考えられる．変形
経路に沿った応力–ひずみ曲線からは横拘束下の [001] 圧縮の理想強度は存在しないのに対し，格子不安定の観点からは同変形条件における臨界ひずみが定義される．その値は Ni が ε33 =–0.140，σ33 =–49.5GPa，Ni3 Al が ε33 =–0.130，σ33 =–39.9GPa である．
したがって，横方向拘束下での [001] 圧縮では理想強度が Case(a) と逆転し，Ni3 Al より Ni の方が高くなる．

64

~

~

2

10

GPa

2

2
12 ,

GPa

2
11 -

B
Magnitude of

Ni
Ni3Al

-0.5

-0.4

-0.3

-0.2

Applied strain

(a)

-0.1

0

100

0
Ni
Ni3Al

-100
-0.5

ε33

-0.4

2
13

B (B +B )-2B
33

11

12

-0.3

-0.2

Applied strain
2
11 -

-0.1

0

ε33

2
12

B

B

-0.3

-0.2

(b)

600

GPa

200

100

66,

400

44,

GPa

200

B

2
13 ,
33

11

12

B (B +B )-2B

Magnitude of

0

-500

0

Magnitude ofMagnitude of

300

2

500

3

10

200

0

Ni

-100

-200
Ni

Ni3Al

Ni3Al

-200
-0.5

-300
-0.4

-0.3

-0.2

-0.1

0

-0.5

-0.4

Applied strain

Applied strain

(c)

(d)

44

-0.1

66

Fig.6.5 Changes in the minor determinants of detBIJ under [001] compression(Case(b):deformation constraint condition)

Table 6.2 Elastic stiﬀness coeﬃcients under [001] compression(Case(b):deformation constraint condition)
ε33
B11 (GPa)
B12 (GPa)
B33 (GPa)
B13 (GPa)ε33
B66 (GPa)

Ni
Ni3 Al-0.140 -0.150 -0.130 -0.140
547.1 573.0 459.1 482.4
248.3 248.4 202.1 201.7
497.0 511.3 404.1 420.4
440.0 472.3 361.7 384.6-0.270 -0.300 -0.270 -0.300
19.8
-51.4
23.2
-58.6

65

0

6.2.3

圧縮変形における結晶構造関係

Case(a) : 横方向応力 0
図 6.6(a) に示すように，Ni は fcc(face centered cubic) 構造であるが，結晶軸方向の異なる bct(body centered tertagonal) 格子 (TypeI) を含むと考えることができる．一方，Ni3 Al では，図 6.6(b) に示すように Ni のみで構成される bct 格子 (TypeI) と，Niと Al が 1:1 で構成される bct 格子 (TypeII) を含むと考えることができる．図 6.6 中に示
すように，fc 格子の 3 辺の長さを λf1 ，λf2 ，λf3 と表すと，λf1 = λf2 = λf3 = 1 のとき fc 格
子は cubic 構造をとる．[001] 方向に λf3 < 1 となる圧縮ひずみを与え，横方向に等方的√な膨張 λf1 = λf2 > 1 を考える．格子長さが λf3 /λf1 = 1/ 2 の関係になるとき，bct 格子はより対称性の高い bcc 格子となる．これを Bain の関係 (Bain 相変態) という (17)(58) ．
本研究における Ni の場合，ε33 =–0.200 のときの横方向ひずみは ε11 =0.122 であった．
√
これより λf3 と λf1 の比を計算すると λf3 /λf1 = (1 + ε33 )/(1 + ε11 ) = 0.713 · · · ≈ 1/ 2 となりほぼ Bain の関係を満たす．図 6.2(a)(b) からわかるように ε33 =–0.200 近傍において，自由エネルギーは極大値をとり σ33 も零，σ11 =σ22 は零に制御しているので，系に働く応力の全成分が零の bcc 構造になっていることがわかる．ただし，エネルギーの極大点にあることから不安定な構造である．また，ε33 =–0.260 近傍では，自由エネル
ギーが極小で応力が零の “dead ”(14) で準安定な bct 構造となる．したがって，圧縮変形における結晶構造の関係，および，前節で得られた不安定条件の関係を整理すると，fcc → Spinodal 不安定→ B66 不安定→ bcc → bct となる．一方，Ni3 Al では ε33 =–0.210のときに σ33 が最も零に近づき，このときの横方向ひずみは ε11 =0.126 である．λf3 /λf1√
は 0.701≈ 1/ 2 となり，fc 格子内の TypeI 格子は bcc 構造に，TypeII 格子は TypeI の中心の Ni が Al で置換された B2 構造に近くなる．ひずみ ε33 =–0.210 は，Ni で TypeIが bcc 構造になるひずみ ε33 =–0.200 とほぼ一致する．しかしながら，Ni3 Al では TypeIと TypeII が混在したするため，結晶全体としては応力は零にならない．なお，ひずみ ε33 =–0.210 は B66 不安定の臨界ひずみと一致している．fcc 構造の B66 不安定によるせん断変形は，内部 bc 格子の非等方変形とも考えられる (17) ．これより，Ni3 Al における B66 不安定は bcc 構造に変化することと密接に関係しているものと考えられる．

66

[001]
compressive direction[010]

compressive direction

[100]

f

f

λ3

λ3

Type I

Type II

f

λ2

2

f

λ1

1

(b) Ni3Al (L12 structure)(a) Ni (fcc structure)

Fig.6.6 Schematic illustration of Bain’s relationship in fcc and L12 ordered structure

Ni3 Al の圧縮変形における結晶構造の関係，および，前節で得られた不安定条件の関係は，fcc(L12 ) → Spinodal 不安定→ [bcc(B2)，B66 不安定] となる．

Case(b) : 横方向変形拘束
Case(a) と同様，図 6.6 の格子変形について考える．横方向変形拘束の [001] 圧縮は，
λf3 < 1，λf1 = λf2 = 1 で表される．図 6.3(b) で示したように，Ni，Ni3 Al ともに ε33 =–√0.300 のときに全応力成分がほぼ等値となる．このとき，λf3 /λf1 は 0.700≈ 1/ 2 で Bainの関係にあり，TypeI，TypeII はそれぞれ bcc，B2 構造となっている．Ni では結晶全体が等方応力下の bcc となる．一方，Ni3 Al は前節と同様，bcc と B2 構造が混在するため等方応力とはならず，わずかに差を生じている．また，Case(a) と同様，ε33 =–0.300
は Ni，Ni3 Al とも B66 不安定を生じるひずみと一致しており，fcc → bcc もしくは B2構造への変化と密接に関係していることが示唆される．

67

6.3

結言

[001] 方向に圧縮変形を受ける Ni および Ni3 Al 単結晶について，横方向の応力が 0 となるよう等方膨張を仮定した場合 (Case(a))，および，横方向ひずみが 0 となるように変形を拘束した場合 (Case(b))，の 2 通りの解析を行った．それぞれの変形経路における応力–ひずみ曲線より示される理想圧縮強度について検討するとともに，各ひずみ下の弾性剛性係数より格子不安定性について検討した．また，圧縮時に生じる fcc → bcc
構造変化 (Bain の関係) の観点からも，応力–ひずみ曲線および格子不安定点について考察した．得られた結果を要約して以下に示す．

(1) Case(a) の解析では，Ni，Ni3 Al ともに応力–ひずみ関係が極小 · 極大を有する 3 次
関数的な曲線となった．特に Ni では圧縮軸の応力が零となる点が極大点前後に現れ，無応力下の別の結晶構造となっている可能性が示唆される．最初に現れる応力–
ひずみ関係の極小点を理想圧縮強度とすると，Ni が ε33 =–0.100，σ33 =–5.44GPa，

Ni3 Al が ε33 =–0.120，σ33 =–5.94GPa となり，降伏ひずみ，降伏応力とも Ni3 Al がNi より高くなった．これは，第 5 章の引張変形下での理想強度の大小関係と逆転している．

(2) Case(b) の解析では，応力–ひずみ関係は上に凸の 2 次関数的な曲線となり，理想圧縮強度は定義できなかった．

(3) 各ひずみ下における系の安定性を評価した結果，Case(a) では Ni は ε33 =–0.110 とε33 =–0.180，Ni3 Al は ε33 =–0.120 と ε33 =–0.210 でそれぞれ Spinodal 不安定，B66
不安定となった．Spinodal 不安定となるひずみは (1) の応力–ひずみ曲線の極小点と一致しており，横方向応力 0 の [001] 圧縮では最初に現れる格子不安定と理想圧
縮強度が一致することが示された．

(4) Case(b) における系の安定性評価では，Spinodal 不安定，B66 不安定がそれぞれ現68

れることが示された．したがって，(2) で応力–ひずみ関係では理想強度は定義されなかったのに対し，格子不安定の観点からは Ni が ε33 =–0.140，σ33 =–49.5GPa，

Ni3 Al が ε33 =–0.130，σ33 =–39.9GPa が弾性限界に対する理想強度となる．Case(a)の場合と異なり，Ni，Ni3 Al の強度の大小関係は再び逆転し，引張強度と同じ傾向となった．

(5) Case(a) における 3 次曲線的な応力–ひずみ関係は，fcc または L12 構造中の bct 格子が bcc または B2 格子となるためであることが示された．特に Ni は ε33 =–0.200のときに不安定な bcc，ε33 =–0.260 のとき準安定な bct となり，すべての応力成分が零となる．Ni3 Al は Al の存在により bcc，B2 格子が混在するため無応力状態の別の結晶構造をとることはない．

(6) Case(b) における結晶構造関係を調べた結果，Ni，Ni3 Al とも B66 条件が負となる
ε33 =–0.300 で fcc，L12 構造中の bct が bcc，B2 構造となる点であることが示された．これより，B66 の不安定条件が fcc → bcc または B2 構造への相変態に密接に関係することが示唆される．

69

第7章
静水圧変形下の理想格子不安定性解析
前章と同様，単位格子の変形で表される Ni および Ni3 Al の理想均一結晶について，[100],[010],[001] の 3 軸に等方的な静水圧引張および圧縮を与える解析を行い，応力–ひずみ関係より得られる理想強度，ならびに各静水圧ひずみ下の格子不安定性について議論する．静水圧引張または圧縮は，き裂先端や押し込み試験における圧子直下の変形に対応し (18) ，転位の発生やへき開の進展などに有益な知見を与えるものと期待さ
れる．

7.1

解析条件

第 4 章と同様，Ni および Ni3 Al の単位格子 (図 4.1) を用いて解析を行った．第 4 章で求めた平衡格子定数 a0 を基準状態とし，静水圧条件下での変形として [100],[010],[001]方向に引張または圧縮のひずみを ∆ε11 =∆ε22 =∆ε33 =0.03 または 0.01 刻みで増加させた (模式図 7.1)．その後，第一原理分子動力学計算により電子状態の収束計算 · 内部構
造緩和を行い，各ひずみ点における安定構造を求めた．カットオフエネルギー Ecut ，k
点数の条件は第 4 章と同じとした．
第 3 章で述べたように，静水圧変形における安定性評価に必要な独立な弾性剛性係70Hydrostatic tension

Hydrostatic compression

Fig.7.1 Simulation condition under hydrostatic tension/compression

数は 3 つ (B11 ，B12 ，B44 ) である．前章までと同様，微小ひずみ摂動 ∆ηij (=0.01 ある
いは 0.005) を与えたときの応力変化 ∆σij から数値的に求めた．得られた弾性剛性係数を用いて式 (3.40) を評価し，その正値性から各ひずみ点における系の安定性を評価
した．

7.2

解析結果及び考察

7.2.1

自由エネルギー–ひずみ関係，応力–ひずみ関係

静水圧ひずみを変化させたときの Ni と Ni3 Al の自由エネルギー–ひずみ関係，および，応力–ひずみ関係をそれぞれ図 7.2(a)(b) に示す．自由エネルギー–ひずみ曲線は，

Ni と Ni3 Al とも平衡状態を極小としてほぼ同じ形となった．引張側は第 5 章の Case(b)の [001] 引張，圧縮側も第 6 章の Case(b) の [001] 圧縮のものに近い．応力–ひずみ関係も同様の傾向を示し，引張側では Ni，Ni3 Al のいずれも ε33 =0.150 で極大をとり，圧縮側では極値を示さない．引張におけるピークを静水圧引張における理想強度とすると，その時の応力は Ni が σ33 =30.0GPa，Ni3 Al が σ33 =28.3GPa となり，Ni と Ni3 Al の強度関係は第 5 章の Case(b) の [001] 引張と同じとなる．圧縮側では応力–ひずみ関係からは静水圧圧縮における理想強度は定義できない．

71

50
Ni

Ni

Ni3Al

Ni3Al

σ33 , GPa

-6

-12StressFree energy

free , eV

0

-18

-24

0

-50

-100

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

Applied strain

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

Applied strain(a) Free energy - strain curves(b) Stress - strain curves

Fig.7.2 Relationships between free energy/stress and applied strain underhydrostatic tension/compression7.2.2

格子不安定性図 7.3(a)∼(c) は，静水圧変形の場合において系の安定性を評価する式 (3.41)∼(3.43)の値の変化である．[001] 引張の場合と異なり，Ni，Ni3 Al とも Spinodal 条件が最初の不安定分岐点となった．Ni は Spinodal 条件，Born 条件，B44 条件の順で，それぞれ

ε33 =0.150，0.160，0.240 で負となった．Ni3 Al は ε33 =0.160 において Spinodal 条件とBorn 条件が同時に負となり，B44 条件が ε33 =0.240 で負となった．Spinodal 条件の臨界ひずみは応力–ひずみのピークとほぼ一致している．各条件に達する臨界ひずみ前後における弾性剛性係数を表 7.1 に示した．Ni および Ni3 Al のいずれも，Spinodal 条件
が負となるひずみでも引張 (膨張) に対する抵抗を表す B11 は負になっていない．このことから，[001] 引張の場合と異なり，静水圧引張における応力–ひずみ関係のピークは結晶が急激に膨張してバラバラになる点ではなく，[001] 圧縮変形時と同様に別の結晶構造への相変態点を表すものと考えられる．Born 不安定となるひずみでは B11 と B12の絶対値の大小関係が逆転しており，[001] 引張と同様に非等方変形に対する分岐点と理解できる．Ni3 Al では Spinodal 不安定と Born 不安定が同じひずみであり，Ni も両者

72

の臨界ひずみは極めて近く，対称性が崩れる方向に変形が進むことが予想される．なお，圧縮側では Ni，Ni3 Al とも (3.41)∼(3.43) 式いずれの値も正の値をとり，不安定となる点は存在しない．すなわち理想的な静水圧圧縮を行った場合，変形抵抗は常に増加し，系は常に安定である．一方，第 6 章の Case(b) の [001] 圧縮で観察したように，系の変形の等方性が崩れると，格子不安定となるひずみが現れるものと思われる．

73

~

1

2

20

2
11 -

B

100

Magnitude ofB

11

Magnitude of

30

3

10

GPa

200

2
12 ,

10

12

B +2B , GPa

~

0

Ni
Ni3Al

-100
-0.3

-0.2

-0.1

0

0.1

Applied strain

0.2

0.3

ε33

11

12

0

0.1

0
Ni
Ni3Al

-10
-0.3

-0.2

-0.1

0

0.1

Applied strain

B +2B

(a)

10

(b)

2
11 -

B

0.2

0.3

ε33

2
12

B

200

Magnitude of44

B , GPa

300

100

0
Ni
Ni3Al

-100
-0.3

-0.2

-0.1

Applied strain

(c)

0.2

0.3

ε33

B

44

Fig.7.3 Changes in the minor determinants of detBIJ under hydrostatic
tension/compression

Table 7.1 Elastic stiﬀness coeﬃcients under hydrostatic tension/compression
ε33
B11 (GPa)
B12 (GPa)ε33
B11 (GPa)
B12 (GPa)ε33
B44 (GPa)

Ni
Ni3 Al0.140 0.150 0.150 0.160
15.0 12.1 16.5
9.2
-4.0
-7.1
-6.1 -11.1
0.150 0.160 0.150 0.160
12.1
3.4
16.5
9.2
-7.1 -11.8 -6.1 -11.1
0.210 0.240 0.210 0.240
0.6
-2.4
0.0
-3.6

74

7.3

結言

静水圧引張または圧縮を受ける Ni および Ni3 Al 単結晶について解析を行い，応力–
ひずみ曲線のピークから理想強度を評価するとともに，各ひずみ下の格子安定性を評価した．得られた結果を要約して以下に示す．

(1) 応力–ひずみ曲線は，Ni，Ni3 Al のいずれも引張側にのみ極大値を示し，圧縮側は単調減少して極値を示さなかった．引張のピークより得られる静水圧引張における理想強度は，Ni が ε33 =0.150，σ33 =30.0GPa，Ni3 Al も同じひずみ ε33 =0.150 で応力は σ33 =28.3GPa となる．圧縮側では極値を示さないため，応力–ひずみ関係
からは静水圧圧縮における理想強度は定義できない．

(2) 各ひずみ下における格子不安定性を評価した結果，引張側では Ni，Ni3 Al いずれ
も Spinodal 条件が最初に現れる不安定となった．ただし，Ni3 Al では同時に Born
条件も負となった．Ni は Spinodal 条件，Born 条件，B44 条件の順で，それぞれ

ε33 =0.150，0.160，0.240 で負となり，Ni3 Al は ε33 =0.160 において Spinodal 条件と Born 条件が同時に負となり，B44 条件が ε33 =0.240 で負となった．圧縮側では系の安定性は常に増加する結果が得られた．

(3) (2) の引張における Spinodal 不安定は Ni，Ni3 Al のいずれも応力–ひずみ曲線のピークにほぼ一致した．しかしながら，Spinodal 不安定前後の弾性剛性係数変化を調べた結果，引張 (膨張) に対する抵抗を表す B11 は負になっていなかった．このことは，[001] 引張の場合と異なり，静水圧引張における応力–ひずみ関係のピー
クは不安定膨張を生じてバラバラになる点ではなく，別の結晶構造への相変態点を表すものと考えられる．

(4) 静水圧圧縮下では系の安定性は常に増加する傾向が得られた．このことは，理想的な静水圧圧縮下では，応力–ひずみ関係および格子安定性からも変形限界は存在

75

せず，変形抵抗は常に増加し，系は常に安定側となることを示している．しかしながら，第 6 章の横方向拘束下の [001] 圧縮で示したように，系の変形の等方性が崩れると格子不安定となるひずみが現れる．

76

第8章
結論
本研究では，Ni 基超合金の γ/γ ′ 界面の力学特性解明につながる新たな知見を得ることを目的として，Ni および Ni3 Al それぞれ単結晶について，[001] 方向引張 · 圧縮，静水圧引張 · 圧縮を第一原理分子動力学法により行い，応力–ひずみ関係を求めるとともに，各ひずみ下の格子不安定性について，弾性剛性係数の正値性から議論した．本論文で得られた結果を以下に総括する．
第 4 章では，まず，Ni および Ni3 Al 単結晶の解析に必要なカットオフエネルギー Ecutおよび k 点数を，系のエネルギーの収束性から決定する予備解析を行い，Ni，Ni3 Al
とも第 5 章以降に行う大変形領域の解析でも高い精度を保てるような Ecut および k 点数を決定した．次に，決定した Ecut および k 点数を用いて Ni，Ni3 Al それぞれの平衡状態における格子定数等の物性値を求める解析を行い，理論値とのずれが誤差範囲内に十分収まることを示した．
第 5 章では，[001] 方向に引張変形を受ける Ni，Ni3 Al の単結晶について，横方向に等方的な Poisson 収縮を仮定した横方向応力 0 の場合 (Case(a))，および，横方向変形を拘束した場合 (Case(b))，の 2 通りの解析を行った．応力–ひずみ曲線のピークより評価した理想引張強度，および，最初に格子不安定となる臨界ひずみを弾性限界とした降伏強度をまとめて表 8.1 に示した．[001] 引張では最初に現れる不安定は Born 条
件に対するものであり，横方向変形に関する弾性剛性係数 B11 ，B12 の大小関係が逆転
し，非等方変形を開始する変形分岐点であった．なお，横方向変形の拘束によって，格子不安定となるひずみが上昇し，応力–ひずみ関係のピーク点が低ひずみ側にシフトす77

るが，Ni と Ni3 Al の強度の大小関係は変わらない．
第 6 章では，[001] 方向に圧縮変形を受けるそれぞれの単結晶について，第 5 章と同様に，横方向に等方膨張を仮定した場合 (Case(a))，および，横方向変形を拘束した
場合 (Case(b))，の 2 通りの解析を行った．応力–ひずみ関係より評価した理想圧縮強度と最初の格子不安定ひずみに基づく強度をまとめて表 8.2 に示した．[001] 圧縮で最初に現れる不安定は Spinodal 条件であり，Case(a) では応力–ひずみ関係のピークと一致する．Case(b) では応力–ひずみ曲線は極値を示さず，理想圧縮強度は定義できないが，Case(a) と同様，Spinodal 条件が負となるひずみが存在した．これより得られる

Ni，Ni3 Al の強度は Case(a) の場合のみ大小関係が逆転して Ni<Ni3 Al となることが示された．また，Case(a) では応力–ひずみ関係が 3 次曲線的形状となったが，これについて結晶構造変化における Bain の関係から考察した．
第 7 章では，静水圧引張または圧縮を受ける Ni，Ni3 Al の単結晶について，応力–ひずみ関係のピークから理想強度を評価するとともに，各ひずみ下の格子不安定性について検討した．応力–ひずみ関係より評価した理想圧縮強度と最初の格子不安定ひずみに基づく強度を表 8.3 に示す．引張側は，応力–ひずみ関係のピークと格子不安定となる点が一致した. 格子不安定は Spinodal 条件に対するものであったが，弾性剛性係数の変化から，膨張に対する不安定ではなく相変態に対するものであることが示唆された．一方，圧縮側は，応力–ひずみ関係は単調に減少して極値を持たず，またいずれの格子不安定条件も負とならなかった．このことは，理想的に等方な圧縮条件下では系の変形抵抗は常に増加し安定側となることを示している．しかしながら，第 6 章のCase(b) の [001] 圧縮で示したように，系の変形の等方性が崩れると格子不安定となる
ひずみが現れる．

78

Table 8.1 Ideal strength and elastic limit evaluated with the lattice instability
under [001] tension

Ni
Ni3 AlNi
Ni3 AlIdeal Strength
Lattice Instability
Born
44
ε33
σ33 (GPa)
ε33
εSpinodalεB
σ33 (GPa)33
33
Case(a) : stress free condition (σ11 =σ22 =0)
0.390
35.3
0.110
0.390
0.240
13.5
0.360
34.0
0.090
0.360
0.300
10.8
Case(b) : deformation constraint condition (ε11 =ε22 =0)
0.330
35.1
0.190
0.360
0.270
28.9
0.330
33.6
0.130
0.360
0.360
21.3

Table 8.2 Ideal strength and elastic limit evaluated with the lattice instability
under [001] compression

Ni
Ni3 AlNi
Ni3 AlIdeal Strength
Lattice Instability
Spinodal66
ε33
σ33 (GPa) ε33
εB
σ33 (GPa)33
Case(a) : stress free condition (σ11 =σ22 =0)-0.100-5.44-0.110-0.180-5.44-0.120-5.94-0.120-0.210-5.94
Case(b) : deformation constraint condition (ε11 =ε22 =0)––-0.150-0.300-49.5––-0.140-0.300-33.9

Table 8.3 Ideal strength and elastic limit evaluated with the lattice instability
under hydrostatic tension/compression

Ni
Ni3 AlIdeal Strength
Lattice Instability
Born
44
ε33
σ33 (GPa)
ε33
εSpinodalεB
σ33 (GPa)33
33hydrostatic expansion/compression (ε11 =ε22 =ε33 )
0.150
30.0
0.150
0.160
0.240
30.0
0.150
28.3
0.160
0.160
0.240
28.3

79

参考文献
(1) Jhi,S.H., Louie,S.G., Cohen,M.L., and Morris,Jr.,J.W., Phys. Rev. lett, 87(2001), 075503.
(2) 太田裕之，三浦英生，北野誠，材料，45, (2002), 1322.
(3) Abrahan,F.F., Materials Science for the 21st Century; The Society of Meterials
Science, A (2001), 195.
(4) Vashishta,P., Bachilechner,M.E., Nakano,A., Campbell,T.J., Kalia, R. K.,
Kodiyalam, S., Ogata, S., and Shimojo, F., Materials Science for the 21st Century; The Society of Meterials Science, A (2001), 213.
(5) Li,J., Cai,W., Chang,J., and Yip,S. , Materials Science for the 21st Century;
The Society of Meterials Science, A (2001), 220.
(6) 香山正憲，田中真梧，材研連講演論文集, (2000)，342.
(7) 尾方成信，北川浩，金属学会誌, 60, (1996), 1079.
(8) 尾方成信，学位論文，材料特性評価のための非経験的分子動力学法の基礎的研究，(1998)(9) Friák,M., Šob,M., and Vitek,V., Phys. Rev. B, 68 (2003), 184101.
(10) Born,M., and Huang,K., Dynamical Theory of Crystal Lattices, (1954) Oxford
UP.
(11) Milstein,F., Phys. Rev. B, 3 (1971), 1130.
(12) Wang,J., Yip,S., Phillpot,S.R., and Wolf,D., Phys. Rev. B, 20 (1993), 4182.
(13) Wang,J., Li,J., Yip,S., Phillpot,S.R., and Wolf,D., Phys. Rev. B, 52 (1995),
12627.
(14) Hill,R., Math. Proc. Camb. Philos. Soc, 77 (1975), 225.
(15) Thompson,J.M., and Shorrock,P.A., Jour. Mech. Phys. Solid., 23 (1975), 21.
80

(16) Milstein,F., and Farber,B., Phys. Rev. lett., 44 (1980), 277.
(17) Milstein,F., Marschall,J., and Fang,H.E., Phys. Rev. lett., 74, (1995), 2977.
(18) Krenn,C.R., Roundy,D., Cohen,M.L., Chrzan,D.C., and Morris,Jr.,J.W., Phys.
Rev. B, 65 (2002), 134111.
(19) 屋代如月，大穂正史，山上勝也，冨田佳宏，材料，52, (2002), 241(20) Beardmore,P., Dias,R.G., and Johnston,T.L., Trans.Metall.Soc., 245, (1969),
1537.
(21) Tien,J.K., and Copley,S.M., Metall.Trans., 3, (1972), 2157.
(22) Giamei,A.F., and Anton,D.L., Metall.Trans., A16, (1985), 1997.
(23) Aihara,T.Jr., Sho,T., and Kawazoe,Y., J.Japan.Inst.Metals, 62, (1998), 978.
(24) 屋代如月，内藤正登，冨田佳宏，日本機械学会論文集，A 編，69(677)，(2002),
210．
(25) Hohenberg,P., and Kohn,W., Phys. Rev., 136, (1964), 864.
(26) Kohn,W., and Sham,L.J., Phys. Rev., 140, (1965), 1133.
(27) Perdew,J.P., and Zunger,A., Phys. Rev.,B, 23, (1981), 5048.
(28) キッテル，固体物理学入門 上 第 7 版 (1998), 丸善.
(29) Bachelet,G.B., Hamann,D.R. and Schlüter,M., Phys. Rev. B, 26, (1982), 4199.
(30) Troullier,N. et al., Phys. Rev. B, 43, (1991), 1993.
(31) Vanderbilt,D., Phys. Rev. B, 4, (1990), 7892.
(32) L.D. ランダゥ，E.M. リフシッツ著，佐々木健，好村磁洋訳，量子力学 1，(1992),
東京図書.

(33) 上田顕，コンピュータシミュレーション, (1990), 朝倉書店.
(34) 香山正憲，固体材料の電子状態の基礎, (1993), 材料学会勉強会資料.
(35) Nielsen,O.H., and Martin,R.M., Phys. Rev. lett., 50, (1983), 697.
81

(36) Nielsen,O.H., and Martin,R.M., Phys. Rev. B, 32, (1985), 3780.
(37) Dacosta,P.G., Nielsen,O.H., and Kunc,K., J. Phys. C : Solid State Phys., 19,(1986), 3163.
(38) Kleinman,L., and Bylander,D.M., Phys. Rev. Lett., 48, (1982), 1425.
(39) 尾口多美夫，佐々木泰造，固体物理, 25, (1990), 857.
(40) Fu,C.L., and Ho,K.M., Phys. Rev. B, 28, (1983), 5480.
(41) Teter,M.P. et al., Phys. Rev. B, 40, (1989), 12255.
(42) King-Smith,R.D., and Vanderbilt,D., Phys. Rev. B, 49, (1994), 5828.
(43) Gillan,M.J., J. Phys. Condens. Matter 1, (1989), 689.
(44) Stich,I. et al. , Phys. Rev. B, 39, (1989), 4997.
(45) Bylander,D.M. et al., Phys. Rev. B, 42, (1990), 1394.
(46) Wallace, D.C., Thermodynamics of Crystals, (1972), Wiley, Newyork(47) 屋代如月，学位論文，金属結晶の変形と破壊の分子動力学シミュレーション並びに局所格子不安定性解析，(1998)(48) Li, J., and Wang, T., J. Phys. Condens. Matter 10, (1998), 9889.
(49) Kresse,G., and Hafner,J., Phys. Rev. B, 47 (1993), 558.
(50) Monkhorst,H.J., and Pack,J.D., Phys. Rev. B, 13, (1976), 5188.
(51) Langreth,D.C., and Perdew,J.P., Phys. Rev. B, 21, (1980), 5469.
(52) Kresse,G., and Hafner,J., Phys. Rev. B, 54 (1996), 11169.
(53) Voter, A., and Chen, S.P., Mat. Res. Soc. Proc., Vol.82, (1987), 175.
(54) 香山正憲:私信「プログラム使用に関するノート Ver 2.1」による.
(55) Kresse,G., and Furthmüller,J., VASP the GUIDE, (2003).
(56) Asato,M, and Hoshino,T., J. Japan Inst. Metals, Vol. 63, No. 6 (1999), pp. 676
82

(57) Milstein,F., and Chantasiriwan,S., Phys. Rev. B, 58, (1995), 6006.
(58) 高木節雄，津崎兼彰，マテリアル工学シリーズ 2 材料組織学，朝倉書店，(2000)

83

付録 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 e284(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] などがある．

85

付録 B
関連発表論文 · 関連講演論文
関連発表論文(1) 屋代如月，大穂正史，山上勝也，冨田佳宏，[001] 方向単軸引張を受ける Si およびAl 単結晶の第一原理分子動力学法による格子不安定性解析，材料，52(3), (2002),
241–246
関連講演論文(1) 山上勝也，屋代如月，冨田佳宏，第一原理計算を用いた Ni 及び Ni3 Al の理想格子
不安定性解析，日本機械学会第 16 回計算力学講演会講演論文集，(2003)86謝辞
本研究を進めていくにあたって終始懇切丁寧な御指導を賜った冨田佳宏教授に深く感謝致します．研究内容について様々な御助言を頂いた安達泰治助教授に感謝致します．研究全般にわたり直接御指導頂いた屋代如月助手に感謝致します．研究室の環境作りにご尽力くださった古宇田由夫技官に感謝致します．
さらに，研究室に配属された当初から，倦むことなく基礎から御指導いただいた大穂正史氏 (現 · 住友電気工業株式会社) をはじめとする諸先輩方に感謝致します．また，励ましあい助け合いながら共に研究生活を過ごした，同輩の吉田十義氏，河東万智子
さん，金井雅之君，中村達彦君，松尾守君，山本裕介君，ならびに，同じ研究分野で有益な議論を交わした久保圭佑君ら講座内諸氏に心から御礼申し上げます．
固体力学研究室で研究生活を共に過ごした，斉藤賢司氏 (現 · 株式会社神戸製鋼所)，
佐藤克也氏，澤田幸秀氏 (現 · ダイハツ工業株式会社)，伊藤智啓氏 (現 · デンソー株式会社)，酒井誠君 (現 · アイシン精機株式会社)，宮崎健治君 (現 · 株式会社タナベ)，黒瀬史治君，内之浦哲也君らにも深く感謝致します．
最後に，予備校や大学 · 大学院の長きに渡る学生生活を経済的かつ精神的に支えてくれた，両親と祖父母に最大限の感謝を致します．ありがとうございました．

87