局所エネルギー・局所応力の第一原理計算法開発と材料界面への適用
香 山 正 憲
田 中 真 悟
椎 原 良 典1)
2)

.

は じめ

fe(r… )d…r

Etot＝に(1)Vsab＝密度汎関数理論に基づく第一原理計算(1)は，物質の結晶構

1V&Etot
1＝&eabV

ftab(r… )d…r(2)V

造から欠陥，表面，粒界，界面等の安定構造や原子・電子挙

詳細は文献( 12 )( 13 )にあるが，エネルギー密度は，全エネ動を高精度に解明する強力な手段である．ノルム保存擬ポテルギーの各項( kinetic 項， Hartree 項，交換相関エネルギーンシャル法(NormConserving PseudoPotential; NCPP 法)項，電子イオン相互作用項，イオンイオン項など)を式やウルトラソフト擬ポテンシャル法( Ultra Soft

Pseudo ( 1 )の形に書き下して構築される．応力密度は，エネルギ

Potential; USPP 法)，射影補強波法(Projector Augmented

ー密度の各項の積分形に式( 2 )の歪での偏微分を施して導PAW 法)(2)(3)など，価電子のみを扱う効率的で高精出される．最終的に通常の平面波基底法で求めた波動関数，度な平面波基底法を CarParrinello に始まる高速計算技法(1)

価 電 子密 度 分 布を 用 い て e ( r… )， tab ( r… )が 計算 さ れ ， FFTと組み合わせることで大規模構造の扱いが可能となった．

(Fast Fourier Transformation高速フーリエ変換)の等間隔CASTEP ， VASP といった汎用ソフトの普及により広範なメッシュの離散データとして与えられる．

Wave;

応用が行われている．しかし，平面波基底法では，全エネル

筆者らは，最も進化した平面波基底法である PAW 法(2)(3)

ギーや応力(4)がスーパーセル全体の積分値，平均値としてしについてエネルギー密度，応力密度を定式化し，産総研開発

か求まらない．欠陥や粒界・界面を扱う見地からは，エネルの PAW 法 汎 用 コ ー ド QMAS ( Quantum Material Simula-

ギーや応力の局所分布が重要である．筆者らは，平面波基底tor )に実装した(14) ． PAW 法は，擬ポテンシャル法( NCPPの第一原理計算法において，局所エネルギー・局所応力を計法，USPP 法)と異なり，原子近傍の電子密度分布や波動関算する手法・プログラムの開発に取り組み，金属等の表面，数を正確に再現する． PAW 法での各密度の表式は文献( 5)欠陥，粒界，界面への適用を進めている(5)(11)．本稿では，を参照．

手法の概要と最近の適用例を紹介し，今後の課題，展望を論じる．

しかし，エネルギー密度，応力密度には，電子の kinetic
項表式の自由度(対称形か非対称形か)があり，選択により値が異なるため non uniqueness の問題( gauge dependent 問.

局所エネルギー・局所応力計算法

題(12)(13) )が生じ，応用が広がっていない．そこで，エネルギー密度，応力密度の gauge dependent 項( kinetic 項の対エネルギーや応力の局所分布を計算する方策として，エネ称形と非対称形の差)が積分してゼロとなるような局所領域

ルギー密度，応力密度をスーパーセル内の実空間メッシュ点

Vi にセル内を分割することを考える．各領域 Vi 内で e(r… )，データとして計算する方法が，平面波基底法( NCPP 法，

tab ( r… ) を 積 分 す れ ば gauge dependent 項 は 消 え る の で ，USPP 法)について提案されている(12)(13) ．通常の方法で与unique な物理量として局所エネルギー E itot，局所応力 siab がえられるスーパーセルの全エネルギー Etot ，応力テンソル求められる(5)．

sab に 対 し ， 実 空 間 セ ル 周 期 関 数 と し て エ ネ ル ギ ー 密 度e ( r… )，応力密度 tab ( r… )が以下のように定義できる( V はスーパーセル体積)．

エネルギー密度，応力密度の gaugeindependent 条件(局所領域 Vi の満たすべき条件)は，各々以下になる．

ff:2r(r… )d…r 3＝

Vi:r(r… )･ …n⊥dS＝0(3)Si

