― 90 ―稲上竹蔵・上林一彦

不純物及び不規則性を考慮した c-BN の体積弾性率計算
稲上竹蔵 ・ 上林一彦

Bulk modulus calculation of c-BN with eﬀect impurity and disorder
Takezo INAGAMI andKazuhiko UEBAYASHI（平成25年11月29日受理）
Cubic-Boron Nitride, c-BN, shows speciﬁc properties such as hardness next to diamond and chemical
stability over a wide range of temperatures. Therefore, c-BN has been studied in mechanical and electronicengineering field. We performed first principle calculation and obtained cohesive energy of the mono
crystal c-BN as a function of lattice constant. We also estimated bulk modulus of c-BN by two methods: oneis based on Birch‒Murnaghan equation of state with volume vs. calculated cohesive energy and the otheris Cohenʼs method with a distance between nearest neighbor atoms. As the result of our calculation, weveriﬁed our calculated bulk modulus of pure c-BN is in good agreement with previous experimental andcalculated work. Moreover, by using the same methods we studied the eﬀects of impurity which is group Ⅱ ,Ⅲ , Ⅳ or Ⅴ element in the periodic table and the eﬀects of disorder in mono crystal c-BN for bulk modulus
within super cell. Our calculate results show the bulk moduli of c-BN with impurity and disorder are lower
than that of pure c-BN by 1.5at% and 3at% respectively.
Keywords: cubic-Boron Nitride, ﬁrst principle calculation, bulk modulus, impurity, disorder, Zinc-Blende1.

緒言

中の不純物及び不規則性が体積弾性率に及ぼす影響に着目し，第一原理計算を用いて非経験的に調査を立方晶窒化ホウ素 (cubic-Boron Nitride，以下 c-BN)とは，周期表の上で炭素 (C) に隣接するホウ素 (B) と窒素 (N) から成る化合物である。c-BN はホウ素 (B) と行った。

2.

第一原理計算について

窒素 (N) が sp3 混成軌道を形成しイオン性を含む共有結合で結び付き，ダイヤモンド構造によく似た閃亜

第一原理計算とは基礎的な物理量 ( 原子番号，価電鉛鉱構造をとる。c-BN の代表的な特徴としてはダイ子，結晶構造 ) のみを使い，それ以外には経験的なパヤモンドに次ぐ硬度を有すること，熱的，化学的にラメータを使わずに量子力学に基づいて電子状態をも優れた安定性を示すことが知られている。このた計算する数値計算手法である。我々は第一原理計算にめ研磨材料，研削材料として工業的に用いられ，そおいて密度汎関数理論に基づくプログラムパッケーの合成法等について実験を主導とし研究が行われて

ジ [Vienna Ab-initio Simulation Package( 以 下 VASP)]いる [1]。一方，バンドギャップが大きく，合成時にを利用した。このプログラム群はウルトラソフト擬

種々の元素をドープすることで p 型，n 型半導体とな

ポ テ ン シ ャ ル，Projector Augmented Wave (PAW)る等の電子工学的な特徴も有した材料であり，半導法と平面波基底を用いて電子状態の計算を実行する

体材料の視点から実験と理論の双方について研究が[12]。交換相関項は後述するように局所密度近似 (LDA)行われている [2] [3]。

や，一般化勾配近似 (GGA) の範囲で取り扱った。カッこのように工学的に注目されている c-BN についトオフエネルギーは 300eV，収束条件はエネルギーて，純粋な c-BN を仮定し体積弾性率を計算した先行変化量が 10-4 eV 以下として計算を行った。また，波研究に対し，我々は機械材料の視点から単結晶 c-BN

数空間の分割数については凝集エネルギーへの影響
平成 26 年２月

― 91 ―
不純物及び不規則性を考慮した c-BN の体積弾性率計算

が 10-3eV/ f.u. となるように，k 点は最小単位胞の計
算では 15625(25 × 25 × 25)，Super Cell の計算では 1000(10 × 10 × 10) とした。

3.

