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. Introduction

In controlling to maximize the functions of various materials, the role which phase diagram has played is very large. Since the microstructure of metals was examined by optical microscopy in the mid-19th century, phase diagrams of alloys and ceramics have been studied by utilizing experimental techniques such as microscopic observation and thermal analysis. However, with the reconstruction of the industry after the Second World War, the metal materials that are familiar to the public also became multicomponent, and therefore such time-consuming experimental phase diagrams failed to cope with the actual material design in terms of labor and speed. In the CALPHAD (Calculation of Phase Diagrams) method,¹⁾ which was newly proposed to incorporate with such problems, the measured values on phase boundaries and thermodynamic properties are analyzed using appropriate thermodynamic models, and the thermodynamic properties of the target system are parameterized. Since the physical properties of the multicomponent system can be reasonably predicted by the thermodynamic extrapolation method, it has the advantage of accurately predicting the actual phase equilibria. Then this method has been greatly developed in the field of material development by utilizing the advantages. However, it is very difficult to apply this technique to the system which has not been clarified in experiments, metastable regions, and non-equilibrium phases. Therefore, referring to the Fe–Mo–B ternary system in this study, the thermodynamic properties of this alloy system are evaluated by the electronic theory calculations based on the first-principles. In addition, using the results, this paper proposes a new technique to theoretically construct alloy phase diagrams at finite temperature, and aims at arranging the present problems.

2. Calculation Methods

In this chapter, we first describe the conditions of the first-principles calculations used in this study, and then explain the calculation principle of the formation energy of a perfectly ordered structure at absolute zero. In addition, we explain a method for calculating free energy based on lattice vibration calculation of atoms and quasi-harmonic approximation as well as the procedure by cluster expansion and variational method as a technique for extending this result to a thermodynamic property at finite temperature. Then, the technique which searches the phase equilibria of the ground state by the genetic algorithm is also described.

2.1. First-principles Calculation Method

In this study, the first-principles calculations were performed by means of the calculation codes Vienna Ab initio Simulation Package (VASP),^2,3,4) using the GGA-PBE for exchange-correlation functional.^5,6) The cutoff energy of the plane wave was assumed to be 400 eV. In the cluster expansion and genetic algorithms described in Sections 2.2.2., and 2.2.3., the k-point sampling uses 17×17×17 grid points for primitive BCC-Fe structures. On the other hand, the number of k-point samplings of the supercell used in the calculation of the quasi-harmonic approximation method is described in 2.2.1. In all calculations, the spin-polarity was considered, and the initial values of magnetic moments of Fe, Mo, and B were set to +2.5 μ_B, 0 μ_B, and 0 μ_B, respectively.

From the total energy obtained from the first-principles calculations, we can evaluate the enthalpy of formation of compounds at absolute zero. As an example, the formation enthalpy (H_AnBm) of a binary compound A_nB_m is represented by Eq. (1).   

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

Here, E_AnBm represents the total energy of the binary compound A_nB_m, and E_A, E_B are the total energy of the elements A and B as the reference states obtained from the first-principles.

2.2. Calculation of Thermodynamic Properties at Finite Temperature

Since the thermodynamic properties evaluated by first-principles calculations are values at absolute zero degrees, it is necessary to evaluate the free energy in order to calculate the phase equilibria at finite temperatures. Components with large contributions at finite temperatures are roughly divided into two types, one is due to lattice vibration and the other is due to randomness of atomic sites due to solid solution. For the former, by calculating the density of states of phonons using first-principles calculations and further incorporating the effect of volume expansion, it is possible to take into account the change of the formation energy with temperature. For the latter, a method of expanding the calculated cohesive energy in clusters for a large number of ordered structures and determining the free energy of solid solutions at finite temperatures by the cluster variational method is used. This section outlines these methods.

2.2.1. Free Energy Calculation of Compounds Considering Lattice Vibration

In order to evaluate the free energy of finite temperature, it is necessary to investigate the effects of lattice vibrations and disordering of atomic configurations. The free energy and specific heat at finite temperatures of pure elements and ordered compounds have large contributions from phonons. By obtaining the force acting between atoms from first-principles calculations, it is possible to calculate the contribution of this phonon.⁷⁾