産業技術総合研究所ユビキタスエネルギー研究部門1)首席研究員，2)主任研究員(〒5638577 池田市緑ヶ丘 1
831)
東京大学生産技術研究所助教
Development of AbInitio LocalEnergy and LocalStress Schemes and Applications to Materials Interfaces; Masanori Kohyama,Shingo Tanakaand Yoshinori Shiihara(Research Institute for Ubiquitous Energy Devices, National Institute of Advanced Industrial
Science and Technology, Ikeda. Institute of Industrial Science, The University of Tokyo, Tokyo)
Keywords: firstprinciples calculation, local energy, local stress, energy density, stress density, grain boundary, impurity segregation
2014年 5 月25日受理[doi:10.2320/materia.53.405]

ま て り あ 第53巻 第 9 号(2014)
Materia Japan



ff

:a:b r(r… )d…r3＝

Vi:b r(r… )e… a･ …n⊥dS＝0(4)Siband の二次モーメント強結合近似解析(18) ( Friedel model )で，各金属の表面応力値がうまく再現できることが判明しr(r… )は価電子密度分布である．各式は，エネルギー密度，応た．従って，論争されてきた遷移金属表面の引張応力の起源

力密度内の電子の kinetic 項の対称形と非対称形の差が，各は，表面原子の配位数減少による dband 幅減少であると結々 :2r(r… )，:a:b r(r… )に比例するためである．導出の詳細は論づけられる．

文献(5)を参照．各式は，セル全体での積分は厳密にゼロであるが，局所領域 Vi 内の積分でもゼロとならねばならな
い．各式の最終形は，局所領域境界全体の表面積分で表されている．

Al, Cu の結晶粒界(8)図， は， Al と Cu の∑ 9( 221 )粒界の計算結果で，価電子密度分布( Al 粒界のみ)と各原子の局所エネルギーを示式( 3 )，( 4 )を満たす gauge independent な領域にスーす．原子埋め込み法( EmbeddedAtom Method; EAM)ポテ

パーセルを分割する具体的方法として，layerbylayer 法とンシャル(19)(20)での値も載せている．Al 粒界では，配位数のBader 積分法を確立した(5)(7)(8)．layerbylayer 法では，層減った原子同士でボンド長の短い再構成ボンドが形成され，状構造の表面 slab や積層欠陥について，表面(積層欠陥面)価電子が集中している．そのボンドの一方を担う原子( No.に平行な面で空間を区切って式( 3 )，( 4 )を満たす原子層11, 12 )は，局所エネルギーが完全結晶比でマイナスであ領域を設定する．表面(積層欠陥面)に平行な二次元周期性かる．相手の原子( No. 9, 14 )もエネルギーは低い．一方，対ら条件式は一次元問題に還元でき，r(r… )を面上で積分した一面に近接原子のない原子(No. 10, 13)は，最も高エネルギー次元 profile šr ( z )の勾配ゼロの条件から，境界面位置 z が決である．粒界全体のエネルギー上昇は，局所エネルギーの和まる．と EAM ポテンシャルの結果で大きな差はない．しかし，一般の欠陥や粒界では，各原子領域への分割を考える．式

EAM ポテンシャルでは電子の集中した再構成ボンドの負の( 3 )の面積分形は，電子密度を勾配に基づいて各原子領域エネルギーは再現できない．価電子が構造乱れに敏感に応答に割り当てる Bader 積分(15) と同じ条件式であり，また式する Al では EAM ポテンシャルの精度は落ちる． Cu 粒界( 4 )の対角和が式( 3 )と同値である．従って，Bader 積分の方法により原子領域(Bader 領域)を設定し，その領域でエネルギー密度，応力密度を積分して，局所エネルギー，局所応力テンソルの対角和(静水圧)を求めることができる．しかし，FFT メッシュの離散データで与えられる r ( r… )に対して，Bader 領域を高精度に決定することは容易ではない．Bader積分の各種アルゴリズムのうちで，YuTrinkle 法(16)が最も信頼がおけることを確認した(8)．これは，確率流の方程式からメッシュ点での「重み」を導出し，離散メッシュの弱点を補う方法である．

.


適用例金属の表面Al(111)表面に layerbylayer 法を適用し，原子層毎の表面平行応力が振動することを見出した(5)．表面第一層に強い引張応力，第二層に圧縮応力，第三，四層に小さい引張応力など，数原子層に渡り表面平行応力が振動する．引張応力と圧縮応力は各原子層の in plane 電子密度値の増加と減少に対応している． Al の価電子は，配位数や原子構造乱れに顕著に応答し再分布する(17) ．表面第一層の配位数減少により ， 電 子 密 度 に 顕 著 な 振 動 が 生 じ ， layer by layer のFriedellike な振動として数原子層に渡り存在し，応力振動が発現すると言える．
fcc 構造の各種の遷移金属(late transition metal)の表面応

