Electron Theory Calculation of Thermodynamic Properties of Steels and Its Application to Theoretical Phase Diagram of the Fe-Mo-B Ternary System

Abstract

In this study, the ground structures of the Fe-Mo-B ternary systems were estimated by first-principles calculations based on genetic algorithm and the free energies of their structures are evaluated by electronic calculations and statistical thermodynamic techniques. In addition, the phase diagram at finite temperature was theoretically constructed using the calculated free energies and the result was compared with the experimental knowledge. The space groups and compositions of many ground structures obtained by the calculations correspond well with the experimental findings, but the agreement is not perfect. However, by including metastable structures by only a few kJ/mol than the ground state, it becomes clear that the appearance of almost all structures can be predicted based on this technique. The new calculation technique of such theoretical phase diagrams suggested in the present study is expected to open up the possibility of estimation of unknown phase diagram, reexamination of experimental phase diagrams and discovery of new phases. On the other hand, examining the calculation conditions for improving the accuracy of energy calculation, consideration of the anharmonicity of atomic vibration, magnetic entropy effect, handling of solid solution, etc. are mentioned as problems requiring further consideration.

1. 緒言

様々な材料の機能を最大限に引き出すための組織制御において，状態図が果たしてきた役割は非常に大きい。19世紀半ばに金属の組織が光学顕微鏡によって調べられ始めて以来，合金やセラミックスなどの状態図は，顕微鏡観察や熱分析などの実験的手法を駆使して研究されてきた。しかし，第2次大戦後の産業の復興とともに身近な金属材料も多成分化したことから，このような時間のかかる実験状態図では労力的にも迅速性の面でも現実の材料設計に対応できなくなった。そこでこのような問題に対処するために新たに提唱された状態図の作成方法がCALPHAD(Calculation of Phase Diagrams)法である¹⁾。CALPHAD法では，適切な熱力学モデルを用いて相境界や熱力学物性値に関する実測値を解析し，対象とする系の熱力学的性質をパラメータ化する。また熱力学的外挿法によって多元系の物性も合理的に予測できるので，実際の相平衡を正確に予測できる特長がある。その利点を生かして材料開発の場において大きく発展してきた状態図計算法である。しかし実験で明らかにされていない系や，準安定領域，非平衡相に対して，この手法を適用することはきわめて難しい。そこで本論文では，Fe-Mo-B三元系を例に第一原理に立脚した電子論計算によってこの合金系の熱力学的性質を評価した。さらにその結果を用いて，有限温度における合金状態図を理論的に構成する新しい手法を提案し，現状の問題点を整理することを研究目的とした。

2. 計算方法

本章では，はじめに本研究で用いた第一原理計算の条件を述べ，絶対零度における完全規則構造の生成エネルギーの計算原理を説明する。さらに，この結果を有限温度の熱力学物性値として拡張する手法として，原子の格子振動計算と準調和近似に基づいて自由エネルギーを計算する方法とクラスター展開・変分法により固溶体の自由エネルギーを計算する手法を説明する。さらに遺伝的アルゴリズムによって基底状態の相平衡を探索する手法について説明する。

2・1　第一原理計算法

本研究では第一原理計算は密度汎関数法による計算コードVienna Ab initio Simulation Package(VASP)^2–4)を使用し，交換相関汎関数にはGGA-PBE^5,6)を用いた。平面波のカットオフエネルギーは400 eVとした。2・2・2，2・2・3節で述べるクラスター展開法および遺伝的アルゴリズム法では規則構造のエネルギー計算においてk点サンプリングは primitive BCC-Fe構造に対し17×17×17のグリッド点を用いた。一方，準調和近似法の計算に用いたスーパーセルにおけるk点サンプリング数は2・2・1で述べる。全ての計算においてスピン分極を含めた計算を行い，磁気モーメントの初期値はFeを+2.5 μ_Bとし，MoとBを0 μ_Bとした。

第一原理計算から得られた全エネルギーから絶対零度における化合物の生成エンタルピーを評価することができる。例として，二元化合物A_nB_mの生成エンタルピー(H_AnBm)は式(1)で表される。

  

$$H_{{\text{A}}_{\text{n}}{\text{B}}_{\text{m}}}=E_{{\text{A}}_{\text{n}}{\text{B}}_{\text{m}}}-\left(\text{n}{E}_{\text{A}}+\text{m}{E}_{\text{B}}\right)$$(1)

ここでE_AnBmは第一原理計算から求めた二元化合物A_nB_mの全エネルギー，E_A,E_Bは基準状態とした元素A，Bにおける単体の全エネルギーである。

2・2　有限温度における熱力学物性値の計算

第一原理計算によって評価した物性値は絶対零度における値であるため，有限温度の相平衡を計算するためには自由エネルギーを評価する必要がある。有限温度で大きな寄与を持つ成分として，原子振動によるものと固溶による原子配置のランダムネスにより生じるものの二種類が挙げられる。前者については，第一原理計算を用いてフォノンの状態密度を計算し，さらに体積膨張の影響を取り入れることで，温度による生成エネルギーの変化を考慮することができる。後者については，多数の規則構造について計算された凝集エネルギーをクラスターで展開し，クラスター変分法で有限温度の固溶体の自由エネルギーを決定する方法がある。本節ではこれらの手法の概要を述べる。

2・2・1　格子振動を考慮した化合物の自由エネルギー計算