As a more detailed calculation procedure, the restoring force generated by the minute displacement of atomic position from the equilibrium position is evaluated using the supercell which extended the crystal lattice. From these restoring forces, the second order force constants for atomic displacements are calculated, and the wave number dependence of the frequency (ω) of the phonons is obtained. This technique is called the harmonic oscillation approximation. From the resulting phonon density of states, the Helmholtz free energy, including the effect of lattice vibration at finite temperature, is given by Eq. (2).   

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

U₀ represents the energy in the equilibrium position, and q_j is the wave vector in the j-th phonon mode. The second term is called the zero‐point energy and originates from the effect that atomic vibrations exist even at absolute zero because of the uncertainty principle. Gibbs free energy is calculated by calculating the volume-energy curve from Helmholtz free energy (F_V) in volume V, and by obtaining the volume V that becomes the energy minimum at a constant pressure condition.   

$$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 code^8,9) was used for quasi-harmonic approximation calculations. In Section 3.2., we show the free-energy calculations for 19 ordered structures, including BCC-Fe. The sizes of supercells and mesh sizes of k-points used in first-principles calculations of their ordered structures are summarized in Table 1. The quasi-harmonic approximation calculation was carried out by preparing a cell with isotropically changing volume after relaxing all atomic configurations including shape and volume of the cell for each structure.

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. Cluster Expansion and Variation Method

In order to calculate the phase equilibrium at finite temperature of a system consisting of many components, it is necessary to consider the mixing free energy in addition to the contribution of phonons. Cluster Expansion and Variation Method (CE-CVM) is a technique for calculating mixing free energies from the effective cluster interaction (ECI) of individual clusters and the entropy of their configurations.¹⁰⁾

Specifically, by introducing spin operators σ and fitting different values of spin operators for each element type, each ordered structure R is represented by a spin arrangement R = (σ₁, σ₂, σ₃, ... σ_i, ...). Here, i is an exponent representing the atomic site position of the ordered structure. Also, clusters can be represented by combinations of spin operators extracted from different sites, for example, (σ_i, σ_j) is a pair of sites (i, j), and (σ_i, σ_j, σ_k) is a three-body cluster of sites (i, j, k). Furthermore, a correlation function φ_α represented by Eq. (4) is introduced by taking the average value of spin products for clusters such as points, pairs, and triplets.   

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

α is a subscript representing the type of cluster, and N_point, N_pair, N_tri are the total numbers of points, pairs, and three-body clusters contained in the ordered structure. Since the value of this correlation function has a one-to-one correspondence with the concentration of clusters contained in the ordered structure, it is possible to express the energy of the ordered structure E_R by the sum of the products of the correlation function and ECI (e_null, e_point, e_pair, e_tri,…).   

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

Although the Eq. (5) can theoretically reproduce the energy of the ordered structure strictly by using an infinite type of cluster, it is necessary to terminate the cluster to be used to a finite number in the actual calculation. Since the interaction of clusters tends to be stronger for short distances and weaker for long distances, the type of clusters considered generally employs a technique of providing a threshold for the coupling distance of clusters and selecting only clusters of smaller sizes. The ECI is summed with the correlation function for the clusters (sub-clusters) contained within the maximum size of the clusters α_max defined therein.   

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

Since the correlation function φ_α is obtained from the ordered structure and the total energy of the left side is obtained from the first-principles calculation, the unknown number is only ECI e_α. However ECI is determined by the least-squares method by providing a large number of relations between different ordered structures and energies. Since there is an arbitrariness in the choice of the sub-cluster combinations used in this approach, a method of determining the combinations of basis clusters by performing cross-correlation error tests to minimize forecast errors is employed.¹¹⁾

Once ECI is determined, the energy of any atomic arrangement can be determined within the accuracy of cluster expansion without first-principles calculations. Also, by considering the number in the case of the arrangement of clusters, the free energy including the entropy of the arrangement can be expressed as follows.   

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

The second term in Eq. (7) is an entropy term, and γ_αS_α shows the entropy contribution from the cluster α using Kikuchi-Barker coefficients.¹²⁾ The free energy is calculated by applying the variational method to Eq. (7) to obtain a correlation function that minimizes the energy.

iCVM codes¹¹⁾ were used for CE-CVM calculations. In BCC Fe–Mo calculations in 2.3.4, the energies of 149 ordered structures with different configurations were used for the analysis. And, the energy for each ordered structure was evaluated after relaxing all atomic configurations including shape and volume of the cell once. Four-body clusters included up to the fifth-nearest-neighbor distance were chosen as the max clusters, and ECI was assigned to 35 kinds of sub-clusters included in them. The cross-correlation error test yielded a forecast error of 13 meV/atom.