力 の layer by layer 法 解 析(6) で は ， 顕 著 な 振 動 は 見 ら れず，表面第一層に強い引張応力が確認された．その応力値と表面の dband 幅や dband center shift に相関が見られた．
一方，表面での配位数や原子間距離の変化を取り入れた d 図1Al の∑ 9 ( 221 )粒界の価電子密度分布と各 Bader領域の局所エネルギー値(完全結晶比， eV )(8) ．
各々二つの{220}原子面のもの．一番下は各原子の局所エネルギーの棒グラフ．緑が第一原理計
算，青が EAM ポテンシャルの結果．

ミ ニ 特 集では d 電子が支配するため，価電子の顕著な再分布や負のエネルギーは見られない．EAM ポテンシャルの結果との違いは Al ほど顕著ではない．
図に粒界の各原子の応力の結果を示す． Al も Cu も，界面で原子が近接して出会うところ( No. 11, 12 )で圧縮応力，原子間が離れているところ(No. 10, 13)で引張応力が生じている． EAM ポテンシャルとの比較では， Al 粒界の再構成ボンドの所で差が顕著であるが， Cu 粒界での一致はかなり良い．
両方の粒界で，原子間が接近しすぎる場所と離れてしまう場所が必然的に生じている．前者(No. 11, 12)では，原子体積や原子間距離が小さく，低エネルギーで圧縮応力が生じ(tighter site(22))，後者(No. 10, 13)では，原子体積や原子間距離が大きく，高エネルギーで引張応力が生じる( loosersite(22) ) ． 両 サ イ ト は 不 純 物 偏 析 に 対 照 的 な 効 果 を 持つ(10)(11)．


図2

Cu の∑9(221)粒界の各 Bader 領域の局所エネルギー値(完全結晶比，eV)(8)．二つの{220}原子面のもの．一番下は各原子の局所エネルギーの棒グ
ラフ．緑が第一原理計算，青が EAM ポテンシャルの結果．

図に bcc Fe の ∑ 11 ( 332 )粒界の第一原理計算結果を示す．古典分子動力学法による構造(23) を出発点に緩和させたもので，〈110〉方向から見て五員環，三員環からなる 53 ユニットと完全結晶ユニットの配列構造を持つ．図に示す通図4

図3Al(上)と Cu(下)の∑9(221)粒界の各原子の応力(対角成分の和)(8) ．正が引張，負が圧縮．第一原理計算(緑)と EAM ポテンシャル(青)の結果を示 す ． EAM ポ テ ン シ ャ ル に よ る 原 子 応 力 はEgami らの定式化(21)に基づく．

ま て り あ 第53巻 第 9 号(2014)
Materia JapanFe の結晶粒界(7)(9)Fe の ∑ 11 ( 332 )粒界の原子配列，原子領域毎の局所エネルギー，局所応力(対角和)，局所Young 率(パネルの上から順に)(9) ．原子配列図の赤線は界面の構造ユニットを示す．横軸は原子位置で，応力は正が引張，負が圧縮，局所エネル
ギー値は完全結晶の値との差．局所 Young 率のパネルの青，赤の横線は，各々バルク結晶値と粒界セルの平均値．界面第一層の 1b 原子の値は，黒丸で示す．

り，そこでのエネルギー上昇，局所応力変化は構造ユニットした場合の不純物原子の局所エネルギー変化， T3 は不純物の原子に集中している．他の回転角の粒界においても，同様偏析に伴う粒界の周囲の金属原子の局所エネルギー変化の構造ユニットの混合配列で構成され，そこでのエネルギ(pure 粒界の場合との差)，T4 は不純物原子が完全結晶に単ー・応力の分布がやはりユニット配列に依存することが観察独で存在する系から不純物原子がなくなる際の周囲の金属原された．Al，Cu と同様，各粒界で高エネルギーで引張応力子の局所エネルギー変化である．計算結果を表に示す．

を持つ原子(looser site，図 4 の 1a， No. 3 の原子)と低エネ

1asite 偏析では，T1 の絶対値が大きいため Eseg が負にな

ルギーで圧縮応力を持つ原子( tighter site ，図 4 の 1b ， No.

る．元の粒界で局所エネルギーの高い 1a 原子が置換される2，No. 4 の原子)が存在している．また，界面の構造乱れとことで偏析利得が生まれており， InitialState Effect と呼ぶ