単体や規則化した化合物の有限温度における自由エネルギーや比熱は，フォノンからの寄与が大きい。第一原理計算から原子間に働く力を求めることで，このフォノンの寄与を計算することが可能である⁷⁾。具体的な計算として，結晶格子を拡張したスーパーセルを用いて，原子位置を平衡位置から微小変位させることで生じる復元力を評価する。これらの復元力から原子変位に対する二次の力定数を計算し，フォノンの振動数(ω)の波数依存性を求める。この手法を調和振動近似と呼ぶ。これによって得られるフォノンの状態密度から，有限温度における格子振動の効果を含むヘルムホルツの自由エネルギーが式(2)のように与えられる。

  

$$F\left(T\right)=U_0+\frac{\hslash }{2}{\displaystyle \sum _{q_{\text{j}}}\omega \left(q_{\text{j}}\right)+k_{B}T}{\displaystyle \sum }_{q_{\text{j}}}\text{ln}\left[1-\text{exp}\left(\frac{-\hslash \omega \left(q_{\text{j}}\right)}{k_{B}T}\right)\right]$$(2)

ここで，U₀は原子位置平衡状態におけるエネルギー，またq_jはj番目のフォノンモードにおける波数ベクトルを表す。 第二項はゼロ点エネルギーと呼ばれるもので，不確定性原理のために絶対零度においても原子振動が存在する効果に由来する。ギブスの自由エネルギーは体積Vにおけるヘルムホルツの自由エネルギー(F_V)から体積-エネルギー曲線を計算し，圧力一定の条件でエネルギー最小となる体積Vを求めることで計算される。

  

$$G\left(T\right)=\underset{V}{\min }\left[F_V\left(T\right)+PV\right]\approx \underset{V}{\min }\left[F_V\left(T\right)\right]$$(3)

準調和近似計算にはPhonopyコード^8,9)を用いた。3・2節ではBCC-Feをはじめとする19の規則構造の自由エネルギーの計算結果を示すが，それらの規則構造の第一原理計算に用いたスーパーセルのサイズおよびk点のメッシュサイズをTable 1にまとめた。各構造についてセルの形状，体積を含む全ての原子配置を一旦緩和させた後に，等方的に体積を変えたセルを作成することで準調和近似計算を行った。

Table 1. Supercell size and mesh size of k-points of structural models used in quasi-harmonic approximation calculations.

Formula	Space group	Size of super cell	Total num. atoms	k-points
BCC-Fe	Im3m	5×5×5	125	2×2×2
BCC-Mo	Im3m	5×5×5	125	2×2×2
β-B	R3m	2×2×2	96	2×2×2
Fe2B	I4/mcm	2×2×3	72	4×4×2
FeB	I41/amd	2×2×2	64	4×4×2
FeB	Pnma	2×2×2	64	4×4×2
FeB2	Pnma	2×2×2	96	4×2×2
Mo2B	I4/mcm	2×2×3	72	4×4×2
Mo3B2	P4/mbm	2×2×2	72	2×2×4
MoB	I41/amd	2×2×2	64	4×4×2
MoB	Cmcm	3×3×3	108	2×2×2
MoB2	R3m	2×2×1	72	4×4×2
Mo2B3	Imm2	2×2×2	80	4×2×2
MoB4	P63/mmc	2×2×2	80	4×4×2
Fe2Mo	P63/mmc	2×2×2	96	2×2×2
Fe7Mo6	R3m	2×2×1	52	2×2×2
Fe3B	Pnma	2×1×2	64	2×4×2
FeMo2B2	P4/mbm	2×2×2	80	2×2×4
Fe2MoB4	Immm	2×2×1	56	4×4×2

2・2・2　クラスター展開-変分法(CE-CVM)による自由エネルギー計算

多成分系の有限温度における相平衡を計算するには，フォノンの寄与に加えて，混合自由エネルギーを考慮する必要がある。クラスター展開・変分法(CE-CVM)は，個々のクラスターの有効クラスター相互作用(ECI)と，その配置のエントロピーから混合自由エネルギーを計算する手法である¹⁰⁾。

具体的には，スピン演算子σを導入し，元素種ごとに異なる値のスピン演算子をあてはめることで，各規則構造Rをスピン配列R=(σ₁，σ₂，σ₃，... σ_i，...)で表記する。ここで は規則構造の原子サイト位置を表す指数である。また，異なるサイトから抽出したスピン演算子の組み合わせによってクラスターを表現することが可能であり，例えば，(σ_i，σ_j)はサイトi，jからなる対クラスター，(σ_i，σ_j，σ_k)はサイトi，j，kからなる三体のクラスターである。さらに点，対，三体などのクラスターについてスピン積の平均値をとることによって式(4)で表される相関関数φ_αを導入する。

  

$$\begin{array}{l}{\phi }_{\text{point}}=\frac{1}{{\text{N}}_{\text{point}}}{\displaystyle \sum }_{\text{i}}{\sigma }_{\text{i}}\\ {\phi }_{\text{pair}}=\frac{1}{{\text{N}}_{\text{pair}}}{\displaystyle \sum }_{\text{i},\text{j}}{\sigma }_{\text{i}}\cdot {\sigma }_{\text{j}}\\ {\phi }_{\text{tri}}=\frac{1}{{\text{N}}_{\text{tri}}}{\displaystyle \sum }_{\text{i},\text{j},\text{k}}{\sigma }_{\text{i}}\cdot {\sigma }_{\text{j}}\cdot {\sigma }_{\text{k}}\end{array}$$(4)

αはクラスターの種類を表す添字であり，N_point，N_pair，N_triは規則構造中に含まれる点，対，三体のクラスターの総数である。この相関関数の値は規則構造に含まれるクラスターの濃度と一対一の対応関係を持つことから，規則構造のエネルギーE_Rを相関関数とECI(e_null，e_point，e_pair，e_tri,…)との積の総和で表すことが可能である。

  