2.2.3. Ground Structure Search

In order to make a phase diagram of a certain system, the identification of the ordered structure which constitutes the phase region is most important. However, in the first-principles calculation, it is difficult to predict a stable ordered structure even if the composition of a substance is determined. This problem is mainly caused by two factors. One is that the number of dimensions describing the structure is 3N+3 (N: number of atoms contained in the unit cell) and large degrees of freedom must be handled, and it is impossible to calculate the energy of all structures comprehensively. The other is that even when the structural relaxation is carried out by calculating the force on the atoms in the first-principles, energy barriers prevent reaching the most stable state because there are many local stable states.

Recently, a computational method to predict the most stable structures from first-principles calculations has been adopted by using genetic algorithm (GA) for these problems.^13,14,15) Genetic algorithms mimic Darwin’s evolutionary theory and are computational methods employing natural selection, genetic and mutating. The calculation is carried out in the following flow.

(1) Energy calculations of randomly created plural ordered structures are carried out, and among them, those with low energy are selected preferentially.

(2) Components such as the local atomic arrangement of those structures are regarded as the genetic, and the structure group of the next generation in which they are crossed and mutated is made.

(3) New generations of energy calculations are performed and low energy ones are re-selected.

(4) Search for the most stable structure by repeatedly executing (2) and (3) and updating the low-energy structure group.

Thus, by imitatively incorporating matings and mutants, it is possible to efficiently explore the stable structure beyond the problem of energy barriers. In addition, the lowest-energy convex-hull consisting of the compositional-energy points of the group of structures is calculated. The structure of the vertex of the convex hull is the ground structure at absolute zero, and the phase diagram based on the first-principles calculation is created.

USPEX^13,14) was used to search the ground state structures of the Fe–Mo–B ternary system by GA. As a specific calculation condition, the number of structures in each generation was made to be 400, and a total of 280 structures obtained by the stable structure to be passed over to the next generation and the structure obtained by the crossing were prepared, and in addition, 80 structures created by the mutation and 40 random structures newly created were added to create a new generation. Qhull code¹⁶⁾ was used to calculate convex-hull, which was generated repeatedly up to 10 generations. A three-dimensional coordinate convex-hull consisting of x_Fe, x_Mo and E of all structures was obtained.

3. Calculation Results and Discussion

3.1. Ground State Diagram of Fe–Mo–B Ternary System

The result of the ground structure search of the Fe–Mo–B ternary system is shown in Fig. 1. The black circles in the figure are the structures obtained in this study, and the square marks are the compounds that have been confirmed experimentally. The solid line is the phase boundary of the two phases composed of the stable compounds searched for by GA, and the triangular region composed of these lines is the invariant reaction in the ground state of this ternary system.

Fig. 1.

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

3.1.1. Compounds Observed in the Experimental Phase Diagram

The crystalline structures of these compounds are summarized in Table 2. This subsection describes the outline of the compounds which constitute the Fe–Mo–B ternary system confirmed by the experiments. In the Fe–Mo–B ternary system, 16 kinds of binary and ternary compounds have been confirmed. Their chemical formulae are shown in Table 2. In this table, the structure names used in the literature¹⁷⁾ were adopted. And, for the structure in which the expression is different from the composition formula, the structure name is also written in parentheses. Mo₂B₃ (MoB_2-x) is referred to as the MoB₂ (MoB_2-x) in the literature,¹⁷⁾ however, this phase appears near the composition of Mo:B=2:3. In addition, the structure search of GA confirmed the corresponding stable structure in the composition of Mo₂B₃. Thus, in Table 2, the name of this structure was changed to Mo₂B₃.

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

Compounds in the Fe–B binary system and Mo–B binary system have both small solid solution widths. On the other hand, in the Fe–Mo binary system, four types of compound have been identified: Fe₂Mo (Laves C14), Fe₇Mo₆ (μ), Fe₃Mo₂ (R) and FeMo (σ). Fe₂Mo (Laves C14) is a stoichiometric compound. In contrast, Fe₇Mo₆ (μ) is a nonstoichiometric compound that appears in the compositional range of 39–44% Mo, Fe₃Mo₂ (R) is 34–39% Mo, and FeMo (σ) is 43–57% Mo, respectively. The ternary compounds Fe₃B (τ₁), FeMo₂B₂ (τ₂), FeMo₈B₁₁ (τ₃), Fe₂MoB₄ (τ₄) show the homogeneity regions by replacing Fe and Mo atoms.