局所 Young 率が相関することが見出された．局所 Youngことができる． 1b site 偏析では，元の粒界の 1b 原子は安率(9)は，セル全体の界面垂直方向の微小引張と微小圧縮での定なため T1 の絶対値は小さい．しかし， 1b site に Si が来応力変化を各原子の引張方向の歪変化で割ったもので，応力ると周囲の Fe 原子と安定な SiFe ボンドを作り， T2 と T3変化はセル全体の平均応力で与え，原子の歪変化は引張と圧が大きく下がり，大きな偏析利得が得られる．典型的な縮下のボロノイ体積変化で近似した．局所応力等は使用してFinalState Effect(不純物ホストの安定な結合形成による利いないが，粒界の力学特性の解明に有用である．
局所磁気モーメントは，looser site で増加，tighter site で得)と言える．一方，1a＋1b 構造では，偏析利得がなくなる(Eseg が正)．Si 同士の隣接で T2 の符号が正になるためであ減少する傾向がある(24) ． majority spin と minority spin のる．以上のような機構解明は，Eseg の値しか得られない従来局所状態密度(LDOS)の占有部分の差で局所磁気モーメント法では不可能であった．こうした解析は，表面の原子・分子が決まる．bcc Fe では majorityspin の dband DOS はほぼ吸着の機構解明にも有用である．

埋まり， minority spin DOS が谷の部分まで埋まる(25) ．原.

子体積や原子間距離の大きい looser site では dband 幅が狭

今 後 の 課 題

く，minorityspin のフェルミレベルの DOS が増え，占有部分が減り，スピン分極が増す．原子体積や原子間距離の小さ手法の概要と最近の適用例を紹介した．今後，異相界面やい tighter site では， d d 混成の増加で minority spin DOS

格子間型の欠陥・不純物，水素やリチウムの吸収，さらに変の占有部分のピークが増え，分極が減り，磁気モーメントが形や反応過程への適用が興味深い．化合物や半導体について減る．minorityspin DOS の形と占有状況が局所原子構造に依存することで，局所磁気モーメントの局所構造依存性が生じる．
Fe 粒界への Si 偏析(10)Si を添加した Fe は，Si 濃度に依存した特異な弾性定数変化や塑性挙動が注目されている．粒界偏析も起きると考えられ，偏析の様子や粒界特性への効果が興味深い．図のモデルに基づき ∑ 11 ( 332 )粒界の Si 偏析の第一原理計算を行った．
偏析エネルギーは，Si 偏析した Fe 粒界セルと pure な Fe結晶セルの全エネルギーの和から， pure な Fe 粒界と Si 原子が単独で存在する Fe 結晶セルの全エネルギーの和を差し引いて算出する．各セルの全エネルギーを原子毎の局所エネルギーに分解して整理すると，四項の和 T1 ＋ T2 ＋ T3 ＋ T4で表される(図)． T1 は粒界にあった金属原子が不純物と

図6

粒界偏析エネルギーの局所エネルギー分解(10)(11)．

表1Fe の ∑11(332)粒界の Si 偏析エネルギーの局所エ ネ ル ギ ー 分 解(10) ． 各 項 の エ ネ ル ギ ー は eV /cell ，偏析エネルギー Eseg は eV / Si atom ． m はセル内の Si 原子数．

入れ替わり完全結晶に移動する場合の金属原子の局所エネル
ギー変化， T2 は結晶中に単独に存在する不純物が粒界偏析

図5

Fe の∑ 11( 332)粒界の Si 原子(赤，紫)偏析モデル(置換型)(10)．1asite，1bsite，及び両方のサイト(1a＋1b)への偏析を扱う．

偏析モデル

mT1T2T3T4Eseg∑11(1a)2－1.02－0.10－1.91＋2.81－0.11∑11(1b)2－0.13－0.20－3.75＋2.81－0.63∑11(1a＋1b)

4－1.15＋0.75－4.45＋5.63＋0.19

ミ ニ 特 集は，複雑な電子密度分布や電子移動に対応した局所領域決定法の検討も必要かもしれない．局所応力について，その物理的な意味，巨視的応力との関係，原子に働く力との関係等の検討が重要である．静電エネルギーの gauge dependent 問
題( Coulomb 形か Maxwell 形かの選択)について， Maxwell形を選択しているが(4)(12)(13) ，別の処理法(26) もあり，検討の余地があろう．いずれにしても，第一原理計算からの情報を最大限引き出すための手法開発は極めて重要である．また，エネルギー・応力の空間分布の粗視化の観点から，ミクロとマクロ(連続体)をつなぐ手段としての役割も期待される．
共同研究者の Somesh Kr. Bhattacharya(前産総研ポスドク，現京都大学), Hao Wang, Vikas Sharma(産総研ポスドク)，石橋章司(産総研)の各氏に感謝します．新学術領域「バルクナノメタル」 22102003 ，科研費 23710107 ， JST 産学共創「ハミルトニアンによる材料強度設計」，構造材料元素戦略研究拠点 ESISM ，高性能汎用計算機高度利用事業「HPCI 戦略プログラム」の支援を受けた．