$$E_R=e_{\text{null}}+e_{\text{point}}{\phi }_{\text{point}}+e_{\text{pair}}{\phi }_{\text{pair}}+e_{\text{tri}}{\phi }_{\text{tri}}+\mathrm{...}$$(5)

式(5)は理論上では無限の種類のクラスターを用いることで厳密に規則構造のエネルギーを再現することが可能であるが，現実の計算においては使用するクラスターを有限の数に打ち切る必要がある。そこでクラスターの相互作用が短距離のものほど強く，長距離になると弱くなる傾向にあるため，クラスターの結合距離に閾値を設けて，それ以下のサイズのクラスターのみを選択する手法が一般的に採用される。そこで定められる最大サイズのクラスターをα_maxとして，これに内包されるクラスター群(サブクラスター)についての相関関数とECIの総和をとる。

  

$$E_R={\displaystyle \sum }_{\alpha }^{{\alpha }_{max}}{e}_{\alpha }{\phi }_{\alpha }$$(6)

相関関数φ_αは規則構造から求まり，左辺の全エネルギーは第一原理計算から得られるために，未知数はECIのe_αのみであるが，多数の異なる規則構造とエネルギーの関係式を用意することで，最小二乗法によりECIを決定することができる。この手法においては，使用するサブクラスターの組み合わせの選択に任意性があるため，交差相関誤差検定を行い予測誤差が最小となるよう基底クラスターの組み合わせを決定する方法が採用される¹¹⁾。

一旦ECIが決定されれば，任意の原子配列のエネルギーは第一原理計算を行わずに，クラスター展開の精度内で決定することができる。また，クラスターの配置の場合の数を考慮することで，配置のエントロピーを含む自由エネルギーを以下のように表すことができる。

  

$$F\left(T\right)={\displaystyle \sum }_{\alpha }^{{\alpha }_{max}}{e}_{\alpha }{\phi }_{\alpha }-T{\displaystyle \sum }_{\alpha }^{{\alpha }_{max}}{\gamma }_{\alpha }{S}_{\alpha }$$(7)

式(7)における第二項のγ_αS_αはKikuchi-Barker係数を用いたクラスターαからの配置エントロピーの寄与である¹²⁾。式(7)に変分法を適用しエネルギーが最小となる相関関数を求めることで自由エネルギーを計算する。

CE-CVM計算にはiCVMコード¹¹⁾を用いた。2・3・4におけるBCC Fe-Moの計算では，配置の異なる149個の規則構造のエネルギーを解析に用いた。またそれぞれの規則構造はセルの形状，体積を含む全ての原子配置を一旦緩和させた後にエネルギーを評価した。第5近接距離までに含まれる4体のクラスターを最大クラスターとし，それに含まれる35種類のサブクラスターにECIを割り当てた。交差相関誤差検定による予測誤差は13 meV/atomとなった。

2・2・3　基底構造探索

任意の系の状態図を作成するにあたり，その状態図を構成する規則相を同定することは最も基本的で重要な情報と言える。しかしながら，第一原理計算では物質の組成が与えられたとしても，その組成における安定な規則構造を予測することが困難である。この問題は主に2つの要因に起因する。一つは構造を記述する変数の次元の数が3N+3(N:単位胞に含まれる原子数)と，大きな自由度を取り扱う必要があり，すべての構造のエネルギーを網羅的に計算することが不可能であるためである。もう一つは，原子核にかかる力を第一原理的に計算することで構造緩和を行う場合においても，局所的安定状態が無数に存在するためにエネルギー障壁に阻まれて最安定状態に到達できないことにある。

近年，この問題に対して遺伝的アルゴリズム(GA)を用いることで，第一原理計算からの最安定構造の予測を行う計算手法が考案されている^13–15)。遺伝的アルゴリズムはダーウィンの進化論を模倣し，自然淘汰と遺伝と突然変異を採用した計算手法であり，具体的な計算は以下の(1)から(4)のプロセスを経て行われる。

(1)ランダムに作成した複数の規則構造のエネルギー計算を行い，それらの中からエネルギーの低いものを優先的に選択する。

(2)それらの構造の局所的な原子配置に関するパラメータを遺伝子として見立て，それらを交配・突然変異させた次世代の構造群を作成する。

(3)新たな世代のエネルギー計算を行い，エネルギーの低いものを再選択する。

(4)(2)および(3)の過程を繰り返し実行し，低エネルギー構造群を更新していくことで最安定構造を探索する。

このように，交配や突然変異を模倣的に取り入れることでエネルギー障壁の問題を越えて，安定構造の探索を効率的に行うことができる。さらに，このようにして得られた構造群の組成-エネルギー点から構成される最低エネルギーの凸包体(convex-hull)を計算する。凸包体の頂点の構造が絶対零度における基底構造であり，第一原理計算をもとにした状態図が作成される。

GAによるFe-Mo-B三元系の基底構造探索にはUSPEX^13,14)を利用した。具体的な計算条件として，各世代における構造の数を400個とし，次世代に引き継ぐ安定構造とその交配により得た構造を合計280個用意し，さらに突然変異により作成した構造を80個，新たに作成したランダム構造を40個を加えることで，新たな世代を作成した。10世代までの世代交代を繰り返し，安定・準安定構造の選別を行った。またconvex-hullの計算には，Qhullコード¹⁶⁾を利用し，得られたすべての構造群についてFeおよびMoの組成比(x_Fe，x_Mo)，およびエネルギー(E)からなる三次元座標からconvex-hullを求めた。