体積弾性率の計算手法について材料の機械的性質の一つである体積弾性率を算出するために 2 通りの手法を用いた。体積弾性率の計算手法である Birch の手法と Cohen の手法について紹介する。

図1(a) 最小単位胞 [6] と (b) 第一原理計算による体積に対する凝集エネルギー変化

密度近似を考慮して c-BN の体積に対する凝集エネルギー変化を計算し，それぞれを比較した ( 図 1 (b)) [6]。

3.1. Birch の手法による体積弾性率の計算
第一原理計算の結果を基に体積弾性率を算出する際に用いた手法の 1 つは Birch の関係式によるもの
で，連続体のひずみ弾性論を基にしている [4]。以下に示す式を用いて，体積と凝集エネルギーの関係から体積弾性率 O を算出する。

この図 1 (b) より，GGA(PW91) と GGA(PBE) の体積に対する凝集エネルギー変化は同程度の値，同様の変化傾向を示した。しかし 2 つの GGA(PW91,PBE) と LDAを比較すると凝集エネルギー変化が約 2(eV/f.u.) 異なる結果となった。また最安定状態の格子定数を比較する と LDA が 3.58(Å)， GGA(PW91) と GGA(PBE) が 共に 3.62(Å) となった。実験的研究による c-BN の格子定数は 3.615(Å) であるので [7]，LDA の近似よりも GGA(3.1)の近似を考慮した方がより実験に近い格子定数を得られる。これは GGA による一般的な改善によると考えられる [8]。

ここで， O はそれぞれ最安定状態の凝集エネル

O

ギーとその体積を表し， O′は体積弾性率の圧力微分である。詳細については付録を参照されたい。

次に，これらの結果から Birch，Cohen それぞれの計算手法を用いて体積弾性率を算出し，計算値 [9]，実験値 [7] と比較した。以下にこれらの比較結果と LDA とGGA(PBE) を用いた体積弾性率の計算結果を示す。

3.2. Cohen の手法による体積弾性率の計算
体積弾性率を求めるもう一方の手法は Cohen の関係式によるものである。この手法は自由電子モデルのエネルギーから体積弾性率を求める手法を拡張し物質の共有結合によるエネルギーギャップを考慮したもので，最安定構造の最近接原子間距離 から体積弾性率 を算出する方法である。この手法はダイヤモンド構造または閃亜鉛鉱構造の結晶構造をとり，周期表上でⅣ族，Ⅲ - Ⅴ族，Ⅱ - Ⅵ族に分類される物質にのみ適用可能とされる。このため，ダイヤモンドや GaN といった物質の体積弾性率を良い精度で再現する [5]。今回は次に示す式 (3.2) を既知のものとして使用した。

表1

電子密度近似による c-BN 最小単位胞の格子定数・体積弾性率の計算結果

図2

GGA (PBE) と LDA による体積弾性率と実験値 [7]，計算値 [9] の比較

(3.2)

4.

第一原理計算による計算結果

4.1. 最小単位胞の計算
純粋な c-BN の最小単位胞 ( 図 1 (a)) について第一原
理計算で LDA，GGA(PW91)，GGA(PBE) の 3 つの電子
秋田高専研究紀要第 49 号

― 92 ―
稲上竹蔵・上林一彦