文献(11) H. Wang, M. Kohyama, S. Tanaka and Y. Shiihara: submitted.
(12) N. Chetty and R. M. Martin: Phys. Rev. B, 45(1992), 6074.
(13) A. Filippetti and V. Fiorentini: Phys. Rev. B, 61(2000), 8433.
(14) S. Ishibashi, T. Tamura, S. Tanaka, M. Kohyama and K.
Terakura: Phys. Rev. B, 76(2007), 153310.
(15) R. F. W. Bader: Atoms in Molecules A Quantum Theory, Oxford UP, Oxford, (1990).
(16) M. Yu and D. R. Trinkle: J. Chem. Phys., 134(2011), 064111.
(17) T. Uesugi, M. Kohyama and K. Higashi: Phys. Rev. B, 68(2003), 184103.
(18) M. Lannoo and P. Friedel: Atomic and Electronic Structure ofSurfaces: Theoretical Foundations, SpringerVerlag, Berlin,(1991).
(19) Y. Mishin, D. Farkas, M. J. Mehl and D. A. Papaconstantopoulos: Phys. Rev. B, 59(1999), 3393.
(20) Y. Mishin, M. J. Mehl, D. A. Papaconstantopoulos, A. F. Voter
and J. D. Kress: Phys. Rev. B, 63(2001), 224106.
(21) T. Egami: Prog. Mater. Sci., 56(2011), 637.
(22) D. I. Thomson, V. Heine, M. C. Payne, N. Marzari and M. W.
Finnis: Acta Mater., 48(2000), 3623.
(23) H. Nakashima and M. Takeuchi: ISIJ, 86(2000), 357.
(24) M. Cak, M. Sob and J. Hafner: Phys. Rev. B, 78(2010),054418.
(25) D. G. Pettifor: Bonding and Structure of Molecules and Solids,Oxford UP, Oxford, (1995), Sec. 8.6.
(26) M. Yu, D. R. Trinkle and R. M. Martin: Phys. Rev. B, 83(2011), 115113.
★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★

( 1 ) R. M. Martin: Electronic Structure: Basic Theory and PracticalMethods, Cambridge UP, Cambridge, (2004).
( 2 ) P. E. Bl äochl: Phys. Rev. B, 50(1994), 17953.
( 3 ) G. Kresse and D. Joubert: Phys. Rev. B, 59(1999), 1758.
( 4 ) O. H. Nielsen and R. M. Martin: Phys. Rev. B, 32(1985), 3780.
( 5 ) Y. Shiihara, M. Kohyama and S. Ishibashi: Phys. Rev. B, 81(2010), 075441.
( 6 ) Y. Shiihara, M. Kohyama and S. Ishibashi: Phys. Rev. B, 87(2013), 125430.
( 7 ) S. Kr. Bhattacharya, S. Tanaka, Y. Shiihara and M. Kohyama:
J. Phys.: Condens. Matter, 25(2013), 13500.
( 8 ) H. Wang, M. Kohyama, S. Tanaka and Y. Shiihara: J. Phys.
Condens. Matter, 25(2013), 305006.
( 9 ) S. Kr. Bhattacharya, S. Tanaka, Y. Shiihara and M. Kohyama:J. Mater. Sci., 49(2014), 3980.
(10) S. Kr. Bhattacharya, M. Kohyama, S. Tanaka and Y. Shiihara:
J. Phys. Condens. Matter, 26(2014), 355005.

ま て り あ 第53巻 第 9 号(2014)
Materia Japan

香山正憲
1985年 東京大学大学院工学系研究科博士課程中退，工学博士
1985年 工業技術院大阪工業技術試験所
2001年 産業技術総合研究所生活環境系特別研究体
2004年 同ユビキタスエネルギー研究部門，現在に至る
専門分野計算材料科学，ナノ界面科学，エネルギー材料◎粒界・界面の第一原理計算，手法・プログラムの開発，電顕との連携研究に従事．https://staff.aist.go.jp/mkohyama/★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★

香山正憲

田中真悟

椎原良典