3. 計算結果と考察

3・1　Fe-Mo-B三元系の基底状態図

Fe-Mo-B三元系の基底探索の結果をFig.1に示した。図中の黒丸は本研究で得られた構造の組成であり，四角のシンボルは実験的に確認されている化合物の組成である。また実線は，GAによって探索された安定化合物から構成される二相の相境界であり，これらから構成される三角形の領域は，この三元系では不変系反応である。

Fig. 1.

Comparison of the ground structure searched by GA and the structure confirmed by experiments.

3・1・1　実験状態図において観察される化合物

これらの化合物の結晶構造をTable 2にまとめた。まず実験から確認されているFe-Mo-B三元系を構成する化合物の概要を述べる。Fe-Mo-B三元系では16種類の二元系，三元系化合物が確認されている。それらの化学式をTable 2に示した。この表では文献¹⁷⁾で使用されている構造名の表記を採用した。また表記が組成式とは異なる構造については，その構造名をカッコ書きで併記した。但し，下記の理由により文献¹⁷⁾ではMoB₂(MoB_2-x)と表記されている化合物をMo₂B₃(MoB_2-x)に改めた。すなわち，この相はMo:B=2:3の組成近傍で現れ，その組成幅が小さいことが状態図に示されている。さらにGAの構造探索からは，Mo₂B₃の組成においてこの相に対応すると考えられる規則構造が確認された。この点については3・1・2節で改めて議論する。

Table 2. Comparison of structures confirmed by experiments and those obtained by GA.

Formula	Experimental data	GA
	Space group	Space group	ΔEfrom hull (kJ/mol)
Fe2B	I4/mcm	I4/mcm	0
FeB	Pnma	I41/amd	0
		Pnma	0.46
Mo2B	I4/mcm	I4/mcm	1.83
Mo3B2	P4/mbm	P4/mbm	2.93
MoB	Cmcm	I41/amd	0
		Cmcm	1.29
MoB2 (Mo2B4)	R3m	R3m	1.83
Mo2B3 (MoB2-x)	P6/mmm	Imm2	1.50
MoB4	P63/mmc	P63/mmc	0.92
Fe2Mo (Laves C14)	P63/mmc	P63/mmc	0
Fe7Mo6 (μ)	R3m	R3m	0
Fe3Mo2 (R)	R3	not confirmed	-
FeMo (σ)	P42/mnm	not confirmed	-
Fe3B (τ1)	Pnma	Pnma	1.57
FeMo2B2 (τ2)	P4/mbm	P4/mbm	0
FeMo8B11 (τ3)	Cmcm	not confirmed	-
Fe2MoB4 (τ4)	Immm	Immm	0
FeB2	not confirmed	Pnma	0
Fe13Mo	not confirmed	P1	0

Fe-B二元系およびMo-B二元系における化合物は，いずれも固溶幅が小さい。一方，Fe-Mo二元系では，Fe₂Mo(Laves C14)，Fe₇Mo₆(μ)，Fe₃Mo₂(R)，FeMo(σ)の4種類の化合物が確認されており，この内Fe₂Mo(Laves C14)は化学量論化合物であるが，Fe₇Mo₆(μ)はMo濃度39-44%，Fe₃Mo₂(R)はMo濃度34-39%，FeMo(σ)はMo濃度43-57%の組成域に現れる非化学量論化合物である。また，三元系化合物であるFe₃B(τ₁)，FeMo₂B₂(τ₂)，FeMo₈B₁₁(τ₃)，Fe₂MoB₄(τ₄)はいずれもFeとMoが置換することによって固溶領域を示す。

3・1・2　基底構造探索によって得られた化合物との比較

GAによる基底構造探索によって得られた二元系，三元系化合物の空間群を判定し，Table 2において実験による安定構造と比較した。この表には，GAでは準安定と判断された構造もいくつか含められている。各構造については，convex-hull からのエネルギー差(ΔE_{from hull})の計算値も記載した。ΔE_{from hull}が0である構造は基底構造を，正の値を持つ構造は準安定構造であることを意味する。またFeB₂(Pnma)とFe₁₃Mo(P1)は，実験的には見出されていない構造である。Table 2中に太字で示した化合物は，GAの基底構造が実験で確認されている化合物と一致した。しかしそれら以外の化合物は，基底構造だけを比較した場合には実験的知見とは対応しない。そこで，この相違について熱力学的見地から考察した。

(1)準安定構造を含めた基底構造が実験結果と一致した化合物

はじめに，FeBとMoBでは，GAではどちらもI4₁/amdが基底構造であるが，これらよりわずかに準安定となる化合物は実験による構造と一致する。これらの準安定構造のconvex-hullからのエネルギー差は，どちらも1 kJ/mol前後であり，有限温度において構造相変態が起こる可能性がある。この点に関する議論は次節の3・2・2と3・3・3で行う。

Mo₂B，Mo₃B₂，MoB₂(Mo₂B₄)，MoB₄，Fe₃B(τ₁)については，基底構造ではないものの，それらよりもわずかにエネルギーの高い構造の中に実験的知見と一致する構造が確認された。ただし，Mo₂B₃(MoB_2-x)の空間群Imm2は実験的知見のP6/mmmと異なっているが，Fig.2に示すように，Imm2はP6/mmmにおけるBと空孔が3：1の比で規則化した構造に対応する。また，Imm2の構造中のB原子をすべて排除したMo格子の空間群はP6/mmmとなる。従って，X線回折実験では原子番号の小さいBの回折強度が弱いため，Imm2であってもP6/mmmと区別できない可能性が指摘される。さらにMoB₂の組成を有するP6/mmmの規則構造では，convex-hullからのエネルギー差が+14.3 kJ/molであるのに対し，Imm2のそれは+1.5 kJ/molであった。これらの考察から，本研究ではMo₂B₃(MoB_2-x)の構造はImm2であると結論づけた。Mo₂B，Mo₃B₂，MoB₂(Mo₂B₄)，Mo₂B₃(MoB_2-x)，MoB₄，Fe₃B(τ₁)の六種類の化合物はいずれもconvex-hull からのエネルギー差が3 kJ/mol以下であり，有限温度の平衡相として出現する可能性がある。この点に関する議論は3・3節において行う。