3.1.2. Comparison with the Compounds Obtained by Ground Structure Search

Space groups of binary and ternary compounds obtained by the ground structure search with GAs were defined and compared with the experimental stable structures in Table 2. The table also includes some of the structures determined to be metastable in the GA. For the structures, calculations for the differences in energy from convex-hull (ΔE_{from hull}) are also provided. A structure with a ΔE_{from hull} of 0 means the ground structure, and a structure with finite values means the metastable structure. Also, FeB₂ (Pnma) and Fe₁₃Mo (P1) are structures not found experimentally.

The compounds shown in bold in Table 2 were consistent with compounds in which the ground structures of GA were experimentally confirmed. However, other compounds did not correspond to the experimental findings when only the ground structure is compared. Then, this difference was examined from the thermodynamic viewpoint.

(1) Compounds of Which Ground Structures Including the Metastable Structures Agree with the Experimental Results

First, for FeB and MoB, the I4₁/amd is the ground structure according to GA, but the compound which becomes slightly less stable than the structure is consistent with the experimental results. The energy differences from convex-hull of these metastable structures are both around 1 kJ/mol, and the structural phase transformation may occur at finite temperatures. Discussions on this point will be made in 3.2.2. and 3.3.3. in the next section.

Mo₂B, Mo₃B₂, MoB₂ (Mo₂B₄), MoB₄, Fe₃B (τ₁) are not the ground structures based on GA analysis, however, the structures to correspond to the experimental results were found in the slightly higher energy structures. The space group of Mo₂B₃ (MoB_2-x) differs from the experimental knowledge. However, as shown in Fig. 2, Imm2 corresponds to structures where B and vacancies in P6/mmm are ordered in a ratio of 3:1. In X-ray diffraction, since the scattering cross section of B with a small atomic number is small, the atomic position of B may not be sufficiently detected. Therefore, it was proven that it became P6/mmm when the space group was judged by excluding B atoms in the structures of Imm2. Furthermore, for the fully ordered structures of P6/mmm with MoB₂ compositions, the energy-difference from convex-hull was +14.3 kJ/mol, whereas that of Imm2 was +1.5 kJ/mol. From these considerations, we conclude that the structure of Mo₂B₃ (MoB_2-x) in this study is Imm2. All six types of Mo₂B, Mo₃B₂, MoB₂ (Mo₂B₄), Mo₂B₃ (MoB_2-x), MoB₄, Fe₃B (τ₁) have energy-differentials from convex-hull less than 3 kJ/mol and may appear as an equilibrium phase at finite temperatures. This discussion will be carried out in Section 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) Compounds of Which Ground Structures Did Not Agree with the Experimental Results

For Fe₃Mo₂ (R), FeMo (σ), FeMo₈B₁₁ (τ₃) there was no corresponding structure in the set of structures searched by GA. Among them, R is composed of 159 atoms,¹⁸⁾ and σ is a disordered structure of a large unit cell composed of 30 atoms.¹⁹⁾ In this study, R was not found because the structure search was conducted by limiting the total number of atoms in the unit cell to 32. In addition, since σ is also a complex structure in which the unit cell consists of 30 atoms, it is considered necessary to further increase the number of generations in order to search for this structure by GA. On the other hand, the crystal structure of FeMo₈B₁₁ (τ₃) is unknown, and there are reports that it is a MoB (R3m)-type structure.^20,21) Although the structure obtained by GA can be regarded as R3m, thermodynamic analysis using the first-principles in the literature¹⁷⁾ has confirmed that the (Mo,Fe)B type solid solution with R3m structure is thermodynamically less stable. For these reasons, this study does not discuss the structures of Fe₃Mo₂(R), FeMo (σ), FeMo₈B₁₁ (τ₃), and they have been excluded from the free-energy evaluation in the next section.

(3) Compounds Not Confirmed by Experiments

