第一原理計算による局所応力解析
椎原良典 1, 香山正憲 2, 田中真悟 2, 田村友幸 1, 石橋章司 1
1 産業技術総合研究所計算科学研究部門, 2 産業技術総合研究所ユビキタスエネルギー研究部門

y.shiihara@aist.go.jp
概要：ペタフロップス級大型計算機上で動作する第一原理計算ソフトウェアの主な対象である表面，界面，粒界，欠陥，アモルファス等を含む大規模原子系においては，局所的な応力状態に関する情報が様々な原子論的メカニズムを議論する上で有用である．本稿では，第一原理計算法による局所応力解析の具体的方法とその適用例について述べる．

1

はじめに近年，表面原子の再構築，界面接合，表面上での自己組織化，粒界破壊等の原子論的メカニズムを議論するための一方法として，量子力学に基づく第一原理応力シミュレーションが注目されている．しかしながら，従来の第一原理計算から得られる応力はスーパーセル上で平均化された巨視的な物理量であり，局所的な応力分布を得ることはできない．ペタフロップス級大型計算機上で動作する第一原理計算ソフトウェアの主な目的は表面，界面，粒界，欠陥，アモルファス等を含む，原子を多数擁する系の解析であり，そのような系においては，セル体積上で平均化された巨視的応力から得られる情報は極めて少ない．従って，原子の不均一構造を反映した局所的応力分布を得ることが可能な，新しい計算法の確立が必要である．
局所応力計算における難点の一つは，局所応力の計算に用いる応力密度が一意に決まる物理量でないことにある[1]．このことは密度汎関数法におけるエネルギー密度の非一意性の問題に由来する．
従来の研究では，単位格子長等の特性長を用いて局所的に平均化することにより，この非一意性を回避していた．しかしながら，この手法を特性長の決まらない原子緩和された系に対して適用することは難しい．
本研究では運動エネルギー密度の非一意性に着目した手法を用いてこの非一意性の問題を回避する．運動エネルギー密度内の非一意応力密度項の積分値が消える部分領域を探索し，その領域内で応力密度を積分することで一意な局所応力を得る．
本稿では手法の詳細を示し，基礎的解析例として
diamond/cBN(110)界面の一次元局所応力分布の計算結果を示す．

2計算手法
応力密度の定義について述べる[2]．巨視的応力

σ ij は全エネルギーの歪微分として式(1)で定義される．

σ ij =1 ∂etot (r )dVV ∂ε ij ∫(1)

V はセル体積であり，etot は密度汎関数法に基づく全エネルギー密度である．式(1)から，σ ij を応力密度 τ ij の汎関数として式(2)で表す．

σ ij (r ) = ∫ τ ij (r )dV(2)すなわち，応力密度はエネルギー汎関数を歪微分したものの被積分項として定義される．
本研究では，エネルギー密度の非一意性の問題は運動エネルギーにのみ存在すると仮定する．密度汎関数法における一電子近似の下で，運動エネ
ルギー密度 t を式(3)で表す．
1t (r ) = fi ∇φi* (r ) ⋅∇φi (r ) + α∇ 2 ρ (r )(3)2
ここで φi は一電子波動関数， fi は軌道の占有数，ρ (r ) は電子密度，αは任意の定数である．t をセル全体で積分すると，任意のαについて運動エネルギーは同じ値となる．すなわち，運動エネルギー密度にはαの任意性があり，そこから得られる応力密度は一意に決まらない．この非一意性の問題を除くため，∇ 2 ρ (r ) の応力密度の積分値がゼロとなる部分領域を求める．任意の部分領域 V ′ 上における∇ 2 ρ (r ) の応力密度の積分値 σ ij L を計算すると式(4)となる．
∂ ∂ρ (r )σ ij L = ∫ ′dVV ∂r∂rji(4)式(4)の体積分を面積分に変換すると，部分領域が満たすべき条件として式(5)が得られる．
∂ρ (r )∫S ′ ∂rj ei ⋅ ndS = 0 (i, j = 1, 2,3) (5)
ここで ei は x, y, z 軸上の単位ベクトルである．式(5)を満たす部分領域上で応力密度を積分することによって，一意に局所応力が得られる．また，その部分領域は電子密度分布によって決まる．

3

計算結果

上述の手法を用いて， diamond/cBN(110), (7+7)
界面について局所応力を計算した．エネルギー密度および応力密度は PAW 法に基づく第一原理計算ソフトウェア QMAS によって計算した[3]．原子構造緩和後およびセル最適化後に得られた格子定数を diamond 結晶および cBN 結晶のものとあわせて Table 1 に示す．表中 a,b,c はそれぞれ x,y,z 軸に平行な辺の長さである．簡単のため，Fig.1 に示すように周期境界条件下の直方体セルを z 方向に一次元分割することを考える．このとき，xy 方向については常に式(5)を満たすことから， σ 33L のみ考慮すればよい．ここでは x-y 面における電子密度の流束で表される式(6)を分割の条件とした．
∂ρ (r )∫ ∂z dxdy = 0 (6)Fig.2 に電子密度流束の分布を示す．破線は式(6)に基づき設定した分割面である．ただし，PAW 法を用いていることから，核近傍の分割面は除いた．
Fig.2 に示されている通り，原子層ごとに部分領域が設定されることから，局所応力は原子層ごとに計算される．Fig.3 に得られた局所応力を成分ごとに示す．破線は界面の位置を表す．面内応力 σ 11 ，Table 1 Lattice constants (Å) of diamond/cBN(110)(7+7), diamond bulk and cubic BN bulk[3]7+7DiamondcBN

a2.541
3.569
3.621

b3.598BN layers

C layers

c17.831

C layersyxz

Fig.1 Schematic view ofseparation in diamond/cBN (110)1-dimensional

σ 22 は，diamond 層で正となっていることに対して，cBN 層では負となっている．これは Table 1 に示した格子ミスマッチを反映したものである．すなわ
ち，diamond 層には引っ張り応力，cBN 層には圧縮応力が作用している．また，セル最適化の過程で原子の z 方向変位は拘束されていないため，σ 33はほぼ 0 となっている．

4

まとめ第一原理計算による局所応力計算法を示した．
本手法により一意な局所応力を計算することができる．diamond/cBN(110), (7+7)界面について計算例を示し，局所応力から界面の不均一原子構造を反映する結果が得られることを示した．今後は，この分割モデルを 3 次元に拡張し，粒界の強度問題等に適用することでその実用性を検証する．

参考文献[1] R. Martin, Electronic structure, Basic theory and517,
CAMBRIDGE
practicalmethods,
UNIVERSITY PRESS, (2005). (See the referencestherein.)[2] A. Filippetti and V. Fiorentini, “Theory andapplications of the stress density”, Phys. Rev. B, 61,8433, (2000).
[3] S. Ishibashi, T. Tamura, S. Tanaka, M. Kohyama
and K. Terakura, “Ab initio calculations ofelectric-field-induced
stressprofilesfordiamond/c-BN (110) superlattices”, Phys. Rev. B,
76, 153310, (2007); http://www.qmas.jp.

Fig.2 Flux distribution of electron density indiamond/cBN (110)

σ22

σ11
σ23
σ33

Fig.3 Local stress components for each atomiclayer in diamond/cBN (110)