Fig. 2.

Comparing (a) MoB₂ P6/mmm structure with (b) Mo₂B₃ Imm2. The dotted circles in figure (b) represent the vacancy positions in P6/mmm structure.

(2)基底構造が実験結果と一致しない化合物

Fe₃Mo₂(R)，FeMo(σ)，FeMo₈B₁₁(τ₃)については，GAによって探索された構造群の中に対応する構造が存在しなかった。この中でRは159原子¹⁸⁾，σは30原子¹⁹⁾から構成される大きな単位胞の不規則構造である。本研究では単位胞の全原子数を32個に制限して構造探索を行ったために，Rが探索されなかったと考えられる。またσについても単位胞が30個の原子によって構成される複雑な構造であるため，GAでこの構造を探索するためには，世代数をさらに大きくする必要があると考えられる。一方，FeMo₈B₁₁(τ₃)の結晶構造は不明で，MoB(R3m)を終端構造とする固溶体であるという報告もある^20,21)。しかし，文献¹⁷⁾における第一原理を用いた熱力学解析では，R3m構造の(Mo，Fe)B型固溶体は熱力学的に不安定であることが確認されている。以上の理由から，本研究ではFe₃Mo₂(R)，FeMo(σ)，FeMo₈B₁₁(τ₃)の構造に関する議論は行わず，次節における自由エネルギー評価の対象から除外した。

(3)実験によって確認されていない化合物

Table 1中の末尾に示したFeB₂，Fe₁₃Moは，対応する規則構造が実験によって確認されていない基底構造である。その内，FeB₂については，VoroshninらによってP6/mmmを持つ構造が実験的に報告されているが²²⁾，その後の実験ではこの相の生成は再現されず，安定性は低いと結論づけられている²³⁾。一方，GAからはVoroshninらの報告とは異なる空間群Pmnaの構造が安定であることを示す結果が得られた。またPmna型のFeB₂は，第一原理計算によるFe-B二元系の基底状態探索においても安定であることが報告されており²⁴⁾，本研究の結果はこれと整合している。有限温度におけるこの構造の熱力学的安定性はこれまで議論されていないので，有限温度の自由エネルギーを計算し，実験事実と比較しながら3・2・2および3・3節で考察を行う。一方，Fe₁₃MoはFig.3に示したようにBCC-Feの規則構造中のFe原子をMoで置換した構造に対応する。Fig.3ではFe₁₃Moを2×2×1に拡張したスーパーセルを示しており，構造内に示したBCC格子で全ての格子点が構成されている。そこでこの構造は3・2・4のBCC固溶体に相当するものとして取り扱った。

Fig. 3.

Crystal structures of Fe₁₃Mo. This diagram shows a supercell with Fe₁₃Mo extended to 2×2×1.

3・2　各相の自由エネルギーの計算結果

本節では，有限温度の相平衡を検討するために，3・1節で得られた規則構造群に対して準調和近似計算による自由エネルギー計算を行った。また固溶体に対しては，2・2・2節のCE-CVM法により一次固溶体の自由エネルギーを評価した。本研究においてはFe，MoへのBの溶解度が小さいことを考慮して，Fe-Mo二元系の固溶体の自由エネルギーのみを計算した。

3・2・1　純物質の自由エネルギー

BCC-Fe，BCC-Mo，β-Bについて，準調和近似計算による自由エネルギー計算を行った。その計算結果をFig.4に示す。横軸に温度，縦軸に単位原子あたりのギブス自由エネルギーを示した。温度の上昇に伴って自由エネルギーは単調に低下するが，温度増加に対するエネルギーの低下量を比較すると BCC-Fe，BCC-Moと比較してβ-Bが最も低下量が小さい。これはB原子の質量がFe，Moに比べて小さいことから格子振動が速く，低エネルギーのフォノンの状態密度が小さいために，振動エントロピーがBCC-Feや BCC-Mo に比べて相対的に小さくなることが原因であると考えられる。また，この図にはDinsdaleによる熱力学データベース²⁵⁾の自由エネルギーの温度依存性を白抜きの丸で示した。熱力学データベースの自由エネルギーについては，T=273 Kにおける第一原理計算の自由エネルギーの値を基準として，それ以上の温度領域における自由エネルギーの温度変化の値をプロットした。いずれの純元素の計算結果もT=1600 Kまでの温度範囲ではデータベースの温度依存性をよく再現する。このことから，本研究で考慮していない格子振動の非調和項の影響はこの温度範囲では十分小さいと考えられる。ただし，BCC-Feに関しては，およそT=1000 Kを超える温度領域からデータベースと計算値がわずかに乖離する傾向が観察される。これはBCC-FeはT=1043 Kにおいて強磁性-常磁性の磁気転移を生じるが，データベースの値にはその磁気変態のエントロピー効果が含まれているのに対し，本硏究結果はこの効果を含まないため実験値との乖離があるものと考えられる。しかしデータベースとの差は小さいことから，本研究ではこの効果を取り入れずにFe-Mo-B三元系の相平衡を計算した。また，FCC-Feは1185~1667 Kの温度範囲で安定相となるが，この相についても常磁性状態における磁気エントロピーを取り扱う必要があるため，本研究ではこの相を取り扱わなかった。磁気エントロピーを考慮するためには，フォノン計算に加えてスピン偏極の配置パターンを変えた様々な計算が必要であり²⁶⁾，大量の計算資源を必要とすることから今後の検討課題とする。