表 1 より，Birch の手法による体積弾性率の計算結付与する際にドープされる元素を主としている [2]。こ果は実験値よりも過大に評価する傾向が見られた。ここで，純粋な (c-BN 32 に対して不純物を仮定した c-BN

れに対し Cohen の手法による体積弾性率の計算結果を (c-BN 32ʼ とする。

は実験値よりも過小に評価する傾向が見られた。また

以下に純粋な (c-BN 32 と不純物元素を仮定した (c-GGA (PW91) と GGA (PBE) による c-BN の最安定状態のBN 32ʼ の体積弾性率の計算結果を示す。 図 4 において格子定数と体積弾性率は同程度の値を示した。

縦軸は体積弾性率，横軸は仮定した不純物元素を示す。

図 2 から LDA を用いた体積弾性率は，実験値に対また黒印 ( ●，▲ ) は N site 置換，白印 ( ○，△ ) は B

して Cohen の手法では良い一致が得られたが，Birchsite 置換を表す。さらに丸印 ( ●，○ ) は Birch の手法の手法では約 12.4% 過大評価している。一方 GGA (PBE)による計算結果，三角印 ( ▲，△ ) は Cohen の手法にを用いた体積弾性率は，実験値に対して Cohen の手法よる計算結果，また，一点鎖線は純粋な c-BN の体積弾

では約 -4.2% の過小評価，Birch の手法では約 2.9% の性率の実験値 [7] を示す。

過大評価となった。第一原理計算による他の先行研究[9] では，実験値に対して Cohen の手法では約 4.4% の過大評価，Birch の手法では約 3.9% の過大評価となっている。このため Cohen，Birch の両手法において体積弾性率が実験値，計算値に比較的良い再現性を示し，先行研究より計算精度が改善され実験値により近づいたGGA (PBE)による近似を以下の計算に適用した。

4.2. 単位胞を Super Cell とした計算
図 1 に示した純粋な c-BN の最小単位胞を基にこれの 32 個分を新たに単位胞とし，Super Cell(c-BN)32 と定義する。この Super Cell(c-BN 32 を基に，不純物元素一個が純粋な (c-BN 32 単位胞内の窒素原子サイト一個と置換する場合 ( 以下 N site 置換，図 3 (a) ) とホウ素原子サイト一個と置換する場合 ( 以下 B site 置換，図3 (b) ) の 2 通りを仮定した。不純物元素一個と c-BNのサイト一個を置換した場合を，不純物を約 1.5at% 含んでいるとする。

(a)N site 置換
図3(b) B site 置換(c-BN)32 を基に不純物を仮定した場合 (a),(b)(a)Impurity X = Be,B,C,N(b)Impurity X = Mg,Al,Si,P
図中の一点鎖線は純粋な c-BN の実験値 [7]図4(c-BN)32ʼ の体積弾性率

本計算では (c-BN 32 単位胞内の最も簡単な場合のみを考慮し，実際の c-BN に考えられる格子間原子や原子空孔が存在する場合は取り扱わない。
不純物として仮定する元素は周期表でⅡ〜Ⅵ族に属する 4Be，5B，6C，7N，12Mg，13Al，14Si，15P とした。これらの元素は c-BN を製造する際に使用する触媒に含まれている元素や，c-BN に p 型，n 型の半導体特性を
平成 26 年２月

― 93 ―
不純物及び不規則性を考慮した c-BN の体積弾性率計算

4.3.

不規則性を含む Super Cell の計算結果

5.

考察

5.1.

不純物を仮定した c-BN の計算結果と純粋なc-BN の実験値との差について
図 4 に示した不純物元素を仮定した場合の c-BN の体積弾性率と純粋な c-BN の実験値 [7] について，両者の差の割合を以下に示す。
表2
図5

Birch の手法による不純物を仮定した c-BN の計算結果と実験値 [7] の体積弾性率の差 (%)(c-BN)32 に不規則性を仮定した場合

表 2 に示した計算結果から c-BN の実験値 [7] に対して，Cohen の手法では純粋な場合の計算結果を含め体積弾性率を小さく評価されている．一方 Birch の手
図6(c-BN)32” の体積弾性率法では実験値に対し純粋な c-BN の体積弾性率の計算結果は 1.9% 過大に評価している．また，この Birchの手法では B site 置換の際に Be,C,Si,P を仮定した場

純粋な (c-BN 32 を基に不規則性を仮定する場合，合，N site 置換の際に C を仮定した場合には体積弾性単位胞内の近接する 1 つの窒素サイトと 1 つのホウ率の低下は見られなかった．特に C を B site 置換した

素サイトを互いに置換し，これを (c-BN 32” とした ( 図場合の体積弾性率の計算結果は純粋な c-BN の実験値5)。図 6 に不規則性を仮定した (c-BN 32” の計算結果に対し 2.3% 大きく評価された．

を示す。この場合，不規則性を約 3at% 含んでいるとする。
図 6 の計算結果において縦軸は体積弾性率，横軸は c-BN の純粋な場合 ( ○，● ) と不規則性を仮定した場合 ( □，■ ) を示す。また白印 ( ○，□ ) は Birchの手法，黒印 ( ●，■ ) は Cohen の手法，一点鎖線は純粋な c-BN の体積弾性率の実験値 [7] を示す。
図 6 から純粋な (c-BN 32 の体積弾性率に比べ不規則性を仮定した (c-BN 32” の体積弾性率が低下することが計算により示された。この結果は Birch の手法，Cohen の手法の両方で見られた。よって計算において純粋な場合に対し不規則性を仮定した場合に体積弾性率が低下する傾向は計算手法に依らないといえる。

秋田高専研究紀要第 49 号

5.2.

純粋な結晶における体積弾性率の実験値と計算結果の差について
c-BN の純粋な場合の体積弾性率について，実験値と本研究の計算結果に差が生じた理由を検討する。
図 4 か ら， 純 粋 な c-BN の 体 積 弾 性 率 の Birch の手 法 に よ る 計 算 結 果 は 実 験 値 [7] に 対 し 約 1.9%(+7.4GPa) 過大評価している。このように計算と実験の体積弾性率の間に差が生じた一因として，理想系と現実系の結晶構造の規則性による違いが考えられる。
本研究の計算では，周期境界条件を用いて理想形として完全に規則化された結晶を用いている。一方で試料を合成する現実系では，一部に不規則性を含むことが考えられる。我々が不規則性を考慮した計算 (4.3 節 )において，3at% の不規則性を含んだ c-BN では現実系

― 94 ―
稲上竹蔵・上林一彦の実験値に対して体積弾性率が約 0.7% (-2.5GPa) 低下する計算結果が得られた。この様に c-BN に不規則性を仮定した場合に計算結果が実験値に近づくことから，現実系での c-BN は不規則性による体積弾性率への影響を受けた可能性が示された。この他に現実系の再現性をより向上させて計算を行うためには，原子空孔や格子間原子をはじめとする格子欠陥等の影響についても計算に考慮する必要がある。

5.3.

Birch の手法と Cohen の手法による体積弾性率の差の改善
前述の Birch の手法による計算結果と実験値の差に対し，Cohen の手法による体積弾性率の計算結果と実験値 [7] との差は約 4.2% となった。また Birch の手法による体積弾性率に対して Cohen の手法による体積弾性率は純粋な場合を比較すると約 6.0% 小さな値を示している。この Birch の手法と Cohen の手法との体積弾性率の差の一因として，式 (3.2) に示した Cohenの手法にイオン性及び共有結合性による効果が明確に含まれていないことが考えられる。イオン性の寄与を経験的に考慮した式は，Cohen [5] によれば(5.1)
で 示 さ れ る。 式 (5.1) は Ⅳ 族 (λ=0)， Ⅲ - Ⅴ 族 ( λ=1)，Ⅱ - Ⅵ族 (λ=2 の半導体にのみ適用できる。また c-BN の原子結合はイオン性を含んだ共有結合となっているため，結合のイオン性 ‐ 共有結合性の度合いを表す係数として b を導入した。 b は 0 から 1までの値をとることとし，完全な共有結合の場合にはb=0 ，完全なイオン結合の場合にはb=1 とする。

この b を用いて式 (5.1) にイオン性と共有結合性の双図中の一点鎖線は純粋な c-BN の実験値 [7]図7Cohen の手法を Birch の手法に近づけた場合の体積弾性率

図 7 より，不純物を仮定した計算においてを用いた Cohen の手法による体積弾性率は Birch の手法による体積弾性率に大まかな部分で近づいた。僅かに異なる部分については各不純物元素自体のイオン性まで考慮することでより Birch の手法に近づくと考えられる。また，第一原理計算による Birch と Cohen の両手法での体積弾性率計算の改善については，本研究と実験の双方を含めたより詳しい検討が必要である。

5.4.

計算で得られた体積弾性率の共有結合半径による考察
ここでは図 4 の (a),(b) に示した，本研究で Birch，Cohen の両手法の計算から得た体積弾性率のみについて，他の文献で得られる元素の共有結合半径を用いて考察を行う。
表3元素の共有結合半径 (Å) [10]方を考慮し，体積弾性率を(5.2)とした。純粋な c-BN の体積弾性率について b を変数として Birch の結果に Cohen の結果を近づけた場合bは 0.44 となった。以下に b=0.44 とし，c-BN に不純物として 4Be，5B，6C，7N を仮定した場合について
Cohen の手法と Birch の手法による体積弾性率を示す。

表 3 から N と B の共有結合半径は本研究では 0.74Åと 0.80Å である。 また図 4 から不純物に N と B より共有結合半径が大きな Be, Mg, Al, Si, P が N,B のいずれのサイトに入った場合でも，本研究での純粋な c-BNの計算結果に対し体積弾性率の低下が見られた。このことから，不純物として置換する元素と母体となる B,Nの共有結合半径の差から生じる結晶格子全体の歪みが本研究で得られた計算結果と一致しており，体積弾性率低下の一因だと考えられる。

平成 26 年２月

― 95 ―
不純物及び不規則性を考慮した c-BN の体積弾性率計算

6.

結言
本研究から以下の 3 点が示された。先ず，純粋な9.

参考文献

[1]

谷口尚 , " 非結晶 BN からの微粒 cBN 焼結体の合c-BN の格子定数と Birch の手法による体積弾性率の成 ," 高圧力の科学と技術 , 13, no.1, pp. 16-23

計算結果は，先行研究による c-BN の実験値 [7] 及び(2003).

計算値 [9] に対し比較的良い再現性が得られた。次[2] 遠藤忠 , "Cubic BN の単結晶 ," 固体物理 , 16, no. 5,pp. 53-58(1981).に，Birch の手法による純粋な c-BN の体積弾性率の計算結果は，不規則性を仮定することで先行研究によ[3]

Rev. Lett., 110, no. 065504, pp. 1-5(2013.)の不純物による体積弾性率への影響について，Birch[4]

F. Birch, "Finite Elastic Strain of Cubic Crystals,"
Phys. Rev., 71, pp. 809-824(1947).

示され，その影響は両手法で同様の傾向が得られた。
[5]

7.

"Functional Complex Point- DefectStructure in a Huge-Size-Mismatch System," Phys.

る実験値に近づくことが示された。最後に，c-BN 内の手法と Cohen の手法を用いて数値計算の視点から

R. Ishikawa,今後の展望

M. L. Cohen, "Calculation of bulk moduli of
diamond and zinc-blende solids," Phys. Rev. B, 32,
no.12, pp. 7988-7990(1985).

本研究では今回，物質の格子定数に対する凝集エ[6]Crystallogr. Comput,” IUCr Newslett.No.7,

ネルギー変化を基にした体積弾性率の計算方法のみ

pp. 106-119(2006).

を用いた。体積弾性率のより詳細な議論を行うためには，電子構造計算による電子状態の詳細な調査や，

K. Momma and F. Izumi, “Commission on[7]

J.S.Zhang,"Elasticity of cubic boron nitride

under ambient conditions," J. Appl. Phys., 109,

今回の計算で仮定した個々の不純物を含む c-BN の体

no. 063521(2011).

積弾性率の実験的測定が望まれる。また今回の計算
では周期性の良い c-BN の単結晶を扱い，不純物につ[8] ナノシミュレーション技術ハンドブック委員会 ,いてもサイト置換のみを仮定した。実際に利用されナノシミュレーション技術ハンドブック , 共立出

る c-BN は焼結体 ( 多結晶体 ) が多いことや，格子間版 (2006).

原子をはじめとする格子欠陥等が存在する可能性も[9]

A. Mahmood,"BULK MODULUS

考慮しなければならない。今回用いた計算手法とは

CALCULAIONS FOR GROUP- Ⅳ CARBIDES AND

異なる方法も導入し，これらの要因を計算に反映さ

GROUP- Ⅲ NITRIDES," Mod. Phys. Lett. B, 18, no.

せることで現実系に近づくと考えられる。以上の改24, p. 1251(2004).

善点を考慮しながら研究を進める事で今回扱った周[10]体や複合材料についても非経験的な取り扱いを可能

井口洋夫 , 元素と周期律 , 裳華房 ,改訂第 16 版 (1978).

期的な単結晶の計算を発展させ，焼結体等の多結晶[11]F. D. Stacey and P. M. Davis, Physics of theEarth Fourth Edition, Cambridge Universityにすることが今後の課題である。

Press (2008).

8.

謝辞

[12]

J. Hafner,“Ab-initio simulations of materials
using VASP,” Journal of computinal chemstry,

本研究を進めるにあたり，熱心な御指導ならびに御教授賜りました上林一彦講師，多くの御意見ならびに御指導賜りました成田章教授，磯部浩一教授，若生昌光教授，金田保則教授に深く感謝致します。
また本研究を進める上で多くの意見を下さりました
上林研究室の皆様，本当にありがとうございました。
本研究は秋田高専学科横断プロジェクトの助成を受けたものです。

秋田高専研究紀要第 49 号

29, 2044(2008).

― 96 ―
稲上竹蔵・上林一彦

10.

付録 :Birch-Murnaghan の状態方程式(10.6)

ここでは連続体のひずみについて考える。連続体の 2 点 に 着 目 し， 変 形 前 の 点 を =( 1, 2, 3)， 変 形
後 の 点 を =( 1, 2, 3) と お く。 こ の 2 点 間 の 変 位 を=( 1, 2, 3) とおき，変形係数をαとする。ここで=α とすると， は以上より，ひずみεは変形係数及びひずみテンソルを用いて(10.7)(10.1)の関係で表せる(図 8)と表すことができる。
一方，連続体の質量を，変位前の基準状態の体積，密度をρ0，変位後の体積を ，密度をρとする。

を0

密度変化の比を (10.7) を用いて表すと(10.8)

図 8 連続体の変形

また，微小領域においてで表し， の微小量の微小量を全微分の形の二乗をとり，各項を整理す

ここで，ひずみεによる系の内部エネルギー変化

るとをひずみεで多項式展開すると(10.2)となる。この(10.9)を二次の微小量まで考慮したオイラー座標のひずみテンソルという。となる。ここではひずみεによる系の内部エネルギー変化(10.3)を考えているので，ε =0 のとき=0。このため 0=0 となる。
また基底状態を仮定し，熱力学第一法則を用いて圧力

ここでは連続体の等方的な変位のみを考え（ ）→（を表し，ひずみεで展開すると）とし，(10.2) の左辺を展開し整理すると(10.4)(10.10)となる。また，両辺の成分について解くと(10.5)

ひずみ ε がゼロのとき，系に圧力は働いていないと考えられるので =0 。よって 1=0 となる。
=0 付近のとき圧力 = 0 ，体積弾性率 = 0 となとなる。これをオイラー座標でのひずみεとおく。

り変形前の基準状態での体積弾性率

また (10.3) に (10.1) を代入すると定義より

0は，一般的な平成 26 年２月

― 97 ―
不純物及び不規則性を考慮した c-BN の体積弾性率計算(10.11)(10.11) と同様に変形前の基準状態で0の圧力微分をとって極限操作を行い， 0′とおく。 = 0 より(10.12)

以上より，ある体積における圧力は (10.10) を展開し，(10.13)と示される。(10.13) は Birch-Murnaghan の状態方程式と呼ばれる [11] 。第 3.1 節で示した (3.1) は，(10.13) を体積について積分することで表せる。

秋田高専研究紀要第 49 号