FeB₂ and Fe₁₃Mo shown at the end of Table 1 are the ground structures for which the corresponding stable structure has not been confirmed by experiments. Concerning FeB₂, formation of the structure with P6/mmm was reported experimentally by Voroshnin et al.²²⁾ However, subsequent experiments did not reproduce the formation of this phase, and it was concluded that the stability is low.²³⁾ On the other hand, the results from GA show that the structure of the space group Pmna is stable, which is different from Voroshnin et al. The Pmna type FeB₂ was confirmed stable also by the ground state search of the Fe–B binary system by first-principles calculation,²⁴⁾ and the results of this study are consistent with this paper. Since the thermodynamic stability of this structure at finite temperatures has not been discussed so far, the free energies of finite temperatures are calculated and discussed in Sections 3.2.2 and 3.3, comparing with experimental results. On the other hand, as shown in Fig. 3, Fe₁₃Mo corresponds to the structure in which Fe atoms in the ordered structure of BCC-Fe are substituted with Mo. In Fig. 3, Fe₁₃Mo is expanded to a 2×2×1 supercell, where all atomic positions are composed of the BCC-lattice illustrated in this figure. Therefore, this structure was treated as being included in the BCC solid solution in 3.2.4.

Fig. 3.

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

3.2. Calculated Free Energy of Each Phase

In this section, to investigate the phase equilibrium at finite temperature, free energy calculations based on quasi-harmonic approximation are performed for the ordered structure groups obtained in the Section 3.1. For the solid solution, the free energies of the primary solid solutions were evaluated by CE-CVM method given in Section 2.2.2. Considering that the solubility of B in Fe and Mo is small, only the free energies of solid solutions in the Fe–Mo binary system are calculated.

3.2.1. Free Energy of Pure Element

Free energy calculations by quasi-harmonic approximation are performed for BCC-Fe, BCC-Mo and β‐B. The results are shown in Fig. 4. Temperature is shown on the horizontal axis and Gibbs free energy per unit atom on the vertical axis. The free energy decreases monotonically with increasing temperature. β-B has the smallest reduction in energies with increasing temperatures compared to BCC-Fe, BCC-Mo. This is probably because the vibrational entropy becomes relatively small compared with BCC-Fe and BCC-Mo, because the lattice vibrations are faster than BCC-Fe, BCC-Mo due to the small mass of B atoms, and hence the density of states of low-energy phonons becomes small. In this diagram, the temperature dependency of the free energies of the thermodynamic database by Dinsdale²⁵⁾ are shown as a white circle. For the free energy of the thermodynamic database, the temperature change of the free energy was plotted on the basis of the value of the free energy calculated by the first-principles method at T = 273 K. The calculated results for both pure elements well reproduce the temperature dependence of the database. From this result, the anharmonic term of lattice vibration that is not considered in the present study is small enough. However, for BCC-Fe, there is a tendency for the databases and calculated values to deviate slightly from the temperature range above T = 1000 K. This is due to the entropy effects of the ferromagnetic-paramagnetic transition of BCC-Fe at T = 1043 K, resulting in a slight stabilization of the database values. Since this effect is not included in this study, a deviation from the experimental values is recognized. However, as the difference from the database is small, the phase equilibria of the Fe–Mo–B ternary system are calculated without taking this effect into account in this study. FCC-Fe is stable in the temperature range of 1185 to 1667 K, however this phase was not considered in the present study because of its difficulty in the treatment of entropy effect by magnetism. To handle this entropy effect, various calculations with different arrangement patterns of the spin polarization are necessary in addition to the phonon calculations, which requires large amount of computational resources.²⁶⁾

Fig. 4.

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

3.2.2. Free Unergies of Fe–B Binary Compounds

We evaluated free energies of four compounds FeB₂, FeB (I4₁/amd), Fe₂B obtained as the ground structures of GA analysis as well as the FeB (Pnma) which are metastable compound at absolute zero degrees. The free energies of the Fe–B binary compounds based on the free energies of pure elements are shown in Fig. 5.

Fig. 5.

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

The free energy of formation of FeB₂, which was not confirmed in the experimental phase diagram, increased slightly with increasing temperatures. On the other hand, in other compounds, thermodynamic stability increases with increasing temperatures. According to the analysis of GA, two types of structures of I4₁/amd and Pnma were obtained for FeB. The calculated free energy of each structure in Fig. 3 shows that the phase transformation takes place between these two structures at about T = 400 K. In other words, under T = 400 K, I4₁/amd is stable, whereas above T = 400 K, the experimental Pnma structures are stable. This result is consistent with the experimental findings.

3.2.3. Free Energies of Mo–B Binary Compounds