Fig. 4.

Temperature dependence of free energy of BCC-Fe, BCC-Mo, and β-B.

3・2・2　Fe-B 二元系化合物の自由エネルギー

GAの基底として得られたFeB₂，FeB(I4₁/amd)，Fe₂B，および絶対零度では準安定化合物であったFeB(Pnma)の四つの化合物の自由エネルギーを評価した。純物質の自由エネルギーを基準としたそれぞれの化合物の自由エネルギーをFig.5に示した。

Fig. 5.

Temperature dependence of free energy of the Fe-B binary compounds.

実験状態図では確認されていないFeB₂の生成自由エネルギーは，温度増加に伴い僅かに増加する。一方，これ以外の化合物においては，いずれも温度の増加とともに熱力学的安定性が増加する。GAの解析によれば，FeBについては，I4₁/amdおよびPnmaの二種類の結構構造が得られた。Fig.5に示した各構造の自由エネルギーの計算結果は，およそT=400 Kでこの二つの構造の間で相変態が起こることを示している。すなわち，T=400 K以下ではI4₁/amdが安定構造であるのに対し，T=400 K以上では実験で報告されているPnma構造が安定化する。この結果は実験による知見と整合する。

3・2・3　Mo-B二元系化合物の自由エネルギー

Mo-B二元系化合物の自由エネルギーの温度変化をFig.6に示した。GAの基底として得られたMo₂B，MoB(I4₁/amd)，MoB₂，および絶対零度では準安定構造であったMoB₄，Mo₂B₃，MoB(Cmcm)，Mo₃B₂の七種類の化合物の自由エネルギーを評価した。この図から，いずれの化合物の生成自由エネルギーも温度によって大きな変化を示さず，高温までそれぞれの熱力学的安定性を保持していることが確認される。

Fig. 6.

Temperature dependence of free energy of the Mo-B binary compounds.

実験によるMoBの安定構造はCmcmであると言われるが，本研究で基底状態として得られたI4₁/amdは計算温度範囲でCmcmよりもわずかに安定性を保っていた。しかしその差は1 kJ/mol程度であり，きわめて小さい値である。従って，計算条件や用いるポテンシャル，固溶によるエントロピー効果や格子振動の非調和性などによって容易に安定性が入れ替わる可能性があると考えられる。

3・2・4　Fe-Mo二元系化合物の自由エネルギー

Fe-Mo二元系ではFe₂Mo，Fe₇Mo₆，Fe₃Mo₂(R)，FeMo(σ)の四種類の化合物が実験から確認されているが，3・1節で述べたようにRおよびσについてはGAで対応する構造が得られなかったため，Fe₂Mo，Fe₇Mo₆のみの自由エネルギーを計算した。また GAによって基底構造として新たに得られた Fe₁₃Moは，BCC-FeにMoが固溶した固溶体とみなせるため，この構造に対しては準調和近似計算は実施しなかった。Fig.7にFe₂MoおよびFe₇Mo₆の生成自由エネルギーの温度依存性を示す。二つの化合物の生成自由エネルギーはいずれも-3 kJ/molから-6 kJ/molであり，Fe-B二元系化合物およびMo-B二元系化合物の生成自由エネルギーに比べて絶対値が小さい。

Fig. 7.

Temperature dependence of free energy of the Fe-Mo binary compounds.

さらにFe-Mo二元系のBCC相は相互溶解度が大きいため，本研究ではこの固溶体の自由エネルギーをCE-CVM法によって計算した。Fig.8は，573 Kおよび1273 Kにおける混合自由エネルギーの計算結果を，Andersson²⁷⁾による熱力学的解析結果とともに示したものである。熱力学的解析とCVM計算の自由エネルギーの差は数 kJ/mol以内であり，いずれの結果もBCC相の二相分離傾向を示している。

Fig. 8.

Mixing free energies of BCC solid solution of Fe-Mo at 573 K and 1273 K by CVM calculation and the thermodynamic analysis²⁷⁾.

3・2・5　Fe-Mo-B三元系化合物の自由エネルギー

Fe-Mo-B三元系では，GAによる基底構造としてFe₃B((Fe,Mo)₃B)，Fe₂MoB₄，FeMo₂B₂が得られたが，それらの化合物の自由エネルギーの温度変化をFig.9に示した。FeMo₂B₂は温度の上昇とともにやや不安定化するのに対し，Fe₂MoB₄は安定化する傾向を示している。

Fig. 9.

Temperature dependence of free energy of the Fe-Mo-B ternary compounds.

3・3　有限温度におけるFe-Mo-B三元系状態図の計算結果と実験状態図との比較

3・2節において得られた化合物と固溶体の自由エネルギーを用いて，有限温度におけるFe-Mo-B三元系の相平衡を計算した。有限温度の相平衡計算は2・2・3節と同様にQhullコード¹⁶⁾を用いることで，準調和近似計算およびCE-CVMにより得た有限温度のエネルギーと組成の三次元座標からなるデータセットからconvex-hullを計算した。なお，BCC Fe-MoにおけるCE-CVMの計算結果については Mo濃度を0.01間隔で変えた組成-エネルギー点をconvex-hull の計算に導入した。計算結果の一例として，実験により相平衡が調べられているT=1273 K，1323 Kにおける理論状態図をFig.10(a)およびFig.11(a)に示した。Fig.10(b)，Fig.11(b)には，比較のために，著者らがCALPHAD法を用いて解析した計算状態図を実験値とともに示した¹⁷⁾。理論状態図は熱力学的解析による計算状態図と実験的知見をよく再現することから，本研究で提案する電子論ベースの新しい状態図計算手法の有効性が確認できる。しかし一部の相平衡には実験的知見との相違が見られるので，この点について指摘していきたい。

Fig. 10.

(a) Theoretical phase diagram for Fe-Mo-B ternary system at T = 1273 K and (b) Calculated phase diagram by thermodynamic analysis¹⁷⁾ with experimental data²⁰⁾.

Fig. 11.

(a) Theoretical phase diagram for Fe-Mo-B ternary system at T = 1323 K and (b) Calculated phase diagram by thermodynamic analysis¹⁷⁾ with experimental data²¹⁾.

Fig.10，Fig.11において太字で記したのは，理論状態図と計算状態図の間で，結晶構造や存在領域に違いがあった化合物を示している。これをまとめると次のようになる。

1)Fe-B二元系では，実験的に確認されていないFeB₂が出現する。

2)Mo-B二元系では，MoB₄の生成が実験的に確認されている。しかしこの構造は理論状態図には出現しない。またMoBの構造が異なっている。

3)Fe-Mo二元系では，Fe₂Moは実験温度範囲では生成しないが，理論状態図では出現する。

4)Fe-Mo-B三元系におけるMo近傍の組成領域における平衡については，実験からは1273 KにおいてMo₂BとFe₇Mo₆が平衡し²⁰⁾，1323 KにおいてMoとFeMo₂B₂が平衡すること²¹⁾が報告されている。一方，理論状態図ではいずれの温度においてもMoとFeMo₂B₂が平衡する。

5)Fe-Mo-B三元系では，FeMo₈B₁₁(τ₃)が実験的に確認されているが，理論状態図では出現しない。これは，3・1・2(2)で述べた理由からここでは議論の対象としない。

6)本計算結果では，化合物は全て化学量論化合物として扱い，原子同士の固溶を考慮しなかった。従ってこの点における実験結果との違いが認められるが，ここでは議論の対象としない。

続いて，1)から4)に関して，FeB₂，MoB₄，MoB(I4₁/amd)，Fe₂Moの自由エネルギーをどの程度変えることで現実の相平衡を再現できるのかを検討した。具体的には，これらの四種類の構造の自由エネルギーを0.1 kJ/mol間隔で変化させて相平衡を計算した。Fig.12はこの操作によって再現された1273 Kにおける理論状態図である。FeB₂は+2.4 kJ/mol，Fe₂Moは+0.8 kJ/mol，Mo₂Bは+0.6 kJ/mol，MoBは+0.8 kJ/mol，MoB₄は-0.9 kJ/mol分のエネルギーを変化させることによって実験状態図が再現される。この検討より，実験状態図と理論状態図の間の自由エネルギー差は数 kJ/mol以内であり，本研究から得られた自由エネルギーは現実のものに非常に近いことが示唆される。

Fig. 12.

Theoretical phase diagram with modified energies of FeB₂, MoB₄, MoB (I4₁/amd), Mo₂B and Fe₂Mo structures at 1273 K. The numerals in the diagram represent the amount of change (in kJ/mol) in the energies given to the structures.

そこで次に，このわずかなエネルギーの乖離の原因とその解決方法を検討する。まず乖離の原因の一つとして，化合物の固溶の効果が指摘される。例えば，Fe₇Mo₆には固溶領域が報告されているが，この固溶を考慮すると配置エントロピー効果によりこの構造が熱力学的に安定化する。このFe₇Mo₆の安定化によってFeとFe₇Mo₆が平衡することでFe₂Moの生成が遮られる可能性がある。また，磁気変態によるエントロピー効果や格子振動の非調和項の検討も必要である。磁気変態によるエントロピーは磁気変態温度が重要になるが，これを第一原理計算から求める手法^28–30)が提案されている。一方，非調和項は高次の力定数を求める方法^31,32)や第一原理分子動力学法^33,34)から評価する手法がある。これらの手法を導入することで自由エネルギー並びに理論状態図の計算精度を上げることができると期待される。

4. 結言

本研究では遺伝的アルゴリズムに基づく第一原理計算によってFe-Mo-B三元系の基底構造を推定し，それらの構造の自由エネルギーを電子論計算と統計熱力学的手法によって評価した。さらにその結果を用いて有限温度における状態図を理論的に構築した。その結果，次のような結論が得られた。

(1)計算により得られた多くの基底構造の空間群や組成は実験的知見とよく対応したが，その一致は完全ではなかった。しかし不一致が見られた構造においても，基底状態よりもわずか数 kJ/molだけ不安定な構造まで含めると，この手法によりほぼ全ての構造の出現を予測することが可能であることが明らかになった。

(2)基底構造と基底に近い不安定構造の自由エネルギーを計算したところ，基底状態での不安定構造が高温で安定化する，いわゆる構造相変態を起こす場合があることがわかった。有限温度で計算された理論状態図を実験結果と比較すると，両者は概ね一致していた。一致の見られなかった相領域でも，数 kJ/mol以内のエネルギー差で実験結果と整合することがわかった。このことから，有限温度の平衡相として現れる可能性がある構造をconvex-hullからのエネルギー差から選別し，理論状態図を作成する方法が有用であることが示唆された。