Figure 6 shows the change in free energies of the Mo–B binary compounds. We evaluated the free energies of seven types of compound: Mo₂B, MoB (I4₁/amd), MoB₂ as the ground structures of GA analysis, and MoB₄, Mo₂B₃, MoB (Cmcm), Mo₃B₂, which were metastable structures at absolute zero degrees. From this figure, it can be confirmed that the free energies of formation of all compounds do not change significantly with temperature and retain their respective thermodynamic stability up to high temperatures.

Fig. 6.

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

Although the stable structure of MoB by experiments is Cmcm, the I4₁/amd obtained as the ground state in this study remains slightly more stable than Cmcm in the calculated temperature range. However, the difference is about 1 kJ/mol, and is extremely small. Therefore, there is a possibility that the stability is easily replaced by the calculation conditions, the potentials to be used in first-principles calculations, the anharmonicity of lattice vibration, and the entropy effect at higher temperatures.

3.2.4. Free Energies of Fe–Mo Binary Compounds

In the Fe–Mo binary system, four types of compound have been identified from the experiments: Fe₂Mo, Fe₇Mo₆, Fe₃Mo (R) and FeMo (σ). However, as discussed in Section 3.1, no corresponding structures have been obtained for R and σ in GA analysis, so the free-energy of Fe₂Mo, Fe₇Mo₆ has been computed. Since Fe₁₃Mo can be regarded as a solid solution of Mo in BCC-Fe, quasi-harmonic approximation calculations were not carried out for this structure. Figure 7 shows temperature dependence of free energy of formation for Fe₂Mo and Fe₇Mo₆. The free energies of both compounds are between −3 kJ/mol to −6 kJ/mol, which are smaller in absolute values than the free energies for the formation of the Fe–B and Mo–B binary compounds.

Fig. 7.

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

Since the BCC phase of the Fe–Mo binary system has a large mutual solubility, the free energies of this solid solution are calculated by CE-CVM method in this study. Figure 8 shows the calculated results of the mixing free energies at 573 K and 1273 K compared with those in the thermodynamic analysis by Andersson.²⁷⁾ The difference in free energies between the thermodynamic analysis and the CVM calculations is within a few kJ/mol. Both these results show a characteristic feature of two-phase separation tendency in the BCC phase.

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. Free Energies of Fe–Mo–B Ternary Compounds

In the Fe–Mo–B ternary system, Fe₃B ((Fe,Mo)₃B), Fe₂MoB₄, FeMo₂B₂ are obtained as the ground structures by GA, and the temperature change of free energies of these compound is shown in Fig. 9. FeMo₂B₂ becomes somewhat less stable with increasing temperatures, while Fe₂MoB₄ and Fe₃B tend to be stabilized.

Fig. 9.

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

3.3. Calculated Results of Fe–Mo–B Ternary Phase Diagram at Finite Temperatures and Comparison with Experimental Phase Diagram

Using the free energies of compounds and solid solutions obtained in the Section 3.2., the phase equilibria of the Fe–Mo–B ternary systems at finite temperatures are obtained. The phase equilibrium calculations employ the method of calculating convex-hull by Qhull code¹⁶⁾ as noted in Section 2.2.3. By inputting the three-dimensional coordinates of the energies and compositions obtained by the quasi-harmonic approximation calculation as well as CE-CVM into Qhull, the convex-hull of this ternary system was calculated. The composition-energy point with the Mo concentration varied at 0.01 intervals was introduced into convex-hull calculation for CE-CVM calculation for BCC Fe–Mo solid solution. As an example of the calculation results, the theoretical phase diagrams for T = 1273 K, 1323 K are shown in Figs. 10(a) and 11(a). For comparisons, calculated phase diagrams by the authors group using CALPHAD method are shown in Figs. 10(b) and 11(b), together with experimental values. The theoretical phase diagram well reproduces the calculated phase diagram and experimental findings by thermodynamic analysis, thus confirming the effectiveness of the new electronic theory-based phase diagram calculation method proposed in this study. However, some phase equilibria differ from experimental findings, so we would like to point out this point.

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.²¹⁾

The bold letters in Figs. 10 and 11 indicate compounds with different crystal structures and compositions between theoretical and calculated phase diagrams. This is summarized as follows.

1) In the Fe–B binary system, FeB₂ which is not confirmed experimentally appears.

2) In the Mo–B binary system, the formation of MoB₄ is experimentally confirmed. However, this structure does not appear in the theoretical phase diagram. MoB is of different structure.