(3)本研究で提案する有限温度の理論状態図の計算精度を向上させるために，1)単位胞の原子数，遺伝的アルゴリズムの世代数などの計算条件の精査，2)原子振動の非調和性，磁気エントロピー効果の考慮などの計算手法の選択，3)固溶体や液相などの不規則構造の自由エネルギーの考慮，などが今後の課題として挙げられる。

謝辞

この成果は，JSPS 科研費16H02387の助成を受けて得られたものです。ここに謝意を表します。

文献

• 1)  N.Saunders and A.P.Miodownik: Calphad (Calculation of Phase Diagrams): A Comprehensive Guide, Pergamon Materials Series, Pergamon, Oxford, UK, (1998), 7.
• 2)   G. Kresse and  J. Furthmüller: Phys. Rev. B, 54(1996), 11169.
• 3)   G. Kresse and  D. Joubert: Phys. Rev. B, 59(1999), 1758.
• 4)   P.E. Blöchl: Phys. Rev. B, 50(1994), 17953.
• 5)   J.P. Perdew,  A. Ruzsinsky,  G.I. Csonka,  O.A. Vydrov,  G.E. Scuseria,  L.A. Constantin,  X. Zhou and  K. Burke: Phys. Rev. Lett., 100(2008), 136406.
• 6)   J.P. Perdew,  K. Burke and  M. Ernzerhof: Phys. Rev. Lett., 77(1996), 3865.
• 7)  A.A.Maradudin, E.W.Montroll, G.H.Weiss and I.P.Ipatova: Theory of Lattice Dynamics in the Harmonic Approximation, 2nd ed., Academic Press, New York, (1971), 5.
• 8)   A. Togo,  F. Oba and  I. Tanaka: Phys. Rev. B, 78(2008), 134106.
• 9)   A. Togo,  L. Chaput,  I. Tanaka and  G. Hug: Phys. Rev. B, 81(2010), 174301.
• 10)   D. De Fontaine: Solid State Phys., 34(1979), 73.
• 11)   M.H.F. Sluiter,  C. Colinet and  A. Pasturel: Phys. Rev. B, 73(2006), 174204.
• 12)   R. Kikuchi: Phys. Rev., 81(1951), 988.
• 13)   A.R. Oganov and  C.W. Glass: J. Chem. Phys., 124(2006), 244704.
• 14)   A.R. Oganov,  A.O. Lyahov and  M. Valle: Acc. Chem. Res., 44(2011), 227.
• 15)   A.O. Lyakhov,  A.R. Oganov,  H.T. Stokes and  Q. Zhu: Comput. Phys. Commun., 184(2013), 1172.
• 16)   C.B. Barber,  D.P. Dobkin and  H. Huhdanpaa: ACM Trans. Math. Softw., 22(1996), 469.
• 17)   K. Takahashi,  K. Ishikawa,  M. Fujioka,  M. Enoki and  H. Ohtani: Tetsu-to-Hagané, 106(2020), 310 (in Japanese).
• 18)   T.A. Velikanova,  M.V. Karpets,  S.U. Artyukh,  S.O. Balanetskii,  V.M. Petyukh,  P.G. Agraval and  M.A. Turchanin: Powder Metall. Met. Ceram., 50(2011), 442.
• 19)   J. Cieslak,  S.M. Dubiel,  J. Przewoznik and  J. Tobola: Intermetallics, 31(2012), 132.
• 20)   E.I. Gladyshevskii,  T.F. Feclotov,  Yu.B. Kuz’ma and  R.V. Slolozdra: Poroshk. Metall., 4(1966), 55 (in Russian).
• 21)   A. Leithe-Jasper,  H. Klesnar,  P. Rogl,  M. Komai and  K. Takagi: J. Jpn Inst. Met., 64(2000), 154 (in Japanese).
• 22)   L.G. Voroshnin,  L.S. Lyakhovich,  G.G. Panich and  G.F. Protasevich: Met. Sci. Heat Treat., 12(1970), 732.
• 23)   T. Van Rompaey,  K.C. Hari Kumar and  P. Wollants: J. Alloy. Compd., 334(2002), 173.
• 24)   A.G. Van Der Geest and  A.N. Kolmogorov: Calphad, 46(2014), 184.
• 25)   A.T. Dinsdale: Calphad, 15(1991), 317.
• 26)   M.Y. Lavrentiev,  D. Nguyen-Manh and  S.L. Dudarev: Phys. Rev. B, 81(2010), 184202.
• 27)   J.-O. Andersson: Calphad, 12(1988), 9.
• 28)   W. Nolting,  A. Vega and  Th. Fauster: Z. Phys. B, 96(1995), 357.
• 29)   N. Rosengaard and  B. Johansson: Phys. Rev. B, 55(1997), 14975.
• 30)   F. Körmann,  A. Dick,  B. Grabowski,  B. Hallstedt,  T. Hickel and  J. Neugebauer: Phys. Rev. B, 78(2008), 033102.
• 31)   J.C. Thomas and  A. Ven: Phys. Rev. B, 88(2013), 214111.
• 32)   T. Tadano and  S. Tsuneyuki: Phys. Rev. B, 92(2015), 054301.
• 33)   A. Glensk,  B. Grabowski,  T. Hickel and  J. Neugebauer: Phys. Rev. Lett., 114(2015), 195901.
• 34)   B. Grabowski,  Y. Ikeda,  P. Srinivasan,  F. Körmann,  C. Freysoldt,  A.I. Duff,  A. Shapeev and  J. Neugebauer: NPJ Comput. Mater., 5(2019), 80.