3) In the Fe–Mo binary system, Fe₂Mo is not formed in the experimental temperature range, but it appears in the theoretical phase diagram.

4) Regarding the equilibrium on the Mo rich side, it was reported that Mo₂B is in equilibrium with Fe₇Mo₆ at 1273 K²⁰⁾ while Mo and FeMo₂B₂ equilibrate at 1323 K²¹⁾ in experimental examination. On the other hand, the theoretical phase diagrams show that Mo and FeMo₂B₂ in equilibrium at both temperatures.

5) In the Fe–Mo–B ternary system, FeMo₈B₁₁ (τ₃) has been confirmed experimentally, but it does not appear in the theoretical phase diagram. For the reason described in 3.1.2.2., it is not discussed here.

6) In this calculation, all compounds were treated as stoichiometric compounds, and solid solution of atoms was not considered. Therefore, although there is a difference from the experimental results in this respect, it is not the subject of discussion here.

As for (1) to (4), it was examined whether the real phase equilibria can be reproduced by changing the free energies of FeB₂, MoB₄, MoB (I4₁/amd), and Fe₂Mo. Specifically, the phase equilibria were calculated by changing the free energies of these four structures at intervals of 0.1 kJ/mol. Figure 12 is a phase diagram at 1273 K which is reproduced by this operation. The change of free energy of FeB₂ is +2.4 kJ/mol, Fe₂Mo is +0.8 kJ/mol, Mo₂B is +0.6 kJ/mol, MoB is +0.8 kJ/mol, and MoB₄ is −0.9 kJ/mol to reproduce the experimental phase diagram. This study shows that the free energy difference between the experimental and theoretical phase diagrams is within a few kJ/mol, and that the free energies obtained from this study are very close to that of the real ones.

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.

Then, we will examine the causes of this slight energy deviation and how to solve them. First, as one of the causes of deviation, the effect of solid solution of the compound is pointed out. For example, homogeneous region has been reported in Fe₇Mo₆, and thus the configuration entropy effect makes the structures thermodynamically stable when this solid solution of Fe and Mo atoms is considered. This stabilization can prevent Fe₂Mo from appearing because of the equilibrium between Fe and Fe₇Mo₆. It is also necessary to investigate the entropy effect due to the magnetic transition and the anharmonic term of the lattice vibration. The entropy due to the magnetic transition is related to the magnetic transition temperature, and there are techniques^28,29,30) which obtains this temperature from the first-principles calculations. On the other hand, the anharmonic term can be evaluated from methods for obtaining higher-order force constants^31,32) or first-principles molecular dynamics methods.^33,34) It is expected that the accuracy of free energy and theoretical phase diagram can be improved by introducing these techniques.

4. Conclusions

In this study, the ground structures of the Fe–Mo–B ternary systems were estimated by first-principles calculations based on genetic algorithms, and the free energies of their structures are evaluated by means of electronic calculations and statistical thermodynamic techniques. In addition, the phase diagram at finite temperature was theoretically constructed using the results. The following conclusions were obtained.

(1) The space groups and compositions of many ground structures obtained by the calculations correspond well with the experimental findings, but the agreement was not perfect. However, it was clarified that the appearance of almost all structures could be predicted by this technique, by including up for the metastable structure by only a few kJ/mol than the ground state.

(2) Calculations of the free energies of the ground structures and nearly ground structures show that there are some cases of so‐called structural phase transformations in which the metastable structures in the ground state are stabilized at high temperatures. When the theoretical phase diagram calculated at finite temperature is compared with the experimental results, both agree almost well each other. Even in the phase region where no agreement was found, it was clarified that the results agreed with the experimental results by changing the energies within a few kJ/mol. This suggests the usefulness of the theoretical phase diagram construction, in which structures appearing as equilibrium phases at finite temperatures are sorted out from energy differences in convex-hull.

(3) In order to improve the calculation accuracy of the theoretical phase diagram at finite temperature proposed in this study, the following future subjects are pointed out: (1) the examination of calculation conditions such as the number of atoms in the unit cell and the number of generations in genetic algorithms, (2) the selection of calculation methods such as the anharmonicity of atomic vibration and the consideration of the magnetic entropy effect, (3) the consideration of free energy of disordered structures such as solid solutions and liquid phase.

Acknowledgments

The authors gratefully acknowledge the financial support by JSPS KAKENHI (grant number, 16H02387).

References

• 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.  Lyakhov 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.
• 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.